Stacking of Canopy Spectral Reflectance from Multiple Growth Stages Improves Grain Yield Prediction under Full and Limited Irrigation in Wheat

Abstract

Grain yield (GY) prediction for wheat based on canopy spectral reflectance can improve selection efficiency in breeding programs. Time-series spectral information from different growth stages such as flowering to maturity is considered to have high accuracy in predicting GY and combining this information from multiple growth stages could effectively improve prediction accuracy. For this, 207 wheat cultivars and breeding lines were grown in full and limited irrigation treatments, and their canopy spectral reflectance was measured at the flowering, early, middle, and late grain fill stages. The potential of temporal spectral information at multiple growth stages for GY prediction was evaluated by a new method based on stacking the multiple growth stages data. Twenty VIs derived from spectral reflectance were used as the input feature of a support vector regression (SVR) to predict GY at each growth stage. The predicted GY values at multiple growth stages were trained by multiple linear regression (MLR) to establish a second-level prediction model. Results suggested that the prediction accuracy (R²) of VIs data from single growth stages ranged from 0.60 to 0.66 and 0.35 to 0.42 in the full and limited irrigation treatments, respectively. The prediction accuracy was increased by an average of 0.06, 0.07, and 0.07 after stacking the VIs of two, three, and four growth stages, respectively, under full irrigation. Similarly, under limited irrigation, the prediction accuracy was increased by 0.03, 0.04, and 0.04 by stacking the VIs of two, three, and four growth stages, respectively. Stacking of VIs of multiple important growth stages can increase the accuracy of GY prediction and application of a stable stacking model could increase the usefulness of data obtained from different phenotyping platforms.

1. Introduction

Wheat production is restricted by the occurrence of climate extremes such as heat and drought, bringing food security challenges for the increasing global population [1]. Drought and heat stress were estimated to reduce worldwide yield by 9 to 10% annually [2], requiring breeders to develop varieties with yield stability under different abiotic stress conditions. Plant breeding programs regularly conduct field evaluations to select candidate genotypes for breeding programs [3], and selection accuracy for high grain-yielding genotypes will enable breeders to expedite the above process. However, conventional methods of measuring secondary traits are complex, with a bottleneck of precise phenotyping because of labor intensiveness and time consumption [4]. Rapid and early evaluation of large experimental plots based on secondary traits through remote sensing is now a routine activity in many plant breeding programs [5]. Proximal sensing has attracted particular attention for predicating traits related to GY over the years. It offers rapid and nondestructive collection of time-series data at low cost for a number of traits from breeding trials grown in different environments [6].

