Entropy Weight Ensemble Framework for Yield Prediction of Winter Wheat Under Different Water Stress Treatments Using Unmanned Aerial Vehicle-Based Multispectral and Thermal Data.

Crop breeding programs generally perform early field assessments of candidate selection based on primary traits such as grain yield (GY). The traditional methods of yield assessment are costly, inefficient, and considered a bottleneck in modern precision agriculture. Recent advances in an unmanned aerial vehicle (UAV) and development of sensors have opened a new avenue for data acquisition cost-effectively and rapidly. We evaluated UAV-based multispectral and thermal images for in-season GY prediction using 30 winter wheat genotypes under 3 water treatments. For this, multispectral vegetation indices (VIs) and normalized relative canopy temperature (NRCT) were calculated and selected by the gray relational analysis (GRA) at each growth stage, i.e., jointing, booting, heading, flowering, grain filling, and maturity to reduce the data dimension. The elastic net regression (ENR) was developed by using selected features as input variables for yield prediction, whereas the entropy weight fusion (EWF) method was used to combine the predicted GY values from multiple growth stages. In our results, the fusion of dual-sensor data showed high yield prediction accuracy [coefficient of determination (R 2) = 0.527-0.667] compared to using a single multispectral sensor (R 2 = 0.130-0.461). Results showed that the grain filling stage was the optimal stage to predict GY with R 2 = 0.667, root mean square error (RMSE) = 0.881 t ha-1, relative root-mean-square error (RRMSE) = 15.2%, and mean absolute error (MAE) = 0.721 t ha-1. The EWF model outperformed at all the individual growth stages with R 2 varying from 0.677 to 0.729. The best prediction result (R 2 = 0.729, RMSE = 0.831 t ha-1, RRMSE = 14.3%, and MAE = 0.684 t ha-1) was achieved through combining the predicted values of all growth stages. This study suggests that the fusion of UAV-based multispectral and thermal IR data within an ENR-EWF framework can provide a precise and robust prediction of wheat yield.

Introduction

Bread wheat is one of the most important food crops that feed 40% of the world population (Liu et al., 2020). The timely and accurate evaluation of the grain yield (GY) before harvest can aid the selection of elite genotypes in large breeding programs (Mcbratney et al., 2005; Panda et al., 2010). Yield advocating traits, such as green biomass, leaf area index (LAI), and chlorophyll contents, have been used for within-season yield prediction (Hassan et al., 2018, 2019a). The canopy temperature is another important indicator of transpiration and leaf water potential under drought and heat stress and can help facilitate the selection of resilient genotypes (Zubler and Yoon, 2020). However, phenotyping most of these traits is destructive, time-consuming, and is associated with a high error probability. Therefore, the nondestructive measurements of the above proxy traits of the GY have been employed to increase the prediction accuracy of crop yield cost-effectively (Yu et al., 2016; Elsayed et al., 2017; Hassan et al., 2019a).

In the past few years, low-altitude remote sensing has attracted interest for its application in high-throughput crop phenotyping (Hassan et al., 2019b; Maimaitijiang et al., 2020). The advances in sensor technology have significantly accelerated the use of unmanned aerial vehicles (UAVs) for data collection with high spectral resolution as compared to satellite platforms (Colomina and Molina, 2014; Sidike et al., 2018). Various types of sensors mounted on UAV platforms, such as multispectral, hyperspectral, RGB, and thermal, are being widely used in the phenotypic evaluation of crops, with satisfactory data accuracy. The UAV-based nondestructive multispectral assessments of the LAI (Comba et al., 2020), biomass (Yue et al., 2019), chlorophyll content (Qiao et al., 2020), nitrogen use efficiency (Yang et al., 2020), senescence (Hassan et al., 2021), and GY (Hassan et al., 2019a) have been reported for several crops. These assessments are based on the spectral reflectance from the canopy of plants in the form of light bands with different wavelengths (Li et al., 2014). Thermal remote sensing is also being applied in precision agriculture to detect water stress (Suyoung et al., 2017) and plant resistance (Ludovisi et al., 2017). Recently, the focus has been increased on combining the data from multiple sources, where a group of datasets from multiple sensors is utilized obtained for plant trait estimation. Multi-source data models have the capability to improve crop trait estimations (Maimaitijiang et al., 2020). The use of canopy temperature and spectral information have been demonstrated to improve the model performance in estimating important plant traits for assessment of biotic/abiotic stress (Appeltans et al., 2020; Zubler and Yoon, 2020) and predicting the yield of soybean (Elmetwalli et al., 2020), barely (Rischbeck et al., 2014), and maize (Zhang et al., 2020). For crop yield prediction, flowering to grain filling stages are highly reliable, with good accuracy and repeatability (Hassan et al., 2019a; Hernandez et al., 2015). The predictions made in most studies have been based on the spectral information of an individual growth stage. The accumulation of VIs from jointing to the grain filling stage using a multiple linear regression algorithm has shown good prediction results in rice (Zhou et al., 2017). Since UAV-based temporal information of multispectral vegetation indices (VIs) and temperature can be obtained cost-effectively from multiple growth stages, combining data across the growth stages could help to achieve higher yield prediction accuracy. Machine learning algorithms have been employed with the canopy spectral features as input to construct models for crop trait evaluation, showing high prediction accuracy and adaptability (Wang et al., 2016; Wang J. et al., 2018). The commonly used machine learning algorithms are the random forest (RF) (Breiman, 2001), support vector machine (SVM) (Sain, 1997), and artificial neural network (ANN) (Bradley, 1995), and these have been successfully used for estimating biomass (Wang et al., 2016), LAI (Wang L. et al., 2018), chlorophyll content (Shah et al., 2019), and water content (Tavakoli and Gebbers, 2019). Among the machine learning algorithms, the emerging elastic net regression (ENR) algorithm has been considered one of the most precise prediction method for regression problems (Hui and Hastie, 2005). The ENR algorithm combines the advantages of ridge regression and least absolute shrinkage and selection operator (LASSO) regression to obtain better prediction results (Ogutu et al., 2012). At present, relatively few studies have been conducted on utilizing information obtained by UAV-based sensors as input to the ENR algorithm for the yield prediction of winter wheat.

