Multispectral Vegetation Indices and Machine Learning Approaches for Durum Wheat (Triticum durum Desf.) Yield Prediction across Different Varieties

Abstract

Durum wheat (Triticum durum Desf.) is one of the most widely cultivated cereal species in the Mediterranean basin, supporting pasta, bread and other typical food productions. Considering its importance for the nutrition of a large population and production of high economic value, its supply is of strategic significance. Therefore, an early and accurate crop yield estimation may be fundamental to planning the purchase, storage, and sale of this commodity on a large scale. Multispectral (MS) remote sensing (RS) of crops using unpiloted aerial vehicles (UAVs) is a powerful tool to assess crop status and productivity with a high spatial–temporal resolution and accuracy level. The object of this study was to monitor the behaviour of thirty different durum wheat varieties commonly cultivated in Italy, taking into account their spectral response to different vegetation indices (VIs) and assessing the reliability of this information to estimate their yields by Pearson’s correlation and different machine learning (ML) approaches. VIs allowed us to separate the tested wheat varieties into different groups, especially when surveyed in April. Pearson’s correlations between VIs and grain yield were good (R² > 0.7) for a third of the varieties tested; the VIs that best correlated with grain yield were CVI, GNDVI, MTVI, MTVI2, NDRE, and SR RE. Implementing ML approaches with VIs data highlighted higher performance than Pearson’s correlations, with the best results observed by random forest (RF) and support vector machine (SVM) models.

1. Introduction

Agriculture in the third millennium faces a number of different challenges, ranging from increasing food production to the provision of new organic products to replace petroleum derivatives, in a context dominated by the effects of global warming and climate change, with variations in seasonal weather patterns and an increase in extreme catastrophic natural events [1,2,3,4]. Concerning these phenomena, agriculture must adapt with appropriate adaptation and mitigation strategies, hence, resilience and sustainability (i.e., carbon sequestration, GHG emission reductions, new soil and crop management techniques). With this regard, the EU Commission in 2020, within the European Green Deal, launched the farm-to-fork strategy to promote a “food systems fair, healthy and environmentally-friendly”, setting a series of ambitious targets for EU agriculture by 2030 [5].

Among arable land crops, durum wheat (Triticum durum Desf.) represents an important food crop in Southern Europe and Mediterranean basin countries with a production of ~20 Mt, representing the main production area of this species [6,7]. Durum wheat and the product derived from its milling, called semolina, is the basic raw material for preparing typical foods of the area, such as pasta, traditional bread, couscous and burghul. Among these products, semolina and pasta production, as well as raw grain, in Italy also play an important role in commercial exports [8,9]. In recent years, the world durum wheat market has shown considerable volatility compared to the past due to a series of conjunctural events (e.g., the war in Ukraine, Canada and USA drought; cloudbursts in late spring in Italy), the sudden increase in the international demand and a critical reduction in yields by leading producing nations. This circumstance brought out potential risks related to the supply of mills and pasta factories, promoting procurement strategies such as supply chain conventions with local farmers or wheat future contracts [10]. Even greater are the risks for the supply of ancient wheat varieties characterised by an important role in biodiversity and promising nutritional aspects [11,12]. Therefore, accurately estimating the yield in a given area is essential to adequately plan industry storage and make the appropriate purchases in time.

Remote sensing (RS) techniques, within the precision agriculture (PA) approach, can provide several solutions to collect data and estimate production [13,14]. Unpiloted aerial vehicles (UAVs) equipped with multispectral (MS) sensors could offer several benefits, such as greater efficiency and accuracy through increased spatial and temporal resolution, providing images at very high resolution (VHR) on a daily basis [15,16,17]. Combining the different MS bands, the obtained vegetation indices (VIs) allow monitoring, analysis and mapping of vegetation temporal and spatial variations [18,19,20]. Numerous studies have highlighted the potential of VIs to predict plant biomass and crop yields on several crops like wheat [21,22], onion [18,23], oat [24,25], barley [26,27,28], tomato [29,30,31] etc. However, the simple linear application of VIs does not always allow reliable yield estimation, especially when spatial and spectral differences are slight [32,33]. A more advanced approach is the application of machine learning (ML) techniques to field and RS data, which allows for better interpretation of patterns and more robust estimations [34,35]. Different ML regression models have been successfully applied to agricultural data like linear, polynomial and logistical regressions, random forest (RF), support vector machines (SVM), neural networks (NN), k-nearest neighbours (k-NN), and stochastic gradient boosting [34]. Few studies have yet been conducted concerning durum wheat, especially in Italy, and numerous aspects like cultivar response to VIs, the best time for UAV monitoring, and the performance of ML approaches need to be fully investigated. Therefore, the objectives of the present study were: (i) evaluate the different spectral responses of the most common durum wheat cultivars grown in Italy; (ii) analyse the correlations between the observed VIs and grain yield for each cultivar; (iii) evaluate the potentials of ML techniques for durum wheat grain yield prediction; (iv) identify for the three previous objectives the best crop survey time between the late tillering (April) and the earing phase (May). The main novelty aspects of the proposed research mainly concern the simultaneous investigation of numerous varieties (30 different durum wheat varieties investigated), the number of indices considered (13 VIs) and their combined use to feed ML models.

2. Materials and Methods

2.1. Experimental Site

The field experiment was set up at the agricultural experimental centre “Casello” of the regional agency ARSAC (Azienda Regionale per lo Sviluppo dell’Agricoltura Calabrese) located in San Marco Argentano, Calabria, Italy (39°38′ N, 16°13′ E, 100 m a.s.l.) (Figure 1). The soil of the experimental site is classified as Fluventic Haploxerepts [36] with a sandy-clay-loam texture, neutral pH (7.7), total N 0.95 g kg^–1 and organic matter 16.9 g kg^–1. The climate is Mediterranean with mild and rainy (70% of total) falls and winters and warm and dry summers (Warm Mediterranean Climate, Csa); mean annual rainfall is 706 mm while the mean yearly air temperature is 14.7 °C (30-year averages). Further information regarding the experimental site characteristics is available in Badagliacca et al. [37].

2.2. Experimental Design and Crop Management

The field experiment was set up as a completely randomised block design (RCBD) with three replications during the 2021–2022 cropping season. Thirty different wheat cultivars were tested: Antalis, Aureo, Beltorax, Bering, Brigante, Claudio, Diogene, Federico II, Fuego, Furio Camillo, Incanto, Iride, LG Fructis, Maciste, Mameli, Marakas, Marco Aurelio, Monastir, Nuraghe, Panoramix, RGT Aventadur, RGT Jasdur, RGT Natur, RGT Voltadur, SY Leonardo, SY Prodigio, Tancredi, Telemaco, Tito Flavio, and Verace (Figure 1). The soil was prepared by mouldboard ploughing to a depth of 30 cm in October, followed by two shallow harrowing operations at 15 cm soil depth in November. Sowing was performed with a plot seeder (Vignoli Tartaro) in the first half of December, in plots of 10 m² (1.44 × 7.00 m) in rows (No. 8) 18 cm apart at 350 viable seeds m^–2 density. The field was fertilised at the sowing (BBCH-scale phase 00) and tillering (BBCH-scale phase 23) stages by broadcasting 200 kg ha^–1 of di-ammonium phosphate (DAP, 18-46-0) at the first stage and 130 kg ha^–1 of urea (46-0-0), at the second stage. Durum wheat grain was harvested at full maturity stage in mid-July 2022, using a plot combine (Wintersteiger Nursery Master).

