Combined Spectral and Spatial Modeling of Corn Yield Based on Aerial Images and Crop Surface Models Acquired with an Unmanned Aircraft System

Abstract

Precision Farming (PF) management strategies are commonly based on estimations of within-field yield potential, often derived from remotely-sensed products, e.g., Vegetation Index (VI) maps. These well-established means, however, lack important information, like crop height. Combinations of VI-maps and detailed 3D Crop Surface Models (CSMs) enable advanced methods for crop yield prediction. This work utilizes an Unmanned Aircraft System (UAS) to capture standard RGB imagery datasets for corn grain yield prediction at three early- to mid-season growth stages. The imagery is processed into simple VI-orthoimages for crop/non-crop classification and 3D CSMs for crop height determination at different spatial resolutions. Three linear regression models are tested on their prediction ability using site-specific (i) unclassified mean heights, (ii) crop-classified mean heights and (iii) a combination of crop-classified mean heights with according crop coverages. The models show determination coefficients R² of up to 0.74, whereas model (iii) performs best with imagery captured at the end of stem elongation and intermediate spatial resolution (0.04 m·px⁻¹). Following these results, combined spectral and spatial modeling, based on aerial images and CSMs, proves to be a suitable method for mid-season corn yield prediction.

1. Introduction

Corn (Zea mays L.) biomass and grain yields vary depending on site, climatic conditions and management decisions. Moreover, variation is likely to occur within fields. Following the idea of Precision Farming (PF), the identification of within-field spatial and temporal variability shows potential to support crop management concepts to meet much of the increasing environmental, economic, market and public pressures on arable agriculture [1]. Management strategies account for (i) environmental issues by adapting the input factors to the demand of the crop and, thus, avoid over- or under-application [2,3], (ii) economic issues by calculating within-field net returns [4] and (iii) possibilities to improve the control and influence of the quality of the product [5].

Yield estimations prior to harvest play a key role in the determination of input factors, like nutrients, pesticides and water, as well as for the planning of upcoming labor- and cost-intensive actions, like harvesting, drying and storage. In addition, bioenergy- and other corn-related industries benefit from these estimations, too [6]. Commonly, farmers use different methods for prediction. Coarse estimations are built on the farmer’s expert knowledge. Better estimations can be drawn from destructive sampling procedures in representative areas [7]. Unfortunately, destructive sampling is very labor-and cost-intensive work. Another approach is using yield maps, providing information about spatial and temporal variability of yields in previous years [8]. Although yield maps give some hints at within-field yield potential, they have limitations in explaining current growing conditions. Thus, reliable information about actual within-field yield estimations is usually drawn from more promising methods. Besides using linear regression models with additional information on crop management [6] or weather and soil attributes [9], several studies demonstrate the power of crop growth models to predict yield [10,11]. Although crop growth models return good estimates, their practical applicability may be limited due to the need of extensive input data for implementation. On a local and regional level, remote sensing products are quite common for estimating corn yield [12,13,14]. For a further increase in accuracy, some authors also combine actual remote sensing data and crop growth models [15,16]. Consequently, PF data has potential to improve crop development and yield prediction with smart management strategies and yield models.

With the upcoming of cheap and handy Unmanned Aircraft Systems (UASs), remotely-sensed data at high spatial and temporal resolutions have become more and more affordable [17]. Many researchers focus on RGB, multi-, hyper-spectral and thermal imaging techniques for crop monitoring [18,19,20], crop and weed discrimination [21,22] or on the generation of Digital Elevation Models (DEMs) [23,24,25,26]. Despite that, less research has been conducted on 3D Crop Surface Models (CSMs) [27,28,29] or on the possibilities of a combined analysis of both 3D and spectral information [30,31,32].

This study focuses on modeling of corn grain yield with a combined spectral and spatial analysis of aerial imagery. Standard imagery, which has been captured by a RGB consumer camera, the most common sensor used on UASs, serves as the data basis. Although RGB imagery carries limited spectral information compared to more sophisticated types, like multi-, hyper-spectral and thermal ones, its high spatial resolution allows one to create detailed CSMs for further crop investigation [28]. In addition to that, spectral information from RGB imagery can be used to determine positions of crops and estimate site-specific crop coverage factors by applying basic methods for crop/non-crop separation [33,34,35].

Recent studies found a high correlation of corn plant height and corn grain yield at early- to mid-season growth stages [36,37,38]. Yin et al. [37] also showed that linear regression models for the prediction of corn grain yield may be the preferred ones, because of their simplicity. Based on these findings, this study’s objective was to assess the potential of CSMs to predict corn grain yield at early- to mid-season growth stages by using mean crop heights and different linear regression models. The underlying hypotheses were to predict corn grain yield with simple linear regression models, building on plot-wise mean crop height as the predictor variable. The mean crop heights were generated in two ways, with and without respect to previously classified crop/non-crop pixels. Additionally, a multiple linear regression model was set up, including the crop coverage factor as a second predictor variable to improve prediction accuracy.

2. Materials and Methods

2.1. Experimental Setup

Ihinger Hof (48.74°N, 8.92°E), a research station of the University of Hohenheim, was chosen to serve as an experimental site for a field trial to predict corn grain yield by aerial imagery and crop surface models. The regional climate is categorized as a temperate climate with an annual average temperature of 7.9 °C and an average precipitation of 690 mm.

Figure 1. Overview of the two-factorial field trial in corn with 64 plots of a size of 36 × 6 m each. Four sowing densities (8–11 seeds·m⁻²) were tested at four different levels of nitrogen fertilization (50, 100, 150 and 200 kg·N·ha⁻¹) in a setup with four replicates.

Figure 1. Overview of the two-factorial field trial in corn with 64 plots of a size of 36 × 6 m each. Four sowing densities (8–11 seeds·m⁻²) were tested at four different levels of nitrogen fertilization (50, 100, 150 and 200 kg·N·ha⁻¹) in a setup with four replicates.

A two-factorial field trial was laid out in a common randomized split-plot design on 27 May 2013, with the corn cultivar “NK Ravello”. Four sowing densities (8–11 seeds·m⁻²) were tested at four different levels of nitrogen fertilization (50, 100, 150 and 200 kg·N·ha⁻¹) in a setup with four replicates. This resulted in 64 plots of a size of 36 × 6 m each and a total trial size of 1.38 ha (see Figure 1). Row spacing was set to 0.75 m, whereas seed spacing was adjusted according to the desired density level (0.115–0.158 m). Harvest and determination of corn grain yield with a moisture content of 14% took place on 28 October 2013, with a Global Navigation Satellite System (GNSS)-assisted combine harvester.

2.2. UAS and Sensor Setup

In this field experiment, a modified MikroKopter (MK) Hexa XL served as the aerial carrier platform to conduct sensor measurements [39]. Equipped with standard MK navigation sensors (Inertial Measurement Unit (IMU) and differential GNSS receiver), it is able to perform user-defined waypoint flights. Assembled with a payload of 1 kg and a lithium polymer battery with a capacity of 5000 mAh, this UAS operates approximately 10 min at an altitude level of 50 m above ground. With an additionally integrated Raspberry Pi Model B computer, it merges its navigation information with observations from attached sensor devices on-the-fly [40,41].

As the imaging sensor, a Canon Ixus 110 IS RGB consumer camera was attached to the UAS [42]. The camera’s sensor resolution was set to a maximum of 4000 × 3000 pixels to achieve a ground resolution of approximately 0.02 m·px⁻¹ at a flight altitude of 50 m. The camera was configured to predefined focal length (5.0 mm), aperture (f/2.8) and exposure time (1/500, 1/800 or 1/1000 s), whereas image triggering was software-controlled via a USB connection with the Raspberry Pi.

2.3. Measurements

Flight missions were performed on three dates during early- and mid-season crop development (beginning of stem elongation, end of stem elongation and end of emergence of inflorescence), referring to Zadoks’ scale’s Z32, Z39 and Z58 [43]. In each mission, aerial images were captured at a scheduled flight altitude of 50 m with an intended overlap of 80% in-track and 60% cross-track to ensure image redundancy. All images have been captured with a nadir view of direction, in clear skies and around noon. Each flight mission produced about 400 images covering all experimental plots with a ground resolution of approximately 0.02 m·px⁻¹. An overview of the flight missions is given in Table 1.

Table 1. Overview of performed flight missions at Zadoks’ scale’s crop growth stages Z32, Z39 and Z58 and the number of images for subsequent processing, flight altitude, approximate image ground resolution, mission time, illumination and wind speed.