The entropy weight algorithm is an emerging method in agricultural studies. It works by allocating the weight-based information entropy of the trait in the model (Li et al., 2011). It has been typically used for feature selection and model combination for combining datasets to assess ecosystem health (Cheng et al., 2020), monitor land-use change (Lu et al., 2014), and evaluate the coverage effectiveness of remote sensing satellites (Li et al., 2018). To the best of our knowledge, the entropy weight method has not been used to predict the yield values from multiple growth stages using UAV datasets. The objectives of this study were (1) to evaluate the potential of UAV-based multispectral and thermal sensors for the yield prediction of wheat using the ENR algorithm, (2) to identify the appropriate wheat growth stage for data collection to maximize the yield prediction accuracy, and (3) to investigate the potential of the entropy weight method in combining the predicted GY values from multiple growth stages.

Materials and Methods

Germplasm and Experimental Design

Field trials were conducted at the experimental station of the Institute of Farmland Irrigation of Chinese Academy of Agricultural Sciences in Xinxiang (113.8°E, 35.2°N) during the 2019–2020 cropping season (Figure 1). In total, 30 winter wheat varieties widely cultivated in the Yellow and Huai Valleys Winter Wheat Zone of China were used in this experiment. Germplasm was planted under three water stress treatments, namely, mild irrigation, moderate irrigation, and high irrigation, to obtain the UAV-based multispectral, thermal, and ground-truth GY data. Irrigation for each treatment was performed in the tillering, wintering, reviving, jointing, heading, and grain filling stages using a laterally moving sprinkler irrigation machine. The irrigation volume was calculated by the flow rate of the sprinkler nozzle and the duration of irrigation. The total irrigation volume for the mild, moderate, and high irrigation treatments were 145, 190, and 240 mm, respectively (Table 1). A completely randomized block design with two replications was adopted for the experiment. The size of each plot was maintained at 11.2 m2 with the dimensions of 8 m × 1.4 m, representing one cultivar with six rows at a spacing of 0.20 m. Field management (e.g., disease and pest control, fertilizer) was maintained at optimal levels depending on the local conditions. In the 2019–2020 growing season, the total precipitation was 115 mm, and the monthly average temperature was highest (23°C) in July and lowest (−6°C) in January. Wheat was harvested using a plot combine harvester in June 2020. The GY of each plot was weighed at a moisture content of approximately 12.5%.

FIGURE 1

Experimental location, design, and management.

TABLE 1

Irrigation strategy for each treatment.

Growth stage	Mild irrigation (mm)	Moderate irrigation (mm)	High irrigation (mm)
Tillering	35	35	35
Wintering	35	35	35
Turning green	20	25	35
Jointing	20	35	50
Heading	20	35	50
Grain filling	15	25	35
Total	145	190	240

Data Acquisition and Processing

