UAV Remote Sensing for High-Throughput Phenotyping and for Yield Prediction of Miscanthus by Machine Learning Techniques

Abstract

Miscanthus holds a great potential in the frame of the bioeconomy, and yield prediction can help improve Miscanthus’ logistic supply chain. Breeding programs in several countries are attempting to produce high-yielding Miscanthus hybrids better adapted to different climates and end-uses. Multispectral images acquired from unmanned aerial vehicles (UAVs) in Italy and in the UK in 2021 and 2022 were used to investigate the feasibility of high-throughput phenotyping (HTP) of novel Miscanthus hybrids for yield prediction and crop traits estimation. An intercalibration procedure was performed using simulated data from the PROSAIL model to link vegetation indices (VIs) derived from two different multispectral sensors. The random forest algorithm estimated with good accuracy yield traits (light interception, plant height, green leaf biomass, and standing biomass) using 15 VIs time series, and predicted yield using peak descriptors derived from these VIs time series with root mean square error of 2.3 Mg DM ha⁻¹. The study demonstrates the potential of UAVs’ multispectral images in HTP applications and in yield prediction, providing important information needed to increase sustainable biomass production.

1. Introduction

Miscanthus is a high yielding perennial biomass crop. Yield is one of the most important traits of Miscanthus [1] and has been the primary focus of the research portfolio on Miscanthus in the last ten years [2,3,4]. Independent and collaborative efforts to breed high-yielding Miscanthus hybrids to produce sustainable biomass for the growing bio-based European market are ongoing in several countries [5,6,7]. In the EU-BBI demo-project GRACE, novel and highly upscalable seed-based Miscanthus hybrids are evaluated [5,7,8,9,10] in seven European countries. Most yield prediction to date has relied on crop growth models driven by climate and soil data with crop specific parameters [11,12,13]. For example, MISCANFOR is a crop growth model specifically developed to predict Miscanthus × giganteus yields in a wide range of environments [14]. It has been widely used and validated at the European [6] and national level [15] for Miscanthus and other perennial biomass crops [16], but new parameterization data to predict yield production of the novel Miscanthus hybrids [5] is required [14].

Yield trait screening and prediction using remote sensing with unmanned aerial vehicles (UAVs) can help both in breeding activities and in obtaining spatial and temporal information for optimizing Miscanthus biomass supply chain logistics, from field to facilities creating bioproducts or biopower [17,18]. Impollonia et al. [19] recently demonstrated the feasibility of moisture content estimation in Miscanthus hybrids using vegetation indices (VIs) derived from UAV-based remote sensing. Remote sensing approaches can also be used to (i) estimate yield-related traits for high-throughput phenotyping (HTP) [20,21,22], (ii) to calibrate crop growth models [23,24], and (iii) to predict the yield of many crops for commercial purposes [25]. To date, few studies have focused on the estimation and prediction of perennial biomass crops traits using remote sensing technologies from satellite or UAV [26,27,28], and only two on Miscanthus [18,19]. Crop traits such as the plant height [29], the fraction of absorbed photosynthetically active radiation (fAPAR) [30], and the aboveground biomass [31] can be estimated from the VIs in combination with machine learning (ML) algorithms. One of the most used ML algorithms in remote sensing analyses for crop traits estimation is random forest (RF) [32,33,34,35]. RF proved to be robust to outliers and noise, does not suffer from overfitting, and can manage a high training size [36]. The use of ML algorithms demonstrates great potential in crop yield prediction [37,38,39,40]. In particular, the use of a VIs’ time series helped to derive descriptors of land surface phenology (LSP, i.e., the spatial and temporal development of the vegetated land surface) [41,42,43,44,45] such as the start of season (SOS), the peak of growing season, the stay-green duration (onset of senescence), the end of the season (EOS), and growing season length [46]. Among the available descriptors of LSP, the peak of a VI is one of the most important descriptors for crop yield prediction, such as the peak of NDVI [47] and EVI2 [17] for grain yield and the peak of GNDVI for biomass yields of perennial grass [26]. Yield prediction is, in many studies, theoretical, as it requires to fit the whole seasonal curve for deriving the peak of the time series of VIs. This renders impossible the yield prediction before having obtained the data of the whole seasonal time series. Being able to perform yield prediction months before harvest with a partial time series should permit optimization of Miscanthus biomass supply chain logistics for practical applications.

In addition, the VIs’ values used for estimation and prediction of crop traits are influenced by many factors, such as sensor characteristics, atmospheric conditions during acquisition, viewing angle, field of view, and sun elevation [48]. In the context of crop phenotyping, where the field trials are often carried out in multi-locations and with different UAV sensors, these factors could have a relevant effect on the compatibility of VIs. Among different sensors characteristics, the full width at half maximum (FWHM) is the main factor that influences the comparability of VIs’ values among different sensors [49]. Indeed, due to the different spectral characteristics of the UAV multispectral sensors available on the market, differences among VIs derived from multiple UAV sensors for the same target can be found [48]. For this reason, there is a need to increase the interoperability of the sensors using equations that are able to overcome these differences through advanced linking procedures between the VIs of sensors [50,51]. The multi-sensor interoperability is an important topic in remote sensing science [52,53,54] when multi-location monitoring is conducted. Overcoming this problematic can be realised by intercalibrating the VIs obtained by different sensors and evaluating the correlation through linear regressions [55,56,57]. This intercalibration can be realised in two main ways: a direct approach where VIs are measured by sensors [58] and an indirect approach where VIs data are simulated by radiative transfer models [59]. A major limitation of the direct regression approach is that it is not transferable because it is site-specific [50]. On the contrary, the indirect approach permits retrieval of VIs’ data from radiative transfer models, such as the PROSAIL model [60,61]. The PROSAIL model is commonly used for the retrieval of biochemical and biophysical crop traits from remotely sensed images for vegetation studies in agricultural applications [62]. PROSAIL has also been used to assess the performances of different satellite sensors for multiple VIs [63], but a similar application is currently lacking in UAV science. Hence, the use of such an approach for VIs’ intercalibration in UAV science might allow for overcoming the problematic of multi-sensor interoperability across different sites.