Date	Growth Stage	Images	Scheduled Altitude (m)	Ground Resolution (m·px−1)	Time	Illumination	Wind (m·s−1)
17/07/2013	Z32	253	50	0.02	11–12 am	clear sky	1
01/08/2013	Z39	198	50	0.02	10–11 am	clear sky	2
15/08/2013	Z58	268	50	0.02	10–11 am	clear sky	2

Table 1. Overview of performed flight missions at Zadoks’ scale’s crop growth stages Z32, Z39 and Z58 and the number of images for subsequent processing, flight altitude, approximate image ground resolution, mission time, illumination and wind speed.

Date	Growth Stage	Images	Scheduled Altitude (m)	Ground Resolution (m·px−1)	Time	Illumination	Wind (m·s−1)
17/07/2013	Z32	253	50	0.02	11–12 am	clear sky	1
01/08/2013	Z39	198	50	0.02	10–11 am	clear sky	2
15/08/2013	Z58	268	50	0.02	10–11 am	clear sky	2

2.4. Image Processing

Prior to processing, the selected original images were reduced in resolution to create four additional datasets of imagery at ground resolutions of 0.04, 0.06, 0.08 and 0.10 m·px⁻¹. These artificial datasets were used to simulate corn grain yield prediction performance at different spatial resolution levels of aerial imagery. Regarding the shape and structure of corn, as well as the applied plant spacing of 0.115 to 0.158 m and a row spacing of 0.75 m, the computed ground resolutions lie somewhere within the leaf and canopy level. As a consequence, high ground resolutions are expected to cover fine structures (leaf level), whereas low resolutions are expected to cover coarse structures (canopy level). The following image processing routine was performed for each dataset and crop growth stage individually.

2.4.1. Orthoimage and Digital Elevation Model

Imagery and corresponding UAS navigation information were used to generate orthoimages and DEMs with the help of the 3D reconstruction software Agisoft PhotoScan 1.0.1 [44]. In a first step of processing, all selected images were aligned, mosaicked and geo-referenced by the software’s feature matching and Structure from Motion (SfM) algorithms. In a similar way as the popular Scale-Invariant Feature Transform (SIFT) approach from Lowe [45], feature detection was performed on each image to generate descriptors for image correspondence detection. Based on the correspondences and initial GNSS image locations, the SfM algorithm reconstructed the 3D scene, camera positions and orientations [46]. In a second step, a DEM was extracted from the 3D scene by applying a natural neighbor interpolation [47]. This DEM represents the geo-referenced surface of the experimental site and is based on altitude values relative to the GNSS’ reference ellipsoid. Generally, absolute crop heights are calculated by subtracting a second DEM, a so-called Digital Terrain Model (DTM), representing the surface of the ground relative to the same reference ellipsoid as the DEM (see Figure 2).

Figure 2. Visualization of DEM and DTM altitudes relative to a commonly shared GNSS reference ellipsoid (red surface). While the DEM represents a surface model of the experimental site (green surface), the DTM represents the surface of the ground. The DTM was approximated by interpolation of ground classified DEM pixels (yellow surface). Absolute crop heights are derived by subtraction of the two surface representations.

Figure 2. Visualization of DEM and DTM altitudes relative to a commonly shared GNSS reference ellipsoid (red surface). While the DEM represents a surface model of the experimental site (green surface), the DTM represents the surface of the ground. The DTM was approximated by interpolation of ground classified DEM pixels (yellow surface). Absolute crop heights are derived by subtraction of the two surface representations.

Therefore, in a third step, a DTM was inferred from the 3D scene by excluding non-ground pixels, which have been previously classified using the software’s automatic classification routine. To ensure the classification of real ground points, the point cloud was subdivided into cells of 7 × 7 m, and each cell’s lowest point was used for triangulation of a coarse initial DTM. After that, the initial DTM was densified by checking whether each remaining point meets the following two requirements: the vertical distance to the DTM-surface lies within a predefined buffer of 0.03 m, and at least one of the vectors to a ground-classified point intersects the DTM-plane with less than a predefined angle of 15°. In a last step, a mosaicked orthoimage, DEM and DTM were exported to three individual GeoTiff raster files for subsequent processing.

2.4.2. Crop Surface Model and Vegetation Indices

Further processing was performed with the statistical computation software, R [48,49,50]. The exported GeoTiff raster files were combined to a single raster stack object containing red, green, blue, DEM and DTM information as individual raster layers. A CSM raster layer was generated by pixel-wise subtraction of DTM layer altitudes from DEM layer altitudes and was added to the raster stack object.

In addition to that, three simple Vegetation Indices (VIs) were derived from the RGB bands containing the pixels’ greenness information in relation to their redness and/or blueness. The Excess Green Index (ExG), Vegetation Index Green (VIg), which is sometimes also referred to as the Normalized Green-Red Difference Index (NGRDI), and an adapted broadband variant of the Plant Pigment Ratio (PPRb) were selected as appropriate VIs to approach a detailed separation of crop and soil pixels [33,34,51,52]. Table 2 lists these VIs’ calculation formulas, which were performed on the raster stack object individually.

Table 2. Vegetation indices applied on the RGB images for pixel based crop/soil separation. The Excess Green Index (ExG) accounts for a combination of green and red, as well as green and blue reflection differences. The Vegetation Index Green (VIg) (sometimes also referred to as the Normalized Green-Red Difference Index (NGRDI)) represents a normalized green and red difference, whereas the adapted broadband variant of the Plant Pigment Ratio (PPRb) makes use of a normalized green and blue difference.

Index Reference	Explanation	Formula
ExG
Woebbecke et al. [33] & Meyer et al. [51]	Excess Green Index	2 × Rgreen − Rred − Rblue
VIg
Gitelson et al. [34]	Vegetation Index Green	$\frac{R_{green}-R_{red}}{R_{green}+R_{red}}$
PPRb
based on Metternicht [52]	Plant Pigment Ratio	$\frac{R_{green}-R_{blue}}{R_{green}+R_{blue}}$

Table 2. Vegetation indices applied on the RGB images for pixel based crop/soil separation. The Excess Green Index (ExG) accounts for a combination of green and red, as well as green and blue reflection differences. The Vegetation Index Green (VIg) (sometimes also referred to as the Normalized Green-Red Difference Index (NGRDI)) represents a normalized green and red difference, whereas the adapted broadband variant of the Plant Pigment Ratio (PPRb) makes use of a normalized green and blue difference.

Index Reference	Explanation	Formula
ExG
Woebbecke et al. [33] & Meyer et al. [51]	Excess Green Index	2 × Rgreen − Rred − Rblue
VIg
Gitelson et al. [34]	Vegetation Index Green	$\frac{R_{green}-R_{red}}{R_{green}+R_{red}}$
PPRb
based on Metternicht [52]	Plant Pigment Ratio	$\frac{R_{green}-R_{blue}}{R_{green}+R_{blue}}$

2.4.3. Plot-Wise Feature Extraction

Features were extracted by a self-developed automatic routine. First, field trial plot information was imported as a polygonal shapefile. For this analysis, plot size was reduced to rectangles of 9 × 6 m around the original plots’ centers to account for plot boundary effects, e.g., sowing or fertilization inaccuracies. Second, a shapefile containing harvested corn yield information was imported, and mean corn yields were determined for each individual plot. Third, mean plot heights were calculated using height information from the CSM layer. Fourth, for each VI layer, all pixels that fall inside a plot were extracted, and five different thresholds were computed on the selected pixels’ aggregated histogram based on the method of Ridler and Calvard [53] and Kort [50]. Consequently, VI layer pixels were classified as non-crop pixels in the case that the pixels’ values were below the defined thresholds and as crop pixels in the case that they were above the defined thresholds, respectively (see Figure 3).

Figure 3. VI-based Ridler thresholding by the example of a 4 × 4 m sub-sample of plot 413 with a sowing density of 11 seeds·m⁻² and nitrogen application of 50 kg·N·ha₋₁. The upper left corner shows the RGB orthoimage, which is displayed at a ground resolution of 0.04 m and at crop growth stage Z39. The second image in the upper row shows the ExG layer, which was derived from the RGB orthoimage. Based on the ExG layer’s histogram, five different thresholds were computed. Threshold r₃ is the original Ridler threshold, whereas the other thresholds represent four variations on the Ridler method (upper right corner). The remaining images show the ExG layer’s classification (green = crop, yellow = soil) based on the five thresholds. In this example, threshold r₃ and r₄ seem to classify best. Thresholds r₁ and r₂ seem to overestimate crop coverage, while r₅ seems to underestimate crop coverage.