2.3. UAV Surveys and Image Processing

The adopted workflow is synthesised in Figure 2.

The multispectral images were acquired in two different wheat stages: late tillering (25 BBCH-scale phase) and heading (55 BBCH-scale phase) on 4 April 2022 and 5 May 2022, respectively. Surveys were performed using the UAV multirotor DJI Phantom 4 Multispectral (DJI Ltd., Shenzhen, China) equipped with a camera generating 1600 × 1300 pixels images at 2 MP resolution. The equipped camera shoots simultaneously six images, one in RGB composition and five monochrome spectral images (i.e., Blue, Green, Red, Red Edge, Near Infrared (NIR)—

Table 1. Spectral characteristics of the multispectral sensors investigated bands.

Band	Central Wavelength [nm]	Bandwidth [nm]
Blue	450	±16
Green	560
Red	650
Red Edge	730
Near Infrared (NIR)	840	±26

).

All flights were performed in a cloud-free condition at the constant speed of 2 ms⁻¹ and 30 m a.g.l. altitude. An automated surveying mission was set using the DJI Ground Station Pro application, ensuring 80% overlap and sidelap. The same mission was used for both flights. Moreover, a set of black and white square (50 cm × 50 cm) ground control points (GCPs) were placed on the field to achieve a better geolocation of the output. Their coordinates were acquired by means of the RTK (real-time kinematic) GNSS (global navigation satellite system) Leica GS12 with a planimetric accuracy of ±2.5 cm and an altimetric one of ±5 cm. The photogrammetric process, consisting of image alignment, stack and radio-metric correction, was performed for both surveys (i.e., April and May 2022) using the Pix4D mapper Pro v.4.3 (Pix4D SA, Lausanne, Switzerland) software. At the end of this process, a multispectral orthomosaic with a geometric resolution of 1 cm has been obtained as output. In order to calibrate the orthomosaic, converting digital number values in reflectance, we used a field spectroradiometer to acquire, in correspondence with the surveyed spectral bands, reflectance values of a grey calibration panel put on the field during the flights. More details about the photogrammetric process can be found in other research previously published by this research group [18,19,38]. To test the correlation between different VIs and crop production data for each analysed variety, a mask was created to eliminate all parts of the images not covered by the crop, limiting the analysis to the designed plots (Figure 1 and Figure 2). A set of VIs (

Table 2. The analysed vegetation indices (Vis) of this study.

Vegetation Index	Acronym	Formula	Reference
Chlorophyll Vegetation Index	CVI	$\frac{\rho NIR}{\rho Green}\ast \frac{\rho Red}{\rho Green}$	[39]
Green Normalised Difference Red Edge Index	GNDRE	$\frac{\rho Red Edge-\rho Green}{\rho Red Edge+\rho Green}$	[40]
Green Normalised Difference Vegetation Index	GNDVI	$\frac{\rho NIR-\rho Green}{\rho NIR+\rho Green}$	[41]
Modified Chlorophyll Absorption Ratio Index	MCARI2	$\frac{1.5\ast [2.5 \left(\rho NIR-\rho Red\right)-1.3 \left(\rho NIR-\rho Green\right)]}{\sqrt{{\left(2\rho NIR+1\right)}^2-\left(6\rho NIR-5\rho Red\right)-0.5}}$	[42]
Modified Triangular Vegetation Index	MTVI	$1.2\ast [1.2 (\rho NIR-\rho Green)-2.5 (\rho Red-\rho Green\left)\right]$
Modified Triangular Vegetation Index 2	MTVI2	$\frac{1.5\ast [1.2 \left(\rho NIR-\rho Green\right)-2.5 \left(\rho Red-\rho Green\right)]}{\sqrt{{\left(2\rho NIR+1\right)}^2-\left(6\rho NIR-5\rho Red\right)-0.5}}$
Normalised Difference Red Edge Index	NDRE	$\frac{\rho NIR-\rho Red Edge}{\rho NIR+\rho Red Edge}$	[43]
Normalised Difference Vegetation Index	NDVI	$\frac{\rho NIR-\rho Red}{\rho NIR+\rho Red}$	[44]
Optimised Soil Adjusted Vegetation Index	OSAVI	$\frac{1.16 (\rho NIR-\rho Red)}{\rho NIR+\rho Red+0.16}$	[45]
Renormalised Difference Vegetation Index	RDVI	$\frac{\rho NIR-\rho Red}{\sqrt{\rho NIR+\rho Red}}$	[46]
Red Edge Triangulated Vegetation Index	RTVI	$100 (\rho NIR-\rho Red Edge)-10 (\rho NIR-\rho Green\left)\right]$	[47]
Simple Ratio	SR	$\frac{\rho NIR}{\rho Red}$	[48]
Simple Ratio Red Edge	SRRE	$\frac{\rho NIR}{\rho Red Edge}$	[49]

) was tested, and for each one, the mean value for each plot was considered.

2.4. Statistics and Machine Learning (ML) Approaches

Field and VIs data were collected and organised in MS Excel^(TM). VIs data and related elaboration were performed separately for the two UAV flight times (April and May) to assess the most useful flying epoch for cultivar separation and grain yield prediction. Statistical and ML analyses were carried out in the RStudio environment. Differences in VIs responses among wheat cultivars were evaluated by one-way analysis of variance (ANOVA) followed by means comparisons using Tukey’s HSD test at the 5% probability level (p-value < 0.05). In addition, a PCA was carried out from data of all calculated VIs to highlight distances in spectral response among the varieties tested. For these two analyses were used “agricolae” [50] and “FactoMineR v1.41” [51] packages. Relationships between wheat yield and VIs were investigated by Pearson’s correlation and different ML approaches. Pearson’s correlation was performed by using “cor()” function [52]. ML approaches were implanted and evaluated by using “Caret” package [53]. Five different ML models, commonly used in agricultural studies, were tested: (1) linear model (LM), (2) random forest (RF), (3) support vector machines (SVM), (4) k-nearest neighbours (k-NN), and (5) neural networks (NN). ML approaches were tested only on the 10 cultivars that showed significant Pearson’s correlation between grain yields and VIs, namely Antalis, Aureo, Beltorax, Bering, Brigante, Federico II, LG Fructis, Marco Aurelio, Panoramix, and RGT Aventadur. For each model, experimental data (i.e., yield and VIs values) were divided into training (70%) and validation (30%) datasets to train and validate their performance. This procedure was repeated 100 times (bootstrap sampling) in order to achieve a comprehensive and reliable model evaluation that can account for the widest variability in the data used for training and validation. The model’s performances were evaluated, comparing the calculated dataset with the validation dataset, using the coefficient of determination (R²), root mean square error (RMSE), and mean absolute error (MAE) observed in 50% of the repetitions performed (50th percentile).