Figure 2 shows the workflow for the data acquisition. A DJI M210 (DJI Technology Co., Shenzhen, China) carrying a RedEdge MX (MicaSense Parrot, France) multispectral camera and a Zenmuse XT2 (DJI Technology Co., Shenzhen, China) thermal sensor was used to collect high-resolution multispectral and thermal images simultaneously. The RedEdge MX featured five spectral sensors, namely, blue (475 nm), green (560 nm), red (668 nm), red-edge (717 nm), and near-IR (842 nm). The RedEdge MX camera automatically adjusts the ambient light effects through the sunshine sensor, thereby minimizing the error in the multispectral images. Zenmuse XT2 contains an 8-mm lens with a 57.12° × 42.44° field of view to record temperature measurements in the 7.5–13.5-μm spectral range with a measurement accuracy of ± 5°C. The DJI ground station was used as an automated flight control system, allowing the user to define the air route and customize the mission plan. Flight mission was executed for all the six growth stages from 11 a.m. to 1 p.m. under a cloudless sky. To avoid the effect of the phenological differences between treatments, the flight missions for each treatment were collected according to the growth stages. To obtain high-resolution images, each flight was set at an altitude of 30 m with 85% front and 80% side image overlapping. Before and after each flight, the calibration board was photographed to convert the digital number (DN) value of the multispectral image into reflectance during subsequent data processing. During the flights, the surface temperature of 12 boards was measured using a handheld thermometer for the radiometric calibration of the thermal images. To obtain the geographic reference of the multisensor UAV image, 18 ground control points (GCPs) were evenly arranged in the field, and their coordinates were measured with a millimeter-level accuracy using a differential global positioning system.

FIGURE 2

Schematic workflow of the methodology used in this study. P denotes the predicted grain yield (GY) value, and C1–C8 indicate the combinations of the values predicted from multiple growth stages. CV, cross-validation; VIs, vegetation indices; NRCT, normalized relative canopy temperature.

The Pix4Dmapper software (Pix4D SA, Lausanne, Switzerland) was employed for the orthomosaic generation using the UAV-based multispectral and thermal IR images. The geographic coordinates (World Geodetic System, 1984) of the GCPs were used in the photogrammetric workflow of Pix4Dmapper to improve the accuracy of the composite orthomosaics. Dense point clouds were generated using the structure-from-motion (SfM) method in Pix4Dmapper along with the photogrammetric workflow. After radiometric correction, the DN values of the multispectral and thermal IR images were converted to reflectance and temperature (°C). To extract the reflectance and temperature information for each plot, the orthomosaic images were segmented into 180 polygon shapes with assigned IDs defining the cultivars under different irrigation treatments. Polygon shape generation and information extraction are completed in QGIS 3.1.0.[1] A total of 22 indices were used in this study, of which 21 VIs were estimated from multispectral reflectance, and 1 index was calculated from the canopy temperature across the irrigation treatments (Table 2).

TABLE 2

Formulae of multispectral vegetation indices and normalized relative canopy temperature.