Figure 3. VI-based Ridler thresholding by the example of a 4 × 4 m sub-sample of plot 413 with a sowing density of 11 seeds·m⁻² and nitrogen application of 50 kg·N·ha₋₁. The upper left corner shows the RGB orthoimage, which is displayed at a ground resolution of 0.04 m and at crop growth stage Z39. The second image in the upper row shows the ExG layer, which was derived from the RGB orthoimage. Based on the ExG layer’s histogram, five different thresholds were computed. Threshold r₃ is the original Ridler threshold, whereas the other thresholds represent four variations on the Ridler method (upper right corner). The remaining images show the ExG layer’s classification (green = crop, yellow = soil) based on the five thresholds. In this example, threshold r₃ and r₄ seem to classify best. Thresholds r₁ and r₂ seem to overestimate crop coverage, while r₅ seems to underestimate crop coverage.

Fifth, for each VI layer and its five identified thresholds, mean plot heights were calculated using the CSM layer height information solely from crop-classified pixels (see Figure 4). Sixth, for each VI layer and its five identified thresholds, plot crop coverage was computed by dividing the number of crop-classified pixels by the total number of pixels in each plot.

Figure 4. Mean crop height computation using the example of a 4 × 4 m sub-sample of plot 413 with a sowing density of 11 seeds·m⁻² and nitrogen application of 50 kg·N·ha⁻¹. The lower part of the figure shows a stack of the RGB orthoimage and the ExG layer classification based on threshold r₄ at a ground resolution of 0.04 m and at crop growth stage Z39. The upper part shows the corresponding CSM layer height information as a 3D representation, colored by the ExG-classification. Mean crop height was calculated by the crop-classified CSM layer heights only and is displayed as a semi-transparent plane.

Figure 4. Mean crop height computation using the example of a 4 × 4 m sub-sample of plot 413 with a sowing density of 11 seeds·m⁻² and nitrogen application of 50 kg·N·ha⁻¹. The lower part of the figure shows a stack of the RGB orthoimage and the ExG layer classification based on threshold r₄ at a ground resolution of 0.04 m and at crop growth stage Z39. The upper part shows the corresponding CSM layer height information as a 3D representation, colored by the ExG-classification. Mean crop height was calculated by the crop-classified CSM layer heights only and is displayed as a semi-transparent plane.

2.5. Modeling Strategy

In the last step of processing, the extracted features were used to model corn grain yield with three different strategies. Based on the findings of Yin et al. [37] that all investigated regression models predict sufficiently well, standard linear regression models were set up for prediction. Assuming that Y_iis the harvested corn grain yield, H_irs is the i-th mean plot height, regardless of any pixel classification, at the r-th ground resolution level and the s-th growth stage, whereas b₀ and b₁ are the regression coefficients. Equation (1) shows a simple linear regression model for corn grain yield prediction, forming strategy S₁.

$$Y_i=(b_0+b_1\times {\overline{H}}_{irs})+{\epsilon }_{irs}$$

(1)

Strategy S₂ was laid out in the same way as strategy S₁, except H_irstv representing the i-th mean plot height calculated from pixels, which were classified as crop by using the v-th VI layer and the t-th Ridler threshold estimate at the r-th ground resolution and s-th growth stage (Equation (2)).

$$Y_i=(b_0+b_1\times {\overline{H}}_{irstv})+{\epsilon }_{irstv}$$

(2)

The third strategy S₃ is a multiple linear regression approach, extending strategy S₂. This approach accounts for a second predictor variable C_irstvrepresenting the i-th plot crop coverage factor, which was computed by the v-th VI layer, and the t-th Ridler threshold estimate at the r-th ground resolution and the s-th growth stage (Equation (3)).

$$Y_i=(b_0+b_1\times {\overline{H}}_{irstv}+b_2\times C_{irstv})+{\epsilon }_{irstv}$$

(3)

While the first two strategies follow the approach of Yin et al. [37], strategy S₃ also considers the crop coverage factor as an additional predictor for expected corn grain yield.

2.6. Statistical Analysis

Statistical analysis was conducted with the statistical computation software, R. The field trial was analyzed as a mixed model using a standard two-way analysis of variance (ANOVA) approach. All modeling strategies for corn grain yield prediction were tested with and without classification-based mean crop heights at all crop growth stages, ground resolutions, Ridler threshold estimates and deduced crop coverage factors. The prediction accuracy of the different modeling strategies was assessed by using R² determination coefficient values as quality indicators. Spatial visualization of predicted and harvested corn grain yield was carried out using the geographical information system QGIS [54].

3. Results and Discussion

3.1. Field Trial

The ANOVA showed a significant influence of nitrogen fertilization on corn grain yield. Significant influences of sowing density, as well as of the interaction of both factors were not detected. The non-significant influence of sowing density was not expected, but might have been caused by the small variability in the range of sowing density levels of 8–11 seeds·m⁻². Detailed results are not presented in the following.

3.2. Image Processing

The 3D reconstruction software, Agisoft PhotoScan 1.0.1, was able to perform image alignment and 3D scene reconstruction for all imagery datasets. Geo-referencing was based on camera location information, derived from GNSS and IMU data. Orthoimage, DEM and DTM computation succeeded for all imagery datasets. Resulting orthoimage ground resolution was at the level of input image ground resolution. As dense point cloud reconstruction is a very hardware-demanding task, the imagery used for DEM and DTM generation was downscaled by a factor of two to save processing time. Although DEM and DTM were exported with the corresponding orthophoto’s ground resolution, the underlying dense point cloud was built with less detail than theoretically possible.

As the produced DTMs are based on the interpolation of previously classified ground points, this method is generally prone to misclassification at dense crop stands and canopy closure. In these situations, only a small amount of ground points will be visible at all, weakening the reliability of the interpolation results. Moreover, some of the classified points may not represent the “real” ground, leading to an underestimation of crop heights. In a homogeneous field, a correction factor could compensate for this underestimation. In an inhomogeneous field, the correction factor would not be constant anymore. To avoid these problems, it is recommended to produce DTMs at sowing stage, without the need for classification and interpolation of large gaps.

Geo-referencing accuracy was assessed by the help of 24 Ground Control Points (GCPs), which were installed permanently and measured with RTK-GNSS equipment. Heavy rainfalls in July silted many of the GCPs. In addition, others have been destroyed by intensive mechanical weed control in between the corn strips. Unfortunately, the GCPs were not renewed before performing flight missions at Z39 and Z58. As a consequence, imagery from these stages lack accurate GCP information. Thus, accuracy assessment was performed on Z32 imagery, only.

Table 3. Resulting root mean squared errors (m) (RMSE) at ground control point (GCP) locations for indirectly (GCP-based) and directly (GNSS- and IMU-based) geo-referenced imagery at Z32 for all image ground resolutions.

Geo-Reference	Coordinate Component	Ground Resolution (m·px−1)
		0.02	0.04	0.06	0.08	0.10
GCPs	Horizontal	0.058	0.063	0.084	0.089	0.082
	Vertical	0.068	0.059	0.051	0.046	0.075
GNSS & IMU	Horizontal	0.430	0.375	0.399	0.409	0.376
	Vertical	0.303	0.273	0.283	0.320	0.379

Table 3. Resulting root mean squared errors (m) (RMSE) at ground control point (GCP) locations for indirectly (GCP-based) and directly (GNSS- and IMU-based) geo-referenced imagery at Z32 for all image ground resolutions.

Geo-Reference	Coordinate Component	Ground Resolution (m·px−1)
		0.02	0.04	0.06	0.08	0.10
GCPs	Horizontal	0.058	0.063	0.084	0.089	0.082
	Vertical	0.068	0.059	0.051	0.046	0.075
GNSS & IMU	Horizontal	0.430	0.375	0.399	0.409	0.376
	Vertical	0.303	0.273	0.283	0.320	0.379