3. Results and Discussion

3.1. Indices (VIs) Responses of Durum Wheat Cultivars

The data and the results of the ANOVA analysis performed on the VIs surveyed in April on the 30 wheat varieties studied are presented below in

Table 3. Vegetation indices (VIs) data observed among the 30 wheat cultivars investigated in April, related statistics (minimum, maximum, difference and mean values) and results of the ANOVA analysis. Different letters indicate significant differences between cultivars (Tukey’s HSD test at p < 0.05). For each VI in bold, the highest value is highlighted, while in italics, the lowest value.

		Vegetation Indices—April
		CVI	GNDRE	GNDVI	MCARI2	MTVI	MTVI2	NDRE	NDVI	OSAVI	RDVI	RTVI	SR	SR RE
Cultivars	Antalis	2.04 c	0.47 b	0.69 b	0.38 b	0.19 c	0.18 b	0.33 c	0.88 a	0.46 c	0.32 b	5.34 a	16.46 d	2.02 e
	Aureo	1.94 e	0.46 c	0.68 b	0.35 c	0.21 a	0.20 a	0.33 c	0.87 a	0.49 a	0.34 a	5.85 a	15.74 d	2.01 e
	Beltorax	1.84 g	0.45 d	0.66 c	0.43 b	0.19 d	0.17 b	0.30 d	0.85 b	0.44 c	0.31 b	4.74 b	14.13 f	1.89 g
	Bering	1.87 f	0.45 c	0.66 c	0.43 b	0.19 d	0.18 b	0.31 d	0.86 b	0.45 c	0.31 b	4.78 b	14.25 f	1.90 g
	Brigante	1.96 e	0.47 b	0.69 b	0.36 c	0.20 b	0.19 a	0.33 c	0.88 a	0.48 c	0.33 a	5.56 a	16.43 d	2.00 e
	Claudio	1.87 f	0.46 c	0.67 b	0.40 b	0.19 c	0.18 b	0.31 c	0.86 a	0.46 c	0.32 b	5.00 b	15.01 e	1.92 g
	Diogene	1.86 f	0.46 b	0.67 b	0.41 b	0.19 c	0.18 b	0.30 d	0.86 b	0.45 c	0.31 b	4.78 b	15.05 e	1.88 g
	Federico II	2.05 c	0.48 b	0.70 a	0.39 b	0.18 e	0.17 b	0.34 c	0.88 a	0.45 c	0.31 b	5.12 b	16.83 d	2.03 d
	Fuego	1.99 d	0.47 b	0.69 b	0.37 c	0.20 c	0.19 b	0.33 c	0.88 a	0.47 c	0.32 b	5.37 a	16.55 d	2.01 e
	Furio Camillo	2.27 a	0.48 a	0.72 a	0.35 c	0.19 c	0.18 b	0.36 a	0.89 a	0.47 c	0.32 b	5.77 a	18.40 a	2.18 a
	Incanto	1.80 h	0.45 d	0.66 c	0.41 b	0.20 b	0.19 b	0.3 d	0.86 b	0.46 c	0.32 b	4.94 b	13.97 f	1.86 h
	Iride	2.00 d	0.46 b	0.69 b	0.39 b	0.19 c	0.18 b	0.33 c	0.87 a	0.46 c	0.32 b	5.30 a	16.03 d	2.01 e
	LG Fructis	1.93 e	0.47 b	0.68 b	0.36 c	0.21 b	0.19 a	0.32 c	0.87 a	0.48 c	0.33 a	5.44 a	15.79 d	1.96 f
	Maciste	2.10 b	0.48 a	0.71 a	0.37 c	0.19 d	0.18 b	0.34 c	0.88 a	0.46 c	0.32 b	5.32 a	17.93 b	2.07 c
	Mameli	1.94 e	0.46 c	0.68 b	0.41 b	0.18 e	0.17 c	0.33 c	0.87 a	0.45 c	0.31 b	4.99 b	15.75 d	1.99 e
	Marakas	1.87 f	0.45 c	0.67 b	0.40 b	0.19 c	0.18 b	0.32 c	0.86 a	0.45 c	0.32 b	5.02 b	15.09 e	1.94 f
	Marco Aurelio	2.00 d	0.47 b	0.69 b	0.38 b	0.20 c	0.18 b	0.33 c	0.87 a	0.46 c	0.32 b	5.30 a	16.17 d	2.00 e
	Monastir	1.97 d	0.46 b	0.68 b	0.38 b	0.20 b	0.19 b	0.32 c	0.87 a	0.47 c	0.33 b	5.38 a	15.18 e	1.96 f
	Nuraghe	1.73 i	0.44 d	0.64 c	0.42 b	0.20 b	0.19 b	0.28 d	0.85 c	0.46 c	0.32 b	4.77 b	13.30 g	1.80 i
	Panoramix	1.59 j	0.42 e	0.61 d	0.48 a	0.18 e	0.17 c	0.25 e	0.83 d	0.43 d	0.30 c	3.97 c	11.86 g	1.70 j
	RGT Aventadur	2.03 c	0.48 a	0.71 a	0.36 c	0.19 c	0.18 b	0.35 c	0.89 a	0.47 c	0.32 b	5.49 a	18.17 a	2.08 c
	RGT Jasdur	1.79 h	0.46 b	0.67 b	0.37 c	0.21 b	0.19 a	0.30 d	0.87 a	0.47 c	0.33 b	5.14 b	15.61 e	1.89 g
	RGT Natur	1.93 e	0.45 d	0.66 b	0.42 b	0.19 d	0.18 b	0.31 c	0.86 b	0.45 c	0.31 b	4.93 b	14.02 f	1.93 f
	RGT Voltadur	1.82 g	0.45 d	0.66 c	0.40 b	0.20 b	0.19 b	0.30 d	0.86 b	0.46 c	0.32 b	5.11 b	14.26 f	1.89 g
	SY Leonardo	1.88 f	0.46 b	0.67 b	0.39 b	0.20 c	0.18 b	0.31 c	0.87 a	0.46 c	0.32 b	5.03 b	15.35 e	1.91 g
	SY Prodigio	1.95 e	0.48 b	0.70 a	0.35 c	0.21 b	0.19 a	0.33 c	0.88 a	0.48 b	0.33 a	5.57 a	17.17 c	2.01 e
	Tancredi	2.25 a	0.48 b	0.71 a	0.37 c	0.19 c	0.18 b	0.36 b	0.88 a	0.46 c	0.32 b	5.58 a	17.22 c	2.14 b
	Telemaco	1.86 f	0.46 c	0.67 b	0.41 b	0.19 d	0.18 b	0.31 c	0.86 b	0.45 c	0.31 b	4.81 b	15.00 e	1.91 g
	Tito Flavio	1.92 e	0.46 b	0.68 b	0.39 b	0.20 c	0.18 b	0.32 c	0.87 a	0.46 c	0.32 b	5.12 b	15.56 e	1.95 f
	Verace	2.02 c	0.46 b	0.69 b	0.39 b	0.19 d	0.18 b	0.34 c	0.87 a	0.45 c	0.32 b	5.25 a	16.20 d	2.05 d
	Min	1.59	0.42	0.61	0.35	0.18	0.17	0.25	0.83	0.43	0.30	3.97	11.86	1.70
	Max	2.27	0.48	0.72	0.48	0.21	0.20	0.36	0.89	0.49	0.34	5.85	18.40	2.18
	Mean	1.94	0.46	0.68	0.39	0.19	0.18	0.32	0.87	0.46	0.32	5.16	15.62	1.96
	Δ VI	0.68	0.06	0.11	0.13	0.03	0.03	0.11	0.06	0.06	0.04	1.88	6.54	0.48
	p-value	<0.001	<0.001	<0.001	0.003	0.003	0.009	<0.001	<0.001	0.019	0.017	0.001	<0.001	<0.001