Proximal remote sensing deals with spectral reflectance of light from the crop canopy and provides information about biophysical, bio-chemical, absorbance, and transmission of electromagnetic energy of plant tissues [7,8]. Based on different wavelengths of spectral reflectance, some vegetation indices (VIs) have been derived from various mathematical combinations, which can significantly reduce the noise interference caused by soil background, atmospheric conditions, and variation in sunlight intensity during spectral measurements for precise estimation of plant properties [9]. VIs are frequently used as a proxy for plant physiological traits such as chlorophyll content, biomass, and leaf area index (LAI) [10,11,12]. Moreover, VIs have also been validated to provide accurate wheat yield prediction results [3,13,14,15]. Some studies demonstrated that VIs from the flowering or grain fill stages provide higher prediction accuracy than other growth stages when estimating wheat yield [13,15]. However, the optimal time point for prediction of GY can vary with genotype and environmental conditions. Thus, the approach of combining multiple VIs was proposed to improve the accuracy of GY prediction. For example, the accumulated VIs ∑PRVI (Nir, Red) and ∑(R_Nir/(R_Red + R_Green)) derived from satellite data from jointing to grain fill stages predicted GY with high accuracy compared to VIs at individual growth stages [16]. Combining VIs extracted from UAV multispectral images for two random growth stages of rice by multiple linear regression has achieved higher yield prediction accuracy than that of an individual growth stage [17]. However, relatively few studies have investigated the use of multi-temporal VIs derived from hyperspectral data for the prediction of GY. The hyperspectral data contain nearly continuous reflectance information from various parts of the canopy with a wide range of wavelengths that can be noisy and cause big-data problems [3,10]. For this, advanced algorithms are required to extract valuable information that can be interpreted into useful plant traits by establishing an empirical model [10]. Different machine learning algorithms with improved and high prediction accuracy have been applied to build predictive models for plant properties [10,18,19,20,21]. A common machine learning algorithm, namely support vector regression (SVR), was widely used for the estimation of physiological parameters in crop species in several studies [22,23,24]. With this, above-ground biomass was estimated using hyperspectral traits as input variables by building the SVM model, and the best hyperparameter combinations were obtained through the gray wolf optimization algorithm [25]. Various studies predicted leaf area index, green leaf area, and green leaf chlorophyll contents with high accuracy in rice, and protein content, moisture, starch and ash in wheat using the SVM regression model [26,27]. These results indicated that SVM could be used in prediction various plant traits. In addition to individual machine learning algorithms to construct crop parameter evaluation models, ensemble methods are gaining more and more attention in the precision agriculture community. The stacking ensemble method was first proposed by Wolpertin [28]. It is an ensemble method that integrates multiple base learners to improve prediction performance. It has been used in different remote-sensing-based applications such as map composites [29], forest cover monitoring [30], and alfalfa yield prediction [31]. Stacking tends to employ heterogeneous learners and utilize differences between them to enhance the final accuracy [31]. Analogously, there is also heterogeneity in canopy spectral information at different crop growth stages, such that integration of hyperspectral data from multiple growth stages through stacking may provide better results for pre-harvest GY prediction in a high throughput manner. Therefore, the objectives of the present study were to (1) compare the prediction accuracy of four individual growth stages using the SVR model for 20 VIs as input features, (2) integrate the predictive ability of multiple growth stages using stacking methods, and (3) identify the best combination of time points for prediction accuracy.

2. Materials and Methods

2.1. Experimental Population and Field Trials

Field trials were carried out at the Chinese Academy of Agricultural Sciences (CAAS) experimental station at Xinxiang in Henan province (35°18′N, 113°52′E), during the 2018–2019 and 2019–2020 cropping seasons. The panel of 207 wheat varieties and breeding lines were planted under two water treatments, full and limited irrigation, in randomized complete blocks with two replications. Seeds was planted in a plot of 4.2 m² (3 × 1.4 m) with six rows spaced 0.2 m apart. Both the full and limited irrigation treatments were irrigated equally with 250 mm of water at the tillering (GS25) stages, while irrigation was applied to the full irrigation treatment at the early jointing (GS35), heading (GS55), and early grain fill (GS70) stages with 250 mm of water. Some natural precipitation was also received in both seasons on average at 5.73 mm. Fertilizer and management for both treatments were optimized according to the local soil conditions in a similar manner. Harvesting was done using a combine harvester.

2.2. Measurement of Canopy Spectral Reflectance

Canopy spectral reflectance measurements were performed using an ASD FieldSpec spectroradiometer (FieldSpec 3, Analytical Spectral Devices ASD, Boulder, CO, USA), which was equipped with a probe containing a 25° field of view. The FieldSpec collects data in the 350–2500 nm spectral range. It had a sampling interval of 1.4 nm between 350 and 1000 nm and of 2 nm between 1000 and 2500 nm, with a resampled spectral resolution of 1 nm. Spectral resolution was 3 nm at 700 nm and 10 nm at 1400 nm. Spectral measurements were carried out at optimum light conditions between 11:00 am and 1:00 pm to avoid confounding effects from variation in solar radiation levels. A BaSO₄ calibration panel was used for spectroradiometer optimization and white referencing before crop reflectance measurement. Optimization and white referencing were repeated every 10 plots. The probe was placed vertically 1 m over the canopy during measurements performed at four different locations per plot and averaged to represent the canopy reflectance of that plot. View Spec software (ASD Inc., Boulder, CO, USA) was used to eliminate noise from spectral curves, calculate the average of numerous spectral curves, and generate a reflectance file. To overcome the noise probability during the spectrum collecting process, the adaptive degree polynomial filter (ADPF) was used. ADPF adds a statistical heuristic to the Savitzky−Golay method to improve signal fidelity while reducing statistical noise. In both 2018–2019 and 2019–2020, the spectral measurements were conducted at the flowering (GS65), early grain fill (EGF, GS73), mid grain fill (MGF, GS85), and late grain fill (LGF, GS90) stages in both irrigation treatments. To avoid the effect of phenological differences between genotypes, spectral data for genotype was collected according to its growth stages.