In addition to direct (GNSS- and IMU-based) geo-referencing, indirect (GCP-based) geo-referencing was conducted on Z32 imagery for enhanced CSM quality assessment. Table 3 lists the resulting root mean squared errors of a comparison of measured and computed GCP coordinates for both methods and all image ground resolutions at Z32. As expected, indirectly geo-referenced imagery showed smaller residuals than the directly geo-referenced one. Horizontal RMSEs for indirectly geo-referenced imagery ranged from 0.058 to 0.089 m, whereas vertical RMSEs ranged from 0.046 to 0.075 m. In contrast to that, horizontal RMSEs for directly geo-referenced imagery ranged from 0.375 to 0.430 m, whereas vertical RMSEs ranged from 0.273 to 0.379 m. The accuracies of both methods are in accordance with the findings of Turner et al. [55] and Ruiz et al. [56], although vertical accuracy performs slightly better than expected. GCP-based accuracy assessment for directly geo-referenced imagery at Z39 and Z58 was not performed. Nevertheless, comparison of identifiable field boundaries with those of Z32 did not show excessive horizontal accuracy errors for all resolutions.

The developed R-routine managed to calculate CSMs, VIs and all threshold variants for every imagery dataset. CSM quality was assessed by comparison of mean plot heights at Z32, derived from accurate and indirectly geo-referenced imagery, with those derived from less accurate and directly geo-referenced imagery. Table 4 shows the resulting root mean squared errors for plot height comparisons, ranging from 0.024 m for high resolution imagery to 0.008 m for low resolution imagery. With a difference of 0.20 m in between the highest and lowest mean plot height at Z32, direct geo-referencing shows little influence on mean plot height computation. Unfortunately, independent reference measurements, e.g., manual height measurements, 3D laser scanning datasets or CSMs, derived by other SfM software packages, were not available to assess absolute CSM accuracy. Therefore, subsequent analyses and results are proven for this dataset, only.

Table 4. Resulting root mean squared errors (m) (RMSE) of comparing mean plot heights calculated from indirectly (GCP-based) and directly (GNSS- and IMU-based) geo-referenced imagery at Z32 for all image ground resolutions.

Value	Coordinate Component	Ground Resolution (m·px−1)
		0.02	0.04	0.06	0.08	0.10
Plot Height	Vertical	0.024	0.010	0.009	0.010	0.008

Table 4. Resulting root mean squared errors (m) (RMSE) of comparing mean plot heights calculated from indirectly (GCP-based) and directly (GNSS- and IMU-based) geo-referenced imagery at Z32 for all image ground resolutions.

Value	Coordinate Component	Ground Resolution (m·px−1)
		0.02	0.04	0.06	0.08	0.10
Plot Height	Vertical	0.024	0.010	0.009	0.010	0.008

Horizontal alignment errors of directly geo-referenced imagery strongly influence the results of automatic feature extraction. To account for misalignment, the polygonal shapefile, containing this field trial’s plot information, was realigned individually for all imagery at all growth stages and image ground resolutions.

The computed original Ridler thresholds r₃ were regarded as suitable for automatic separation of crop and soil, as well as most of the threshold variants r₂ and r₄. In contrast to that, threshold variants r₁ and r₅ showed results of crop overestimation at threshold r₁ and underestimation at threshold r₅, respectively (see e.g., Figure 3). However, mean plot heights H_irstv and crop coverage factors C_irstv were computed for all strategies at every threshold level r₁₋₅ for subsequent comparison of prediction performance.

3.3. Modeling Strategy

All results of the applied corn grain yield prediction strategies are summarized in Table 5, whereas Figure 5 visualizes the most important findings. Strategy S₃ was evaluated for collinearity of its predictor variables, mean crop height and crop coverage. Critical collinearity at any crop growth stage was not found. As all strategies built on data from one growing period, leave-one-out cross-validation was conducted to evaluate each model’s predictive quality. Table 6 shows the resulting root mean squared errors of prediction (RMSEP), ranging from 0.67 to 1.28 t·ha⁻¹ (8.8% to 16.9%).

Crop growth stage Z32 was neglected in Figure 5, as none of the strategies resulted in R² determination coefficient values higher than 0.56. As crops were still small and stems were beginning to elongate, crops’ leaves were not overlapping at this point in time. Lacking canopy closure, the prediction models had to account for information contained at the leaf level. Therefore, imagery with highest resolutions of 0.02 and 0.04 m·px⁻¹ performed best and showed significant R² values. In contrast, lower resolution datasets did not provide much detail, resulting in low R² values. Strategy S₃ was generally able to significantly improve prediction accuracies of strategies S₁ and S₂ for all VIs by adding the crop coverage factor as the second predictor variable. Although the highest resolution imagery of 0.02 m·px⁻¹ performed best at this stage, even higher resolutions may be more appropriate for CSM and, thus, mean plot height generation. Reaching maximum R² values of 0.56 and considering additional environmental impacts on crop growth during the growing season, none of the applied strategies was assessed to be reliable for early-season corn grain yield prediction.

Figure 5. Resulting determination coefficients R² of modeling strategies S_1­3 for all VIs and aerial image ground resolutions at crop growth stages Z39 and Z58. Grey values represent R₂ values for strategy S₁, whereas black values represent strategy S₂ at Ridler threshold r₃ and colored values represent strategy S₃ at Ridler threshold r₃, respectively. In addition to the R² values of strategies S₂ and S₃ at Ridler threshold r₃, minimum and maximum R² values of the four remaining threshold variants are indicated as range bars for every aerial image ground resolution individually.

Figure 5. Resulting determination coefficients R² of modeling strategies S_1­3 for all VIs and aerial image ground resolutions at crop growth stages Z39 and Z58. Grey values represent R₂ values for strategy S₁, whereas black values represent strategy S₂ at Ridler threshold r₃ and colored values represent strategy S₃ at Ridler threshold r₃, respectively. In addition to the R² values of strategies S₂ and S₃ at Ridler threshold r₃, minimum and maximum R² values of the four remaining threshold variants are indicated as range bars for every aerial image ground resolution individually.

Table 5. The resulting determination coefficients R² of the prediction of corn grain yield by applying strategies S₁₋₃ for all combinations of VIs, aerial image ground resolutions, crop growth stages and computed Ridler thresholds. Significance codes for predictor variable crop height are represented as *in superscript, whereas significance codes for predictor variable crop coverage factor are represented as *in subscript (appearing only in strategy S₃).

