Field trials were conducted at the research station of the Technical University of Munich in Dürnast, Freising (48.40630° N, 11.69535° E) in the growing seasons of 2020/2021 and 2021/2022 further referred to as seasons 2021 and 2022. The soil at this location can be characterized by a homogeneous Cambisol with 20.8% clay, 61.5% silt and 16.6% sand. Precipitation during this period was 595 mm and 415 mm in the first and the second season, respectively. The average temperature was 7.2°C in the first and 8.0°C in the second season ( Supplementary Figure 1 ). A lot of precipitation around flowering characterized season 2021 whereas the season 2022 suffered too little precipitation at the end of the tillering stage. Climate data was collected from a weather station (Station id 5404) operated by the Climate Data Center of the German Weather Service located a few hundred meters from the trials. The temperature was aggregated to phenologically meaningful growing degree-days (GDD) ( Equation 1 ) (Bonhomme, 2000):

where Tmean_d is the mean temperature for day d after sowing as determined by Equation 2 , maxT_d,h and minT_d,h are hourly maximum and minimum temperatures for day d and baseT is the base temperature, which was set to 0°C.

A panel consisting of 19 diverse European winter wheat elite varieties (Triticum aestivum) in 2021 and 18 varieties in 2022 was grown in plots with a size of 10 m x 1.85 m. All varieties grown in 2022 were grown in 2021 as well ( Table 1 ). The plots were arranged in a randomized complete block design with four replicates, resulting in 76 plots in 2021. In 2022, the 72 plots were part of a bigger trial, which was arranged as a randomized strip-plot design with four replicates as well. Orthophotos of the trials can be found in Supplementary Figure 2 in the appendix. All plots used for this study were fertilized by applying 180 kg N ha^-1 in three equal splits at BBCH 25, 32 and 65. Plant protection was carried out according to local practice. Sowing took place on the 10.11.2020 and the 20.10.2021 and all plots were harvested at full maturity on the 03.08.2021 and the 26.07.2022, respectively.

The entire plots were harvested using a combine harvester. The water content of the grains was determined by weighing the grains after harvest, drying them at 65°C until constant weight was reached and weighing them again. The final yield was normalized to a moisture content of 14%. In each year, the three varieties with the lowest average yield were classified as low yielding and the three varieties with the highest average yield as high yielding. The phenology of each plot was visually rated using the BBCH scale (Meier et al., 2009) on a plot level. Leaf area index (LAI) was measured using a Licor 2000 leaf area meter (LI-COR Biosciences, Lincoln, USA.) with a 45° view cap to minimize operator influence. One measurements was taken at the top of the canopy and four measurements were taken under the canopy at three different locations per plot, which were then averaged.