2.3. Vegetation Indices

Vegetation indices (VIs) are usually used to predict crop yield. In this study, 20 VIs (

Table 1. Vegetation indices used in this study.

Vegetation Index	Name	Formula	Reference
CAI	Cellulose absorption index	$0.5\cdot \left({\mathrm{R}}_{2000}+{\mathrm{R}}_{2200}\right)-{\mathrm{R}}_{2100}$	[32]
CARI	Chlorophyll absorption ratio index	${\mathrm{R}}_{700}\cdot \mathrm{abs}\left(\mathrm{a}\cdot 670+{\mathrm{R}}_{670}+\mathrm{b}\right)/{\mathrm{R}}_{670}\cdot {\left({\mathrm{a}}^2+1\right)}^{0.5}$ $\mathrm{a}=\left({\mathrm{R}}_{700}-{\mathrm{R}}_{550}\right)/150$ $\mathrm{b}={\mathrm{R}}_{550}-\left(\mathrm{a}\cdot 550\right)$	[33]
CI	Curvature index	${\mathrm{R}}_{675}\cdot {\mathrm{R}}_{690}/{\mathrm{R}}_{683}^2$	[34]
DWSI5	Disease-water stress indices	$\left({\mathrm{R}}_{800}+{\mathrm{R}}_{550}\right)/\left({\mathrm{R}}_{1660}+{\mathrm{R}}_{680}\right)$	[35]
Datt5	Datt	${\mathrm{R}}_{672}/{\mathrm{R}}_{570}$	[36]
GNDVI	Green normalized Difference vegetation index	$\left({\mathrm{R}}_{800}-{\mathrm{R}}_{550}\right)/\left({\mathrm{R}}_{800}+{\mathrm{R}}_{550}\right)$	[37]
MPRI	Modified photochemical reflectance index	$\left({\mathrm{R}}_{515}-{\mathrm{R}}_{530}\right)/\left({\mathrm{R}}_{515}+{\mathrm{R}}_{530}\right)$	[38]
MSAVI	Modified soil adjusted vegetation index	$0.5\cdot (2\cdot {\mathrm{R}}_{800}+1-{\left({\left(2\cdot {\mathrm{R}}_{800}+1\right)}^2-8\cdot \left({\mathrm{R}}_{800}-{\mathrm{R}}_{670}\right)\right)}^{0.5})$	[39]
MSI	Moisture stress index	${\mathrm{R}}_{1600}/{\mathrm{R}}_{817}$	[40]
MTCI	Meris terrestrial chlorophyll index	$\left({\mathrm{R}}_{742}-{\mathrm{R}}_{702}\right)/\left({\mathrm{R}}_{742}+{\mathrm{R}}_{702}\right)$	[41]
NDLI	Normalized difference lignin index	$\left(\log \left(1/{\mathrm{R}}_{1754}\right)-\log \left(1/{\mathrm{R}}_{1680}\right)\right)/\left(\log \left(1/{\mathrm{R}}_{1754}\right)+\log \left(1/{\mathrm{R}}_{1680}\right)\right)$	[42]
NDMI	Normalized Difference Moisture Index	$\left({\mathrm{R}}_{820}-{\mathrm{R}}_{1600}\right)/\left({\mathrm{R}}_{820}+{\mathrm{R}}_{1600}\right)$	[43]
NDVI	Normalized difference vegetation index	$\left({\mathrm{R}}_{750}-{\mathrm{R}}_{705}\right)/\left({\mathrm{R}}_{750}+{\mathrm{R}}_{705}\right)$	[44]
NDWI	Normalized difference water index	$\left({\mathrm{R}}_{872}-{\mathrm{R}}_{1245}\right)/\left({\mathrm{R}}_{872}+{\mathrm{R}}_{1245}\right)$	[45]
OSAVI	Optimized soil-adjusted vegetation index	$\left(1+0.16\right)\cdot \left({\mathrm{R}}_{800}-{\mathrm{R}}_{670}\right)/\left({\mathrm{R}}_{800}+{\mathrm{R}}_{670}+0.16\right)$	[46]
PRI	Photochemical reflectance index	$\left({\mathrm{R}}_{531}-{\mathrm{R}}_{570}\right)/\left({\mathrm{R}}_{531}+{\mathrm{R}}_{570}\right)$	[47]
PWI	Plant water index	${\mathrm{R}}_{900}/{\mathrm{R}}_{970}$	[48]
SRPI	Simple ratio pigment index	${\mathrm{R}}_{430}/{\mathrm{R}}_{680}$	[49]
SWIR LI	Short wave infrared litter index	$3.87\cdot \left({\mathrm{R}}_{2210}-{\mathrm{R}}_{2090}\right)-27.51\cdot \left({\mathrm{R}}_{2280}-{\mathrm{R}}_{2090}\right)-0.2$	[50]
VREI4	Vogelmann red edge index 4	$\left({\mathrm{R}}_{734}-{\mathrm{R}}_{747}\right)/\left({\mathrm{R}}_{715}+{\mathrm{R}}_{720}\right)$	[51]

Note: R denotes reflectance.

) consisting of multiple spectral regions were used as input features of SVR for yield prediction. To remove the design effect, the best linear unbiased estimates (BLUEs) of genotypes for GY and for each of the 20 VIs were calculated in each of the four growth stages.