Ground Res. (m·px−1)			ExG					VIg					PPRb
Z	Sx	rx	0.02	0.04	0.06	0.08	0.10	0.02	0.04	0.06	0.08	0.10	0.02	0.04	0.06	0.08	0.10
Z32	S1		0.48***	0.26***	0.12**	0.05	0.09*	0.48***	0.26***	0.12**	0.05	0.09*	0.48***	0.26***	0.12**	0.05	0.09*
Z32	S2	r1	0.54***	0.27***	0.12**	0.05	0.09*	0.53***	0.25***	0.11**	0.05	0.08*	0.46***	0.25***	0.11**	0.05	0.09*
Z32	S2	r2	0.55***	0.24***	0.11**	0.04	0.08*	0.55***	0.21***	0.09*	0.03	0.07*	0.47***	0.21***	0.08	0.03	0.08*
Z32	S2	r3	0.55***	0.23***	0.08*	0.03	0.08*	0.46***	0.16***	0.05	0.02	0.05	0.46***	0.18***	0.05	0.02	0.07*
Z32	S2	r4	0.53***	0.20***	0.06*	0.02	0.06*	0.36***	0.11**	0.02	0.01	0.04	0.42***	0.14**	0.03	0.01	0.04
Z32	S2	r5	0.51***	0.17***	0.04	0.01	0.05	0.25***	0.09*	0.01	0.01	0.05	0.37***	0.09*	0.02	0.00	0.03
Z32	S3	r1	0.54***	0.30*	0.21*	0.20**	0.21**	0.53***	0.25*	0.13	0.07	0.10	${0.52}_{**}^{***}$	${0.30}_{*}^{*}$	0.19*	0.13*	0.18*
Z32	S3	r2	0.55***	${0.36}_{**}^{*}$	0.32***	0.31***	0.33***	0.56***	${0.31}_{**}^{*}$	0.23**	0.19***	0.21**	${0.53}_{**}^{***}$	${0.35}_{***}^{*}$	0.27***	0.23***	0.27***
Z32	S3	r3	0.55***	${0.38}_{***}^{*}$	0.35***	0.34***	0.38***	${0.52}_{**}^{***}$	${0.37}_{***}^{*}$	0.32***	0.31***	0.34***	${0.52}_{**}^{***}$	${0.37}_{***}^{*}$	0.33***	0.30***	0.33***
Z32	S3	r4	0.55***	${0.37}_{***}^{*}$	0.33***	0.32***	0.38***	${0.47}_{***}^{***}$	${0.35}_{***}^{*}$	0.31***	0.32***	0.37***	${0.50}_{**}^{***}$	${0.36}_{***}^{*}$	0.32***	0.30***	0.33***
Z32	S3	r5	0.52***	${0.33}_{***}^{*}$	0.29***	0.29***	0.34***	${0.41}_{***}^{***}$	${0.31}_{***}^{*}$	0.28***	0.29***	0.36***	${0.46}_{**}^{***}$	${0.33}_{***}^{*}$	0.30***	0.29***	0.32***
Z39	S1		0.59***	0.68***	0.68***	0.63***	0.59***	0.59***	0.68***	0.68***	0.63***	0.59***	0.59***	0.68***	0.68***	0.63***	0.59***
Z39	S2	r1	0.60***	0.69***	0.68***	0.62***	0.58***	0.59***	0.68***	0.68***	0.62***	0.59***	0.57***	0.69***	0.68***	0.62***	0.58***
Z39	S2	r2	0.62***	0.70***	0.68***	0.62***	0.58***	0.61***	0.70***	0.68***	0.62***	0.58***	0.59***	0.70***	0.68***	0.62***	0.58***
Z39	S2	r3	0.63***	0.71***	0.68***	0.62***	0.58***	0.63***	0.70***	0.68***	0.62***	0.58***	0.60***	0.70***	0.68***	0.62***	0.58***
Z39	S2	r4	0.64***	0.71***	0.68***	0.62***	0.58***	0.63***	0.70***	0.67***	0.61***	0.57***	0.61***	0.70***	0.68***	0.62***	0.58***
Z39	S2	r5	0.64***	0.71***	0.68***	0.62***	0.58***	0.63***	0.70***	0.66***	0.60***	0.56***	0.62***	0.70***	0.68***	0.62***	0.58***
Z39	S3	r1	0.60***	0.69***	0.69***	0.62***	0.58***	0.59***	0.68***	0.68***	0.63***	0.60***	0.59***	0.69***	0.68***	0.63***	0.60***
Z39	S3	r2	${0.65}_{*}^{***}$	0.71***	0.70***	0.64***	0.60***	0.62***	0.70***	0.68***	0.63***	0.59***	0.59***	0.71***	0.69***	0.62***	0.58***
Z39	S3	r3	${0.68}_{**}^{***}$	0.72***	${0.72}_{**}^{***}$	${0.67}_{**}^{***}$	${0.63}_{**}^{***}$	0.63***	0.70***	0.68***	0.62***	0.59***	0.60***	0.71***	0.69***	0.63***	0.59***
Z39	S3	r4	${0.69}_{**}^{***}$	${0.73}_{*}^{***}$	${0.73}_{**}^{***}$	${0.69}_{***}^{***}$	${0.66}_{***}^{***}$	0.64***	0.71***	0.68***	0.62***	0.59***	0.61***	0.72***	0.70***	0.63***	0.59***
Z39	S3	r5	${0.70}_{***}^{***}$	${0.73}_{*}^{***}$	${0.74}_{***}^{***}$	${0.70}_{***}^{***}$	${0.68}_{***}^{***}$	0.64***	0.71***	0.68***	0.62***	0.58***	0.62***	0.71***	${0.70}_{*}^{***}$	0.64***	0.60***
Z58	S1		0.62***	0.68***	0.64***	0.64***	0.67***	0.62***	0.68***	0.64***	0.64***	0.67***	0.62***	0.68***	0.64***	0.64***	0.67***
Z58	S2	r1	0.59***	0.68***	0.65***	0.64***	0.67***	0.64***	0.68***	0.64***	0.64***	0.68***	0.59***	0.69***	0.65***	0.65***	0.68***
Z58	S2	r2	0.55***	0.69***	0.65***	0.65***	0.67***	0.64***	0.68***	0.64***	0.64***	0.67***	0.56***	0.69***	0.65***	0.65***	0.68***
Z58	S2	r3	0.52***	0.69***	0.65***	0.65***	0.68***	0.64***	0.68***	0.64***	0.64***	0.67***	0.53***	0.69***	0.65***	0.65***	0.68***
Z58	S2	r4	0.49***	0.69***	0.65***	0.65***	0.67***	0.63***	0.68***	0.64***	0.64***	0.67***	0.52***	0.69***	0.65***	0.65***	0.68***
Z58	S2	r5	0.46***	0.69***	0.65***	0.65***	0.66***	0.62***	0.68***	0.63***	0.63***	0.66***	0.52***	0.69***	0.65***	0.65***	0.68***
Z58	S3	r1	0.60***	0.69***	0.65***	0.64***	0.67***	0.64***	0.68***	0.65***	0.64***	0.68***	0.59***	0.69***	0.65***	0.65***	0.69***
Z58	S3	r2	0.55***	0.69***	0.65***	0.65***	0.68***	0.65***	0.68***	0.65***	0.64***	0.67***	0.56***	0.69***	0.65***	0.66***	0.69***
Z58	S3	r3	0.53***	0.69***	0.65***	0.65***	0.68***	0.65***	0.69***	0.65***	0.64***	0.67***	0.53***	0.69***	0.65***	0.65***	0.69***
Z58	S3	r4	${0.53}_{*}^{***}$	0.69***	0.65***	0.66***	0.69***	${0.66}_{*}^{***}$	0.69***	0.66***	0.65***	0.67***	0.52***	0.69***	0.65***	0.65***	0.68***
Z58	S3	r5	${0.54}_{**}^{***}$	0.69***	0.66***	0.67***	0.68***	${0.66}_{**}^{***}$	${0.70}_{*}^{***}$	${0.66}_{*}^{***}$	0.65***	0.67***	0.52***	0.69***	0.65***	0.65***	0.68***

Significance Codes R^{2H(sig.height)}R²_{C(sig.coverage)}^***: p< 0.001 ^**: p< 0.01 ^*: p< 0.05

Table 5. The resulting determination coefficients R² of the prediction of corn grain yield by applying strategies S₁₋₃ for all combinations of VIs, aerial image ground resolutions, crop growth stages and computed Ridler thresholds. Significance codes for predictor variable crop height are represented as *in superscript, whereas significance codes for predictor variable crop coverage factor are represented as *in subscript (appearing only in strategy S₃).

