Assessment of Ensemble Learning to Predict Wheat Grain Yield Based on UAV-Multispectral Reflectance

Abstract

Grain yield is increasingly affected by climate factors such as drought and heat. To develop resilient and high-yielding cultivars, high-throughput phenotyping (HTP) techniques are essential for precise decisions in wheat breeding. The ability of unmanned aerial vehicle (UAV)-based multispectral imaging and ensemble learning methods to increase the accuracy of grain yield prediction in practical breeding work is evaluated in this study. For this, 211 winter wheat genotypes were planted under full and limited irrigation treatments, and multispectral data were collected at heading, flowering, early grain filling (EGF), and mid-grain filling (MGF) stages. Twenty multispectral vegetation indices (VIs) were estimated, and VIs with heritability greater than 0.5 were selected to evaluate the models across the growth stages under both irrigation treatments. A framework for ensemble learning was developed by combining multiple base models such as random forest (RF), support vector machine (SVM), Gaussian process (GP), and ridge regression (RR). The R² values between VIs and grain yield for individual base models were ranged from 0.468 to 0.580 and 0.537 to 0.598 for grain yield prediction in full and limited irrigation treatments across growth stages, respectively. The prediction results of ensemble models were ranged from 0.491 to 0.616 and 0.560 to 0.616 under full and limited irrigation treatments respectively, and were higher than that of the corresponding base learners. Moreover, the grain yield prediction results were observed high at mid grain filling stage under both full (R² = 0.625) and limited (R² = 0.628) irrigation treatments through ensemble learning based stacking of four base learners. Further improvements in ensemble learning models can accelerate the use of UAV-based multispectral data for accurate predictions of complex traits like grain yield in wheat.

1. Introduction

Bread wheat is an important food crop that meets the energy needs of more than 1/4th of the world population [1]. Sustainable wheat production is largely determined by the adaptivity of cultivars to various stress environments. Drought and heat stress are major detrimental factors that can cause up to 10% of yield reduction in wheat [2,3]. Drought and heat stress causes early wilting and high rate of senescence by breaking down the chloroplast tissues which lead to yield losses [3]. Due to frequent fluctuations in these environmental stresses, the rate of genetic improvement is reported less than 1% in wheat. This is the major challenge for wheat breeders to develop elite cultivars for future food security. Wheat breeding programs usually perform early field assessments for candidate selection under different growth conditions [4], and secondary traits such as biomass and leaf area index (LAI) are exploited to predict grain yield before harvest [5]. These early assessments greatly help the breeders to make decisions and shorten the evaluation time [6]. However, field-based phenotyping is considered a bottleneck in wheat breeding because it is time-consuming, destructive, and has a high error probability [7,8]. In recent years, the advances in high-throughput phenotyping platforms (HTPP) such as an unmanned aerial vehicle (UAV) carrying multispectral sensors have provided a non-destructive and rapid approach to collect data from multiple sites at low cost [9,10,11,12,13].

The UAV-based multispectral information is measured mainly in the form of vegetation indices (VIs), such as normalized difference vegetation index (NDVI), green normalized difference vegetation index (GNDVI), and ratio vegetation index (RVI). This information can be used to detect the biomass, leaf area index, and chlorophyll level of plants [14,15,16]. Previously, these VIs have been taken as secondary traits to predict grain yield with high heritability [5,14]. In addition, it has been reported that the combination of different VIs to assess the plant physiological traits could contribute to a higher prediction accuracy compared with a single VI [17,18,19,20]. However, it is difficult for traditional models such as linear regression to use multiple traits for grain yield prediction. Recently, some machine-learning algorithms that can perform flexible nonlinear mappings between a large number of VIs have been introduced to predict crop grain yield [21,22,23]. Previously, random forest (RF) [24], support vector machine (SVM) [25], Gaussian process (GP) [26], and ridge regression (RR) [27] were evaluated by combining different spectral traits. These models have greatly improved the prediction accuracy for many traits in various plants [28,29,30,31]. For example, RF was successfully applied to evaluate chlorophyll [32], biomass [28,29], and LAI [33] in winter wheat. While SVM and RR showed high prediction accuracy and robustness for predicting wheat and soybean yield [11,34]. The GP-based machine learning models were exploited to estimate the LAI of forest [35] and nitrogen status in wheat [36].

Although several machine learning methods have achieved good prediction accuracy for several traits, the integration of various algorithms to obtain higher model performance is not explored. Some of the above models are traits- or species-specific, and the integration of different algorithms can improve the ability to use various traits at the same time. This helps the machine learning models to perform large-scale data analysis and predictions. Recently, the ensemble learning method that integrates multiple base learners is attracting great attention [37,38,39]. Stacking regression is an ensemble method that is proposed by Wolpert in 1992 [40] and statistically principled by Breiman in 1995 [41]. It has been applied to map composite [42], forest coverage [43], and PM2.5 monitoring [44]. In addition, in the field of plant phenotyping, stacking regression was combined with hyperspectral data to evaluate the photosynthetic capacity in tobacco and alfalfa yield [37]. To the best of our knowledge, the stacking regression method has not been applied to predict grain yield in wheat.

The main objectives of this study are as follows: (1) The ability of the UAV-based multispectral imagery for grain yield prediction is evaluated; (2) the effect of the ensemble learning method for improving grain yield prediction accuracy is investigated; and (3) the effect of base learners on developing accurate ensemble learning models is analyzed.

2. Materials and Methods

2.1. Plant Materials and Field Trials

Based on a number of UAV-based vegetation indices, a panel of 211 elite winter wheat genotypes was exploited to evaluate the ensemble learning models for grain yield prediction, where 186 accessions were collected from the Yellow and Huai Valleys Winter Wheat Zone (YHVWWZ) of China and 25 from five other countries.