Broad sense heritability was estimated to check the repeatability of all VIs across the growth stages and treatments (

Table A1. Heritability results for vegetation indices.

Vegetation Index	Full Irrigation				Limited Irrigation
	Flowering	EGF	MGF	LGF	Flowering	EGF	MGF	LGF
CAI	0.54	0.52	0.57	0.53	0.50	0.54	0.58	0.49
CARI	0.73	0.58	0.71	0.70	0.56	0.74	0.56	0.35
CI	0.80	0.66	0.65	0.58	0.66	0.63	0.68	0.37
DWSI5	0.80	0.56	0.79	0.76	0.63	0.69	0.72	0.41
Datt5	0.81	0.83	0.84	0.79	0.75	0.70	0.73	0.64
GNDVI	0.82	0.78	0.77	0.84	0.66	0.75	0.71	0.38
MPRI	0.74	0.67	0.76	0.77	0.69	0.66	0.57	0.38
MSAVI	0.46	0.46	0.51	0.71	0.48	0.53	0.55	0.45
MSI	0.78	0.58	0.68	0.72	0.70	0.69	0.66	0.42
MTCI	0.85	0.82	0.83	0.85	0.73	0.80	0.80	0.36
NDLI	0.40	0.39	0.40	0.58	0.41	0.36	0.42	0.39
NDMI	0.83	0.67	0.78	0.77	0.72	0.76	0.73	0.40
NDVI	0.81	0.75	0.77	0.84	0.66	0.75	0.73	0.46
NDWI	0.83	0.56	0.75	0.73	0.71	0.68	0.71	0.46
OSAVI	0.52	0.45	0.57	0.78	0.40	0.60	0.64	0.46
PRI	0.76	0.75	0.74	0.78	0.63	0.66	0.64	0.65
PWI	0.83	0.70	0.81	0.78	0.73	0.76	0.75	0.59
SRPI	0.77	0.76	0.76	0.79	0.66	0.65	0.75	0.50
SWIR LI	0.58	0.53	0.60	0.62	0.59	0.52	0.61	0.44
VREI4	0.85	0.82	0.84	0.86	0.72	0.78	0.79	0.34

Abbreviations: EGF, early grain fill; LGF, late grain fill; MGF, mid grain fill.

). Broad sense heritabilities (H²) were estimated by:

$$H^2={\sigma }_G^2/\left({\sigma }_G^2+{\sigma }_{GE}^2/e+{\sigma }_{\epsilon }^2/\left(er\right)\right)$$

(1)