Acronym	Index	Formulae	Developer(s)
CIRE	Chlorophyll index RedEdge	(RNIR/RRE)−1	Gitelson et al., 2003
DVI	Difference vegetation index	RNIR−RR	Tucker et al., 1979
EVI	Enhanced vegetation index	2.5×(RNIR−RR)/(1RNIR + 6×RR−7.5×RB)	Huete et al., 2002
GNDVI	Green normalized difference vegetation index	(RNIR−RG)/(RNIR + RG)	Gitelson et al., 1996
MCARI1	Modified chlorophyll absorption in reflectance index 1	(RREG−RR−0.2×((RREG−RG))×(RREG/RR)	Daughtry et al., 2000
PSRI	Plant senescence reflectance index	(RR−RB)/RNIR	Merzlyak et al., 1999
MSR	Modified simple ratio index	((R_NIR/R_R)-1)/R_NIR/R_R1	Chen, 1996
MTCI	MERIS terrestrial chlorophyll index	(RNIR−RREG)/(RRE−RR)	Dash and Curran, 2004
MTVI2	Modified triangular vegetation index 2	1.5×[1.2×(R_NIR-R_G)-2.5×(R_R-R_G)]/CPSTABLEENTER[2×(R_NIR+1)^2-6×R_NIR+5×R_R-0.5]	Haboudane et al., 2004
NDVI	Normalized difference vegetation index	(RNIR−RR)/(RNIR + RR)	Rouse, 1972
NDVIRE	Normalized difference vegetation index RedEdge	(RNIR−RREG)/(RNIR + RREG)	Elsayed et al., 2015
NLI	Nonlinear vegetation index	(RNIR×RNIR−RR)/(RNIR×RNIR + RR)	Goel and Qin, 1994
OSAVI	Optimized soil-adjusted vegetation index	(RNIR−RR)/(RNIR + RR + 1.6)×1.16	Rondeaux et al., 1996
PPR	Plant pigment ratio	(RG−RB)/(RG + RB)	Metternicht, 2003
RDVI	Re-normalized difference vegetation index	(R_NIR-R_R)/(R_NIR+R_R)	Wang et al., 1998
RVI	Ratio vegetation index	(RNIR/RR)	Ba Ret and Guyot, 1991
NRI	Nitrogen reflectance index	(RG−RR)/(RG + RR)	Diker and Bausch, 2003
SAVI	Soil-adjusted vegetation index	(RNIR−RR)/(RNIR + RR + 0.5)×1.5	Huete, 1988
SIPI	Structure insensitive pigment index	(RNIR−RB)/(RNIR + RB)	Penuelas et al., 1995
TCARI	Transformed chlorophyll absorption ratio index	3[(RREG−RR)−0.2×(RREG−RG)×RREG/RR]	Eitel et al., 2007
TVI	Triangular vegetation index	0.5×[120×(RNIR−RG)−200×(RR−RG)]	Broge and Leblanc, 2001
NRCT	Normalized relative canopy temperature	T_i-T_minT_max-T_min	Elsayed et al., 2015

R

Gray Relational Analysis

In a gray relational analysis (GRA), a system with incomplete information is called a gray system, meaning that the relationship between the factors is uncertain (Aslan et al., 2012). When the experiment is unclear or when the experimental method cannot be implemented accurately, a gray analysis can help overcome the drawbacks in statistical regression (Jin et al., 2013). For example, there is a close relationship between VIs and yield; however, their detailed relationships remain unclear. Therefore, the main purpose of the GRA is to measure the degree of relationship within this system by analyzing the gray relationships between VIs and GY. The GRA procedure includes the following steps:

The reference series reflect the characteristics of the system behavior, and the comparison series influences the system behavior. In this study, the GY was considered the reference series, and each index was considered a comparison series. The reference sequence is represented by the following formula:

where n represents the number of samples, and n is 180 in this study. Comparison data series can be expressed as follows:

There are m comparison data series, each containing n-values.

Data in each factor column in the system may have different dimensions, making it difficult to compare or obtain a correct conclusion when comparing. Therefore, to ensure the reliability of the results, the following non-dimensional processing of the data is generally required when performing the GRA:

The calculation of the difference data series △ is as follows:

The gray relational coefficient ξ(k) for the kth data point in the ith difference data series can be expressed as follows:

where △min and △max are the global maximum and minimum values in the difference data series, respectively. △ is the kth value in the △ difference data series, and ζ is the distinguishing coefficient: ζ∈ [0, 1]. In this study, the distinguishing coefficient is set to 0.5.

Generally, the average value of the gray relational coefficient is taken as the gray relational degree (GRD), which is expressed as follows:

Elastic Net Regression

To avoid the instability of the LASSO solutions when the input features are highly correlated (e.g., a large number of VIs constructed from limited bands), the ENR has been proposed to analyze the high-dimensional data. The ENR is an extension of the LASSO, which is robust to severe multicollinearity among the input features (Ogutu et al., 2012). The ENR combines the penalties of the ridge regression (ℓ1) and LASSO (ℓ2) and can be expressed as follows:

On setting α = λ2/(λ1 + λ2), the ENR is equivalent to the minimizer of the following:

subject to , where Pα(β) represents the penalty of the ENR. The ENR can be considered LASSO and ridge regression when a = 0 and 1, respectively. The ℓ1 part of the ENR is used for automatic variable selection, while the ℓ2 part encourages grouped selection and stabilizes the solution paths with respect to the random sampling, thereby improving the prediction results. By introducing a grouping effect when selecting the variable, a group of highly correlated input features tends to have similar coefficients. The ENR can choose the groups of correlated features when these groups are unknown in advance. Notably, the ENR selects more than n variables when p>> n, which is different from the LASSO. In this study, there is inevitably a high correlation between the various VIs. Therefore, the ENR will be an ideal choice when using VIs as the input features for yield prediction.

Modeling Framework

In this study, a 10-fold outer cross-validation method was adopted to train and test the model. To avoid contingency, we conducted 50 iterations for the outer cross-validation, resulting in a total of 500 models. The average of the accuracy evaluation index generated from the 500 models was used to evaluate the model performance. In the process of outer cross-validation, the inner cross-validation and the grid search were conducted for parameter tuning of the ENR models (Figure 2). In the outer cross-validation, the VIs and GY data were randomly divided into 10 equal subsets. One of them was used for testing each time, and the remaining nine subsets were used for training. Each training set of the outer cross-validation was evenly divided into 10 sets, similar to the outer cross-validation. One of them was used for testing, and the nine subsets were used for training. During the inner cross-validation process, multiple combinations of the candidate parameters were set in the inner training set for model construction and then tested on the inner test set. Each parameter combination was tested 10 times, and the hyperparameter combination with the lowest average test error was set for the outer cross-validation for model training. This study uses the R package ‘‘caret’’[2] to construct the ENR model for yield prediction. In the “caret” package, the parameters to be adjusted are the fraction and quadratic penalty parameter lambda. Table 3 represents the candidate ranges of these two parameters.

TABLE 3

Candidate hyperparameters for elastic net regression.

Number	Lambda	Fraction	Number	Lambda	Fraction
1	0.050	0.000E+00	16	0.541	3.162E-03
2	0.083	1.000E-04	17	0.574	4.047E-03
3	0.116	1.280E-04	18	0.607	5.179E-03
4	0.148	1.638E-04	19	0.640	6.629E-03
5	0.181	2.096E-04	20	0.672	8.483E-03
6	0.214	2.683E-04	21	0.705	1.086E-02
7	0.247	3.433E-04	22	0.738	1.389E-02
8	0.279	4.394E-04	23	0.771	1.778E-02
9	0.312	5.623E-04	24	0.803	2.276E-02
10	0.345	7.197E-04	25	0.836	2.913E-02
11	0.378	9.211E-04	26	0.869	3.728E-02
12	0.410	1.179E-03	27	0.902	4.771E-02
13	0.443	1.509E-03	28	0.934	6.105E-02
14	0.476	1.931E-03	29	0.967	7.814E-02
15	0.509	2.471E-03	30	1.000	1.000E-01

Moreover, we tested the model performance on the test samples in the cross-validation procedure to test the adaptability of the model. The coefficient of determination (R2), root mean square error (RMSE), relative root-mean-square error (RRMSE), and mean absolute error (MAE) were adopted to evaluate the model performance. The calculation formulae of the parameters are as follows:

where n represents the number of samples, y and are the measured and predicted GY of sample i, and is the average value of the measured GY. The higher the R2-value, the lower the RMSE, RRMSE, and MAE values and the better the performance of the model for GY prediction.

Entropy Weight Method

The ENR algorithm was independently implemented at each growth stage. Instead of using these results to predict the GY individually, we proposed an entropy weight fusion (EWF) model that combines the predicted results from the different growth stages via weights obtained during the model training stage. The basic mechanism of the entropy weight method is to use the entropy to characterize the degree of disorder in the system (Farhadinia, 2017). In this method, the relative error between the predicted and measured values of the GY obtained in an individual growth stage by the selected ith prediction model can be expressed as follows:

where i = (1, 2, 3…m), j = (1, 2, 3…n), and y represents the predicted value of the yield forecast model for the ith individual growth stage on the jth plot. The process for calculating the weights is as follows:

The relative error ratio was calculated between the predicted value of the ith individual growth stage and the measured value at plot j:

where . The entropy value h was calculated for the relative error in the ith individual growth stage prediction:

The relative error variation coefficient was determined based on the principle of the opposite of the entropy value and its degree of variation:

The weight was then obtained for the predicted output value from a single growth stage:

The weights were obtained by combining the output forecast values from the multiple growth stages. The final output forecast value can be expressed as follows:

The higher the entropy of the prediction error sequence of a single growth stage, the lower the degree of variation and the greater the weight. The entropy weight method fully considers the relative error in the output prediction value from the different growth stages. Therefore, the predicted results from the multiple growth stages complement each other to improve the accuracy of yield prediction. In this study, eight combinations were created to evaluate the accuracy of the entropy weight method for GY prediction. Table 4 represents a detailed description of each combination.

TABLE 4

Combination of different growth stages used in the model for grain yield prediction.

Combination	Combination of growth stages
C1	Jointing, booting, heading
C2	Jointing, booting, heading, flowering
C3	Jointing, booting, heading, flowering, grain filling
C4	Jointing, booting, heading, flowering, grain filling, maturity
C5	Booting, heading, flowering
C6	Heading, flowering, grain filling
C7	Flowering, grain filling, maturity
C8	Heading, flowering, maturity

Results

Descriptive Statistics of Grain Yield

The distribution of yield from wheat plots is shown in Figure 3. The GY was normally distributed under all the irrigation treatments as well as across the treatments. The GY was found to be higher under high irrigation treatments than under the moderate and mild irrigation treatments. The mean GY values for the high, moderate, and mild irrigation treatments were 7.09, 5.99, and 4.40 t ha–1, respectively. The highest coefficient of variation (19.51%) was observed in the mild irrigation treatment and the lowest (12.70%) in the high irrigation treatment. The overall range of the GY data across the irrigation treatment was 2.79–8.64 t ha–1, with a data variation of up to 24.35%. Across treatment data with this type of variation can help evaluate the prediction accuracy of the model.

FIGURE 3

Grain yield distribution curves under various irrigation treatments.

Results of Gray Relational Analysis and Feature Selection

A total of 22 indices were ranked using the GRA method. Table 5 lists the results for all the growth stages. The GRD of the normalized relative canopy temperature (NRCT) ranked first for the jointing, booting, and flowering stages and relatively high for the heading (rank = 10) and grain filling (rank = 9) stages. However, the NRCT ranked last at maturity. The rankings for most VIs were unstable across all the growth stages. For example, plant pigment ratio (PPR) and difference vegetation index (DVI) had a high ranking at both jointing and booting but ranked low at flowering and grain filling. In accordance with the GRA mechanism, the higher the GRD between the main and the reference sequence, the more closely the sequences are related, which indicates a close relationship between the NRCT and the yield during the multiple growth stages.

TABLE 5

Ranking of indices using the gray relational analysis (GRA) for six growth stages.

Rank	Jointing	Booting	Heading	Flowering	Grain filling	Maturity
1	NRCT	NRCT	NLI	NRCT	MSR	TVI
2	PPR	PPR	NDVI	RVI	RVI	NLI
3	OSAVI	SAVI	TCARI	MSR	NDVI	SAVI
4	MTVI2	OSAVI	MSR	NDVI	NLI	RDVI
5	SAVI	RDVI	RVI	NLI	NRI	EVI
6	RDVI	MTVI2	MTVI2	OSAVI	OSAVI	DVI
7	MCARI	PSRI	OSAVI	MTVI2	MTVI2	OSAVI
8	PSRI	TCARI	PPR	TCARI	CIRE	MTVI2
9	TCARI	DVI	SAVI	RDVI	NRCT	PSRI
10	DVI	TVI	NRCT	NRI	NDVIRE	MCARI
11	TVI	MCARI	SIPI	SAVI	GNDVI	MTCI
12	EVI	EVI	RDVI	GNDVI	RDVI	NDVI
13	SIPI	SIPI	MCARI	SIPI	SIPI	RVI
14	GNDVI	NLI	PSRI	PSRI	PSRI	MSR
15	NRI	NDVI	EVI	EVI	SAVI	CIRE
16	NDVIRE	MSR	DVI	CIRE	EVI	NDVIRE
17	NLI	RVI	TVI	NDVIRE	TCARI	SIPI
18	CIRE	MTCI	GNDVI	DVI	TVI	NRI
19	MTCI	CIRE	NRI	TVI	DVI	GNDVI
20	NDVI	NRI	NDVIRE	MTCI	MTCI	TCARI
21	MSR	GNDVI	CIRE	PPR	MCARI	PPR
22	RVI	NDVIRE	MTCI	MCARI	PPR	NRCT

To further explore the features with better performance and to reduce the dimensionality of the data, the top feature was iteratively added into the ENR. The performance of the model (i.e., MAE) in the training process was updated until all the features were inputted into the ENR (Figure 4). Among the six developmental stages, the grain filling stage yielded the lowest error, and it tended to be stable when the number of features was 19. The highest error was observed in jointing, and the model showed a stable tendency after inputting 16 features. The appropriate numbers of input features for the booting, heading, flowering, and maturity stages were found to be 18, 18, 22, and 22, respectively.

FIGURE 4

Model training error as a function of the number of features. The order of input of features depends on the gray relational degree (GRD). MAE, mean absolute error.

Performance of Elastic Net Regression Model for Individual Growth Stage

To analyze the improvement of the thermal data for yield prediction accuracy, the model was first built using the features extracted from the multispectral images (Figure 5). The mean prediction values for the grain filling stage was R2 = 0.461, followed by flowering (R2 = 0.432), heading (R2 = 0.422), maturity (R2 = 0.417), booting (R2 = 0.290), and jointing (R2 = 0.130). Figure 6 represents the accuracy assessment results of the ENR model for GY predictions by using both thermal and multispectral features. The results show that the dual-sensor data fusion method achieves higher prediction accuracy at all measured stages compared to using single multispectral sensor-based features. As with using only multispectral features, the ENR showed the highest prediction results with a low error at the grain filling (R2 = 0.667) stage. The mean prediction results at jointing, booting, heading, flowering, and maturity were R2 = 0.544, R2 = 0.571, R2 = 0.602, R2 = 0.640, and R2 = 0.527, respectively.

FIGURE 5

Statistical distributions of (A) coefficient of determination (R2), (B) root mean square error (RMSE), (C) relative root-mean-square error (RRMSE), and (D) mean absolute error (MAE) of the elastic net regression (ENR) algorithm for GY prediction using multispectral features in test phases. JS, jointing stage; BS, booting stage; HS, heading stage; FS, flowering stage; GFS, grain filling stage; MS, maturity stage.

FIGURE 6

Statistical distributions of (A) R2, (B) RMSE, (C) RRMSE, and (D) MAE of the ENR for GY prediction using both multispectral and thermal features in test phases. JS, jointing stage; BS, booting stage; HS, heading stage; FS, flowering stage; GFS, grain filling stage; MS, maturity stage.

After obtaining the predicted GY using thermal IR and multispectral features, the regression between the predicted GY from the various stages was conducted (Figure 7). High correlations ranging from R2 = 0.59 to R2 = 0.89 between adjacent growth stages were observed across the growth stages. Moreover, the greater the interval between the growth stages, the lower the R2-value. For example, the R2 between the jointing stage and the booting, heading, flowering, grain filling, and maturity stages were 0.78, 0.66, 0.64, 0.52, and 0.41, respectively. In comparison, the correlations between the predicted yield in the maturity and other growth stages were lower, with quite weak regression values ranging from R2 = 0.41 to R2 = 0.59. There were differences in the distribution curves of predicted GY values, which provides complementary information.

FIGURE 7

Regression plots, density curve, and R2-values between predicted GY in six developmental stages.

Performance of Entropy Weight Fusion Method

For comparison with the EWF method, multispectral and thermal features from multiple stages were used as the inputs of ENR to the training model. The results indicated that combining the features of multiple stages increases the accuracy of yield prediction than individual stages (Figure 8). The C4 yielded the highest R2-value of 0.725, followed by C3 (R2 = 0.717) and C2 (R2 = 0.691). The remaining combinations achieved similar prediction accuracy (R2 = 0.669–0.681). However, the obvious fluctuations of accuracy parameters (R2, RMSE, RRMSE, and MAE) with wide ranges were observed.

FIGURE 8

Statistical distributions of (A) R2, (B) RMSE, (C) RRMSE, and (D) MAE of the ENR model that uses both multispectral and thermal features from different stages as inputs.

Figure 9 represents the performance of the EWF model in predicting the GY using the combined predicted values from the multiple growth stages. Comparing with the individual growth stages, the EWF model also provides a more accurate result regardless of the number of stages adopted. Among the eight combinations, the optimal test results of the EWF model were observed in C4, with a mean R2 of 0.729. An increase of 0.062 compared with the highest mean R2-value was observed in the grain filling stage (R2 = 0.667). Moreover, the RMSE, RRMSE, and MAE values were reduced to 0.831 t ha–1, 14.3%, and 0.684 t ha–1, respectively. A low prediction was observed in C1 (R2 = 0.681). Compared to C5 (R2 = 0.692), C6 (R2 = 0.678), C7 (R2 = 0.677), and C8 (R2 = 0.688), C2 (R2 = 0.721) and C3 (R2 = 0.719) had a more accurate predictions. The fluctuations in the accuracy parameters (R2, RMSE, RRMSE, and MAE) of the EWF model were more moderate compared to the model that directly used multistage features as inputs (Figures 8, 9), which again demonstrates the stability of the EWF method. A paired t-test was utilized to assess whether the EWF models performed statistically high in terms of the R2-values compared with the other models (Figure 10). The results showed significantly high R2-values for the EWF model in all the growth stage combinations.

FIGURE 9

Statistical distributions of (A) R2, (B) RMSE, (C) RRMSE, and (D) MAE of the entropy weight fusion (EWF) method for GY prediction in the test phases.

FIGURE 10

Results from paired t-test between model R2 obtained from the EWF method and the individual stages. *** significant at P ≤ 0.001; JS, jointing stage; BS, booting stage; HS, heading stage; FS, flowering stage; GFS, grain filling stage; MS, maturity stage.

Discussion

The UAV-based phenotyping is an emerging technique in practical crop breeding. Previous studies have shown that the UAV-based features and the machine learning model can be used together to predict crop yields in breeding work with a large number of crop genotypes (Osval et al., 2017; Fei et al., 2021). In this study, ENR is a relatively new machine learning algorithm being used for yield prediction. ENR combines the properties of ridge regression and LASSO (Ogutu et al., 2012), both of which have been successfully applied to crop yield prediction (Kang et al., 2021; Shafiee et al., 2021). The incorporation of multiple VIs adds collinearity to the models, and the ENR is robust to severe multicollinearity among the input features (Ogutu et al., 2012). Another reason for using ENR was the simplicity of the linear model compared to other machine learning algorithms such as RF or ANNs, which makes the model run less time-consuming and more efficient to train.

Several VIs such as normalized difference vegetation index (NDVI) and green normalized difference vegetation index (GNDVI) have been evaluated to monitor crop health under stress and predict the GY. Most of the multispectral VIs that have been reported were species-specific and easily saturated (Hatfield and Prueger, 2013). Therefore, it is challenging to predict important crop traits using a single VI (Wang et al., 2016). The successful use of multiple VIs to improve the prediction accuracy of important traits in crops has been reported in many studies (Wang et al., 2016; Jin et al., 2020). In this study, 21 multispectral VIs and 1 temperature index (i.e., NRCT) were measured in multiple growth stages to validate the UAV data and ENR and check their accuracy for GY prediction. For accurate yield predictions and to avoid model overfitting, the machine learning algorithms may benefit from using a feature selection algorithm to reduce the dimensionality of the data to an appropriate level (Yoosefzadeh-Najafabadi et al., 2021). The GRA is a widely accepted approach in feature selection (Deris et al., 2013; Lu et al., 2019; Yao et al., 2019; Miswan et al., 2021). The results in this study show that GRA can reduce the dimensionality of the input features to some extent.

The model performed poorly when using multispectral VIs to predict yield. This could be due to the saturation issue associated with the visible-near-infrared (Vis-NIR) sensor for dense vegetation such as wheat, soybean, and rice (Thenkabail et al., 2000; Tilly et al., 2015). The fusion of multispectral and thermal features includes canopy spectral and temperature information outperformed for yield prediction, which was consistent with previous reports (Maimaitijiang et al., 2017, 2020). Temperature is closely related to plant physiological processes such as transpiration, leaf water potential, and photosynthesis (Sagan et al., 2019). Generally, high canopy temperature is negatively correlated with crop yield (Tattaris et al., 2016; Sagan et al., 2019). Previously, the UAV-based thermal IR data has been successfully applied to plant trait evaluation (Gonzalez-Dugo et al., 2013; Ludovisi et al., 2017; Liu et al., 2018; Raeva et al., 2019; Crusiol et al., 2020). Previous studies have also shown that combining thermal data with multispectral data outperformed as compared to the fusion of spectral and structural information from RGB and multispectral images for the prediction of LAI, biomass, chlorophyll, and nitrogen in soybean (Maimaitijiang et al., 2017).

The spatial heterogeneity of the ground changes among the developmental stages of the crop could lead to a significant difference in the prediction accuracy across the growth stages (Juliane et al., 2014). The results of yield prediction based on the individual growth stages are similar to previous reports, i.e., wheat yield prediction accuracy was higher at grain filling stages under different growth environments (Hernandez et al., 2015; Hassan et al., 2019a). During the grain filling stage, the starch, protein, and organic matter produced by photosynthesis are transported to the grain (Guan et al., 2017), and this stage is closely linked to the thousand-grain weight. Therefore, the accuracy of yield estimation was highest at the grain filling stage. In addition, the reduction in the greenness and chlorophyll level after the grain filling stages due to the decrease in the degree of dry matter accumulation in the leaves of plants could influence the detection accuracy of VIs based on the red and near-IR light (Yue et al., 2017). This reduces the model performance in the late developmental stage, which causes a decrease in yield prediction accuracy at maturity. In addition, crop canopy information at varying growth stages is associated with different yield elements, and a combination of remote sensing data from multiple growth stages can effectively improve the yield prediction accuracy.

Another main objective of this study was to use an appropriate method to acquire the prediction values from a combination of temporal remote sensing data across the growth cycle. Although previous studies used temporal VIs for yield prediction, most of them used a single VI (Wang et al., 2014; Zhou et al., 2017), which can be influenced by different degrees of saturability or soil background (Wang et al., 2016). We first directly used the multispectral and thermal features from multiple stages as inputs to ENR, and this method was able to obtain higher yield prediction accuracy than individual stages, but the accuracy parameters fluctuated more compared to EMF and were slightly lower than the prediction accuracy of the EMF method for some combinations. This may be due to the redundant information generated by the accumulation of features from multiple growth stages. In addition, the excessive dimensionality of input features also poses the risk of overfitting the machine learning model (Feng et al., 2017; Coolen et al., 2020). Among the combinations of EMF, the prediction accuracy of C2 was comparable to a combination with the highest prediction accuracy of C4. The data required for C2 can be obtained at the flowering stage, which is appropriate for application in practical management. The results of this study suggest that the fusion of multispectral and thermal features within an entropy weight ensemble framework can provide accurate wheat yield predictions. However, more comprehensive studies, such as studies of different crop varieties in different environments, are needed to determine the most accurate and efficient multistage data for combination.

Conclusion

A rapid and nondestructive method for an accurate GY prediction of wheat is desired in breeding programs. In this study, an ensemble framework was developed to increase the GY prediction accuracy by integrating the predicted values from multiple stages using the UAV-based multispectral and thermal IR imagery. The test results showed that the prediction accuracy of the grain filling stage was the highest among the six growth stages. The ensemble method outperformed the individual stage-based GY prediction in terms of accuracy. Combining the features of the first four growth stages allows for early and accurate yield prediction to aid in decision-making. This study offers a new method for GY prediction through UAV-based remote sensing, and it can help in large breeding activities.

Data Availability Statement

The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding author/s.

Author Contributions

SF analyzed the data and wrote the manuscript. ZC and YX directed the trial and provided the main idea. QC and ZL helped to collect data. MH, YM, and MS provided comments and suggestions to improve the manuscript. All authors read and agreed to the published version of the manuscript.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest. The reviewer AR declared a past co-authorship with the authors MH, YX to the handling editor.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.