Ground Res. (m·px−1)			ExG					VIg					PPRb
Z	Sx	rx	0.02	0.04	0.06	0.08	0.10	0.02	0.04	0.06	0.08	0.10	0.02	0.04	0.06	0.08	0.10
Z32	S1		0.48***	0.26***	0.12**	0.05	0.09*	0.48***	0.26***	0.12**	0.05	0.09*	0.48***	0.26***	0.12**	0.05	0.09*
Z32	S2	r1	0.54***	0.27***	0.12**	0.05	0.09*	0.53***	0.25***	0.11**	0.05	0.08*	0.46***	0.25***	0.11**	0.05	0.09*
Z32	S2	r2	0.55***	0.24***	0.11**	0.04	0.08*	0.55***	0.21***	0.09*	0.03	0.07*	0.47***	0.21***	0.08	0.03	0.08*
Z32	S2	r3	0.55***	0.23***	0.08*	0.03	0.08*	0.46***	0.16***	0.05	0.02	0.05	0.46***	0.18***	0.05	0.02	0.07*
Z32	S2	r4	0.53***	0.20***	0.06*	0.02	0.06*	0.36***	0.11**	0.02	0.01	0.04	0.42***	0.14**	0.03	0.01	0.04
Z32	S2	r5	0.51***	0.17***	0.04	0.01	0.05	0.25***	0.09*	0.01	0.01	0.05	0.37***	0.09*	0.02	0.00	0.03
Z32	S3	r1	0.54***	0.30*	0.21*	0.20**	0.21**	0.53***	0.25*	0.13	0.07	0.10	${0.52}_{**}^{***}$	${0.30}_{*}^{*}$	0.19*	0.13*	0.18*
Z32	S3	r2	0.55***	${0.36}_{**}^{*}$	0.32***	0.31***	0.33***	0.56***	${0.31}_{**}^{*}$	0.23**	0.19***	0.21**	${0.53}_{**}^{***}$	${0.35}_{***}^{*}$	0.27***	0.23***	0.27***
Z32	S3	r3	0.55***	${0.38}_{***}^{*}$	0.35***	0.34***	0.38***	${0.52}_{**}^{***}$	${0.37}_{***}^{*}$	0.32***	0.31***	0.34***	${0.52}_{**}^{***}$	${0.37}_{***}^{*}$	0.33***	0.30***	0.33***
Z32	S3	r4	0.55***	${0.37}_{***}^{*}$	0.33***	0.32***	0.38***	${0.47}_{***}^{***}$	${0.35}_{***}^{*}$	0.31***	0.32***	0.37***	${0.50}_{**}^{***}$	${0.36}_{***}^{*}$	0.32***	0.30***	0.33***
Z32	S3	r5	0.52***	${0.33}_{***}^{*}$	0.29***	0.29***	0.34***	${0.41}_{***}^{***}$	${0.31}_{***}^{*}$	0.28***	0.29***	0.36***	${0.46}_{**}^{***}$	${0.33}_{***}^{*}$	0.30***	0.29***	0.32***
Z39	S1		0.59***	0.68***	0.68***	0.63***	0.59***	0.59***	0.68***	0.68***	0.63***	0.59***	0.59***	0.68***	0.68***	0.63***	0.59***
Z39	S2	r1	0.60***	0.69***	0.68***	0.62***	0.58***	0.59***	0.68***	0.68***	0.62***	0.59***	0.57***	0.69***	0.68***	0.62***	0.58***
Z39	S2	r2	0.62***	0.70***	0.68***	0.62***	0.58***	0.61***	0.70***	0.68***	0.62***	0.58***	0.59***	0.70***	0.68***	0.62***	0.58***
Z39	S2	r3	0.63***	0.71***	0.68***	0.62***	0.58***	0.63***	0.70***	0.68***	0.62***	0.58***	0.60***	0.70***	0.68***	0.62***	0.58***
Z39	S2	r4	0.64***	0.71***	0.68***	0.62***	0.58***	0.63***	0.70***	0.67***	0.61***	0.57***	0.61***	0.70***	0.68***	0.62***	0.58***
Z39	S2	r5	0.64***	0.71***	0.68***	0.62***	0.58***	0.63***	0.70***	0.66***	0.60***	0.56***	0.62***	0.70***	0.68***	0.62***	0.58***
Z39	S3	r1	0.60***	0.69***	0.69***	0.62***	0.58***	0.59***	0.68***	0.68***	0.63***	0.60***	0.59***	0.69***	0.68***	0.63***	0.60***
Z39	S3	r2	${0.65}_{*}^{***}$	0.71***	0.70***	0.64***	0.60***	0.62***	0.70***	0.68***	0.63***	0.59***	0.59***	0.71***	0.69***	0.62***	0.58***
Z39	S3	r3	${0.68}_{**}^{***}$	0.72***	${0.72}_{**}^{***}$	${0.67}_{**}^{***}$	${0.63}_{**}^{***}$	0.63***	0.70***	0.68***	0.62***	0.59***	0.60***	0.71***	0.69***	0.63***	0.59***
Z39	S3	r4	${0.69}_{**}^{***}$	${0.73}_{*}^{***}$	${0.73}_{**}^{***}$	${0.69}_{***}^{***}$	${0.66}_{***}^{***}$	0.64***	0.71***	0.68***	0.62***	0.59***	0.61***	0.72***	0.70***	0.63***	0.59***
Z39	S3	r5	${0.70}_{***}^{***}$	${0.73}_{*}^{***}$	${0.74}_{***}^{***}$	${0.70}_{***}^{***}$	${0.68}_{***}^{***}$	0.64***	0.71***	0.68***	0.62***	0.58***	0.62***	0.71***	${0.70}_{*}^{***}$	0.64***	0.60***
Z58	S1		0.62***	0.68***	0.64***	0.64***	0.67***	0.62***	0.68***	0.64***	0.64***	0.67***	0.62***	0.68***	0.64***	0.64***	0.67***
Z58	S2	r1	0.59***	0.68***	0.65***	0.64***	0.67***	0.64***	0.68***	0.64***	0.64***	0.68***	0.59***	0.69***	0.65***	0.65***	0.68***
Z58	S2	r2	0.55***	0.69***	0.65***	0.65***	0.67***	0.64***	0.68***	0.64***	0.64***	0.67***	0.56***	0.69***	0.65***	0.65***	0.68***
Z58	S2	r3	0.52***	0.69***	0.65***	0.65***	0.68***	0.64***	0.68***	0.64***	0.64***	0.67***	0.53***	0.69***	0.65***	0.65***	0.68***
Z58	S2	r4	0.49***	0.69***	0.65***	0.65***	0.67***	0.63***	0.68***	0.64***	0.64***	0.67***	0.52***	0.69***	0.65***	0.65***	0.68***
Z58	S2	r5	0.46***	0.69***	0.65***	0.65***	0.66***	0.62***	0.68***	0.63***	0.63***	0.66***	0.52***	0.69***	0.65***	0.65***	0.68***
Z58	S3	r1	0.60***	0.69***	0.65***	0.64***	0.67***	0.64***	0.68***	0.65***	0.64***	0.68***	0.59***	0.69***	0.65***	0.65***	0.69***
Z58	S3	r2	0.55***	0.69***	0.65***	0.65***	0.68***	0.65***	0.68***	0.65***	0.64***	0.67***	0.56***	0.69***	0.65***	0.66***	0.69***
Z58	S3	r3	0.53***	0.69***	0.65***	0.65***	0.68***	0.65***	0.69***	0.65***	0.64***	0.67***	0.53***	0.69***	0.65***	0.65***	0.69***
Z58	S3	r4	${0.53}_{*}^{***}$	0.69***	0.65***	0.66***	0.69***	${0.66}_{*}^{***}$	0.69***	0.66***	0.65***	0.67***	0.52***	0.69***	0.65***	0.65***	0.68***
Z58	S3	r5	${0.54}_{**}^{***}$	0.69***	0.66***	0.67***	0.68***	${0.66}_{**}^{***}$	${0.70}_{*}^{***}$	${0.66}_{*}^{***}$	0.65***	0.67***	0.52***	0.69***	0.65***	0.65***	0.68***

Significance Codes R^{2H(sig.height)}R²_{C(sig.coverage)}^***: p< 0.001 ^**: p< 0.01 ^*: p< 0.05

Z39 was identified as the crop growth stage with the best prediction performance. Figure 5 points out the most interesting findings. Generally, all VIs performed well, although the best results were achieved using ExG. High and intermediate ground resolutions of 0.04 and 0.06 m·px⁻¹ showed R² values of up to 0.74 for strategy S₃. However, strategy S₃ improved results for ExG only. VIg and PPRb did not show significant improvements. Strategy S₂ outperformed strategy S₁ for resolutions of 0.02 and 0.04 m·px⁻¹, whereas at intermediate and low ground resolutions, strategies S₁ and S₂ did not differ in prediction accuracy. Coarse VI layer information and the beginning of canopy closure seemed to level out differences of simple plot mean height computation and the classification-based one. Unexpectedly, the highest resolution of 0.02 m·px⁻¹ performed worse than high/intermediate resolutions. Although strategies S₂ and S₃ significantly improved prediction using ExG, highest resolution strategies appeared to be prone to higher noise and a scale effect, as the level of resolution leads to analysis in between leaf and canopy level. As a consequence, CSM and classification results may be biased.

Table 6. The resulting root mean squared errors of prediction (RMSEP) of the leave-one-out cross-validation for evaluation of the predictive quality of applying strategies S₁₋₃ for all combinations of VIs, aerial image ground resolutions, crop growth stages and computed Ridler thresholds.