where ${\sigma }_G^2$ and ${\sigma }_{\epsilon }^2$ represent the genotype and error variances, respectively, ${\sigma }_{GE}^2$ is genotype × environment interaction variance, e and r were the numbers of environments and replicates, respectively.

2.4. Support Vector Regression

Support vector machine (SVM) is a supervised machine learning algorithm widely used in data analysis and pattern recognition (Figure 1) [52]. It is used to find the best separating hyperplane on the feature space to maximize the interval between the positive and negative samples on the training set [53]. The learning strategies of SVM are based on the principle of structural risk minimization by minimizing empirical risk and confidence range. The introduction of the loss function allows SVM to be extended to solve the problem of regression [54], i.e., support vector regression (SVR). The SVR uses either linear or nonlinear kernel functions depending on whether the data relationship is linear or nonlinear. SVR can be defined as follows:

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

(2)

where a is the additional hyperplane alongside the regression line and b represents the bias. $k\left(x_i,x\right)$ represents the kernel function. Radial basis function (RBF) was selected as kernel function in this study because of its superiority to process the nonlinear data with less hyperparameters and smaller generalization errors than other functions [55].

2.5. Stacking Method

In this study, we tested a new idea of stacking the predictive ability using data from multiple growth stages instead of individual growth stages to increase the accuracy for GY prediction. The stacking growth stages method is similar to the stacking regression. The whole process is divided into two levels. SVR and multiple linear regression (MLR) integrates the predictability of VIs from multiple growth stages. In the stacking process, all datasets are initially randomly divided into training datasets and test datasets (3:1), and fivefold cross-validation (CV) on the training set is performed by SVR. The five models trained during the CV process also separately predict the GY on the test sets to produce different predictions for the five groups before averaging them. After implementing the above process for multiple growth stages, the out-of-sample predictions matrix generated from the training data and predictions matrix generated from the test data would be as the new training dataset and new test dataset for second-level regression model (MLR), respectively. Analogously, fivefold CV was performed in the second level regression to obtain final predictions. SVR and MLR models were constructed using the function “svm” in R-package “e1071” and function “lm”, respectively.

To avoid the influence of chance factors, the division of the original data into training and test datasets was repeated 40 times randomly. Each fivefold CV generated five models, resulting in 200 models in the first- and second-level model training stages after 40 random divisions, so that the prediction accuracy was assessed by the mean value of coefficient of determination (R²) and the root mean square error (RMSE) of the 200 tests, and the equations of the evaluation indices are shown in the following equations:

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

(3)

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

(4)

where n represents the number of samples, $y_i$ and ${\hat{y}}_i$ represent the actual and the predicted grain yields of sample i, respectively, $\overline{y}$ represents the mean of the actual grain yield. Models with higher R² and lower values of RMSE suggest better prediction.

When the stacking method involved two growth stages, joint time-points I-II, I-III, I-IV, II-III, II-IV, and III-IV indicated that stacking spectral traits from flowering and early grain fill, flowering and mid grain fill, flowering and late grain fill, early and mid grain fill, early and late grain fill, and mid and late grain fill, respectively. When the stacking method involved three growth stages, joint time-points I-II-III, I-II-IV, II-III-IV, and I-III-IV indicated that stacking spectral traits from flowering, early and mid grain fill; flowering, early and late grain fill; early, mid and late grain fill; and flowering, mid and late grain fill, respectively.

2.6. Hyperparameter Tuning Based on Grid Search and Cross-Validation

At the first level of the stacking growth stages method, the training data and outer CV were used to form an outer sample matrix and to verify the prediction accuracy in the test data. Moreover, the inner CV [56] was used for fine tuning the hyperparameters of SVR by the grid search method [57]. In the outer CV, the original training dataset was randomly and evenly divided into five subsets, and each time one of them was used for validation, and the remaining four subsets for training. Each training set of the outer CV was split, with 10% of the data in the inner validation set and the remaining 90% in the inner training set. The inner CV was implemented with the following values for parameters cost (0.520, 0.525, 0.530, 0.535, 0.540, 0.545, 0.550, 0.555, 0.560, and 0.565) and gamma (0.675, 0.680, 0.685, 0.690, 0.695, 0.700, and 0.705). The hyperparameters combination with the highest average prediction accuracy was treated as the best and transferred to the outer CV to train the model using the outer training set.