Spectral measurements were acquired using a Phantom 4 Multispectral RTK (DJI, Shenzhen, China) unmanned aerial vehicle (UAV) in the first year. The UAV captures reflectance in wavelengths of 450, 560, 650, 730 and 840 nm and measures the incoming sunlight by a sensor on top of the UAV. Flight height was set to 10 m above ground level (AGL) resulting in a ground sampling distance (GSD) of 0.7 cm. In the second year, images were acquired using a MicaSense Dual Camera Kit (AgEagle Aerial Systems Inc., Wichita, USA) capturing reflectance in wavelengths of 444, 560, 650, 717, 842 nm. The camera was mounted to a DJI Matrice M300 RTK UAV (DJI, Shenzhen, China). Flight height was set to 30 m AGL resulting in a GSD of 2.5 cm. In both years, overlap in both directions was set to 90%. Before and after each flight, images of a panel with a known reflectance were taken. Flights were carried out twice per week during heading and flowering stages and once per week at other stages. First flight was carried out on the 25.03.2021 and the 24.02.2022 and the last flight on the 20.07.2021 and the 27.07.2022 when the canopies were fully senescent. This resulted totally in 19 flights in 2021 and 22 flights in 2022. Images were taken around the solar noon and under sunny conditions, if possible. The images from each flight were mosaicked using the Agisoft Metashape Professional 1.8.4 (Agisoft, St. Petersburg, Russia) structure-from-motion software and were radiometrically calibrated using the reflectance panels on the ground and the incident light sensor on the UAV. The processing parameters used for all flight dates were similar ( Figure 1 ). The point cloud was georeferenced using the real-time kinetic global positioning system (RTK-GPS) integrated into the UAS, with the RTK correction signal provided by SAPOS (Deutsche Landesvermessung, 2023) in 2021 and ground control points were used in 2022. Orthomosaics acquired in 2021 were resampled to the same GSD as in 2022 by the average method implemented in gdal (Gdal/Ogr Contributors, 2023). Reflectance of individual bands was extracted by calculating the median of a specific region of interest (ROI) representing a plot using a custom Python 3.7 script (Python Software Foundation, https://www.python.org/).

To predict yield on a plot level and classify yield performance groups, we employed Random Forest (RF) machine learning models in R 4.2 (R Core Team, 2021). We optimized the number of trees per forest to 500 and used the R package caret (Kuhn, 2008). The number of trees per forest was set to 500 and the number of features per node was optimized by minimizing the root mean square error (RMSE) for the regression models and the accuracy for the classification models if more than one feature was available as in the moving window model.

Pearson correlation coefficient between yield and spectral features was calculated using measurements taken during tillering and flowering. At this date, most varieties were in the mid to end flowering and the correlation of VIs and yield was maximal for most VIs. The performances of the regression RF models were assessed by the coefficient of determination (R²) ( Equation 3 ) as well as the RMSE ( Equation 4 ) using a 10-fold cross validation that was repeated 3 times and averaged:

Where x_i and y_i represent the observed and the predicted yield, x ¯ i and y ¯ i represent the mean of the observed and the predicted yield, respectively. n represents the number of samples. The performances of the classification RF models were assessed by the accuracy ( Equation 5 ) of the prediction using a 10-fold cross validation that was repeated 3 times and averaged:

Substantial grain yield variation was observed between experimental plots, with the season 2021 yielding about 128 g m^-2 less than the second season ( Table 1 ). Variety Julius was classified into the low yielding group in both years, while Skyfall and Dagmar were classified as high yielding in both years. The temperature sum to achieve a specific growth stage did not show significant differences between the yield groups. However, in both years a tendency towards an advanced phenology in the low yielding varieties could be observed ( Figure 2 ). The Leaf Area Index (LAI) was notably higher in the first year of the trial compared to the second. In 2021, the high-yielding varieties showed a significantly higher LAI during the stem elongation, the booting and the late grain filling stage than the low-yielding varieties. This difference could not be observed in the second year ( Figure 2 ).

Figures 3A–D show the Pearson correlation coefficient examining the relationship between reflectance, vegetation indices and TFs for the two years of trial at the end of the tillering and end of the flowering stages for the two years of trial. The analysis reveals that most features exhibit high correlations with one another during the tillering stage in both years, with few exceptions. Exceptions are the REDEDGE band in 2021, the REDEDGE and NIR bands in 2022 and the CORRELATION TF in both years. Correlation to yield at the tillering stage is high for all features in 2021 with a maximal R of 0.61 for the NIR band ( Figure 3A ). In 2022, correlations to yield are generally low, the highest correlation was found for the CORRELATION TF (R = 0.38) ( Figure 3B ).

At the flowering stage, correlation coefficients among features tend to decrease, especially between features of distinct feature groups. For instance, VIs are strongly correlated among themselves, while TFs similarly demonstrate robust correlations within their group. However, features belonging to the VI group are weakly with features belonging to the TF group. The absolute correlation of the features to yield generally decreases for TFs, increases for VIs and single band reflectances show few differences. In 2022, the absolute correlations increase towards the flowering stage, except for the REDEGE reflectance. For both years, the NDRE was the feature showing the highest correlation to yield, whereas the REDEDGE reflectance was not correlated to yield.

Reflectance of the GREEN and the RED band decreased with plant growth during the tillering and stem elongation stages followed by an increased with the onset of the grain filling stage. This trend resulted in a minimal reflectance around booting and flowering stages ( Figure 4 ). In 2021, the high yielding group consistently displayed significantly lower GREEN and RED reflectances across various stages ranging from tillering to late grain filling. However, in the subsequent season, these marked differences between yield groups only emerged from the latter half of the stem elongation stage onward.

Both the RVI and NDRE indices showed significant disparities between the low and high yielding groups in 2021 at all recorded dates, except for the first one. In 2022, these distinctions were noticeable from the flowering stage to the initial segment of the stem elongation stage. Particularly, NDRE exhibited differences from the conclusion of the stem elongation stage onward ( Figure 4 ). The TF CORRELATION decreased until booting and increased towards the end of grain filling in 2021.The feature does not show a clear development with time in 2022. The DISSIMILARITY TF decreases from the tillering until the booting stage and slight increases until harvest in both years. Significant differences between the two yield groups were found in 2021 for stages ranging from tillering to the end of flowering. In 2022, differences were found on few dates after flowering only ( Figure 4 ).

The performance of the yield prediction models depends highly on the features and the time point selected. Overall, the 2021 season demonstrated superior results, exhibiting lower average RMSE in contrast to the 2022 season. Generally, predictions improve from the tillering to the booting stage and deteriorate after flowering ( Figure 5 ). In 2021, the most successful feature was the Normalized Difference Red Edge (NDRE) index at the booting stage, displaying an RMSE of 49.1 g m^-2. In 2022, the Difference Vegetation Index (DVI) at the grain filling stage showed the best performance (RMSE = 60.6 g m^-2). Individual bands, such as the RED band, produced good results, particularly in the 2021 season. However, on average they were outperformed by the VIs in both seasons ( Table 4 ). The TFs were the feature group that led to regression models with the highest RMSE. On average, the RMSE was 10.3 g m^-2 higher for the TFs than the VIs in 2021.

Although data smoothing had a marginal influence the seasonal average RMSE, but worsened a few minimal RMSE. Additionally, smoothing altered the phenological stage at which the minimal RMSE was attained, illustrating shifts in the optimal time points for yield prediction. For instance, after smoothing, the NDRE feature in 2021 displayed an optimal time point for yield prediction during the early grain filling stage as opposed to the unsmoothed time points, which indicated the booting stage as optimal ( Table 4 ).

Combining features from three adjacent time points notably enhanced the average RMSE models by approximately 10 g m^-2 and reduced the minimal RMSE by 5 g m^-2 for both years across all feature groups. Employing the moving window method revealed that the lowest RMSE of 45.0 g m^-2 was achieved using the RED reflectance band at the flowering stage, while in 2022, the Vegetation Index (VI) CCII reached an RMSE of 50.8 g m^-2. The combination of TFs from various dates yielded models that were comparable to models constructed with a single reflectance band or VI. Specifically, the DISSIMILARITY feature reached a minimal RMSE of 57.4 g m^-2 at the booting stage in 2021, and the HOMOGENEITY feature achieved an RMSE of 56.9 g m^-2 at the flowering stage in 2022 ( Table 4 ).

The efficacy of classification models depends on the chosen features and time points, paralleling the observation in regression models. On average, the models in the first season exhibited higher accuracy compared to those in the subsequent season. However, the maximum accuracies achieved in both years were similar. Notably, in 2021, an accuracy of 0.938 was attained using the BLUE band during the late grain filling stage, while the HOMOGENEITY TF in 2022 achieved the maximum accuracy of 0.915 during the grain filling stage. In 2021, VIs generally outperformed single-band reflectances and TFs for yield predictions, whereas in 2022, TFs slightly surpassed the other two groups. The model performances in both years exhibited substantial variability across different dates. Although data smoothing mitigated some fluctuations, significant differences between adjacent dates persisted. Similar to the regression models, the average and maximum performance of the models remained largely consistent after the smoothing of data points.

Incorporating more than one date into the models generally yielded a slightly higher average accuracy for all features in both years. On average, the best-performing feature was the RED band with an accuracy of 0.758 in 2021 and the Modified Chlorophyll Absorption Ratio Index (MCARI) in 2022, displaying an average accuracy of 0.615. However, the maximal accuracies remained largely unaltered compared to those achieved using single dates ( Table 5 ). Despite this, fluctuations between individual dates were reduced ( Figure 6 ). A trend towards increased accuracy with time was evident in 2021, particularly notable with the RED band, consistently yielding high accuracies after heading. In the second year, accuracies improved and were higher than 0.5 for all features at the end of the heading stage. Moreover, fluctuations between dates were notably reduced in this subsequent year ( Figure 6 ).

Rededge bands have been widely studied for assessing crop performance and yield in various crops, including wheat (Horler et al., 1983; Pavuluri et al., 2015). Canopy reflectance in the red, rededge and near-infrared (NIR) wavelength range is influenced by two primary optical properties of canopies: chlorophyll absorption in the red region and multiple scattering effects on the NIR due to canopy structural properties. The red-edge region is more sensitive to chlorophyll content than to leaf area (Xie et al., 2018). Hence, the variability in LAI might have been bigger or more important for yield formation compared to the chlorophyll content in our specific panel. Moreover, the correlation of reflectance in the red-edge region with yield is known to change quickly with the exact wavelength measured (Pavuluri et al., 2015), making the selection of the exact wavelength difficult and leading to inconsistent results.

In contrast, visible bands (Blue, Green and Red) can be more sensitive to yield-related variations in chlorophyll, and biomass accumulation during the tillering and the stem elongation stage until the beginning of the booting stage. They are known to be correlated to a certain extent to both, chlorophyll concentration and LAI (Daughtry et al., 2000). Accordingly, our results showed that the RED, GREEN, and BLUE bands were among the most effective spectral features for yield prediction, exhibiting significant differences between high- and low-yielding varieties at almost all measurement dates in the first year of trial and during several in the second year. Their reflectances decrease during the transition from the stem elongation to the beginning of the booting stage when LAI and chlorophyll density is known to be maximal (Hinzman et al., 1984; Hinzman et al., 1986) due to the optical properties of chlorophyll. From heading until harvest the reflectances in these bands increase due to senescence when chlorophyll degradation takes place (Spano et al., 2003).

The NIR region is known to be sensitive to leaf area and especially ground cover (Korobov and Railyan, 1993), making it a useful band for predicting biomass and therefore yield. Our results indicate that the NIR band performed best during the stem elongation stage for yield prediction and at the booting stage when there were significant differences in (LAI) between the two yield groups. In the second season, the differences in the reflectances were significantly different only once before booting in the second season possibly due to a lack of differences in LAI between the two yield groups. This aligns with the findings by Korobov and Railyan (1993), who reported a higher correlation of NIR reflectance with dry matter and ground cover the during booting stage compared to later stages. Thus, normalizing the difference of the NIR and the REDEDGE reflectance in the form of the NDRE index, showed a good performance for chlorophyll estimation (Barnes et al., 2000).

Usually, VIs containing information from the rededge region of the spectrum are considered being more sensitive to chlorophyll absorption in dense canopies (Nguy-Robertson et al., 2012). It is expected that combining the highly LAI-sensitive NIR band with the rededge band that contains more information about leaf pigments in the canopy and therefore improves the performance of our yield prediction model at the flowering to early grain filling stages.

The performance of yield prediction and classification depends highly on the phenological stage of the crop. Our study found that the flowering stage and early grain filling stage allowed for the best predictions of yield and classification of varieties in winter wheat, which is consistent with the findings of several other studies (Hassan et al., 2019; Prey et al., 2020; Prey et al., 2022). From a physiological stand point of view, at the time around flowering the crop has to provide enough assimilates in order to maximize the number of fertile florets per spike, leading to a higher number of kernels per spike and finally a higher yield (Fischer, 1985). Therefore, estimating biomass and chlorophyll content at these stages is optimal for yield estimation. Unfortunately, the spectral signal often saturates at these stages making the estimation challenging, especially in high yielding years (Prey et al., 2020).

Early differences in biomass and LAI dynamics between wheat genotypes are well-documented (Pang et al., 2014; Grieder et al., 2015) and Raun et al. (2001) proposed to follow the biomass formation after dormancy for yield prediction. Marti et al. (2007) hypothesized that a high biomass at the end of the tillering stage detected by NDVI estimates the growth during the stem elongation stage which in turn is crucial for the number of kernels per area produced (González et al., 2005). Yield prediction at tillering therefore already yielded some success especially in the first year of the trial. From tillering to harvest, wheat is known to compensate for, as an example, a low stand count by altering the number of yield components (Holen et al., 2001), therefore predictions at these early stages are prone to changes later in the season, subsequently leading to prediction errors.

After flowering, the RMSE of the models for prediction increases again in the first year whereas they start to increase later during grain filling in the second year. The flowering stage is longer in the first compared to the second year of trial (232 GDD and 134.4 GDD). The reason for this difference is unclear, the end of flowering is generally difficult to rate, since the dry anthers tend to fall to the ground or are washed of by rain. Therefore, the end of flowering might have been later in the second year as well. Differences in stay-green characteristics can influence yield, especially by influencing the weight of single grains (Wu et al., 2012) which were higher in the high yielding varieties compare to the lower yielding ones. A further indication of this difference is the slightly advanced phenology of the lower yielding varieties in both years and the higher LAI measured in the first year of the trial.

The classification models showed a less clear trend during the season and fluctuations between dates are more severe than for the regression models. Garriga et al. (2017) did not find a difference in classification model performance between the anthesis and the grain filling stages. Successful classification in our case can be achieved directly by identifying high yielding plots or indirectly by identifying varieties. The algorithm learns either to recognize high yield or to classify a certain variety, regardless of their yield potential. Wheat variety classification for yield classification for breeding purposes other than by Garriga et al. (2017) is scarce. Works often focus on the classification of kernels (Porker et al., 2017; Khojastehnazhand and Roostaei, 2022), which is done under laboratory environments and therefore less prone to errors induced by the measurement conditions. Still, if the optimal feature at the optimal time point is selected, classification of low and high yielding varieties can be a promising tool for plant breeding applications (Garriga et al., 2017).

Our study found that single-band reflectance, such as the RED band, was as effective as or even more effective than vegetation indices (VIs) for predicting yield, especially in the first year of trial. The RED band is known to be related to leaf area index (LAI), although this relationship is often non-linear (Hinzman et al., 1984) and therefore requires non-linear methods such as RF to perform well for yield prediction. Pavuluri et al. (2015) found a saturation of RED reflectance when predicting yield, which can also be found in our prediction models. In contrast, VIs typically show good linear correlations with grain yield, with NDVI being widely used for yield prediction (Duan et al., 2017; Hassan et al., 2019). Furthermore, the VIs show a more consistent performance between dates compared to the single band reflectances. The NDRE was the best performing feature for yield assessment, which is in accordance with other studies (Prey et al., 2022) possibly due to its strong correlation to biomass (Argento et al., 2021). Many VIs have been screened by Prey et al. (2020) and few have been showing a consistent performance over the years, which makes a general selection difficult, similar to our study. Further, VIs narrow down the information that is accessible and Vatter et al. (2022) found good performances for yield prediction when using 11 wavebands from a multispectral camera that were fed to a deep learning model. This might be good strategy to obtain the optimal prediction model from multispectral cameras without any prior knowledge and the need for feature selection.

TFs are complex in their calculation and they offer a variety of possible ways of calculation, possible combinations with underlying rasters and ways to be calculated. Detailed information on how TFs are calculated is often lacking (Zheng et al., 2019; Wang et al., 2021; Zhang et al., 2021). Therefore, TFs still have to be examined in detail and their parameters optimized under different experimental conditions and scenarios of sensing data collection. We calculated TFs in a standardized way, but still found a high variability between dates. They are further known to be highly dependent on the GSD and therefore, the flight height (Zheng et al., 2019). In our study smoothing aided in enhancing the stability of the yield prediction models, although it did not improve their performances. A novel approach was presented by Herrero-Huerta et al. (2020) who calculated so-called canopy roughness directly on the point cloud from the structure from motion processing and showed its correlation to biomass. Often, the TFs are difficult to interpret and their link to yield relevant canopy traits is often unclear. Originally, the TFs were developed for classification (Haralick et al., 1973) and therefore classifying varieties corresponds more to their intended purpose than yield prediction. TFs are further often used in combination because there might be additional information (Wang et al., 2021; Liu et al., 2022), especially in the later stage, when they are not strongly correlated to single band reflectances and VIs anymore, as indicated by our results.

Models using individual dates showed generally a worse performance than models containing three adjacent dates. Furthermore, the performance between different dates is fluctuating strongly, even if the phenological stage is similar. Prey et al. (2022) found combining data from multiple dates to yield more improved predictions especially if the features used were performing poorly. In our case, the improvement for the yield prediction was in a similar range, regardless of the initial performance of the feature. In rice, Zhou et al. (2017) found that combining data from different growth stages by a multilinear regression model, can improve the estimation accuracy. For practical applications, finding the optimal date might be difficult and requires very close monitoring of the phenology, which can be very diverse among varieties as in our variety-testing panel. Therefore combining multiple dates might be especially performant, if features are used that show a high fluctuation between dates such as the TFs. The RF algorithm however is capable of dealing with different suitability of dates and therefore neglect the ones that do not perform well by attributing different importance to the features. The downside of the method is that, obtaining the additional measurements requires substantial work. Therefore, the number and time points measured should be considered carefully and be optimized in future studies.

The red-edge position and its shape is often used to estimate the stress status of field crops (Guyot et al., 1988; Boochs et al., 1990). However, it is obvious that the dynamics (time series) of the Red-edge band is difficult to interpret compared to the visible bands. During the early stages of tillering, the red-edge reflectance increased, possibly due to an increase of ground cover, whereas later it decreased again, when the canopy height increased during the SE stage. At the beginning of the heading stage, another increase in the red-edge reflectance could be observed, accompanied with the increase of reflectance in the visible bands. However, in contrast to other bands, the Red-edge reflectance decreases with the onset of senescence at the early grain-filling stage, possibly due to a reduction in chlorophyll and a shrinking canopy structure (Wang et al., 2022). However, fluctuation also occurs during the mentioned stable period running from the beginning of booting to the end of flowering. These fluctuations can be of various origins. For instance, the appearance of the canopy might change significantly due to the emergence of the spikes. Although this study was unable to exploit the entire shape of the red-edge reflectance, due to limitations in our multispectral camera having one band in the red-edge region, future work should further advance the understanding of the dynamics of red-edge reflectance and responsible canopy characteristics. Also, features should, in addition to their performance for yield prediction be assessed regarding their heritability (H²) since breeders are interested in knowing the genetic variation underlying a trait or in our case a spectral feature. Generally, this study shows that a trait time series followed by smoothing and a moving window allows for more stable predictions when also not better predictions.