. Among the different VIs, Aureo, Furio Camillo, Panoramix (only for MCARI2), and RGT Aventadur cultivars showed the highest while Beltorax and Panoramix (for most of the indices) showed the lowest values; intermediate values were retrieved in the other cultivars. The differences among cultivars were always significant for all VIs investigated (p < 0.05). The most notable differences among wheat varieties’ VIs responses were observed for SR, RTVI, CVI, SR_RE, MCARI 2, GNDVI, and NDRE.

Table 4. Vegetation Indices (VIs) data observed among the 30 wheat cultivars investigated in May, related statistics (minimum, maximum, difference and mean values) and results of ANOVA analysis. Different letters indicate significant differences between cultivars (Tukey’s HSD test at p < 0.05). For each VI in bold, the highest value is highlighted, while in italics, the lowest value.

		Vegetation Indices—May
		CVI	GNDRE	GNDVI	MCARI2	MTVI	MTVI2	NDRE	NDVI	OSAVI	RDVI	RTVI	SR	SR RE
Cultivars	Antalis	2.89 d	0.48 c	0.74 b	0.17 a	0.19 a	0.18 a	0.41 b	0.81 c	0.32 a	0.47 a	6.37 a	17.25 c	2.41 c
	Aureo	2.35 l	0.45 f	0.70 e	0.18 a	0.21 a	0.20 a	0.37 f	0.79 f	0.33 a	0.48 a	6.35 a	14.87 f	2.20 h
	Beltorax	2.34 l	0.45 f	0.68 f	0.16 a	0.19 a	0.18 a	0.36 g	0.77 g	0.31 a	0.45 a	5.66 a	14.89 f	2.15 i
	Bering	2.67 g	0.47 d	0.72 c	0.16 a	0.18 a	0.17 a	0.39 d	0.80 e	0.31 a	0.45 a	5.78 a	16.02 d	2.29 f
	Brigante	2.76 e	0.49 a	0.73 b	0.17 a	0.19 a	0.18 a	0.39 c	0.81 c	0.32 a	0.46 a	6.13 a	17.79 c	2.33 d
	Claudio	2.49 j	0.46 e	0.71 e	0.18 a	0.2 a	0.19 a	0.37 f	0.79 e	0.33 a	0.47 a	6.13 a	15.25 e	2.19 h
	Diogene	2.73 e	0.48 b	0.73 b	0.17 a	0.19 a	0.18 a	0.40 c	0.81 c	0.32 a	0.46 a	6.03 a	17.73 c	2.34 d
	Federico II	3.07 c	0.49 a	0.74 a	0.15 a	0.17 a	0.16 a	0.41 b	0.81 c	0.31 a	0.44 a	5.76 a	17.16 c	2.42 c
	Fuego	2.52 i	0.47 c	0.72 c	0.18 a	0.20 a	0.19 a	0.39 c	0.81 d	0.33 a	0.48 a	6.36 a	17.74 c	2.32 e
	Furio Camillo	3.14 b	0.48 b	0.75 a	0.16 a	0.18 a	0.17 a	0.43 a	0.82 b	0.32 a	0.46 a	6.36 a	18.47 b	2.56 a
	Incanto	2.55 h	0.46 e	0.70 e	0.16 a	0.18 a	0.17 a	0.37 f	0.79 f	0.31 a	0.45 a	5.69 a	15.06 e	2.21 g
	Iride	2.67 g	0.47 c	0.73 c	0.17 a	0.19 a	0.18 a	0.40 c	0.81 d	0.32 a	0.46 a	6.11 a	16.96 d	2.34 d
	LG Fructis	2.68 g	0.46 e	0.72 d	0.16 a	0.18 a	0.17 a	0.39 c	0.79 e	0.31 a	0.44 a	5.77 a	15.72 e	2.31 e
	Maciste	3.06 c	0.49 a	0.76 a	0.17 a	0.19 a	0.18 a	0.43 a	0.83 a	0.33 a	0.47 a	6.47 a	19.30 a	2.52 a
	Mameli	2.63 g	0.48 b	0.74 a	0.18 a	0.20 a	0.19 a	0.41 b	0.82 b	0.34 a	0.48 a	6.66 a	19.47 a	2.43 c
	Marakas	2.56 h	0.47 d	0.72 c	0.18 a	0.20 a	0.19 a	0.39 c	0.80 d	0.33 a	0.47 a	6.30 a	16.53 d	2.30 e
	Marco Aurelio	2.77 e	0.47 c	0.72 d	0.15 a	0.18 a	0.17 a	0.39 d	0.79 e	0.31 a	0.44 a	5.67 a	16.41 d	2.30 e
	Monastir	2.70 f	0.47 c	0.73 c	0.17 a	0.20 a	0.18 a	0.40 c	0.81 d	0.33 a	0.47 a	6.33 a	16.54 d	2.34 d
	Nuraghe	2.37 k	0.44 f	0.68 f	0.16 a	0.19 a	0.17 a	0.35 h	0.77 h	0.31 a	0.44 a	5.44 a	13.42 h	2.09 j
	Panoramix	2.18 m	0.44 g	0.68 f	0.16 a	0.19 a	0.17 a	0.34 i	0.77 g	0.31 a	0.45 a	5.30 a	13.95 h	2.06 k
	RGT Aventadur	2.88 d	0.48 b	0.75 a	0.17 a	0.19 a	0.18 a	0.42 a	0.82 b	0.33 a	0.47 a	6.57 a	18.59 b	2.48 c
	RGT Jasdur	2.41 j	0.44 f	0.67 f	0.14 a	0.17 a	0.16 a	0.34 i	0.76 i	0.30 a	0.42 a	5.01 a	13.08 i	2.07 k
	RGT Natur	2.64 g	0.45 f	0.71 e	0.16 a	0.19 a	0.18 a	0.39 d	0.79 e	0.32 a	0.46 a	6.02 a	14.60 g	2.28 f
	RGT Voltadur	2.47 j	0.46 e	0.71 e	0.18 a	0.20 a	0.19 a	0.38 e	0.80 e	0.33 a	0.47 a	6.31 a	15.85 e	2.25 g
	SY Leonardo	2.54 h	0.46 e	0.71 e	0.16 a	0.19 a	0.18 a	0.38 e	0.79 e	0.32 a	0.45 a	5.87 a	15.82 e	2.25 g
	SY Prodigio	2.44 j	0.46 d	0.72 d	0.19 a	0.21 a	0.20 a	0.38 d	0.80 d	0.34 a	0.48 a	6.48 a	16.75 d	2.27 f
	Tancredi	3.19 a	0.49 a	0.75 a	0.17 a	0.19 a	0.18 a	0.43 a	0.82 b	0.32 a	0.46 a	6.50 a	17.98 c	2.54 a
	Telemaco	2.65 g	0.48 c	0.73 c	0.17 a	0.2 a	0.18 a	0.39 c	0.81 d	0.33 a	0.47 a	6.21 a	16.79 d	2.31 e
	Tito Flavio	2.46 j	0.46 e	0.71 e	0.16 a	0.18 a	0.17 a	0.37 f	0.79 e	0.31 a	0.45 a	5.62 a	15.49 e	2.20 h
	Verace	2.69 f	0.48 c	0.74 a	0.18 a	0.20 a	0.19 a	0.42 b	0.82 b	0.34 a	0.48 a	6.73 a	18.54 b	2.44 c
	Min	2.18	0.44	0.67	0.14	0.17	0.16	0.34	0.76	0.30	0.42	5.01	13.08	2.06
	Max	3.19	0.49	0.76	0.19	0.21	0.20	0.43	0.83	0.34	0.48	6.73	19.47	2.56
	Mean	2.65	0.47	0.72	0.17	0.19	0.18	0.39	0.80	0.32	0.46	6.07	16.47	2.31
	Δ VI	1.01	0.05	0.09	0.05	0.04	0.04	0.09	0.07	0.04	0.06	1.72	6.39	0.50
	p-value	<0.001	<0.001	<0.001	0.310	0.120	0.220	<0.001	<0.001	0.471	0.416	0.525	<0.001	<0.001