3. Results

3.1. Phenotypic Variation

The descriptive statistics and distribution of measured grain yield (GY) are shown in Figure 2. GY was normally distributed across the irrigation treatments in both growing seasons. Limited irrigation reduced the GY by 17.5% and 9.3% compared with full irrigation in growing seasons 2018–2019 and 2019–2020, respectively. Figure 2 shows that some plots in the limited irrigation treatments have higher yield than in the full irrigation treatments. This is because some varieties in this experiment came from foreign countries and are not suitable for local land conditions. The yield of these varieties under full irrigation treatments are also much lower than the yields of local varieties under limited irrigation treatments. The coefficients of variation (CV) of GY under full and limited irrigation treatments were similar in both growing seasons, whereas the BLUE values of GY in the full irrigation treatment had a higher CV compared to the limited irrigation treatment. Whereas, heritability results for all vegetation indices (VIs) were estimated high (up to 0.86) at most of the growth stages under both irrigation treatments (Table A1).

3.2. Model Performance of Individual Growth Stages

All 20 VIs were used as input features in SVR to predict GY for four individual growth stages under both irrigation treatments (Figure 3). In the full irrigation treatment, the cross-validation of the model testing phase showed slight variations in predictions at different growth stages, with mean R² = 0.60 (mean RMSE = 0.73 t ha⁻¹) at flowering, mean R² = 0.62 (mean RMSE = 0.71 t ha⁻¹) at EGF, mean R² = 0.66 (mean RMSE = 0.69 t ha⁻¹) at MGF, and mean R² = 0.63 (mean RMSE = 0.71 t ha⁻¹) at LGF. The SVR model results in the limited irrigation treatment showed lower mean R² compared to the full irrigation treatment across growth stages. The higher GY prediction accuracy was observed at MGF with mean R² = 0.42 (mean RMSE = 0.59 t ha⁻¹) compared to EGF (mean R² = 0.40; mean RMSE = 0.60 t ha⁻¹), flowering (mean R² = 0.39; mean RMSE = 0.61 t ha⁻¹), and LGF (mean R² = 0.35; mean RMSE = 0.64 t ha⁻¹) stages. The regression between the predicted GY from the various stages were analyzed (Figure 4). There was a significant positive correlation between the predicted values for each growth stage. Under both irrigation treatments, there was high correlation (r = 0.632–0.846) between yield predictions for the first three growth stages. Yield predictions at the LGF were less correlated (r = 0.482–0.808) with those at other growth stages.

3.3. Model Performance of the Stacking Method

The distributions of accuracy parameters of the MLR algorithms in predicting GY by stacking VIs of two growth stages are shown in Figure 5. There were six different combinations when two of the four growth stages were selected. In the full irrigation treatment, joint time-points II-IV yielded the highest mean R² of 0.72 (mean RMSE = 0.63 t ha⁻¹), followed by I-IV with mean R² = 0.71 (mean RMSE = 0.63 t ha⁻¹). I-III, II-III, and III-IV yielded similar mean R² values as 0.69 (mean RMSE = 0.65 t ha⁻¹), 0.70 (mean RMSE = 0.65 t ha⁻¹), and 0.70 (mean RMSE = 0.64 t ha⁻¹). I-II yielded the lowest prediction accuracy with the mean R² of 0.67 (mean RMSE = 0.67 t ha⁻¹). The prediction accuracy of III-IV in limited irrigation was the highest with mean R² of 0.46 and mean RMSE of 0.58 t ha⁻¹, followed by I-III (mean R² = 0.45, mean RMSE = 0.58 t ha⁻¹). The joint time-points I-II, I-IV, II-III, and II-IV yielded similar mean R² (0.42, 0.43, 0.43 and 0.42) with mean RMSE of 0.59 t ha⁻¹, 0.59 t ha⁻¹, 0.59 t ha⁻¹ and 0.60 t ha⁻¹, respectively.