In summary, the overall objectives of this study, based on UAV remote sensing, were: (i) to estimate crop traits (light interception, plant height, green leaf biomass, and standing biomass) for supporting breeding programs and for providing modelling parameters for Miscanthus, (ii) to predict yield to obtain spatial and temporal information for improving the logistics biomass supply chain of Miscanthus, and (iii) to explore the potential impact of the timeliness on the yield prediction, by evaluating the performance of the yield prediction model using the peak derived from a partial time series of VIs. To achieve these overarching objectives, UAV multispectral images and ground phenotypic data were collected at two locations within the same multi-environment trial: one in Italy and one in the UK. These data were analysed using: (i) the PROSAIL model to simulate crop spectral signatures in order to intercalibrate VIs of two different common multispectral sensors, (ii) the random forest (RF) algorithm to estimate crop traits using the VIs time series and to predict yield using the peak descriptor derived from the VIs time series, for the reliability it achieved in previous studies [17,26,47], and (iii) the generalized additive model (GAM) to derive the peak from a complete and partial time series of VIs.

2. Materials and Methods

2.1. Experimental Design

The field trials were conducted in two locations (Figure 1): PAC 1, located in San Bonico in the Italian province of Piacenza (NW Italy) (45°00′11.70″ N, 9°42′35.39″ E), and TWS 1, located in Trawscoed near Aberystwyth in Wales (UK) (52°24′59.8″ N, 4°04′02.6″ W). These sites are two of the seven plot scale (PS) trials conducted within the EU-BBI GRACE demo-project. In PAC 1, the climate is temperate with a mean annual precipitation of 792 mm, while the climate in TWS 1 is oceanic with a mean annual precipitation of 984 mm. The trials were established in April 2018 with 14 Miscanthus hybrids with n = 4 replicates for a total of n = 56 plots. The trials were planted with eight novel intraspecific M. sinensis × M. sinensis hybrids (M. sin × M. sin), five novel interspecific hybrids M. sinensis × M. sacchariflorus (M. sin × M. sac), and M. × giganteus as the control genotype (for more details see Impollonia et al. [19]).

2.2. Phenotypic and Yield Measurements

Phenotypic measurements were taken from the emergence of the crop in spring 2020 until the winter harvest in the early months of 2021. This season will thereafter be referred to as the 2020 growing season. These phenotypic measurements were carried out in the two locations and on four out of the fourteen hybrids in the trial: GRC 3 (M. sin × M. sin), GRC 14 (seeded) and GRC 15 (clonal) (M. sin × M. sac), and GRC 9 (standard clonal M. × giganteus). The measurements of this study were carried out during the third growing season. Five contiguous plants along a central row in each plot were used for monitoring during the growing season. The following list of crop traits were measured during the growing season: plant height (cm) and light interception (%) were measured weekly, and green leaf biomass (Mg DM ha⁻¹) and standing biomass (Mg DM ha⁻¹) were measured fortnightly; 172 and 145 data were collected for light interception, 240 and 204 for plant height, 232 and 316 for green leaf biomass, and 268 and 328 for standing biomass in PAC 1 and TWS 1, respectively. Plant height was measured from the soil to the height of the last ligule of the tallest stem using a graduated pole until the crop reached complete flowering or started to senesce in November. Light interception was measured with a lab-constructed 1 m “line ceptometer” with 10 photodiodes at 10 cm spacings, generating an electric current that is converted with simple circuitry to a voltage linearly proportional to the light intensity. Light intensity was measured above the canopy and at the base of each of the five selected plants. Light interception measurements were carried out from emergence until full canopy closure (around 95% of light is intercepted by the crop canopy) on a weekly basis and later at a lower frequency. Standing aboveground biomass was estimated on a monthly basis, starting after emergence in 2020 until harvest in winter 2021. Ten randomly sampled shoots per plot each fortnight (aka “serial cuts”) were related to the final quadrat yields at spring harvest and used to back calculate the seasonal dynamics of above ground biomass (Mg DM ha⁻¹) from spring emergence until final harvest in the following spring. Each 10 shoot serial cuts were separated into green leaf, brown leaf, and stem fractions. The fresh weight measured in the field and the dry weight measured after oven drying to a constant weight at 80 °C were used to calculate the mass and moisture contents of each fraction and were scaled to Mg DM ha⁻¹. The crop’s phenological stages were estimated using thermal time accumulation, following the method proposed on M. × giganteus by Tejera et al. [64]. Two main phenological stages were estimated: vegetative growth and senescence. The GDD accumulation started at plant emergence, and the accumulation of 1500 GDD was used as the threshold marking the difference between the vegetative growth and senescence. The final harvested yield was measured for all 14 hybrids, using a quadrat area of 6.6 m² (10 plants per plot planted at 1.5 plants m⁻² (M. × giganteus and M. sin × M. sac) or 20 plants per plots at 3 plants m⁻² (M. sin × M. sin)) and a cutting height of 10 cm. In each plot, the fresh weight of all plants in the quadrat was recorded, and a subsample of approximately five stems per plot was used to determine the moisture content and calculate the yield in Mg DM ha⁻¹. Plants were harvested on 02 February 2021 (Days of year (DOY): 33) at PAC 1 and on 08 March 2021 (DOY: 67) at TWS 1.