shows the VIs responses of the 30 wheat cultivars tested and the ANOVA results related to the survey carried out in May. The highest values were recorded in SY Prodigio, Maciste, Federico II, Furio Camillo, Mameli, Tancredi, and Verace, while the lowest were observed in RGT Jasdur, RGT Natur, Panoramix, and Marco Aurelio; intermediate values were retrieved in the other cultivars. Differing from what was noted in April, in May, not all VIs showed statistically significant differences (p < 0.05) among the cultivars. In particular, significant differences were observed for CVI, GNDRE, GNDVI, NDRE, NDVI, SR, and SR_RE, whereas not substantial were the differences for MCARI2, MTVI, MTVI2, OSAVI, RDVI, and RTVI. Among the wheat cultivars, the highest differences in VI responses were observed for SR, RTVI, CVI, and SR_RE.

The PCA analysis carried out with the VI response data of the 30 wheat cultivars separately for the April and May surveys are presented in Figure 3 and Figure 4, respectively. With regard to the first survey (April), PC1 accounted for 78.0% of the total variance, while PC2 accounted for 19.3%. In particular, PC 1 discriminated Panoramix, Nuraghe, Beltorax, Bering from Aureo, SY Prodigio, RGT Aventadur, Maciste, Tancredi, and Furio Camillo. These cultivar groups were separated from the others located in the central part of the chart. Conversely, PC2 separated Panoramix, Monastir, SY Prodigio, Brigante, LG Fructis, Incanto, RGT Voltadur, RGT Jasdur, Nuraghe, and Aureo from RGT Aventadur, Verace, Mameli, Federico II, Maciste, Tancredi, and Furio Camillo and both of these groups from the other varieties sited in the middle of the graph.

The PCA developed from the VIs data of May highlighted that PC1 accounted for 78.0% and PC2 for 19.3%, respectively. Compared with the PCA of April, the differences between cultivars were limited; PC1 separated RGT Jasdur, Beltorax, Nuraghe and Panoramix from Verace, Mameli, RGT Aventadur, Tancredi and Maciste, and both from the other cultivar situated in the middle of the chart. PC2 has discriminated against Marakas, Beltorax, Fuego, Nuraghe, Panoramix, Claudio, RGT Voltadur, SY Prodigio, and Aureo from Marco Aurelo, Maciste, Tancredi, Furio Camillo, and Federico II; both these groups were separated from other cultivars, the most numerous, that were sited in the middle part of the chart.

Overall, the information retrieved highlighted a significantly different VIs response of most of the tested cultivars, especially at the first survey epoch, at late tillering. This evidence agrees with Hassan et al. [54], who compared the NDVI response of 32 different wheat varieties in China, and Marino and Alvino [55], who compared the response of three indices on 10 different varieties of durum wheat in Italy. Therefore, it is crucial to consider the different spectral responses of each cultivar [55] before RS and PA implementation and in the translation of these technologies from one field to another. The reduction in the differences among cultivars from April to May could be ascribed to a different and contrasting phenomenon affecting plant response, such as saturation [56], decline in plant vigour and awn presence [57].

3.2. Pearson’s Correlations Analysis between Grain Yield Data and Vegetation Indices (VIs) Responses