Ground Res. (m·px−1)			ExG					VIg					PPRb
Z	Sx	rx	0.02	0.04	0.06	0.08	0.10	0.02	0.04	0.06	0.08	0.10	0.02	0.04	0.06	0.08	0.10
Z32	S1		0.93	1.11	1.20	1.25	1.21	0.93	1.11	1.20	1.25	1.21	0.93	1.11	1.20	1.25	1.21
Z32	S2	r1	0.88	1.11	1.21	1.25	1.21	0.89	1.13	1.21	1.25	1.22	0.94	1.12	1.21	1.25	1.21
Z32	S2	r2	0.86	1.13	1.21	1.25	1.22	0.87	1.15	1.22	1.26	1.23	0.93	1.14	1.22	1.25	1.22
Z32	S2	r3	0.87	1.14	1.23	1.26	1.22	0.95	1.17	1.24	1.26	1.24	0.94	1.16	1.24	1.26	1.23
Z32	S2	r4	0.88	1.15	1.24	1.26	1.23	1.03	1.20	1.26	1.27	1.24	0.98	1.19	1.26	1.27	1.24
Z32	S2	r5	0.90	1.17	1.25	1.27	1.24	1.11	1.22	1.27	1.27	1.24	1.02	1.21	1.27	1.28	1.26
Z32	S3	r1	0.90	1.10	1.16	1.17	1.15	0.90	1.14	1.23	1.26	1.23	0.91	1.10	1.19	1.23	1.19
Z32	S3	r2	0.88	1.04	1.07	1.08	1.06	0.87	1.08	1.14	1.16	1.15	0.91	1.07	1.13	1.16	1.13
Z32	S3	r3	0.88	1.02	1.05	1.05	1.02	0.90	1.03	1.07	1.08	1.05	0.91	1.04	1.08	1.10	1.07
Z32	S3	r4	0.88	1.03	1.06	1.07	1.03	0.95	1.05	1.08	1.07	1.03	0.93	1.04	1.08	1.09	1.07
Z32	S3	r5	0.90	1.07	1.10	1.09	1.06	1.00	1.08	1.10	1.09	1.04	0.96	1.07	1.09	1.10	1.07
Z39	S1		0.83	0.73	0.74	0.79	0.83	0.83	0.73	0.74	0.79	0.83	0.83	0.73	0.74	0.79	0.83
Z39	S2	r1	0.83	0.71	0.73	0.80	0.83	0.83	0.73	0.74	0.80	0.83	0.85	0.72	0.74	0.80	0.83
Z39	S2	r2	0.80	0.70	0.73	0.80	0.83	0.81	0.70	0.74	0.80	0.83	0.83	0.71	0.74	0.80	0.83
Z39	S2	r3	0.79	0.70	0.73	0.80	0.83	0.79	0.70	0.74	0.80	0.84	0.81	0.70	0.74	0.80	0.83
Z39	S2	r4	0.78	0.69	0.73	0.80	0.84	0.78	0.70	0.75	0.81	0.85	0.81	0.70	0.74	0.80	0.83
Z39	S2	r5	0.77	0.69	0.74	0.80	0.84	0.78	0.70	0.76	0.82	0.86	0.80	0.70	0.74	0.80	0.83
Z39	S3	r1	0.84	0.72	0.74	0.81	0.85	0.84	0.73	0.74	0.80	0.83	0.84	0.73	0.75	0.80	0.83
Z39	S3	r2	0.78	0.71	0.73	0.79	0.83	0.81	0.71	0.75	0.81	0.84	0.84	0.70	0.74	0.81	0.85
Z39	S3	r3	0.75	0.69	0.70	0.76	0.79	0.80	0.71	0.76	0.81	0.84	0.83	0.69	0.73	0.80	0.84
Z39	S3	r4	0.73	0.68	0.68	0.73	0.76	0.79	0.70	0.75	0.81	0.84	0.81	0.69	0.73	0.80	0.83
Z39	S3	r5	0.71	0.67	0.68	0.72	0.74	0.78	0.70	0.75	0.81	0.85	0.81	0.70	0.73	0.79	0.83
Z58	S1		0.82	0.72	0.77	0.77	0.74	0.82	0.72	0.77	0.77	0.74	0.82	0.72	0.77	0.77	0.74
Z58	S2	r1	0.85	0.72	0.76	0.76	0.73	0.79	0.72	0.77	0.77	0.73	0.85	0.71	0.76	0.76	0.73
Z58	S2	r2	0.89	0.71	0.76	0.76	0.73	0.79	0.72	0.77	0.77	0.73	0.89	0.71	0.76	0.76	0.73
Z58	S2	r3	0.93	0.71	0.76	0.76	0.73	0.80	0.72	0.77	0.77	0.74	0.91	0.71	0.76	0.76	0.73
Z58	S2	r4	0.96	0.71	0.76	0.76	0.73	0.80	0.72	0.78	0.77	0.74	0.93	0.71	0.76	0.76	0.73
Z58	S2	r5	0.98	0.71	0.76	0.76	0.74	0.82	0.72	0.78	0.77	0.75	0.93	0.72	0.76	0.76	0.73
Z58	S3	r1	0.85	0.72	0.77	0.77	0.74	0.80	0.73	0.78	0.78	0.74	0.86	0.72	0.77	0.76	0.72
Z58	S3	r2	0.90	0.72	0.77	0.77	0.74	0.80	0.73	0.78	0.78	0.74	0.90	0.72	0.77	0.76	0.72
Z58	S3	r3	0.94	0.72	0.77	0.76	0.73	0.79	0.72	0.77	0.78	0.74	0.92	0.72	0.77	0.76	0.73
Z58	S3	r4	0.93	0.72	0.77	0.75	0.73	0.78	0.71	0.76	0.77	0.74	0.93	0.72	0.77	0.76	0.73
Z58	S3	r5	0.92	0.72	0.76	0.75	0.73	0.78	0.71	0.76	0.76	0.74	0.93	0.72	0.77	0.77	0.73

Table 6. The resulting root mean squared errors of prediction (RMSEP) of the leave-one-out cross-validation for evaluation of the predictive quality of applying strategies S₁₋₃ for all combinations of VIs, aerial image ground resolutions, crop growth stages and computed Ridler thresholds.

Ground Res. (m·px−1)			ExG					VIg					PPRb
Z	Sx	rx	0.02	0.04	0.06	0.08	0.10	0.02	0.04	0.06	0.08	0.10	0.02	0.04	0.06	0.08	0.10
Z32	S1		0.93	1.11	1.20	1.25	1.21	0.93	1.11	1.20	1.25	1.21	0.93	1.11	1.20	1.25	1.21
Z32	S2	r1	0.88	1.11	1.21	1.25	1.21	0.89	1.13	1.21	1.25	1.22	0.94	1.12	1.21	1.25	1.21
Z32	S2	r2	0.86	1.13	1.21	1.25	1.22	0.87	1.15	1.22	1.26	1.23	0.93	1.14	1.22	1.25	1.22
Z32	S2	r3	0.87	1.14	1.23	1.26	1.22	0.95	1.17	1.24	1.26	1.24	0.94	1.16	1.24	1.26	1.23
Z32	S2	r4	0.88	1.15	1.24	1.26	1.23	1.03	1.20	1.26	1.27	1.24	0.98	1.19	1.26	1.27	1.24
Z32	S2	r5	0.90	1.17	1.25	1.27	1.24	1.11	1.22	1.27	1.27	1.24	1.02	1.21	1.27	1.28	1.26
Z32	S3	r1	0.90	1.10	1.16	1.17	1.15	0.90	1.14	1.23	1.26	1.23	0.91	1.10	1.19	1.23	1.19
Z32	S3	r2	0.88	1.04	1.07	1.08	1.06	0.87	1.08	1.14	1.16	1.15	0.91	1.07	1.13	1.16	1.13
Z32	S3	r3	0.88	1.02	1.05	1.05	1.02	0.90	1.03	1.07	1.08	1.05	0.91	1.04	1.08	1.10	1.07
Z32	S3	r4	0.88	1.03	1.06	1.07	1.03	0.95	1.05	1.08	1.07	1.03	0.93	1.04	1.08	1.09	1.07
Z32	S3	r5	0.90	1.07	1.10	1.09	1.06	1.00	1.08	1.10	1.09	1.04	0.96	1.07	1.09	1.10	1.07
Z39	S1		0.83	0.73	0.74	0.79	0.83	0.83	0.73	0.74	0.79	0.83	0.83	0.73	0.74	0.79	0.83
Z39	S2	r1	0.83	0.71	0.73	0.80	0.83	0.83	0.73	0.74	0.80	0.83	0.85	0.72	0.74	0.80	0.83
Z39	S2	r2	0.80	0.70	0.73	0.80	0.83	0.81	0.70	0.74	0.80	0.83	0.83	0.71	0.74	0.80	0.83
Z39	S2	r3	0.79	0.70	0.73	0.80	0.83	0.79	0.70	0.74	0.80	0.84	0.81	0.70	0.74	0.80	0.83
Z39	S2	r4	0.78	0.69	0.73	0.80	0.84	0.78	0.70	0.75	0.81	0.85	0.81	0.70	0.74	0.80	0.83
Z39	S2	r5	0.77	0.69	0.74	0.80	0.84	0.78	0.70	0.76	0.82	0.86	0.80	0.70	0.74	0.80	0.83
Z39	S3	r1	0.84	0.72	0.74	0.81	0.85	0.84	0.73	0.74	0.80	0.83	0.84	0.73	0.75	0.80	0.83
Z39	S3	r2	0.78	0.71	0.73	0.79	0.83	0.81	0.71	0.75	0.81	0.84	0.84	0.70	0.74	0.81	0.85
Z39	S3	r3	0.75	0.69	0.70	0.76	0.79	0.80	0.71	0.76	0.81	0.84	0.83	0.69	0.73	0.80	0.84
Z39	S3	r4	0.73	0.68	0.68	0.73	0.76	0.79	0.70	0.75	0.81	0.84	0.81	0.69	0.73	0.80	0.83
Z39	S3	r5	0.71	0.67	0.68	0.72	0.74	0.78	0.70	0.75	0.81	0.85	0.81	0.70	0.73	0.79	0.83
Z58	S1		0.82	0.72	0.77	0.77	0.74	0.82	0.72	0.77	0.77	0.74	0.82	0.72	0.77	0.77	0.74
Z58	S2	r1	0.85	0.72	0.76	0.76	0.73	0.79	0.72	0.77	0.77	0.73	0.85	0.71	0.76	0.76	0.73
Z58	S2	r2	0.89	0.71	0.76	0.76	0.73	0.79	0.72	0.77	0.77	0.73	0.89	0.71	0.76	0.76	0.73
Z58	S2	r3	0.93	0.71	0.76	0.76	0.73	0.80	0.72	0.77	0.77	0.74	0.91	0.71	0.76	0.76	0.73
Z58	S2	r4	0.96	0.71	0.76	0.76	0.73	0.80	0.72	0.78	0.77	0.74	0.93	0.71	0.76	0.76	0.73
Z58	S2	r5	0.98	0.71	0.76	0.76	0.74	0.82	0.72	0.78	0.77	0.75	0.93	0.72	0.76	0.76	0.73
Z58	S3	r1	0.85	0.72	0.77	0.77	0.74	0.80	0.73	0.78	0.78	0.74	0.86	0.72	0.77	0.76	0.72
Z58	S3	r2	0.90	0.72	0.77	0.77	0.74	0.80	0.73	0.78	0.78	0.74	0.90	0.72	0.77	0.76	0.72
Z58	S3	r3	0.94	0.72	0.77	0.76	0.73	0.79	0.72	0.77	0.78	0.74	0.92	0.72	0.77	0.76	0.73
Z58	S3	r4	0.93	0.72	0.77	0.75	0.73	0.78	0.71	0.76	0.77	0.74	0.93	0.72	0.77	0.76	0.73
Z58	S3	r5	0.92	0.72	0.76	0.75	0.73	0.78	0.71	0.76	0.76	0.74	0.93	0.72	0.77	0.77	0.73