The prediction accuracy was also significantly improved under both irrigation treatments when three growth stages were stacked (Figure 6). In the full irrigation treatment, the prediction accuracy of I-II-IV, II-III-IV, and I-III-IV are similar with mean R² of 0.73 (mean RMSE = 0.61 t ha⁻¹), 0.72 (mean RMSE = 0.61 t ha⁻¹), and 0.72 (mean RMSE = 0.62 t ha⁻¹), respectively. I-II-III yielded the lower prediction accuracy with mean R² of 0.70 and mean RMSE of 0.64 t ha⁻¹. For the limited irrigation treatment, similar results were observed where higher predictive performance was observed by II-III-IV and I-III-IV with a mean R² value of 0.46 (mean RMSE = 0.58 t ha⁻¹) and 0.47 (mean RMSE = 0.57 t ha⁻¹). I-II-III and I-II-IV yielded the same prediction accuracy with mean R² = 0.44 and mean RMSE = 0.59 t ha⁻¹.

Figure 7a presents the modeling performance of the MLR algorithm that stacked VIs of all four growth stages to predict GY in test phases. In the full irrigation treatment, MLR yielded a mean R² value of 0.73 with mean RMSE = 0.61 t ha⁻¹. An increase in mean R² = 0.07 and a reduction of 0.08 t ha⁻¹ in the mean RMSE was observed compared to individual stage with the best prediction accuracy. Similar results were observed for limited irrigation; the mean R² increased (0.46) and the mean RMSE (0.57 t ha⁻¹) decreased compared with the MGF (mean R² = 0.42; mean RMSE = 0.59 t ha⁻¹). Figure 7b shows the distribution of the regression coefficients of each individual growth stage within the MLR model to explain the principle of the higher prediction accuracy of the stacking method. A larger coefficient within the MLR indicated a higher weight in the stacking procedure. Under full irrigation, mean coefficients (standard deviation) at flowering, EGF, MGF, LGF were 0.27 (0.12), 0.40 (0.10), 0.05 (0.11), and 0.54 (0.08), respectively, indicating that the stacking performance depended heavily on LGF. Under the limited irrigation treatment, the highest coefficient was found in the MGF (0.57 ± 0.18), followed by the LGF (0.46 ± 0.19), flowering (0.38 ± 0.18), and EGF (−0.06 ± 0.20).

4. Discussion

Accurate estimations of VIs through proximal remote sensing from the wheat canopy to predict within season grain yield could accelerate crop breeding. Hyperspectral techniques with high spectral resolution and strong spectral continuity can collect continuous and fine spectral curves of objects at visible and near-infrared wavelengths [58,59]. In this study, both irrigation treatments yielded the best prediction accuracy at MGF with high heritability (VIs) up to 0.84, which can be an appropriate stage for collecting spectral data. Our results were quite similar to several studies which have reported the use of VIs derived from hyperspectral data to predict yield in different crops such as alfalfa [31], wheat [60], maize [61], rice [62], and potato [63]. However, these studies only used spectral information from a single growth period to predict yield. Previous studies have shown that the combination of multi-stage vegetation indices can improve crop yield predictions [16,17]. In this study, we used the power of multiple VIs to predict GY based on recent advances in stacking regression algorithms. Results showed that the prediction accuracy of GY can be significantly increased through using data from multiple growth stages in the stacking process. In the full irrigation treatment, both the joint time-points I-II-IV and II-III-IV (stacking spectral traits from all four stages) showed best prediction accuracy in 11 joint time-points, and the joint time-point FE was more practical for predicting GY because it used the first two growth stages that could help breeders in early evaluation of yield with higher prediction accuracy for a large number of genotypes. Joint time-points I-III-IV provided the best prediction accuracy in the limited irrigation treatment. Although the R² value of the LGF was not the best among the four growth stages in both treatments, all joint time-points with the highest prediction accuracy included LGF. In the stacking of four growth stages in both treatments, LGF was assigned a high weight within the MLR (Figure 7b), which was consistent with the results of regression analysis that the predicted yield values at the LGF were less correlated with the first three stages in both irrigation treatments (Figure 4). Therefore, VIs of LGF can provide more additional information related to yield, leading to the fact that adding the hyperspectral traits collected at LGF into the stacking process can effectively improve the prediction accuracy.