Pearson’s correlations analysis between grain yield data and VIs showed different behaviours of the 30 wheat varieties at the different VIs and a specific sensitivity to the survey epoch (

Table 5. Pearson’s correlation coefficients calculated between the grain yield of each wheat cultivar and its vegetation indices (VIs) response surveyed in April. For each VI in bold, the highest value is highlighted, while in italics, the lowest value.

			Vegetation Indices—April
		CVI	GNDRE	GNDVI	MCARI2	MTVI	MTVI2	NDRE	NDVI	OSAVI	RDVI	RTVI	SR	SR RE
Cultivars	Antalis	0.74	0.97	0.94	0.63	0.15	0.23	0.82	0.98	0.51	0.51	0.41	0.90	0.76
	Aureo	0.98	0.62	0.80	0.74	0.75	0.74	0.93	0.56	0.74	0.74	0.87	0.69	0.94
	Beltorax	0.19	0.49	0.36	0.56	0.80	0.76	0.08	0.38	0.66	0.67	0.45	0.19	0.02
	Bering	0.99	0.99	0.99	0.99	0.92	0.95	0.99	0.98	0.98	0.98	0.98	0.99	0.99
	Brigante	0.98	0.38	0.39	0.01	0.03	0.02	0.34	0.15	0.02	0.23	0.02	0.05	0.31
	Claudio	0.19	0.14	0.49	0.28	0.19	0.20	0.99	0.42	0.22	0.23	0.39	0.39	0.99
	Diogene	0.85	0.97	0.96	0.96	0.82	0.88	0.96	0.98	0.90	0.91	0.92	0.99	0.97
	Federico II	0.97	0.91	0.91	0.95	0.99	0.99	0.92	0.87	0.98	0.98	0.98	0.89	0.93
	Fuego	0.25	0.21	0.74	0.77	0.01	0.05	0.88	0.83	0.23	0.25	0.67	0.76	0.88
	Furio Camillo	0.61	0.26	0.31	0.33	0.88	0.70	0.35	0.15	0.53	0.50	0.54	0.17	0.35
	Incanto	0.76	0.96	0.95	0.57	0.20	0.26	0.95	0.66	0.35	0.37	0.60	0.82	0.95
	Iride	0.38	0.08	0.11	0.03	0.14	0.07	0.18	0.03	0.10	0.01	0.01	0.09	0.19
	LG Fructis	0.91	0.96	0.90	0.99	0.98	0.99	0.83	0.84	0.99	0.99	0.99	0.93	0.83
	Maciste	0.87	0.03	0.15	0.54	0.63	0.62	0.22	0.42	0.58	0.58	0.53	0.23	0.15
	Mameli	0.72	0.56	0.49	0.34	0.29	0.30	0.42	0.39	0.35	0.34	0.35	0.36	0.39
	Marakas	0.99	0.06	0.67	0.29	0.13	0.15	0.93	0.49	0.27	0.26	0.62	0.44	0.94
	Marco Aurelio	0.88	0.24	0.54	0.02	0.47	0.33	0.73	0.39	0.04	0.03	0.15	0.38	0.77
	Monastir	0.03	0.19	0.11	0.08	0.09	0.08	0.02	0.22	0.12	0.11	0.02	0.06	0.01
	Nuraghe	0.70	0.57	0.54	0.63	0.68	0.68	0.57	0.48	0.64	0.64	0.68	0.64	0.62
	Panoramix	0.27	0.98	0.76	0.92	0.99	0.99	0.58	0.92	0.98	0.97	0.73	0.81	0.57
	RGT Aventadur	0.72	0.99	0.99	0.99	0.98	0.98	0.99	0.91	0.99	0.99	0.97	0.99	0.97
	RGT Jasdur	0.21	0.84	0.75	0.77	0.77	0.76	0.43	0.98	0.78	0.78	0.62	0.66	0.34
	RGT Natur	0.02	0.03	0.04	0.24	0.83	0.76	0.02	0.05	0.53	0.51	0.14	0.02	0.01
	RGT Voltadur	0.30	0.88	0.78	0.80	0.80	0.80	0.58	0.89	0.81	0.81	0.66	0.80	0.53
	SY Leonardo	0.73	0.61	0.69	0.45	0.15	0.19	0.73	0.67	0.28	0.30	0.49	0.65	0.73
	SY Prodigio	0.92	0.81	0.99	0.89	0.05	0.05	0.95	0.97	0.33	0.37	0.71	0.99	0.92
	Tancredi	0.04	0.07	0.04	0.43	0.88	0.79	0.02	0.19	0.59	0.60	0.27	0.15	0.01
	Telemaco	0.03	0.05	0.13	0.06	0.01	0.02	0.25	0.18	0.04	0.04	0.08	0.17	0.23
	Tito Flavio	0.99	0.23	0.87	0.03	0.41	0.29	0.96	0.42	0.08	0.07	0.56	0.67	0.96
	Verace	0.93	0.92	0.99	0.99	0.99	0.99	0.98	0.99	0.99	0.99	0.99	0.99	0.99
	Min	0.02	0.03	0.04	0.01	0.01	0.02	0.02	0.03	0.02	0.01	0.01	0.02	0.01
	Max	0.99	0.99	0.99	0.99	0.99	0.99	0.99	0.99	0.99	0.99	0.99	0.99	0.99
	Mean	0.61	0.53	0.61	0.54	0.53	0.52	0.62	0.58	0.52	0.53	0.55	0.56	0.61

).