At Z58, results were strongly influenced by the occurrence of canopy closure. Hence, neither strategy S₂ nor strategy S₃ were able to significantly improve the corn grain yield prediction performance of strategy S₁. Moreover, highest resolution strategies showed similar patterns as in Z39. Except using VIg, imagery at a ground resolution of 0.02 m·px⁻¹ seemed to underlay CSM and misclassification as in Z39. All other resolutions performed comparatively well, independent of applied strategy and VI. Although, these resolutions did not reach the maximum R² values of Z39, they were still considered as suitable for prediction.

Figure 6. Spatial illustration of plot-wise distribution of harvested corn grain yield (top), corn grain yield predicted by strategy S₃ at crop growth stage Z39, with ExG at Ridler threshold r₄ and an aerial image ground resolution of 0.04 m·px⁻¹ (middle) and the resulting prediction error of this strategy (bottom). For this strategy, the total root mean squared error of prediction (RMSEP) equals 0.68 t·ha⁻¹ (8.8%).

Figure 6. Spatial illustration of plot-wise distribution of harvested corn grain yield (top), corn grain yield predicted by strategy S₃ at crop growth stage Z39, with ExG at Ridler threshold r₄ and an aerial image ground resolution of 0.04 m·px⁻¹ (middle) and the resulting prediction error of this strategy (bottom). For this strategy, the total root mean squared error of prediction (RMSEP) equals 0.68 t·ha⁻¹ (8.8%).

Table 7 summarizes the key findings. The most suitable resolution and modeling strategy depends on the crop growth stage. Due to row-based cultivation of corn and missing canopy closure, early growth stages require very high resolution imagery for accurate CSM computation and classification-based separation of crop and soil. Therefore, strategies S₂ and S₃ result in higher R₂ values than strategy S₁ (R₂ ≤ 0.56). With ongoing crop development and beginning canopy closure, high resolution imagery and crop/soil classification gets less and less important. Highest resolution imagery showed a significant reduction of prediction accuracy at mid-season growth stages. All other imagery resolutions performed almost equally well (approximately 0.60 ≤ R₂ ≤ 0.70) at all strategies S_1–3 within these stages. Best prediction results were achieved by applying strategy S₂ and especially strategy S₃ at Z39 (R² ≤ 0.74). Although strategy S₃ proved to have good performance at this specific growth stage, further investigation of the influence of crop coverage factor C_irstv on the prediction results of this multiple linear regression strategy seems of great interest.

Table 7. Overview of the best performing parameters for early- to mid-season corn grain yield prediction at different crop growth stages. So far, the increase in prediction performance in strategy S₃ appears to underlay an unknown factor. Therefore, strategy S₃ is listed in brackets.

	Growth Stage
	Z32	Z39	Z58
Ground Resolution	highest/high	high/intermediate	high/intermediate/low
Vegetation Index	ExG	ExG	VIg
Prediction Strategy	S2 / (S3)	S2 / (S3)	S1 / S2 / (S3)

Table 7. Overview of the best performing parameters for early- to mid-season corn grain yield prediction at different crop growth stages. So far, the increase in prediction performance in strategy S₃ appears to underlay an unknown factor. Therefore, strategy S₃ is listed in brackets.

	Growth Stage
	Z32	Z39	Z58
Ground Resolution	highest/high	high/intermediate	high/intermediate/low
Vegetation Index	ExG	ExG	VIg
Prediction Strategy	S2 / (S3)	S2 / (S3)	S1 / S2 / (S3)

These findings indicate the best corn grain yield prediction at mid-season crop growth stages Z39 and Z58. They are in accordance with the findings of Yin et al. [37]. Nevertheless, none of the strategies showed results comparable to the best predictions of Yin et al. [37]. Depending on the growth stage and crop rotation system, Yin et al. [37] stated significant determination coefficients of 0.25 ≤ R² ≤ 0.89, whereas low R² values were achieved at early-season growth stages, only.

Applying strategy S₃ at Z39, Figure 6 visualizes plot-wise prediction results and compares them to the harvested corn grain yield. Using ExG at Ridler threshold r₄ and an aerial image ground resolution of 0.04 m·px⁻¹, the total RMSEP equals 0.68 t·ha⁻¹ (8.8%). Although this strategy performed best, the ANOVA of the field trial’s input factors did not show significant influence of sowing density on corn grain yield. As strategy S₃ utilizes computed crop coverage C_irstvas the estimator for sowing/stand density, the increase in prediction performance seems to underlay another factor, correlated with C_irstv. Other combinations of strategy S₃ and VIg/PPRb did not show improved results compared to strategy S₂.

4. Conclusions

This work shows the potential of exploiting spectral and spatial information from UAS-based RGB imagery for predicting corn grain yield in early- to mid-season crop growth stages. RGB imagery was used to compute crop surface models and to extract crop height information. In combination with RGB-based VI information, three different linear regression models were tested for the prediction of corn grain yield with R² determination coefficients of up to 0.74 and RMSEP ranging from 0.67 to 1.28 t·ha⁻¹ (8.8% to 16.9%).

Generally, all tested VIs performed almost equally well at any crop growth stage. The same applies to tested classification thresholds r₂₋₄. Although some of the more extreme thresholds r₁ and r₅ showed satisfying results, these thresholds cannot be recommended, because of potential over- or under-estimation of crop coverage.

The most suitable resolution and modeling strategy depends on the crop growth stage. Due to row-based cultivation of corn and missing canopy closure, early growth stages require very high resolution imagery for accurate CSM computation and classification-based separation of crop and soil. Compared to using simple unclassified mean crop heights (S₁), prediction results significantly improve, when accounting for additional crop/soil classification information (S₂ and S₃). With ongoing crop development and beginning canopy closure, high resolution imagery gets less and less important, sometimes even disadvantageous, due to higher noise. Good prediction results are achieved at intermediate resolutions by considering crop coverage as the second predictor variable (S₃). With the completion of canopy closure, neither high resolution imagery nor crop/soil classification show potential to further improve prediction. Concluding these findings, combined spectral and spatial modeling, based on aerial images and CSMs, proves to be a suitable method for mid-season corn yield prediction.