2.3. UAV Multispectral Data and Vegetation Indices

Unmanned aerial vehicle (UAV) multispectral data acquisition flights were performed approximately fortnightly from 24 April 2020–1 February 2021 at the PAC 1 (25 flights) site and from 9 June 2020–25 February 2021 at the TWS 1 (17 flights) site (Table S1).

Table 1. Unmanned aerial vehicles (UAVs) and multispectral cameras used to perform flights in the two locations, with their respective characteristics of spectral bands, central wavelength, and full width at half maximum (FWHM).

Location	UAV	Multispectral Camera Characteristics
		Model	Band	Centre (nm)	FWHM (nm)
PAC 1	DJI M210 RTK	MicaSense RedEdge-Mx	Blue	475	32
			Green	560	27
			Red	668	14
			Red edge	717	12
			Near-infrared	840	57
TWS 1	DJI M210	SlantRange 4P	Blue	470	100
			Green	550	100
			Red	650	40
			Red edge	710	20
			Near-infrared	850	100

details the specifications of the UAVs and the multispectral cameras used at the two sites. All the flights were performed between 11 a.m. and 3 p.m., with the flight altitude above ground level (AGL) fixed at 50 m and 40 m at PAC 1 and TWS 1, respectively. The forward and lateral overlap was set at 80% and 75% of the images, respectively. Light sensors mounted on top of the UAVs were used for the radiometric calibration of the images. In addition, in PAC 1, images of the reflectance panel provided by MicaSense were also taken for calibrating images, to reduce the effects of the day’s changing lighting conditions during the flight [65]. The radiometric calibration and orthomosaic generation were done using the Pix4D mapper (Pix4D, S.A., Lausanne, Switzerland). The 15 vegetation indices (VIs) were calculated as illustrated in

Table 2. List of the vegetation indices evaluated in random forest models for crop traits estimation and yield prediction of Miscanthus.

VIs	Equation	Reference
Datt1	$\frac{\mathrm{NIR}-\mathrm{RedEdge}}{\mathrm{NIR}+\mathrm{Red}}$	[66]
EVI2	$2.4\frac{\mathrm{NIR}-\mathrm{Red}}{1+\mathrm{NIR}+\mathrm{Red}}$	[67]
GNDVI	$\frac{\mathrm{NIR}-\mathrm{Green}}{\mathrm{NIR}+\mathrm{Green}}$	[68]
GOSAVI	$\frac{\mathrm{NIR}-\mathrm{Green}}{\mathrm{NIR}+\mathrm{Green} +0.16}$	[69]
greenWDRVI	$\frac{0.1\mathrm{NIR}-\mathrm{Green}}{0.1\mathrm{NIR}+\mathrm{Green}}+\frac{1-0.1}{1+0.1}$	[70]
MSAVI	$\frac{2\mathrm{NIR} +1-\sqrt{{\left(2\mathrm{NIR} +1\right)}^2-8\left(\mathrm{NIR} - \mathrm{Red}\right)}}{2}$	[71]
MTVI1	$1.2\left(1.2 \left(\mathrm{NIR}-\mathrm{Green}\right)-2.5\left(\mathrm{Red}-\mathrm{Green}\right)\right)$	[72]
MTVI2	$1.5\frac{1.2\left(\mathrm{NIR} - \mathrm{Green}\right)-2.5\left(\mathrm{Red} - \mathrm{Green}\right)}{\sqrt{{\left(2\mathrm{NIR} +1\right)}^2-6\mathrm{NIR} -5\sqrt{\mathrm{Red}}-0.5 } }$	[72]
NDRE	$\frac{\mathrm{NIR}-\mathrm{RedEdge}}{\mathrm{NIR}+\mathrm{RedEdge}}$	[73]
NDVI	$\frac{\mathrm{NIR}-\mathrm{Red}}{\mathrm{NIR}+\mathrm{Red}}$	[74]
OSAVI	$\left(1+0.16\right)\frac{\mathrm{NIR}-\mathrm{Red}}{\mathrm{NIR}+\mathrm{Red} +0.16}$	[75]
OSAVI2	$\left(1+0.16\right)\frac{\mathrm{NIR}-\mathrm{RedEdge}}{\mathrm{NIR}+\mathrm{RedEdge} +0.16}$	[76]
rededgeWDRVI	$\frac{0.1\mathrm{NIR}-\mathrm{RedEdge}}{0.1\mathrm{NIR}+\mathrm{RedEdge}}+\frac{1-0.1}{1+0.1}$	[70]
SAVI	$\left(1+0.5\right)\frac{\mathrm{NIR}-\mathrm{Red}}{\mathrm{NIR}+\mathrm{Red} +0.5}$	[77]
WDRVI	$\frac{0.1\mathrm{NIR}-\mathrm{Red}}{0.1\mathrm{NIR}+\mathrm{Red}}+\frac{1-0.1}{1+0.1}$	[70]

, using the orthomosaic. The mean of the VIs was extracted for each plot using the polygons of the experimental designs that were drafted in AutoCAD (Autodesk; San Rafael, CA, USA) and georeferenced with QGIS software (QGIS Development Team, 2021).

2.4. Using the PROSAIL Model to Intercalibrate VIs from Different Multispectral Sensors