In particular, with regard to the April survey (Table 5), indices that showed the greatest number of significant correlations (R² > 0.7) with grain yield were CVI (11), GNDVI (10), MTVI (10), MTVI2 (9), NDRE (10), and SR RE (10). In contrast, the index that showed the worst performance was MCARI2. Among the different cultivars, those that showed more good correlations with VIs were Antalis, Aureo, Bering, Diogenes, Federico II, Fuego, LG Fructis, Panoramix, RGT Aventadur, RGT Jasdur, RGT Voltadur, SY Prodigio, Tito Flavio, and Verace (Table 5). With regard to the May survey, the most significant number of relevant correlations (R² > 0.7) were observed for CVI (9), GNDRE (10), NDRE (9), NDVI (9), OSAVI (9), RDVI (9), RTVI (9), SR (10), and SERE (9) (

Table 6. Pearson’s correlation coefficients calculated between the grain yield of each wheat cultivar and its vegetation indices (VIs) response surveyed in May. For each VI in bold, the highest value is highlighted, while in italics, the lowest value.

			Vegetation Indices—May
		CVI	GNDRE	GNDVI	MCARI2	MTVI	MTVI2	NDRE	NDVI	OSAVI	RDVI	RTVI	SR	SR RE
Cultivars	Antalis	0.59	0.66	0.73	0.01	0.05	0.04	0.58	0.79	0.02	0.02	0.02	0.10	0.45
	Aureo	0.03	0.81	0.22	0.39	0.48	0.45	0.04	0.32	0.46	0.45	0.45	0.08	0.07
	Beltorax	0.93	0.97	0.99	0.99	0.99	0.99	0.99	0.99	0.98	0.98	0.98	0.99	0.99
	Bering	0.76	0.09	0.62	0.12	0.06	0.08	0.84	0.39	0.11	0.12	0.12	0.44	0.87
	Brigante	0.99	0.75	0.63	0.97	0.95	0.97	0.02	0.59	0.99	0.99	0.99	0.02	0.21
	Claudio	0.38	0.01	0.02	0.03	0.04	0.01	0.03	0.04	0.02	0.01	0.04	0.05	0.03
	Diogene	0.52	0.49	0.62	0.26	0.12	0.18	0.79	0.64	0.21	0.24	0.24	0.92	0.84
	Federico II	0.99	0.98	0.12	0.34	0.37	0.37	0.24	0.69	0.37	0.38	0.38	0.13	0.34
	Fuego	0.02	0.89	0.58	0.45	0.37	0.41	0.33	0.53	0.43	0.42	0.42	0.93	0.33
	Furio Camillo	0.85	0.16	0.31	0.58	0.48	0.56	0.76	0.21	0.84	0.89	0.89	0.96	0.97
	Incanto	0.27	0.71	0.80	0.15	0.18	0.18	0.73	0.77	0.33	0.34	0.34	0.23	0.21
	Iride	0.90	0.68	0.99	0.50	0.48	0.55	0.45	0.96	0.89	0.92	0.92	0.11	0.35
	LG Fructis	0.99	0.91	0.88	0.90	0.85	0.88	0.80	0.86	0.92	0.92	0.92	0.80	0.79
	Maciste	0.57	0.87	0.64	0.01	0.03	0.02	0.29	0.68	0.18	0.02	0.04	0.02	0.45
	Mameli	0.56	0.48	0.64	0.97	0.99	0.98	0.78	0.63	0.92	0.92	0.92	0.65	0.78
	Marakas	0.99	0.90	0.99	0.67	0.60	0.65	0.95	0.99	0.77	0.77	0.77	0.98	0.94
	Marco Aurelio	0.61	0.56	0.77	0.99	0.99	0.99	0.97	0.75	0.99	0.98	0.98	0.98	0.99
	Monastir	0.34	0.10	0.69	0.62	0.58	0.59	0.98	0.97	0.55	0.55	0.55	0.96	0.71
	Nuraghe	0.04	0.89	0.99	0.99	0.99	0.99	0.98	0.98	0.99	0.99	0.99	0.99	0.93
	Panoramix	0.28	0.88	0.78	0.60	0.61	0.62	0.57	0.82	0.67	0.67	0.67	0.71	0.49
	RGT Aventadur	0.94	0.57	0.97	0.86	0.81	0.83	0.69	0.97	0.89	0.91	0.91	0.96	0.69
	RGT Jasdur	0.01	0.18	0.20	0.16	0.20	0.19	0.17	0.23	0.21	0.22	0.22	0.10	0.15
	RGT Natur	0.54	0.35	0.02	0.01	0.02	0.02	0.03	0.05	0.02	0.02	0.02	0.05	0.02
	RGT Voltadur	0.33	0.32	0.28	0.96	0.98	0.97	0.28	0.24	0.88	0.89	0.89	0.33	0.31
	SY Leonardo	0.02	0.37	0.02	0.17	0.01	0.03	0.02	0.14	0.01	0.02	0.01	0.90	0.02
	SY Prodigio	0.08	0.49	0.45	0.99	0.98	0.98	0.33	0.70	0.98	0.99	0.99	0.84	0.23
	Tancredi	0.03	0.92	0.35	0.66	0.73	0.72	0.21	0.56	0.72	0.72	0.72	0.49	0.25
	Telemaco	0.99	0.97	0.95	0.13	0.48	0.28	0.85	0.97	0.11	0.15	0.15	0.87	0.87
	Tito Flavio	0.97	0.66	0.45	0.46	0.42	0.42	0.44	0.35	0.40	0.40	0.40	0.52	0.51
	Verace	0.99	0.98	0.92	0.32	0.63	0.58	0.03	0.73	0.45	0.37	0.37	0.96	0.01
	Min	0.01	0.01	0.02	0.01	0.01	0.01	0.02	0.04	0.01	0.01	0.01	0.02	0.01
	Max	0.99	0.98	0.99	0.99	0.99	0.99	0.99	0.99	0.99	0.99	0.99	0.99	0.99
	Mean	0.55	0.62	0.59	0.51	0.52	0.52	0.51	0.62	0.54	0.54	0.54	0.57	0.49

). Also, for the survey of May, the worst performance was observed for the correlations with MCARI2. Among the tested cultivars, those that showed more relevant correlations between their yield variations and VIs were Beltorax, Brigante, Furio Camillo, Iride, LG Fructis, Mameli, Marakas, Marco Aurelio, Nuraghe, Panoramix, RGT Aventadur, RGT Voltadur, SY Prodigio, Tancredi, and Verace (Table 6).

Some wheat varieties showed a contrasting Pearson correlation trend between the two survey epochs. Indeed, Aureo, Bering, Diogene, Federico II, Fuego, RGT Jasdur, and Verace highlighted an increased number of relevant correlations with VIs measured in April, whereas Beltorax, Brigante, Furio Camillo, Iride, Mameli, Marakas, Marco Aurelio, Nuraghe, Tancredi, and Telemaco have highlighted good correlations with the May data. LG Fructis, Panoramix, and RGT Aventadur showed a high number of significant correlations in April and May.

The VIs are used as sensitive and reliable indicators to assess crop status like many plants’ growth, health and productivity. Relationships between VIs and crop grain yield based on linear regression are widely used because their easy computational implementation provides a simple and efficient tool for obtaining good predictions within certain limits. Moreover, although hyperspectral sensors are more used today, they have a high cost that is not always affordable for medium and small farms (the most diffused in the Mediterranean area) and provide data that are often redundant [58,59], and which generally require very high hardware resources and computation times that are not always compatible with obtaining timely and reliable information, which is crucial in PA.

Concerning the present experiment, if in April better correlations were observed between VIs and grain yield, in May, a better predictive ability among the different VIs and wheat varieties was retrieved. Therefore, this evidence confirms that the best period of wheat survey for yield prediction begins around two months before presumed harvest time [54,55,60,61] when the yield predictions are stabilised, and the leaves have not started the ageing process, which alters their reflectivity and makes them unsuitable for providing information about the status of the plant [62]. The different correlation performances among the tested wheat varieties confirm, as also observed in VIs responses, their differential spectral behaviour and responsiveness to yield variations. In particular, as reported by several authors, not all varieties responded equally to the different VIs and not always to both survey times [55,63,64]. Therefore, a careful choice of varieties and VIs is needed for a proficient application of PA techniques, considering the varieties’ different responses and the indices’ sensitivities [55,65].

Among VIs, globally on all varieties, CVI, GNDVI, NDRE, and SR_RE highlighted the best correlations (R² > 0.6) in April, while GNDRE and NDVI had the best in May. That is in disagreement with the experiment conducted by Gonzalez-Dugo et al. [33] that reported a non-significant correlation between VIs and wheat grain yield in Mediterranean conditions. In addition, the specific condition of the present experiment, where the in-field variability was minimised (homogeneity of soil characteristics and equal management system) and a limited difference in grain yield between plots (~0.6 t ha⁻¹, on average), highlighted a high sensitivity of responsive indices.