The experiment was conducted at Xinxiang (35°18′N, 113°52′E), Henan province of China during the 2019–2020 cropping season (Figure 1). The genotypes were planted under two irrigation treatments (i.e., full irrigation and limited irrigation) in randomized complete blocks (RCBD) with two replications. A total of 844 plots were phenotype under both irrigation treatments using UAV platform. Each plot consisted of six rows with a length of 3 m, a width of 1.4 m, and inter-row spacing of 0.2 m. The trials were planted on 26 October 2019 with a seeding rate of 270 seedings/m² and harvested on 2 June 2020. Both irrigation treatments irrigated equal water at the tillering stages, while the full irrigation treatment also flooded in heading and early grain filling stages with 2250–2700 m³ ha⁻¹ of water. The fertilizer and management of both treatments were optimized equally. In addition, the seasonal precipitation was 110.2 mm, and harvesting was done at full maturity using a combine harvester.

2.2. UAV Platform and Flight Mission

A RedEge MX multispectral sensor (Micasense Parrot, Seattle, WA, USA) (https://wwwy.micasense.com/parrotsequoia) (Accessed: 30 November 2020) was mounted on a DJI M210 (SZ DJI Technology Co., Shenzhen, China) for multispectral imagery. The built-in GPS of the Camera provided the position and orientation of the images for data geo-referencing. The RedEdge MX 4.0 can simultaneously collect images of five different bands with GPS information at specific intervals and with the same resolution. The resolution of blue, green, red, red-edge, and near-infrared bands was 32 nm, 27 nm, 14 nm, 12 nm, and 57 nm (half maximum bandwidth), respectively. The detailed information of the band wavelengths is listed in

Table 1. Specification of the multispectral camera used in this study.

Band	Bandwidth	Wavelength	Definition	Image Resolution
Blue	475	32	1.4 mp	1280 × 960
Green	560	27	1.4 mp	1280 × 960
Red	668	14	1.4 mp	1280 × 960
Red-edge	717	12	1.4 mp	1280 × 960
Near infrared	842	57	1.4 mp	1280 × 960

. Flights were conducted in clear and cloudless weather conditions between 11:00 a.m. to 2:00 p.m. The UAV flew over the trial area in the fully automatic flight mode set by the DJI GS Pro software (https://www.dji.com/cn/ground-station-pro) (Accessed: 1 December 2020). Images for radiometric calibration of the RedEge MX were captured on the flat ground before and after each flight through a calibrated reflectance panel provided by Micasense. To obtain images with high resolution, the forward and side image overlapping was set to 85% and 80%, respectively. The flight altitude was maintained at 40 m for heading and flowering stages, and at 30 m for early and mid grain filling stages. The ground resolution was around 3 cm for the images taken at the altitude of 40 m and 2.5 cm for the altitude of 30 m. The detailed information of the UAV-based data acquisition is presented in

Table 2. Flight details of the unmanned aerial vehicle imagery system.

Growth Stage	Zadok’s Stage	Flight Altitude (m)	Snap Shoot Interval (s)	Ground Resolution (cm)
Heading	ZS-56	40	1.5	3.0
Flowering	ZS-65	40	1.5	3.0
Early grain filling	ZS-73	30	1.5	2.5
Mid-grain filling	ZS-85	30	1.5	2.5

ZS, Zadok’s stage.

.

In total, around 2000 to 2500 multispectral images were captured from each flight to make mosaic images. Pix4D mapper (Version 1.4, Pix4d, Lausanne, Switzerland) was exploited for orthomosaic generation after each aerial imagery of the experimental field. Meanwhile, plot segmentation for each image was done through QGIS 3.1.0 (https://www.qgis.org/) (Accessed: 30 November 2020). The polygon shapes with a specific ID were designed as masks to extract the spectral data of each plot. The spectral information was extracted using the computer vision algorithms provided by the ENVI software [5].

2.3. Vegetation Indices and Ground Data

As listed in

Table 3. Introduction of the multispectral vegetation indices.

Vegetation Index	Full Name	Equation	Reference
NDVI	Normalized Difference Vegetation Index	(NIR − R)/(NIR + R)	[45]
SAVI	Soil-Adjusted Vegetation Index	(NIR − R)/(NIR + R + 0.5) × 1.5	[46]
OSAVI	Optimized Soil-Adjusted Vegetation Index	(NIR − R)/(NIR + R + 1.6) × 1.16	[47]
NRI	Nitrogen Reflectance Index	(G − R)/(G + R)	[48]
GDNVI	Green Normalized Difference Vegetation Index	(NIR − G)/(NIR + G)	[49]
SIPI	Structure Insensitive Pigment Index	(NIR − B)/(NIR + B)	[50]
PSRI	Plant Senescence Reflectance Index	(R − B)/NIR	[51]
CRI	Carotenoid Reflectance Index	1/G + 1/NIR	[52]
EVI	Enhanced Vegetation Index	2.5 × (NIR − R)/(1 + NIR + 6 × R − 7.5 × B)	[53]
MSR	Modified Simple Ratio Index	((NIR/R) − 1)/((NIR/R) +1) × 0.5	[54]
NLI	Nonlinear Vegetation Index	(NIR × NIR − R)/(NIR × NIR + R)	[55]
RDVI	Re-normalized Difference Vegetation Index	(NIR − R)/(NIR + R) × 0.5	[56]
TVI	Transformational Vegetation Index	(NDVI + 0.5)0.5	[57]
MTVI	Modified Triangular Vegetation Index	1.5 × [1.2 × (NIR − G) − 2.5 × (R − G)]/[(2 × (NIR − G) − 6 × NIR + 5 × R0.5)0.5–0.5]	[58]
NDRE	Red edge Normalized Difference Vegetation Index	(NIR − REG)/(NIR + REG)	[59]
DVIREG	Red-edge Difference Vegetation Index	NIR − REG	[60]
OSAVIREG	Red-edge optimized Soil-Adjusted Vegetation Index	(NIR − REG)/(NIR + REG + 1.6) × 1.16	[60]
RDVIREG	Red-edge Re-normalized Difference Vegetation Index	(NIR − REG)/(NIR + REG)0.5	[60]
MSRREG	Red edge modified Simple Ratio Index	((NIR/REG) − 1)/((NIR/REG) +1)0.5	[60]
MTCI	MERIS Terrestrial Chlorophyll Index	(NIR − REG)/(REG − R)	[61]

, twenty vegetation indices were calculated as secondary traits to predict the grain yield. Most of the VIs contain near-infrared and red bands that have been used to quantify the biomass, pigment content, etc. in many crops. Grain yield was estimated for each plot after it was harvested at full maturity.

2.4. Stacking Regression Models for Ensemble Learning

Ensemble learning is a machine learning paradigm where multiple base learners are integrated to solve regression or classification problems. Stacking regression is an ensemble learning model proposed by Wolpert (1992) [40] to blend the predictors and improve prediction accuracy. It is commonly used to generate ensembles of heterogeneous predictors. In this study, four regression models including random forest (RF) [24], support vector machine (SVM) [25], Gaussian process (GP) [26], and ridge regression (RR) [27] were integrated for stacking regression-based ensemble learning. The R package “caret” (version 6.0−86) in R 4.0.2 (https://CRAN.R-project.org/package=caret) (Accessed: 30 November 2020) is exploited to build the base learners and the stacking regression framework. The fundamental of stacking regression is illustrated in Figure 2. The data pairs, canopy multispectral reflectance, and grain yield were randomly and evenly split into ten parts, and then one of the parts was used for the test. The predictions for each fold were made by training the model and performing tenfold cross-validation. Based on this, an out-of-sample predictions matrix (OSPM) was obtained. During the tenfold cross-validation, grain yield predictions of each regression model were generated separately to check the results of the base learner on the test set before they were averaged. An OSPM with a dimension of m × n (m is the number of base learners, and n is the number of samples in the training set) was obtained after the m base learners completed the above process. Then, the OSPM was used to train the level-2 regression model to make final predictions. The multiple linear regression (MLR) was used as the level-2 model to avoid collinearity among the prediction results for grain yield. The tenfold cross-validations were also conducted for the level-2 model to reduce uncertainty in prediction results. Especially, the same split process (with tenfold cross-validation) was performed in all models to ensure fair comparisons between the methods. To avoid uncertainty in results, the process of dividing data into training and test sets was repeated 20 times with tenfold cross-validation. This process generated 200 models, and the average prediction accuracy of these 200 models in the test set was taken as the final evaluation index.

2.4.1. Random Forest

Random forest (RF) regression can be regarded as a machine learning model that combines a large number of regression trees [24]. The final output prediction results are the averaged value of all the trees. Regression trees represent a set of conditions or restrictions, and these trees are constructed by bootstrap sampling from a training sample set. This construction strategy overcomes the shortcomings of data and is easy to overfit in the case of complex trees. The key step in constructing an RF is splitting regression trees, and the splitting criterion is based on selecting the input variable with the lowest Gini Index:

$${\mathrm{I}}_{\mathrm{G}}\left({\mathrm{t}}_{\mathrm{X}\left({\mathrm{x}}_{\mathrm{i}}\right)}\right)=1-{\displaystyle \sum }_{\mathrm{j}=1}^{\mathrm{m}}\mathrm{f}{\left({\mathrm{t}}_{\mathrm{X}\left({\mathrm{x}}_{\mathrm{i}}\right)},\mathrm{j}\right)}^2$$

(1)

2.4.2. Support Vector Machine

Support vector machine (SVM) is derived from statistical learning theory and the minimum structural risk principle. SVM is widely applied to data analysis and pattern recognition [25]. The purpose of SVM is to minimize the error by adding a hyperplane and maximizing the margin between positive and negative samples in the training set. The introduction of the loss function allows SVM to solve nonlinear regression problems. Support vector regression can be defined as follows:

$$f\left(x\right)={\displaystyle \sum }_{i=1}^n\left({\widehat{a}}_i-a_i\right)k\left(x_i,x\right)+b$$

(2)

where a represents the additional hyperplane alongside the regression line; k (x_i, x) is the kernel function, and b is the bias.

2.4.3. Gaussian Process

Gaussian process (GP) is a type of probabilistic kernel machine based on Bayesian and statistical learning theory [26]. It has been extensively used in the field of machine learning. In probability theory and mathematical statistics, GP stands for the observed values in a continuous domain (such as time or space), and it can be regarded as the distribution and inference in a function space. In the process of prediction, GP maximizes the type-II maximum likelihood through the boundary likelihood of observations. In addition, it adjusts the hyperparameters and calculates the posterior distribution of unknown observations. When the observed variable space of GP belongs to a real field, a regression can be conducted to make predictions.

2.4.4. Ridge Regression

Ridge regression was proposed by Tikhonov in 1943 [27] and generalized by Hoerl and Kennard in 1970 [62]. It shrinks the regression coefficients by penalizing them or constraining their possible values. Specifically, ridge regression minimizes the sum of squared errors including an L2-norm penalty on the size of the parameter estimates. The coefficient estimates are, and the corresponding formula is as follows:

$${\widehat{\beta }}_{RR}={}_{\mathsf{\beta } }{}^{\arg \min } {\displaystyle \sum }_{i=1}^N{\left(y_i-{\widehat{\beta }}_0-{\displaystyle \sum }_{j=1}^p{x}_{ij}{\widehat{B}}_j\right)}^2+\lambda {\displaystyle \sum }_{j=1}^p{\widehat{\beta }}_j^2{}^{ }$$

(3)

where $\lambda$ controls the amount of shrinkage; N represents the number of observations; y is the dependent variable; p and ${\widehat{B}}_{j }$ are the number of independent variables and the value of the jth coefficient, respectively.

2.4.5. Cross-Validation and Hyperparameter Tune

At level-1 of the stacking regression, the tenfold cross-validation was exploited to form a sample matrix, and this procedure is considered as outer cross-validation. Meanwhile, the inner cross-validation and grid search were conducted to fine-tuning the hyperparameters of the base learners shown in Figure 2. In the outer cross-validation, the original data set was randomly divided into 10 equal subsets (Figure 3). Each time, one of the subsets was used for validating, and the remaining nine subsets for training. Each training set used for the outer cross-validation was also randomly and evenly split into 10 sets, with 10% of the data in the inner validating set and the remaining 90% in the inner training set. In the inner cross-validation, different combinations of candidate hyperparameters were set to construct the model on the inner training set. Then, the constructed model was validated on the inner validating set. Each hyperparameter combination was validated 10 times, and the hyperparameter combination with the highest average validation accuracy on the outer training set was transferred to the outer cross-validation to train the ideal model. The detailed hyperparameters of each machine learning are listed in

Table 4. Candidate hyperparameters of each machine learning algorithm.

Number	RF		SVM		RR	GP
	Ntree	Mtry	Cost	Gamma	Lambda	Sigma
1	405	3	0.250	0.450	0.00058	0.41
2	410	4	0.263	0.453	0.00063	0.42
3	415	5	0.275	0.455	0.00067	0.43
4	420	6	0.288	0.458	0.00072	0.44
5	425	7	0.300	0.460	0.00077	0.45
6	430	8	0.313	0.463	0.00083	0.46
7	435	9	0.325	0.465	0.00089	0.47
8	440	10	0.338	0.468	0.00095	0.48
9	445	11	0.350	0.470	0.00102	0.49
10	450	12	0.363	0.473	0.00110	0.50
11	455	13	0.375	0.475	0.00118	0.51
12	460	14	0.388	0.478	0.00126	0.52
13	465	15	0.400	0.480	0.00136	0.53
14	470	16	0.413	0.483	0.00146	0.54
15	475	17	0.425	0.485	0.00156	0.55
16	480	18	0.438	0.488	0.00168	0.56
17	485	19	0.450	0.490	0.00180	0.57
18	490	20	0.463	0.493	0.00193	0.58
19	495	21	0.475	0.495	0.00207	0.59
20	500	22	0.488	0.498	0.00058	0.60

ntree means the number of regression trees, and mtry means the number of input variables per node; cost is the parameter that controls the trade-off between minimization of the model’s complexity and minimization of the training error; gamma is the parameter for radial basis kernel function, and it determines the distribution of the data mapped to the new feature space; lambda is model penalization value; sigma is the parameter that inverses kernel width for the radial basis kernel function; RF, random forest; SVM; support vector machine; GP, Gaussian process; and RR, ridge regression.

.

2.5. Model Performance Evaluation

To evaluate the model performance for grain yield prediction, the coefficients of determination (R²) and root mean square error (RMSE) were calculated through the following equations:

$$R^2=1-\frac{{{\displaystyle \sum }}_{i=1}^n{\left(y_i-{\widehat{y}}_i\right)}^2}{{{\displaystyle \sum }}_{i=1}^n{\left(y_i-\overline{y}\right)}^2}$$

(4)

$$RMSE=\sqrt{\frac{1}{n}{\displaystyle \sum }_{i=1}^n{\left(y_i-{\widehat{y}}_i\right)}^2}$$

(5)

where n is the number of samples; $y_i$ and ${\widehat{y}}_i$ are the measured and the predicted grain yield of sample i, respectively; $\overline{y}$ represents the mean of the measured grain yield. The model with a higher value of R² and lower values of RMSE can predict grain yield better.

2.6. Statistical Analysis

A mixed linear model was exploited to test the significance of variation between genotypes, irrigation treatments, and their interactions for vegetation index and grain yield. The model is presented as follows:

Y = Xβ + Zµ + ε

(6)

where Y represents the response demonstrated by fixed effect (β) and random effect (µ) with random error (ε). X and Z are fixed and random effects, respectively. Heritability was estimated through the following formula:

H² = σ_g²/(σ_g² + σ_ε²/r)

(7)

where r is the number of replicates per treatment; σ_g² and σ_ε² indicate the genotypic and error variances, respectively [63]. The spectral traits with low heritability were considered as noise features in previous studies [4]. In this study, the VIs with heritability less than 0.5 were not used as the input features to construct grain yield prediction models.

3. Results

3.1. Phenotypic Analysis

The distribution of averaged grain yield data of irrigation treatments is given in Figure 4. The grain yield under limited irrigation treatment was 8% less than that under full irrigation treatment. It can be seen from the results listed in

Table A5. Significance test and heritability of the vegetation indexes extracted in the mid-grain filling stage.

Vegetation Index	Genotype (G)	Treatment (T)	G × T	H2
	F-Value	F-Value	F-Value	Full Irrigation	Limited Irrigation
NDVI	9.091 ***	1121.021 ***	1.365 *	0.80	0.81
SAVI	9.091 ***	1121.035 ***	1.365 *	0.80	0.81
OSAVI	9.091 ***	1121.020 ***	1.365 *	0.80	0.81
NRI	4.860 ***	1682.199 ***	1.275 *	0.63	0.69
GNDVI	13.186 ***	614.600 ***	1.223 *	0.86	0.86
SIPI	8.883 ***	765.657 ***	1.357 *	0.78	0.82
PSRI	8.278 ***	779.147 ***	1.544 ***	0.75	0.81
CRI	6.456 ***	2.374	0.928	0.71	0.76
EVI	4.408 ***	576.88 ***	1.096	0.60	0.68
MSR	7.881 ***	1353.925 ***	1.104	0.76	0.79
NLI	5.014 ***	1205.454 ***	1.184	0.62	0.71
RDVI	6.774 ***	785.527 ***	1.080	0.72	0.79
TVI	9.124 ***	1073.981 ***	1.398 *	0.80	0.81
MTVI2	4.487 ***	791.938 ***	1.113	0.61	0.68
NDRE	17.429 ***	528.708 ***	1.089	0.88	0.91
DVIREG	8.626 ***	997.224 ***	1.044	0.78	0.81
OSAVIREG	17.429 ***	528.698 ***	1.089	0.88	0.91
RDVIREG	7.046 ***	761.566 ***	1.144	0.73	0.77
MSRREG	18.188 ***	526.719 ***	1.106	0.88	0.91
MTCI	17.312 ***	636.939 ***	1.111	0.88	0.90

*, **, ***, significant at p < 0.05, p < 0.01 and p < 0.001, respectively.

that the coefficient of variation (CV) under full irrigation treatment was 15.7%, which was slightly higher than that under limited irrigation treatment (14.6%). The analysis of variation (ANOVA) results revealed significant variations (p < 0.001) among the genotypes for all vegetation indices (VIs) across the growth stages and grain yield. In addition, there was a significant difference between the two irrigation treatments for all traits (

Table A1. Main statistical parameters of grain yield under two irrigation treatments.

Treatment	Mean (t ha−1)	CV (%)	F-Value			H2
			Genotype (G)	Treatment (T)	G × T
Full Irrigation	7.59	15.7	10.401 ***	264.432 ***	0.98	0.85
Limited Irrigation	6.92	14.6				0.89

***, significant at p < 0.001; CV means coefficient of variation.

,

Table A2. Significance test and heritability of the vegetation indexes extracted in the heading stage.

Vegetation Index	Genotype (G)	Treatment (T)	G × T	H2
	F-Value	F-Value	F-Value	Full Irrigation	Limited Irrigation
NDVI	4.854 ***	743.179 ***	1.263 *	0.51	0.78
SAVI	4.855 ***	743.178 ***	1.262 *	0.51	0.78
OSAVI	4.856 ***	743.226 ***	1.262 *	0.51	0.78
NRI	4.599 ***	1284.147 ***	1.344 *	0.52	0.80
GNDVI	5.998 ***	235.842 ***	1.104	0.70	0.74
SIPI	3.439 ***	564.833 ***	1.222 *	0.50	0.62
PSRI	1.148	0.771	1.156	0.14	0.13
CRI	7.058 ***	7.053 **	0.968	0.78	0.72
EVI	7.777 ***	480.102 ***	1.482 **	0.82	0.71
MSR	3.952 ***	710.450 ***	1.022	0.47	0.75
NLI	4.184 ***	1446.413 ***	1.261 *	0.44	0.74
RDVI	4.298 ***	325.237 ***	0.988	0.56	0.73
TVI	4.890 ***	741.107 ***	1.278 *	0.51	0.78
MTVI2	8.176 ***	1005.586 ***	1.436 **	0.81	0.76
NDRE	17.346 ***	139.202 ***	1.166	0.88	0.90
DVIREG	7.126 ***	1302.732 ***	0.928	0.65	0.82
OSAVIREG	17.346 ***	139.204 ***	1.166	0.88	0.90
RDVIREG	9.764 ***	1103.257 ***	1.349 *	0.73	0.87
MSRREG	17.958 ***	138.719 ***	1.177	0.88	0.91
MTCI	17.304 ***	184.176 ***	1.187	0.88	0.91

*, **, ***, significant at p < 0.05, p < 0.01 and p < 0.001, respectively.

,

Table A3. Significance test and heritability of the vegetation indexes extracted in the flowering stage.

Vegetation Index	Genotype (G)	Treatment (T)	G × T	H2
	F-Value	F-Value	F-Value	Full Irrigation	Limited Irrigation
NDVI	6.436 ***	948.477 ***	1.408 *	0.71	0.77
SAVI	6.436 ***	948.481 ***	1.408 *	0.71	0.77
OSAVI	6.436 ***	948.440 ***	1.408 *	0.71	0.77
NRI	3.687 ***	502.703 ***	1.132	0.59	0.57
GNDVI	6.622 ***	740.962 ***	1.463 **	0.73	0.78
SIPI	5.619 ***	2053.124 ***	1.565 ***	0.67	0.75
PSRI	1.153	2.017	1.158	0.64	0.13
CRI	1.929 ***	3.624	1.281 *	0.34	0.54
EVI	1.725 ***	154.691 ***	1.481 **	0.35	0.51
MSR	5.995 ***	1008.409 ***	1.371 *	0.70	0.76
NLI	3.189 ***	694.061 ***	1.511 **	0.46	0.72
RDVI	2.950 ***	253.226 ***	1.195	0.44	0.70
TVI	6.458 ***	937.304 ***	1.413 *	0.71	0.77
MTVI2	2.009 ***	292.869 ***	1.487 **	0.39	0.59
NDRE	10.989 ***	35.906 ***	1.044	0.77	0.91
DVIREG	3.525 ***	244.143 ***	1.347 *	0.49	0.80
OSAVIREG	10.989 ***	35.908 ***	1.044	0.77	0.91
RDVIREG	4.918 ***	209.064 ***	1.526 ***	0.56	0.85
MSRREG	11.078 ***	34.183 ***	1.046	0.77	0.91
MTCI	10.682 ***	73.627 ***	1.053	0.77	0.90

*, **, ***, significant at p < 0.05, p < 0.01 and p < 0.001, respectively.

,

Table A4. Significance test and heritability of the vegetation indexes extracted in the early grain filling stage.

Vegetation Index	Genotype (G)	Treatment (T)	G × T	H2
	F-Value	F-Value	F-Value	Full Irrigation	Limited Irrigation
NDVI	6.407 ***	131.876 ***	1.484 **	0.69	0.79
SAVI	6.407 ***	131.873 ***	1.484 **	0.69	0.79
OSAVI	6.407 ***	131.889 ***	1.484 **	0.69	0.79
NRI	5.019 ***	906.706 ***	1.424 **	0.51	0.79
GNDVI	9.13 ***	45.687 ***	1.452 **	0.79	0.83
SIPI	7.784 ***	17.953 ***	1.705 ***	0.75	0.83
PSRI	1.156	0.303	1.162	0.75	0.13
CRI	5.786 ***	442.256 ***	1.082	0.69	0.72
EVI	4.652 ***	751.814 ***	1.129	0.64	0.67
MSR	5.834 ***	139.053 ***	1.327 *	0.65	0.79
NLI	3.265 ***	812.713 ***	1.144	0.45	0.56
RDVI	5.635 ***	2.048	1.104	0.64	0.78
TVI	6.404 ***	131.186 ***	1.49 **	0.69	0.78
MTVI2	4.417 ***	886.16 ***	1.106	0.60	0.66
NDRE	19.785 ***	67.904 ***	1.349 *	0.88	0.93
DVIREG	2.376 ***	450.393 ***	0.804	0.37	0.35
OSAVIREG	19.785 ***	67.9 ***	1.349 *	0.88	0.93
RDVIREG	3.432 ***	342.097 ***	0.746	0.58	0.48
MSRREG	20.303 ***	67.16 ***	1.375 *	0.88	0.93
MTCI	19.822 ***	45.575 ***	1.401 *	0.88	0.93

*, **, ***, significant at p < 0.05, p < 0.01 and p < 0.001, respectively.

and Table A5).

The heritability (H²) was high in the mid-grain filling stage, which ranged from 0.60 to 0.91 for all VIs (Table A1, Table A2, Table A3 and Table A4). The VIs with heritability values greater than 0.5 were selected as the input features for the model to make grain yield prediction for wheat in a particular growth stage. Under full irrigation treatment, the numbers of VIs with a heritability value greater than 0.5 were 17, 14, 18, and 20 for heading flowering EGF and MGF, respectively, and the the numbers of VIs with a heritability value greater than 0.5 under limited irrigation treatment were 19, 19, 18, and 20 across growth stages (Figure 5).

3.2. Performance Base Learners for Grain Yield Prediction

Under full irrigation treatment, the individual base learners achieving the best yield prediction results for the four growth stages are as follows: RR with mean R² of 0.488 and mean RMSE of 0.859 t ha⁻¹ for the heading stage; GP with mean R² of 0.520 and mean RMSE of 0.831 t ha⁻¹ for the flowering stage; RF with mean R² of 0.602 and mean RMSE of 0.711 t ha⁻¹ for the early grain filling stage, and RR with mean R² of 0.604 and mean RMSE of 0.777 t ha⁻¹ for the mid grain filling stage (Figure 6 and

Table 5. Test accuracies (R²) of base models and ensemble models for grain yield prediction.

Model	Full Irrigation				Limited Irrigation
	Coefficient of Determination (R2)				Coefficient of Determination (R2)
	Heading	Flowering	EGF	MGF	Heading	Flowering	EGF	MGF
RF	0.485	0.505	0.602	0.531	0.541	0.604	0.599	0.626
SVM	0.435	0.469	0.540	0.519	0.510	0.580	0.506	0.537
GP	0.465	0.520	0.589	0.563	0.550	0.595	0.512	0.573
RR	0.488	0.515	0.588	0.604	0.561	0.611	0.531	0.612
RF-SVM	0.488	0.522	0.620	0.562	0.543	0.614	0.617	0.628
RF-GP	0.492	0.534	0.620	0.589	0.555	0.612	0.618	0.625
RF-RR	0.502	0.515	0.611	0.601	0.564	0.619	0.609	0.629
SVM-GP	0.467	0.526	0.591	0.567	0.549	0.600	0.517	0.579
SVM-RR	0.472	0.547	0.616	0.610	0.569	0.617	0.541	0.578
GP-RR	0.496	0.551	0.618	0.622	0.570	0.616	0.542	0.613
RF-SVM-GP	0.490	0.528	0.619	0.588	0.549	0.612	0.615	0.628
RF-SVM-RR	0.499	0.533	0.624	0.607	0.563	0.620	0.613	0.627
SVM-GP-RR	0.493	0.546	0.619	0.623	0.564	0.617	0.540	0.613
RF-GP-RR	0.500	0.544	0.621	0.620	0.568	0.616	0.613	0.628
RF-SVM-GP-RR	0.498	0.538	0.622	0.620	0.562	0.617	0.611	0.628

The ensemble models with prediction accuracy lower than base models are highlighted by underline. Abbreviations: RF, random forest; SVM; support vector machine; GP, Gaussian process; and RR, ridge regression.

). SVM did not achieve the best predictions for any of the growth stages. Under limited irrigation treatment, the individual base learners achieving the best yield prediction results for the four growth stages are as follows: RR with mean R² of 0.561 and mean RMSE of 0.699 t ha⁻¹ for the heading stage; RR with mean R² of 0.611 and mean RMSE of 0.660 t ha⁻¹ for the flowering stage; RF with mean R² of 0.599 and mean RMSE of 0.668 t ha⁻¹ for the early grain filling stage, and RF with mean R² of 0.626 and mean RMSE of 0.641 t ha⁻¹ for the mid grain filling stage. In addition, it can be seen that the prediction results under limited irrigation treatment were slightly higher than those under full irrigation treatment.

As shown in Figure A1 and Figure A2, the Pearson’s correlations between the grain yield from predictions generated by the machine learning algorithms across the four growth stages and the two irrigation treatments were high (r = 0.81–0.97). By contrast, the density curves of the grain yield predictions were slightly different from each other. Moreover, the base learners differed significantly in the distribution interval of the prediction accuracy.

3.3. Ensemble Approach for Grain Yield Prediction

To evaluate the performance of the stacking method, the grain yield prediction is obtained by using different combinations of the base learners for the four growth stages and the two irrigation treatments. The grain yield prediction results obtained by combining two of the base learners are presented in Figure 7 and Table 5. The model combinations achieving the best mean prediction results for the four growth stages under full irrigation treatment are as follows: RF-RR with mean R² of 0.502 and mean RMSE of 0.843 t ha⁻¹ for the heading stage, GP-RR with mean R² of 0.551 and mean RMSE of 0.808 t ha⁻¹ for the flowering stage, RF-GP with mean R² of 0.620 and mean RMSE of 0.709 t ha⁻¹ for the early grain filling stage, and GP-RR with mean R² of 0.622 and mean RMSE of 0.752 t ha⁻¹ for the mid grain filling stage. The model combinations achieving the best mean prediction results for the four growth stages under the limit irrigation treatment are as follows: GP-RR with mean R² of 0.570 and mean RMSE of0.691 t ha⁻¹ for the heading stage, RF-RR with mean R² of 0.619 and mean RMSE of 0.643 t ha⁻¹ for the flowering stage, RF-GP with mean R² of 0.618 and mean RMSE of 0.648 t ha⁻¹ for the early grain filling stage, and RF-RR with mean R² of 0.629 and mean RMSE of 0.639 t ha⁻¹ for the mid grain filling stage.

As for the stacking regression combing three base learners, the model combinations with high values of R² and lowest value of RMSE for the four growth stages under full irrigation treatment are as follows: RF-GP-RR with mean R² of 0.500 and mean RMSE of 0.845 t ha⁻¹ for the heading stage, SVM-GP-RR with mean R² of 0.546 and mean RMSE of 0.813 t ha⁻¹ for the flowering stage, RF-SVM-RR with mean R² of 0.624 and mean RMSE of 0.709 t ha⁻¹ for the early grain filling stage, and SVM-GP-RR with mean R² of 0.623 and mean RMSE of 0.752 t ha⁻¹ for the mid grain filling stage (Figure 8). The model combinations under limited irrigation treatment, the model combinations achieving the best mean prediction results for the four growth stages under full irrigation treatment are as follows: RF-GP-RR with mean R² of 0.568 and mean RMSE of 0.693 t ha⁻¹ for the heading stage, RF-SVM-RR with mean R² of 0.620 and mean RMSE of 0.643 t ha⁻¹ for the flowering stage, RF-SVM-GP with mean R² of 0.615 and mean RMSE of 0.650 t ha⁻¹ for the early grain filling stage, and RF-SVM-GP with mean R² of 0.628 and mean RMSE of 0.639 t ha⁻¹ for the mid grain filling stage.

Figure 9 shows the performance of the stacking regression combing all the four base machine learning models for grain yield prediction of the four growth stages. Under full irrigation treatment, the mean values of R² for the stacking regression were improved to 0.498, 0.538, 0.622, and 0.620 for the heading stage, flowering stage, early grain filling stage, and mid-grain filling stage, respectively. Meanwhile, the mean values of RMSE were reduced to 0.851 t ha⁻¹, 0.820 t ha⁻¹, 0.709 t ha⁻¹, and 0.713 t ha⁻¹ for the heading stage, flowering stage, early grain filling stage, and mid-grain filling stage, respectively. Similar findings were also observed under limited irrigation treatment. The mean values of R² were up to 0.562, 0.617, 0.611, and 0.628, and the mean values of RMSE were 0.702 t ha⁻¹, 0.647 t ha⁻¹, 0.651 t ha⁻¹, and 0.642 t ha⁻¹, for the heading stage, flowering stage, early filling stage, and mid-grain filling stage, respectively. We estimated the spatial distribution of the grain yield at the plot scale based on the ensemble model by combining four base learners at mid grain filling. The difference of grain yield under the two irrigation treatments can be observed directly in Figure A3.

In terms of the mean R² and RMSE, most of the ensemble models achieved higher prediction accuracy than the individual best models (Table 5), which confirmed the effectiveness of the ensemble model implemented in this study. In addition, the overall accuracy indicated that the prediction performance of the ensemble model was proportional to the number of base learners (Figure 10).

3.4. Regression Coefficient Results for a Secondary Model

Figure 11 illustrates the distribution of the regression coefficients for the four base learners within the secondary model (MLR). A higher regression coefficient of a particular base learner indicated a larger weight in the stacking procedure. Under full irrigation treatment, the stacking performance for the heading stage was strongly influenced by the RR and GP models with high mean regression coefficients of 0.89 and 0.61, respectively. Similar results were observed for the flowering stage, where the mean regression coefficients of the RR and GP were 0.59 and 0.51, respectively. The impacts of the RF and SVM models on the stacking regression were relatively small with the mean regression coefficients of −0.01 and −0.02, respectively. The RF had the highest impact on the stacking performance for the early grain filling stage with a mean regression coefficient of 0.50, followed by the SVM with a mean regression coefficient of 0.33, RR with a mean regression coefficient of 0.03, and GP with a mean regression coefficient of 0.21. For the mid-grain filling stage, the RF and RR models had similar weights with the mean regression coefficients of 0.46 and 0.45, respectively, followed by the GP and SVM with the mean regression coefficients of 0.08 and 0.03. Under limited irrigation, RF showed a higher influence on the stacking performance with a mean regression coefficient of 0.61, while the GP, RR, and SVM models had mean regression coefficients of 0.31, 0.14, and −0.06, respectively. For the flowering stage, early grain filling stage, and mid-grain filling stage, the impact of the RF-based model was with the regression coefficient of 0.53, 0.55, and 0.05, respectively; the impact of the RR-based model was with the regression coefficient of 0.41, 0.52, and 0.50, respectively; the impact of SVM-based model was with the regression coefficient of 0.24, 0.39, and 0.13, respectively, the impact of GP-based model was with the regression coefficient of −1.14, −0.41, and 0.36, respectively. Overall, the results indicated that the stacking regression achieved higher prediction accuracy by allocating a more reasonable weight of base learners under various modeling conditions.

4. Discussion

The UAV-based multispectral vegetation indices have been increasingly exploited to predict plant physiological traits by crop breeding studies [5,64]. In this study, 20 VIs covering all five light bands captured by the multispectral sensor were taken to evaluate the machine learning models for grain yield prediction. The high heritability of the VIs indicated an excellent accuracy of these traits for grain yield prediction (Table A1, Table A2, Table A3 and Table A4). Among these VIs, GNDVI, NDRE, OSAVIREG, MSRREG, and MTCI had a heritability greater than 0.5 under both the treatments and across the growth stages. These five VIs could be selected as the best ones to predict variations among the genotypes for grain yield. Previously, some studies have reported that GNDVI and NDRE were the best predictors for grain yield and nutrient uptake efficiencies across the growth stages [5,64]. Therefore, these VIs can be used to select high-yielding genotypes with high accuracy in large breeding programs. Concerning water stress condition between different irrigation treatments, it has been found that thermal images show a correlation between minor changes in water stress that are undetectable by the multispectral indices as normalized difference vegetation index (NDVI) [65]. Thermal imagery can help diagnosis of water stress in plants, causes by the stomatal closure, which determines the reduction of the transpiration rate and decreasing evaporative cooling increases leaf temperature. In this regard, canopy thermal imaging represents a fast and practical way to evaluate and estimate crop water status, indicating plant’s water content. In some cases, relationship of canopy temperature has also been reported high with grain yield at early grain filling stages but lower at later grain fulling stages [66]. This might be due to low greenness and photosynthetic activity at late maturation stages. When wilting starts and transpiration decreases, the temperature difference between plots could get small which weakened the discrimination ability of canopy temperature at plot level growth. Therefore, data fusion of multispectral and canopy temperature imagery for GY prediction at a critical growth stage can increase the accuracy of prediction analysis. Meanwhile, the introduction of machine learning models can further improve the prediction accuracy and decisions during selection. To achieve precise prediction results and reduce the risk of overfitting, machine learning algorithms usually adopted feature selection strategies to reduce the data dimension to a suitable level [34]. In this study, VIs with heritability greater than 0.5 were selected across the growth stages, which means that all the base learners and their combinations were evaluated at maximum input data accuracy and repeatability. Previously, feature selection algorithms such as recursive feature elimination (RFE) [37] and Boruta [67] have been applied to prediction analysis. Different from RFE and Boruta, the feature selection based on heritability can be performed without knowing the predictor variables in breeding work. The successful use of heritability to reduce the number of input features has been reported in the previous study [4], and it achieved high prediction accuracy. Therefore, taking VIs with high heritability as input data increased the prediction accuracy and reliability of the models to assess the variations of genotypes.

Previously, several studies have exploited different machine learning models such as RR, SVM, RF, and GP to predict grain yield and physiological attributes of plants on remote sensing data sets [32,33,34,35,36]. This study aims to evaluate the combination of the above base learners to form an ensemble learning approach for grain yield prediction in different growth stages. It is difficult for a single machine learning model to predict several plant attributes using similar algorithms. For example, random forest algorithms are among the most powerful machine learning algorithms because of their proven accuracy, adaptability, and simplicity [32]. Such algorithms have been used in applications ranging from forest growth monitoring to winter wheat leaf chlorophyll content estimation [32,68]. However, the model performance of random forest is not as good as ridge regression in most cases considered by this study. Different characteristics of trait composition made it difficult for an algorithm of any model to maintain the prediction ability and accuracy in the case of several varied traits as input data [34,37]. The RR model has been reported to have high accuracy and robustness under most of the modeling conditions [69,70,71,72]. In our study, better prediction results were also obtained by the RR model compared to other models (Figure 6 and Table 5). It might be attributed to the biased nature of the model for collinear data analysis and comprehensive relationship for most of the traits [69,73]. However, the prediction ability of the RR model could be low in some cases, due to traits, growth stages, and growing conditions [74]. Therefore, it is important to find a method that can combine the advantages of multiple models to achieve an improved prediction accuracy under various growth conditions for grain yield prediction.

Ensemble learning approaches have been reported to increase the diversity of algorithms by combining different base learners, and the combination of more heterogeneities of base learners improves the ensemble learning model with higher prediction ability. When implementing a stacking method, it is necessary to include self-sufficient, independent, and diverse base learners for analysis [34,37]. In this study, four machine learning algorithms with different principles and internal structures were successfully combined, which achieved higher accuracy for predicting grain yield than the single base learner. In addition, the accuracy parameters (R² and RMSE) of the base learners exhibited more fluctuations with the wide ranges than that of the ensemble approaches (Figure 7, Figure 8 and Figure 9), indicating the stability of the ensemble approaches on new data.

Among all the combinations with two, three, and four base learners, the one that contained the RR model performed well in terms of prediction accuracy across growth stages and treatments. Meanwhile, the regression coefficient analysis of combining the four base learners showed a higher weight of RR in most of the cases. Therefore, the RR model had a great influence on the ensemble method for grain yield prediction. Combining two or three base learners can contribute to higher prediction accuracy than combining four base learners in some cases, but these model combinations were unstable and performed poorly in other situations. Conversely, the model combination of four base learners exhibited good accuracy and stability in all the cases, which is of great significance to the practical applications of the prediction model. Overall, the stacking regression method can achieve higher prediction accuracy than individual base learners, and the improvement was proportional to the number of base learners (Figure 10). This is consistent with the conclusion of a previous study for estimating hourly and continuous ground-level PM2.5 concentrations [44]. The stacking regression can be further optimized for grain yield prediction. The results of this study showed that the larger the number of base learners, the higher the accuracy of the final model (Figure 10). It means that more machine learning algorithms should be incorporated by the stacking regression method. Therefore, deep learning-based regression methods such as multilayer perceptron (MLP) [75] can be used as a new algorithm in the stacking procedure, while multiple birth SVM regression [76] and parallel RF regression [77] can be respectively used as variants of the SVM and RF regression for further improvement of the model. In addition, stacking a large number of base learners requires a level-2 model to perform the multicollinearity data analysis in the model [78]. Thus, ridge regression, least absolute shrinkage, selection operator (LASSO), and elastic net regression (ENET) can be used as level-2 models for collinearity analysis [78,79,80,81]. The grain yield predictions made by single models under full and limited irrigation treatments were similar across the growth stages, indicating the adaptability of UAV data for grain yield prediction. The prediction results obtained from ensemble learning were very similar under both irrigation treatments i.e., full irrigation treatment (R² = 0.625) and limited irrigation treatment (R² = 0.628) at the mid-grain filling stage. This result indicated that the mid-grain filling stage is the most appropriate stage to predict the grain yield under full and limited irrigation treatments. This is rational since grain filling is the stage where wheat transfers organic matter such as starch and protein produced by photosynthesis from the vegetative organs to the grains, and the vegetation indexes in this period are closely related to the final quality of thousand-grain weight [82].

5. Conclusions

Recently, there is an increasing focus on using multi-model fusion for data analysis to increase the prediction accuracy in crops. We have successfully used the ensemble learning method to stack the multiple base leaners to increase the grain yield prediction accuracy in wheat. The experimental results showed that the stacking of multi base learners possess the capability to improve the traits estimations effectively as compared to simple regression methods. Since UAV based phenotyping platforms can estimate spectral information from multiple growth stages cost-effectively, the use of ensemble learning approaches is important in increasing the accuracy of within season yield predictions. Our results illustrated the usefulness of ensemble learning approach for yield prediction using multi-stage multispectral data. We also demonstrated the weight of each base leaner such as random forest (RF), support vector machine (SVM), Gaussian process (GP), and ridge regression (RR) in ensemble leaning model development for yield prediction. Grain yield prediction results were high at a mid grain filling stage when a four base learner combination was used in an ensemble learning model as compared to other growth stages and base leaner combinations. To date, relatively few studies have been done to use the information obtained by UAV based sensors as inputs in ensemble learning model prediction of grain yield in winter wheat. Further validations of ensemble learning methods on multiple crops and UAV based data are required to increase its validity and authenticity in crop breeding.