The PROSAIL model was used to simulate crop spectral signatures to intercalibrate the VIs (Table 2), calculated from the two different multispectral sensors used in this study (Table 1). PROSAIL can simulate the canopy reflectance, between 400 and 2500 nm, by combining the PROSPECT and SAIL models. The PROSPECT model simulates the optical properties of the leaves using four input parameters: leaf structure parameter (N), leaf chlorophyll content (LCC), relative leaf equivalent water thickness (C_wr), and leaf dry matter content (C_m). The SAIL model simulates the bidirectional reflectance of a canopy using six input parameters: leaf area index (LAI), leaf inclination distribution (LIDF), hotspot parameter (hot), solar zenith angle (tts), observer zenith angle (tto), and relative azimuth angle (psi). The canopy and leaf parameters for Miscanthus were retrieved from available data in the literature [78,79]. The hsdar R package [80] was used to simulate the canopy reflectance of the PROSAIL model [60,61] using the function PROSAIL which uses the FORTRAN code of the PROSAIL model (Version 5B). The look-up table (LUT) generated included 430,080 parameter combinations, following the min-max ranges of input parameters summarized in

Table 3. Ranges of input parameters for the PROSAIL model for the generation of the LUT.

Parameter		Abbreviation	Unit	Values (Step)
Leaf	Structure parameter	N	Unitless	1–2 (1)
	Chlorophyll content	LCC	µg cm−2	10–80 (10)
	Relative equivalent water thickness	Cwr	%	20–80 (20)
	Dry matter content	Cm	g cm−2	0.01–0.025 (0.005)
Canopy	Leaf area index	LAI	m2 m−2	1–8 (1)
	Leaf inclination distribution	LIDF		Spherical
	Hotspot parameter	hot	m m−1	0.05–0.45 (0.2)
	Solar zenith angle	tts	deg	20–80 (10)
	Observer zenith angle	tto	deg	5–10 (5)
	Relative azimuth angle	psi	deg	180–220 (10)
	Structure parameter	N	Unitless	1–2 (1)

.

The present work considers all potential parameter combinations for LUT generation because the novel Miscanthus hybrids evaluated have not previously been studied in the literature, and because the crop monitoring was performed on the whole growing season (vegetative growth and senescence). The spectral reflectance simulated were resampled based on the UAV sensor characteristics (Table 1) and the 15 VIs used in this study were calculated. For each VI, regression analysis was conducted to intercalibrate the VI values calculated from the two multispectral sensors. Linear regressions were performed using the VIs and not the spectral bands, in order to (i) evaluate the different sensitivity of the VIs to sensor characteristics and (ii) identify which VIs need an intercalibration procedure between sensors. The final database of VIs was built by scaling the SlantRange sensor VIs data toward the MicaSense sensor.

2.5. Time Series of VIs and Peak Derivation

The 15 Vis, calculated from UAV multispectral images of the two sites and intercalibrated using the linear models derived by the PROSAIL simulation, were smoothed using a generalized additive model (GAM) to generate daily VIs time series. The GAM is a non-parametric regression model which allows non-linear fitting of the variables. GAM models were fitted in the R package “mgcv” [81]. GAM fitting of VIs allows for removal of the outliers and regularizes the time series [19,82,83]. The time series of VIs for each plot were fitted against the modified days of the year (DOY_M). The DOY_M was used to overcome the problem of having non-sequential DOY throughout the growing season, as this last season overlaps two different years. DOY_M were calculated as conventional DOY for the first year of the growing season (2020) and as DOY + 365 starting from the 1st of January for the year 2021.

The daily time series of the VIs were used to estimate the crop traits by linking the traits values measured in the field with the VIs’ values of the time series. To predict yield, the peak descriptor was chosen among several land surface phenology (LSP) descriptors, due to the reliability it achieved in previous studies [17,26,47]. The peak descriptor is defined as the maximum value and was derived from the GAM fitting of each VI time series. Two types of VIs time series were evaluated to derive the peak descriptor value: the complete and partial times series. The complete Vis time series were obtained by fitting the GAM models to the VIs’ data acquired throughout the entire growing season. The partial time series were obtained in two steps. Firstly, by fitting the GAM models to the VIs’ values, obtained from a reduced set of initial UAV flights over the crop. The first seven flights (175 DOY in PAC 1 and 266 DOY in TWS 1) were selected to form this reduced set of initial UAV flights, based on the physiology and phenology of Miscanthus, to cover most of the vegetative phase, the seventh flight being close to the peak of biomass accumulation. Secondly, the VIs’ data obtained from the following UAV flights were added one by one to the data fitted in the first step, up to the end of the season. The peak derived from the complete time series of VIs were used for the yield prediction modelling, and those from the partial time series of VIs were used to analyse the variation of peak values and the model performance in order to assess feasibility of “early season” yield prediction during Miscanthus crop growth.

The random forest (RF) algorithm was used to estimate the crop traits (Section 2.2) and to predict the yield of Miscanthus hybrids. The RF models were created using the caret R package [84]. Three steps in RF modelling were adopted. Firstly, RF models were created for (i) traits estimation, using the 15 VIs of the 4 Miscanthus hybrids (see Section 2.2) and phenological stages (encoded as 0 and 1 for the stages of vegetative growth and senescence, respectively) and (ii) yield prediction, using the peak values derived from the complete time series of the 15 investigated VIs of the 14 Miscanthus hybrids. Secondly, the importance of variables for each tested RF model was calculated by the dropout loss of the root mean square error (RMSE), and the variables with a median of RMSE loss greater than 0 were selected. Thirdly, the RF models used to estimate the variable of interest (i.e., crop traits or yield) were created using only the selected variables. The optimal size of the variable subset (“mtry”) for each model was obtained by a grid-searching method, using repeated k-fold cross-validation (ten-fold cross validation repeated five times). The training dataset was created using a stratified random sampling method [31]: data from both locations and genotypes were split into 2/3 of the dataset for training and 1/3 for testing, based on data distribution. The variable importance was calculated by the dropout loss of the RMSE (i.e., the increase of prediction in the RMSE [85]) using the DALEX package [86]. To demonstrate the uncertainty of importance estimation, the variable importance was calculated for 10 permutations [85]. The RF models’ performances were evaluated through the coefficient of determination (R²), the root mean square error, and the normalized root mean square error (NRMSE) [87]. The NRMSE was calculated by normalizing the RMSE using the mean of the measured values [19]. The RMSE and NRMSE metrics were also used to compare the performance of the yield prediction model using test datasets, with the peak derived from the complete and partial time series of VIs.