Typically, NDVI correlates with the fraction of photosynthetically active radiation [66] and its ability to assess plant vigour and yield is widely recognised in wheat (i.e., [54,56,62,63,67,68]). In our study, its potential to predict yields emerged in May, but not in April, on a significant number of varieties, highlighting the best performance among all VIs equal to 0.62. GNDVI is another widely used VI for wheat yield estimation correlated with chlorophyll content [41]. In our study, GNDVI, across all 30 cultivars, achieved a fairly good performance of 0.61 in April, while just one month later (May survey), its performance significantly decreased (for both overall R² value and number of varieties with R² > 0.7). This observation agrees with Yang et al. [69], who observed a better correlation with wheat biomass at an early survey than at the grain filling stage due to an improved ability of this index to predict crop water stress. Other good correlations were observed from NDRE, GNDRE, and SR_RE, which use red edge and green bands, permitting good discrimination among cultivars as well as grain yield productivity estimation, especially at an early survey, detecting water stress and chlorophyll content in plant tissue thus well-characterising crop canopy status [69,70]. These bands allow an improved assessment of vegetation status overcoming saturation and showing higher performance, according to Fu et al. [62]. Among the various VIs tested, the CVI showed good correlation performance in April, being able to assess the chlorophyll content on the canopy [71], which in turn correlated well with grain yield. On the contrary, the correlations observed for OSAVI and RDVI were particularly poor, especially if compared with what has been observed by other authors [55,56].

3.3. Machine Learning (ML) Approaches for Grain Yield Estimation

The validation results, observed in 50% of the 100 repetitions performed of the five models tested, are presented below in

Table 7. Machine learning (ML) algorithms validation performance-related data of the two survey periods: coefficient of determination (R²), root mean square error (RMSE) and mean absolute error (MAE).

		April			May
		R2	RMSE	MAE	R2	RMSE	MAE
Models	Linear model (LM)	0.82	0.37	0.31	0.82	0.49	0.39
	Random forest (RF)	0.88	0.18	0.16	0.84	0.36	0.27
	Support Vector Machine (SVM)	0.87	0.35	0.27	0.81	0.40	0.34
	K-nearest neighbors (k-NN)	0.86	0.27	0.23	0.85	0.44	0.36
	Neural network (NN)	0.71	1.00	1.90	0.68	1.94	1.90

. In general, all tested models showed a good predictive ability, evidenced by the observed coefficient of determination never below 0.68. Among the survey epochs, better performance was observed when models were calibrated with April’s VIs data. Regardless of the epoch of VIs monitoring (the models’ trend was the same among them), lower performance was highlighted by NN. In contrast, the best prediction was achieved by applying the RF model. RF, SVM, and k-NN showed high performances very close to each other, especially when calibrated with spectral response data recorded in April. Moreover, analysing the values of RMSE and MAE, a good level of error was observed for RF, with an RMSE of 0.18 t ha⁻¹ and k-NN (RSME = 0.27 t ha⁻¹), confirming the goodness of predictions by these two models. In contrast, the predictions of the NN model were inaccurate (Table 7).

In general, the performance achieved by the different ML approaches was always higher than Pearson’s correlations between VIs and grain yield, confirming their reliability on crop yield prediction by capturing non-linear relationships and showing robustness against spurious data, according to several authors [35,62,67,72,73,74]. In particular, ML results across the different models showed R² values never below 0.68, while the best Pearson’s correlations R² values have never exceeded a value of 0.62. Considering that ML algorithms were validated in a different dataset (different training and validation datasets), compared to Pearson’s correlations that were tested on the whole dataset, further highlights the predictive capabilities of this approach. The RF at both survey epochs showed the best predictive performance in terms of all the statistics calculated (R², RMSE and MAE), in accordance with several authors [34,75,76,77]. Following was the performance of k-NN. Our results are in accordance with Bebie et al. [78] and Chergui [79], who observed the best performance by RF and k-NN models for the durum wheat grain yield prediction cultivated in a Mediterranean environment, and Yue et al. [73] and Zhou et al. [67], who observed better performance from the RF and SVM models for wheat biomass prediction.

In particular, the potential of RF lies in its structure as an ensemble learning method where many decision trees are trained, validated and mediated to achieve the best prediction by minimising variance [77]. k-NN uses another approach by finding relationships between independent variables and the predicted outcome by averaging the observations in the same neighbourhood [80]. The SVM model is an approach similar to linear regression where the trained function is a straight line, referred to as a hyperplane, that best fits the data points while minimising the errors that stand on key points, called support vectors, that determine the position and orientation of the hyperplane [81]. Based on the results achieved in the present study, all three methods have well interpreted the data with a non-linear correlation by overcoming issues such as data noise (affected by saturation and or soil reflectance interference), collinearity and overfitting handling numerous input variables providing good prediction of durum wheat yields across different cultivars. This can enable more accessible and more direct use of information obtained from drone MS surveys without needing to select the best indices. Surprisingly, although widely used in different predictive approaches in agriculture and other sectors, the NN showed the worst performance, even lower than LM, typically chosen as the reference algorithm. A similar trend was observed by Wang et al. [82] that highlighted the lower wheat biomass prediction ability of ANN compared with RF and SVM. This finding can be ascribed to the power of RF, k-NN and SVM to elaborate a small amount of data instead of the NN model that works with a large amount of data. In addition, the constructed NN with only three layers was probably not the best to analyse the many (even affected by collinearity) and different pieces of information from the VIs. Also, in the case of the ML approaches, using early VIs data from the April survey provides a better grain yield prediction than the May survey, as postulated above and according to Fu et al. [62] and Zhu et al. [83].

4. Conclusions

Reliable techniques and methods to discriminate plant reflectance responses and assess their status and productivity are of fundamental importance for implementing PA. This approach can be very important to improve crop productivity and sustainability by enabling reliable estimation of yields to plan storage, sales, and purchases and, thus, food security. This study aimed to compare the spectral response of thirty durum wheat varieties and to predict their grain yields based directly on VIs or using different ML approaches fed with VIs data. In particular, this study involved the largest number of durum varieties commonly grown in Italy (i.e., thirty), analysing their response at two different and specific phenological stages. Spectral VIs have allowed us to separate the tested wheat varieties into groups, and differently if surveyed in April or May. Although the VIs showed, for most cases, significant differences among the tested varieties, especially in April, CVI, NDRE, RTVI, and SR showed the most remarkable differences between the cultivars. The performance of correlations between grain yield and VIs showed a significant variability among the tested varieties. Correlations were good (R² > 0.7) for ten varieties, while lower performances were observed for others. The VIs that best correlated with grain yield were CVI, GNDVI, MTVI, MTVI2, NDRE, and SR_SR. The ML approaches permitted the improvement of the yield prediction from VIs data, especially when RF and SVM models were used. For both yield prediction approaches, VIs correlations and ML, the April survey allowed us to estimate better than the May survey. The present study involved thirty different cultivars grown in a fairly homogeneous soil with an unique management system, studied during a single crop cycle and monitored at only two points in the crop cycle. Results could be affected by many variables, like different soil properties, management choices, climate conditions, surveying dates and spatial and spectral resolution of multispectral sensor. Therefore, based on this consideration, further study should be carried out to include soil and management variability and implementing different soil types and climate conditions with a more intense UAV monitoring campaign to fully assess the potentiality of VIs and ML approaches to predict durum wheat grain yield under different conditions.