The grain yield is normally reflected in the three components such as thousand grain weight (TGW), spike number (SN), and grain number (GN) per spike [64]. All of these components mainly influence by photosynthesis per unit area during the reviving-heading stage, determined by the water, fertilizer across the growth stages [65]. VIs from different growth stages can capture various factors related to GY, and utilizing the spectral parameters of multiple growth stages could better predict GY. Our results showed higher prediction accuracy than previous studies, which used potential individual growth stages to predict crop yield [16,17]. Moreover, the simultaneous use of multiple VIs of various growth stages as independent variables will lead to a large increase in the number of input features of the machine learning and thus generate a redundant dataset, which is not conducive to the construction of the model. In addition, numerous input features can bring the risk of overfitting. Our results revealed superiority of the stacking method over individual growth stage-based predictions in capturing variation among genotypes for yield related traits.

Water deficiency generally influences GY by affecting the rate of senescence [66]. Rapid senescence is linked to a shorter grain fill period, resulting in lower GY [13]. Therefore, it is necessary to conduct experiments under various irrigation levels to breed drought-resistant varieties. In this study, GY predictions were higher in the full irrigation treatment compared to limited irrigation for a single growth stage, as well as the stacked growth stage method. The results were consistent with the previous findings that drought influences the quality of spectral data and impacts the accuracy of trait prediction [3,13]. The canopy area and coverage of crops were low under water stress, causing canopy spectral reflectance to be easily disturbed by the soil background, thus affecting the data accuracy. Another reason for the lower prediction accuracy in the limited irrigation treatment might be a lower range of variation in GY (coefficient of variation, 9.82%) of BLUE compared to the full irrigation treatment (coefficient of variation, 12.62%), hence decreasing the magnitude of the data variation and reducing the predictive ability of the model [67].

The experimental results presented here were based solely on VIs calculated from a proximal remote sensing instrument, and lead to the inefficiency of collecting canopy spectral information in large-scale breeding experiments. The integration results presented in this paper are worthwhile, particularly for high-throughput phenotyping platforms (HTPP) of GY prediction of various crops. UAV equipped with imaging equipment have been widely used in the evaluation of crop physiological parameters [3,13,19] and can rapidly measure the canopy spectra of experiments covering large areas. They can even provide spectral information for the entire growth period for the stacking process for better yield prediction. Moreover, improved regression models that were used in this study or in totally new machine learning algorithms can be employed when establishing GY prediction models for single growth stages. For instance, deep learning-based regression methods, such as deep neural networks [67] and deep belief networks [68], could be used as a new technology in the stacking procedure. In addition, the prediction performance of different types of regression techniques in each individual growth stage may be varies. Therefore, the optimum prediction accuracy could be obtained by choosing the most suitable machine learning method for each growth stage. For a secondary learner, there will inevitably be overlapping information in the spectral data between adjacent stages during canopy spectral data collection with high time resolution. It resulted in strong collinearity between the predicted yields of different stages and could be applied in collinear data analysis-based approaches such as ridge regression [69] or least absolute shrinkage and selection operator [21] as a secondary learner to get more reliable integration results.

5. Conclusions

In this study, the stacking of VIs from multiple growth stages was tested to increase the prediction accuracy for GY. Pre-harvest prediction of yield can help improve the efficiency of selection for breeding efforts, as well as optimize management practices. In addition, instead of using vegetation indices of individual growth stages to predict wheat yield, we proposed a new stacking method for integrating spectral information from multiple growth stages to increase the yield prediction accuracy. The test results demonstrated that the stacking method provided more stable information with high prediction accuracy than individual growth stage results in both irrigation treatments. The stacking methods could be valuable for HTPP datasets from ground and aerial remote sensing platforms. We suggest that the use of improved algorithms from machine learning and deep learning could increase the accuracy of GY prediction models. Ultimately, this should accelerate the efficiency of plant breeding. This study concluded that there is potential for the method of stacking multiple growth stages to be further improved. First, more growth stages could be incorporated into the stacking procedure. For example, heading time is the key stage of SN formation, and it was not included in this study.