3. Results

3.1. PROSAIL Model for Intercalibration of VIs Derived from Different Multispectral Sensors

The outputs of the linear regressions analysis for each VI of interest (see Table 2), as calculated from PROSAIL model simulation considering the specific spectral response function of the two sensors (y-the MicaSense VIs and x-the SlantRange VIs), are reported in Figure 2. EVI2, MSAVI, and SAVI were the three VIs with the slope values closest to 1 and intercept values closest to 0 (Figure 2). The slope values of these VIs were 0.99 for EVI2, 0.99 for MSAVI, and 0.97 for SAVI, while the intercept was respectively 0.02 for EVI2, 0.02 for MSAVI, and 0.02 for SAVI. OSAVI had a similar relation (slope: 0.91) but demonstrated higher variability at VI values lower than 0.6. Datt1, NDRE, and OSAVI2 had a relationship between the two sensors with a slope close to 1 but they displayed a different intercept, indicating an offset of the SlantRange sensor due to the underestimation of the VI compared to MicaSense (Figure 2). A slope close to 1, associated with a high variability for the whole range of VI values, was observed for MTVI1 (slope: 1) and MTVI2 (0.97). GNDVI, GOSAVI, greenWDRVI, and NDVI illustrated, instead, the highest differences between the two sensors at the lowest values of VI. On the contrary, rededgeWDRVI demonstrated the highest differences at high VI values.

3.2. Importance of Variables in Machine Learning Models

The importance of the random forest (RF) models’ variables for crop trait estimations and for yield prediction are illustrated in Figure 3. The analysis of variables’ importance for the RF models was evaluated by drop-out loss of the RMSE for each variable, compared to the full model [85,86]. The phenological stage (“Stage”) was the most important variable for estimating plant height, green leaf biomass, and standing biomass (Figure 3). For crop trait estimations, the two most important VIs were NDVI and MTVI1 for light interception, rededgeWDRVI and NDVI for plant height, and greenWDRVI and GNDVI for green leaf biomass and standing biomass. For yield prediction using VIs peaks values, the most important VIs (with a median of RMSE loss greater than 0) were greenWDRVI, NDVI, WDRVI, GNDVI, and MTVI2 (Figure 3) and the peak values of these five VIs will therefore be used for RF modelling.

3.3. Machine Learning Model for Crop Traits Estimation

The crop trait distribution of contrasting Miscanthus hybrids, measured at two locations, is illustrated in Figure 4. The frequency distribution of the traits (light interception, plant height, green leaf biomass, and standing biomass), which values were used for training and testing the models, demonstrated that lower values were recorded in TWS 1 than in PAC 1 (Figure 4). For light interception, the range was from 5.2% to 100% in TWS 1 and PAC 1, and the mean was 58% in TWS 1 and 81.5% in PAC 1. The mean of plant height was 211 cm in PAC 1 and 147 cm in TWS 1, with the range from 28 cm to 344 cm and from 22 cm to 280 cm in PAC 1 and TWS 1, respectively. The range and the mean of green leaf biomass were from 0.14 Mg DM ha⁻¹ to 14.5 Mg DM ha⁻¹ and 5 Mg DM ha⁻¹ in PAC 1 and from 0.05 Mg DM ha⁻¹ to 6.3 Mg DM ha⁻¹ and 1.5 Mg DM ha⁻¹ in TWS 1. For standing biomass, the range was from 0.5 Mg DM ha⁻¹ to 46.4 Mg DM ha⁻¹ in PAC 1 and from 0.5 Mg DM ha⁻¹ to 21.1 Mg DM ha⁻¹ in TWS 1, and the mean was 8 Mg DM ha⁻¹ in TWS 1 and 18.9 Mg DM ha⁻¹ in PAC 1.

Overall, the RF model estimated crop traits with good model performance (Figure 5 and Figure 6). Among the crop traits, light interception was estimated with the greatest accuracy (NRMSE of 10.9%; Figure 5). High model accuracy was also achieved for the estimation of plant height (21.8% NRMSE), while the lowest model accuracies were observed for green leaf biomass and standing biomass (42.2% and 45.3% of NRMSE, respectively; Figure 5). For these last parameters, the RF model demonstrated good accuracy from low to intermediate values, while above values of 5 Mg DM ha⁻¹ of green leaf biomass and 20 Mg DM ha⁻¹ of standing biomass the model performances dropped. The NRMSE performance metrics for each location and hybrid are reported in Figure 6. Generally, no relevant differences were observed between the two locations (Figure 6). In PAC 1, the GRC 3 hybrid demonstrated the worst performance for green leaf biomass, for standing biomass, and for plant height, while in TWS 1, the GRC 14 hybrid demonstrated the worst performance for all the crop traits considered, except for plant height.

3.4. Machine Learning Model for Yield Prediction

The frequency distribution of yield measured for the 14 Miscanthus hybrids in PAC 1 and TWS 1 is provided in Figure 7. On average, the yield of Miscanthus was higher in PAC 1 than in TWS 1. In PAC 1, the highest yield was recorded by M. sin × M. sac (a mean of 18.3 Mg DM ha⁻¹). M. sin × M. sin productivity averaged 11.3 Mg DM ha⁻¹, while M. giganteus was less productive (mean of 9.6 Mg DM ha⁻¹). In TWS 1, the highest yield was recorded by M. sin × M. sin (mean of 9.4 Mg DM ha⁻¹). M. sin × M. sac productivity averaged 8.2 Mg DM ha⁻¹ while M. giganteus was less productive (mean of 6.6 Mg DM ha⁻¹). The RF model, trained and tested with the yield values reported in Figure 7, enabled a reliable prediction of Miscanthus yield for all hybrids, using the peak derived from the complete time series of VIs (Figure 8). The RF model obtained a RMSE value of 2.3 Mg DM ha⁻¹ and a NRMSE value of 19.7% (Figure 8a). In PAC 1, M. sin × M. sac demonstrated the lowest NRMSE value, while M. sin × M. sin demonstrated the highest NRMSE value. On the contrary, in TWS 1, M. sin × M. sin demonstrated the lowest NRMSE value while M. sin × M. sac demonstrated the highest NRMSE value (Figure 8b). The modified days of the year (DOY_M) of the peak of the Vis (greenWDRVI, NDVI, WDRVI, GNDVI, and MTVI2) are reported in Figure 8c. On average, the VIs reached the peak earlier in PAC 1 (172 DOY_M–20 June) than in TWS 1 (DOY_M 232–19 August). In PAC 1, all hybrids reached the peak at the same time, while in TWS 1, all M. sin × M. sin hybrids were the earliest to reach the peak (228 DOY_M-15th August), while M. sin × M. sac hybrids and M. × giganteus reached the peak along a wide timespan ranging from the end of September until mid-November (Figure 8c).

3.5. Time Series of VIs and Yield Prediction Analysis

The complete time series of the five VIs identified as the most important for yield prediction (see Section 3.2 and Figure 3) are reported in Figure 9a for PAC 1 and Figure 10a for TWS 1. In PAC 1, all VIs values recorded throughout the growth of Miscanthus were the highest for M. sin × M. sac and the lowest for M. giganteus (Figure 9a), following the same order of the mean yield measured in the field (Figure 7). In TWS 1, a similar time series of all the Vis were recorded for the M. sin × M. sac and M. giganteus (Figure 10a). In particular, the peaks of M. sin × M. sac and of M. giganteus were reached later than that of M. sin × M. sin, as illustrated in Figure 8c. The variation throughout the season of the peak, derived by fitting the VIs via the generalized additive model (GAM), are displayed in Figure 9b for PAC 1 and Figure 10b for TWS 1. In PAC 1, the difference between the value of the peak derived from the complete time series of VIs and the value of the peak derived from the partial time series of VIs is the lowest (close to zero) after the 302 DOY_M (end-October). Before this date, the peaks’ differences are lower for M. sin × M. sin than M. sin × M. sac and M. giganteus (Figure 9b). In TWS 1, the difference between the value of the peak derived from the complete time series of VIs and the value of the peak derived from the partial time series of VIs is lowest after 331 DOY_M (end-November). As for PAC 1, the peaks differences are lower for M. sin × M. sin than M. giganteus and M. sin × M. sac (Figure 10b) in TWS 1 before 331 DOY_M. The timeline of the performance of the RF model tested with the VIs peak from partial time series is reported in Figure 11. In PAC 1, the NRMSE decreased until 302 DOY_M for M. giganteus and M. sin × M. sac while for M. sin × M. sin, it remained stable for all UAV flights performed from 175 DOY_M onward. In TWS 1, no relevant differences in NRMSE were observed between the UAV flights performed from 266 DOY_M onward.

4. Discussion

The use of UAV-based remote sensing provides a great potential for high-throughput phenotyping (HTP) at the plot scale, with applications in both breeding and in estimating the quality and quantity of the biomass needed for optimising downstream management of biomass fluxes. In this study, vegetation indices (VIs) and their peak were derived from unmanned aerial vehicles’ (UAVs) mounted multispectral sensors to estimate crop traits (light interception, plant height, green leaf biomass, and standing biomass), and to predict the final harvestable yield of novel Miscanthus hybrids and common M. × giganteus grown at two sites (Italy and UK).

4.1. The Importance of VIs Intercalibration Procedure for Multi-Sensor Interoperability

Intercalibrating the VIs of multi-sensors is relevant for remote sensing crop monitoring [54], particularly when the objective is to build models to estimate crop traits or to predict yield, and when sensors with different spectral characteristics are used. Indeed, the models might not reach the same accuracy if the VIs are calculated with sensors with different spectral band characteristics. The spectral signatures simulated from the PROSAIL model were used in this study to intercalibrate Vis, calculated from two common multispectral cameras (MicaSense RedEdge-MX and SlantRange 4P). Ideally, the VIs from any sensor can be then intercalibrated against the VIs of a selected reference sensor. This intercalibration approach is commonly adopted for VIs obtained from different satellites, [59] but, to our knowledge, this is the first time that such an approach is applied to UAV multispectral sensors. The PROSAIL model was used to simulate the canopy reflectance based on the value of specific Miscanthus traits, such as LAI, chlorophyll content, dry matter content, relative equivalent water thickness, leaf inclination distribution, and two site-specific values: solar zenith angle and relative azimuth angle. This approach was applied to the values of 15 VIs. EVI2, MSAVI, and SAVI were the three VIs with the closest 1:1 relation, which indicates that these VIs have a very low sensitivity to the spectral characteristics of the two sensors used in the present work. Similar results with no need for intercalibration for the same VIs were reported by Li et al. [88], comparing ETM+ and OLI satellite imageries. However, the other VIs evaluated in this study (Datt1, GNDVI, GOSAVI, greenWDRVI, MTVI1, MTVI2, NDRE, NDVI, OSAVI, OSAVI2, rededgeWDRVI, and WDRVI) demonstrated a higher sensitivity to sensor characteristics, underlining the importance of intercalibrating VIs for multi-sensor interoperability. In particular, the VIs based on green and red-edge bands demonstrated a higher variation than the VIs based on the red band, at all ranges of values. These differences are explained by the differences in the spectral characteristics of MicaSense and SlantRange in the green and red-edge bands. In particular, the SlantRange sensor has a broader green FWHM (100 nm) and a different central wavelength (710 nm) of the red edge compared to the MicaSense sensor (green FWHM: 27 nm and red edge central wavelength: 717 nm). This difference was already reported to cause considerable signal differences in other studies [89,90].

This study highlighted the importance of intercalibrating different multispectral sensors to increase interoperability in remote sensing. Kim et al. [91] and Villaescusa-Nadal et al. [92] reported that the use of linear regression equations to intercalibrate multi-sensors contributes to significantly correct (up to 50%) the effects of different spectral characteristics on VIs. However, it would be interesting to validate this procedure by flying simultaneously on the same field with two UAV sensors, comparing the values of the VIs with and without the intercalibration procedure, in order to evaluate the improvement in terms of VIs compatibility and multi-sensor interoperability. In fact, even if the spectral characteristics of the multispectral sensor are the factor that influences the compatibility of the VIs of different sensors the most [49], other factors cause differences of Vis, such as the atmospheric conditions during acquisition [48]. In the UAV images acquisitions, the changing light and meteorological conditions during the flights can affect the quality of the spectral data [65]. Therefore, this procedure is limited by other factors that cannot be considered by applying a simulated regression coefficient. However, using equations that are able to intercalibrate VIs derived by multiple multispectral sensors can reduce differences in Vis, improving crop monitoring and modelling for estimation of crop traits and prediction of yield.

4.2. Estimating Miscanthus Traits with Machine Learning

This study estimated Miscanthus traits using the random forest (RF) machine learning algorithm. The RF model was trained with the data collected on three novel seed-based Miscanthus hybrids and the common rhizome-based genotype M. × giganteus, at two contrasting locations (North-West Italy and Mid-West Wales). The RF algorithm, using 15 common VIs, successfully estimated crop traits, solving the non-linear responses between VIs and crop traits observed by Li et al. [28] for other perennial crops. The estimation of the crop traits from the time series of VIs acquired by UAV-based remote sensing can generate more data useful to calibrate existing Miscanthus crop models and to re-fine these models for novel Miscanthus hybrids in contrasting environments. The crop trait estimated with the greatest accuracy was light interception, which demonstrated a RMSE of 8.4%, the accuracy being especially good at high values of light interception (Figure 5). This result is in agreement with Guillen-Climent et al. [93], who found that the fraction of intercepted photosynthetically active radiation (fIPAR) was successfully estimated by a ML algorithm. Upreti et al. [30] found similar values of NRMSE (12.06% ) using the RF tree bagger approach for estimating the fraction of absorbed photosynthetically active radiation (fAPAR) of durum wheat. Good model accuracy was also achieved for plant height estimation (RMSE = 42 cm and NRMSE = 21.8%; Figure 5). A similar RMSE value (41 cm) was found by Han et al. [31] for plant height estimation of maize using the crop surface model, and by Tao et al. [94] (NRMSE = 21.2%) in the estimation of plant height of winter wheat using UAV hyperspectral images. The worst model accuracy was found for the green leaf biomass and standing biomass, with 1.3 Mg DM ha⁻¹ and 5.8 Mg DM ha⁻¹ of RMSE, respectively (Figure 5). However, for these parameters, the model demonstrated good accuracy from low to intermediate values, while, above values of 5 Mg DM ha⁻¹ of green leaf biomass and 20 Mg DM ha⁻¹ of standing biomass, the model performance dropped. The model could be affected by errors in the estimation in these intervals due to the fewer data used to train the model [95]; in fact, most of them were collected in PAC 1 (Figure 4). The different levels of accuracy of the models in estimating light interception and plant height compared to green leaf biomass and standing biomass could also be related to the period in which the measurements were taken and to the response of VIs during senescence. Field measurements of light interception and plant height were carried out in each environment from emergence until they peaked, which explains their good model accuracy. Estimation of crop traits with data coming from single UAV flights across the growing seasons is particularly affected after biomass reaches its maximum value in autumn. After biomass peak, with the start of the senescence period, the values of VIs start to decrease (Figure 9a and Figure 10a) [96], while Miscanthus green leaf biomass and standing biomass values remained stable or slightly decreased during this period. The difference of rate of decrease between the VIs and crop traits during senescence is a key aspect to consider in remote sensing estimation of crop traits. The importance of the senescence stage in the crop traits estimation is confirmed by the results on the variables’ importance (Figure 3). Indeed, the phenological variable “Stage” was the most important variable in the estimation of plant height, green leaf biomass, and standing biomass. The plant height of the GRC 3 (a M. sin × M. sin planted at high density) were poorly estimated, and this could be due to its canopy architecture and flowering time. The M. sin × M. sin hybrid has many distinguishable stems flowering (where plant height is measured), but its leaves are particularly curved and attached along the stem at a lower height than M. sin × M. sac. In addition, this genotype was transplanted at higher densities (3 plants m⁻²) and flowered earlier (end of August) compared to other genotypes, which flowered in early autumn. This more “prostrate” canopy architecture (with a higher stem segment between inflorescences and bent leaves) introduced noise in the plant height estimation from UAVs [97], since most of the reflectance comes from bent leaf mass. This noise caused by changes in plant architecture and the onset of flower can be seen in the NRMSE values of the RF models at the PAC 1 site (Figure 6). In fact, the earlier a genotype with prostrate architecture flowers, the worse is the estimation of plant height and biomass from an UAV.

4.3. Yield Prediction Using Machine Learning and Peak of VIs

The random forest (RF) trained with the peak derived from complete time series of five VIs acquired by UAVs was able to predict the yield of the 14 Miscanthus hybrids. The peak value of greenWDRVI, NDVI, WDRVI, GNDVI, and MTVI2 resulted the most important features for RF modelling, as derived from previous analysis (see Section 3.2, Figure 3). The RF model accurately predicted the yield, with 2.3 Mg DM ha⁻¹ of RMSE and 19.7% NRMSE (Figure 8a). The peak for Miscanthus hybrids occurs on average in mid-summer and early autumn in southern/warm (Italy) and northern/cold (UK) locations, respectively. The importance of the peak as land surface phenology (LSP) descriptors for yield prediction was already reported by Prasad et al. [98], who found that the peak had the highest correlation with cotton yield prediction compared to other LSP descriptors. Similar results were reported by Montazeaud et al. [47], who found a high correlation between the peak of NDVI and the yield, and by Liu et al. [17], who found that the EVI2 peak was a good predictor of grain yield. Among the VIs used in the RF model, the peak of the VIs based on the green band as GNDVI and greenWDRVI (Figure 3) were the most important variables for predicting not only Miscanthus yield at harvest, but also in estimating standing biomass and green leaf biomass during the growing season (Figure 4). Similar results for GNDVI were found in switchgrass and other warm-season perennial grasses [26]. In order to assess the capability of the model to predict the yield months before harvest (i.e., using only early season UAV acquisition and not waiting to perform UAV monitoring during the entire crop season), the RF model was calculated using the peak derived from the partial time series of Vis, and performance analysed. In the UK, the RF model accurately predicted the yield five months before harvest for all Miscanthus hybrids. In Italy, the yield of the M. sin × M. sin hybrids can be predicted with good accuracy seven months before harvest while M. giganteus and M. sin × M. sac hybrids required more time, as a good accuracy was obtained 3–4 months before harvest. This capability of the model to predict yield months before harvest is possible if no damage to plants later in the season occurs, thus considering a steady development process. Given the good accuracy achieved by the RF model even when peak value is assessed on a limited number of data (see Section 3.5, Figure 11), this approach can be considered as a suitable method to predict yield several months before harvest. However, it is possible that no time series analyses (i.e., a complete series of dense observation from sowing to harvest) are required, and the 5 VIs values could be used directly as inputs for yield prediction modelling, if the UAV data is acquired in the range of expected peak occurrence. This opportunity was also reported for switchgrass [26] using multispectral images acquired from multiple UAV flights, underlying the possibility of reducing the number of multispectral observations (e.g., 25 UAV flights were performed in PAC 1) to 2–3 UAV flights in proximity of the VIs’ peak period. In this context, the identification of VIs’ peak period, which is generally related to crop physiology and phenology (reaching the maximum production), is important in terms of cost (reduction of the number of observations) and of operationality of the yield prediction model. This model operability, intended as a capability of the model to accurately predict the yield some months before the harvest, is extremely relevant for improving the logistics of the biomass supply chain of Miscanthus and for supporting the improvement of crop modelling with remote sensing data.

5. Conclusions

This study demonstrated that vegetation indices (VIs) derived from unmanned aerial vehicle (UAV) multispectral images acquired in Italy and the UK can be successfully used in random forest (RF) machine learning (ML) algorithms to estimate the light interception, plant height, green leaf biomass, and standing biomass, and to predict the yield of novel Miscanthus hybrids using the peak derived from the VIs’ time series. This study evaluated the timeline of the performance of the model using peaks derived from the partial VIs’ time series, and the RF model demonstrated a good capability to predict the yield months before the harvest, both in Italy and in the UK, by using a limited number of UAV observations. In particular, the results suggested that the VIs’ values acquired during the peak period (without using the complete or partial time series) could be used directly for yield prediction, increasing the model operability. Yield prediction and high-throughput phenotyping, based on ML algorithms and on UAV remote sensing, can improve the logistics of the biomass supply chain, supporting breeding programs, and improving the crop modelling of novel Miscanthus hybrids. UAV platforms are suitable tools for HTP applications, as they enable the monitoring of small plots or field scale trials with numerous genotypes, due to their ability to capture high-resolution images. However, the satellite platforms are more suited for yield prediction, as they can collect data of many fields simultaneously and can develop applications to predict commercial yield at regional and national scales. In addition, this study reported for the first time a methodology to overcome the issue of multi-sensor interoperability among UAV multispectral sensors. The use of intercalibrating equations derived from the PROSAIL model proved to be a powerful tool to intercalibrate VIs from multi-sensors with different spectral characteristics. Although this intercalibration procedure is relevant for the upscaling of models from experimental plots to field by intercalibrating the UAV with satellites sensor characteristics, it is limited because it only considers the spectral sensor characteristics and no other factors, such as light and meteorological conditions during the flights, which may affect the quality of the spectral data and cannot be considered by applying a simulated regression coefficient.