Estimation of Winter Wheat Yield Using Multiple Temporal Vegetation Indices Derived from UAV-Based Multispectral and Hyperspectral Imagery

Abstract

Winter wheat is a major food source for the inhabitants of North China. However, its yield is affected by drought stress during the growing period. Hence, it is necessary to develop drought-resistant winter wheat varieties. For breeding researchers, yield measurement, a crucial breeding indication, is costly, labor-intensive, and time-consuming. Therefore, in order to breed a drought-resistant variety of winter wheat in a short time, field plot scale crop yield estimation is essential. Unmanned aerial vehicles (UAVs) have developed into a reliable method for gathering crop canopy information in a non-destructive and time-efficient manner in recent years. This study aimed to evaluate strategies for estimating crop yield using multispectral (MS) and hyperspectral (HS) imagery derived from a UAV in single and multiple growth stages of winter wheat. To accomplish our objective, we constructed a simple linear regression model based on the single growth stages of booting, heading, flowering, filling, and maturation and a multiple regression model that combined these five growth stages to estimate winter wheat yield using 36 vegetation indices (VIs) calculated from UAV-based MS and HS imagery, respectively. After comparing these regression models, we came to the following conclusions: (1) the flowering stage of winter wheat showed the highest correlation with crop yield for both MS and HS imagery; (2) the VIs derived from the HS imagery performed better in terms of estimation accuracy than the VIs from the MS imagery; (3) the regression model that combined the information of five growth stages presented better accuracy than the one that considered the growth stages individually. The best estimation regression model for winter wheat yield in this study was the multiple linear regression model constructed by the VI of ‘$\left({\mathrm{b}}_1-{\mathrm{b}}_2\right)/\left({\mathrm{b}}_3-{\mathrm{b}}_4\right)$’ derived from HS imagery, incorporating the five growth stages of booting, heading, flowering, filling, and maturation with r of 0.84 and RMSE of 0.69 t/ha. The corresponding central wavelengths were 782 nm, 874 nm, 762 nm, and 890 nm, respectively. Our study indicates that the multiple temporal VIs derived from UAV-based HS imagery are effective tools for breeding researchers to estimate winter wheat yield on a field plot scale.

1. Introduction

Wheat (Triticum aestivum L.), as the third largest cereal crop in the world, plays an important role in world food production and food security strategies. According to the Food and Agriculture Organization of the United Nations (FAO), more than 220 million ha are sown with wheat, with over 770 million tons of wheat being produced in 2021 [1]. China accounts for over 23 million ha of wheat crop and is responsible for nearly 18% (about 137 million tons) of all wheat produced worldwide [1]. Wheat is a major food source for the people in North China. The North China Plain is one of China’s primary wheat-growing areas (mainly winter wheat). However, in this region, winter wheat becomes more vulnerable to drought stress due to the continental monsoon climate, leading to decreased output. Groundwater irrigation is one main source of water supply during the winter wheat growth season of the North China Plain, but it is severely restricted according to the sustainable development policy. Thus, in order to fulfill the dual purpose of food security protection and water conservation, it has become necessary to breed drought-resistant varieties of winter wheat to maintain production. As the ultimate detection target, the crop yield becomes an important selection parameter of winter wheat breeding. However, in the empirical breeding process, yield can only be measured after the crop growth cycle, which is costly, labor-intensive, and time-consuming for the breeding researchers. Therefore, establishing a method that could estimate winter wheat yield on a field plot scale within a short time would help them in selecting drought-resistant varieties.

Satellite remote sensing data have been widely applied for non-destructive crop yield estimation over various large-scale regions (from local to national, continental, and global) since the 1970s [2,3]. The yield estimation models based on these data have demonstrated reasonable crop yield estimation accuracy in the large-scale regions based on satellite imagery and have been widely used due to their convenience and simplicity [2,3,4,5,6,7,8]. However, for breeding researchers, the application of satellite data in yield estimation is often hampered by the high spatial heterogeneity within the small areas of breeding fields, missing data during critical crop growth stages, and high costs due to the coarse spatial resolution and fixed passing time and band setting.

In recent years, the development of sensor technologies in unmanned aerial vehicles (UAVs) has promoted their application for data acquisition [9]. Compared to satellite remote sensing, UAV remote sensing has an improved spatial, spectral, and temporal resolution and is associated with lower costs and greater flexibility and versatility [10], making UAVs increasingly popular in precision agriculture [11,12]. Previous studies have reported the relationship between crop yield and crop phenotypic parameters such as plant height, leaf nitrogen content, leaf area index, above-ground biomass, and so on [13,14,15,16]. UAV platforms have exerted a beneficial effect in the retrieval of a wide array of crop characteristics that are associated with yield [17,18,19,20], and their use can help meet breeding researchers’ requirements regarding crop yield estimation on a small plot area scale within a short amount of time.

Furthermore, vegetation indices (VIs) derived from UAV-mounted multispectral (MS) and hyperspectral (HS) sensors have been widely used to estimate crop yield [21,22]. Duan et al. developed a method based on VIs that were derived from UAV multispectral data to correlate with rice phenotyping and estimate grain yield [13]. García-Martínez et al. estimated corn grain yield by combining vegetation indices, canopy cover, and plant density using multispectral and RGB images acquired via the use of UAVs [23]. Ramos et al. proposed a random forest algorithm that performed well in tests with a ranking-based strategy that focused on predicting maize crop yield using UAV-based multispectral vegetation indices [24]. However, the application of UAV-mounted hyperspectral sensors for agricultural monitoring is limited by the weight of imaging systems and the complexity of image processing [25,26]. As a result, they have rarely employed by breeding researchers in the study of crop breeding. Moreover, the VIs used in previous studies have been mostly derived from limited growth stages within the crop growing season, which may increase the risk of missing critical spectral features in other growth stages [27,28]. Therefore, estimating crop yield through combining UAV-mounted HS sensor-based VIs in several critical crop growing stages could be helpful to integrate critical spectral features throughout the crop growing season and lead to a higher estimation accuracy, which would be beneficial to breeding researchers.

For this study, motivated by the need to boost the efficiency of winter wheat breeding, we aimed to evaluate yield estimation among various winter wheat breeding cultivars using multi-temporal vegetation indices derived from UAV-mounted multispectral and hyperspectral sensors. To accomplish this objective, we (1) investigated the correlation between the winter wheat yield and 19 UAV-based multispectral VIs and 17 band combination types of hyperspectral data at single growth stages and multiple growth stages, respectively, and (2) identified the best VI and the best time for estimating winter wheat yield using UAV data.

2. Materials and Methods

2.1. Experimental Setup

The experimental site was set up at the Dry-Land Farming Institute of Hebei Academy of Agricultural and Forestry Sciences (DFI) at Hengshui City, Hebei Province, China (37°54′15.63″N, 115°42′29.32″E, World Geodetic System 1984) (Figure 1). The area has a semi-arid temperate and monsoonal climate characterized by four distinct seasons, with an average yearly temperature of 13.3 °C and a yearly precipitation of 497.1 mm.

The experimental site design included eleven winter wheat cultivars: C1 (Chang8744), C2 (Shimai22), C3 (Luyuan472), C4 (Shimai15), C5 (HengH1603), C6 (Xinmai28), C7 (Jimai418), C8 (Shannong28), C9 (Nongda212), C10 (Heng4399), and C11 (Jimai22). Each cultivar was then split into seven different irrigation groups and three repeats (1.5 m × 6 m in size) according to a randomized block design. All cultivars were planted on 15 October 2020 at a density of 375 plants/m². A base fertilizer, pure nitrogenous fertilizer (225 kg/ha), P₂O₅ (112.5 kg/ha), and K₂O (112.5 kg/ha) were applied before sowing. No additional fertilizers were used for the growth of the winter wheat discussed in this study. The irrigation date of each irrigated sub-plot is shown in

Table 1. The irrigation date and corresponding growth stages of different irrigated sub-plots in the study.

Irrigation Group	Irrigation Date (d/m/y)	Growth Stage	Total Irrigation Volume (m3/ha)
A	3 April 2021	Jointing stage	1500
	3 May 2021	Flowering stage
B	None	--	0
C	29 November 2020	Overwintering stage	750
D	10 March 2021	Regreen stage	750
E	3 April 2021	Jointing stage	750
F	10 April 2021	Jointing stage	750
G	18 April 2021	Booting stage	750

. The irrigation volume of each time was 750 m³/ha. The total precipitation during the 2020–2021 growing season in the site was 43.9 mm.

2.2. Data Acquisition

2.2.1. Ground Truth Data

All winter wheat cultivars were harvested on 8 July 2021. The yield of each sub-plot was weighted and normalized to a moisture content of 13% and is expressed as ‘t/ha’. According to the experimental design, 231 samples were measured, while 17 measurements were eliminated as outliers, including 6 plots of Repeat 1, 4 plots of Repeat 2, and 7 plots of Repeat 3. The statistics of measured winter wheat yield are outlined in

Table 2. Descriptive statistics of measured winter wheat yield (t/ha).

Parameter	Number of Samples	Minimum	Maximum	Mean	Standard Deviation	Coefficient of Variation
Grain yield	214	46.46	11.24	8.58	1.28	14.97%

. The mean measured grain yield values under different irrigation groups and winter wheat cultivars are shown in Figure 2. Winter wheat yield differences across different irrigation groups and different cultivars were assessed using a one-way analysis of variance after checking the normality assumption at 0.05 probability level (

Table 3. Results of our one-way analysis of variance of winter wheat yield for different irrigation groups and cultivars at 0.05 probability level.

	F-Value	p-Value
Different irrigation groups	13.73	0.00
Different winter wheat cultivars	9.84	0.00

). There were significant differences among the different irrigation groups and different winter wheat cultivars. The grain yield of irrigation group B was lower than other irrigation groups, and group A had the highest yield. The C6 cultivar showed the poorest yield, while C9 performed best.

2.2.2. Multi-Sensor UAV Data

UAV images derived from multispectral (MS) and hyperspectral (HS) sensors were employed in this study (Figure 3). The UAV campaign was conducted under low wind speed and clear sky conditions between 10:00 a.m. and 2:00 p.m. local time to reduce the influence of atmospheric and solar radiation. The overlap percentages in the forward and lateral flying directions of both UAVs were 80% and 70%, respectively. The acquisition dates and details corresponding to the growth stages regarding both UAVs are shown in

Table 4. The details of UAV multispectral and hyperspectral imagery acquisition in the study.

	Acquisition Date (d/m/y)	Growth Stage
UAV imagery	18 April 2021	Booting stage
	28 April 2021	Heading stage
	12 May 2021	Flowering stage
	21 May 2021	Filling stage
	2 June 2021	Maturation stage

.

The DJI P4 Multispectral (DJI Technology Co., Ltd., Shenzhen, China) was used to collect multispectral images, including 5 sensors of blue, green, red, red-edge, and near-infrared with wavelengths of 456 nm (±16 nm), 560 nm (±16 nm), 650 nm (±16 nm), 730 nm (±16 nm), and 840 nm (±26 nm), respectively. The sensors used a 1/2.9 inch complementary metal-oxide-semiconductor (CMOS). The field of view was 62.7°, the focal length was 5.74 mm, the f-number was f/2.2, and focus was kept at the infinite point (∞). The MS images were obtained at a flight height of 50 m with an accuracy of flight altitude of 0.1 m, and the corresponding ground pixel resolution was 1.6 cm. The MS orthomosaic maps were generated using Pix4D mapper (Pix4D SA, Lausanne, Switzerland).

A DJI M600 Pro (DJI Technology Co., Ltd., Shenzhen, China) equipped with a Pika L hyperspectral camera (Resonon, Inc., Bozeman, MT, USA) was used to capture the hyperspectral images discussed in the study. The hyperspectral camera has 150 bands in a spectral range of 400–1000 nm with a spectral resolution of 4 nm. The HS images were acquired at a flight height of 50 m with an accuracy of flight attitude of 0.5 m, and the corresponding ground pixel resolution was 3.0 cm. SpectrononPro (Resonon, Inc., Bozeman, MT, USA) and ENVI 5.3 (Esri Inc., Redlands, CA, USA) were used to generate HS orthomosaic maps.

2.3. Vegetation Indices Calculation

The average of a 0.8 m × 4 m image area was used for band reflectance value calculation for each sub-plot. The area was approximately in the center of each sub-plot to eliminate the effect of the marginal areas of each sub-plot and other neighboring sub-plots. The extracted images are herein referred to as UAV images.

A large number vegetation indices have been proposed for grain yield estimation. A total of 19 MS VIs that were previously used for crop yield and yield-related phenotypic characteristic estimation were calculated, with average reflectance being derived from UAV MS images [3,18], and the 19 indices are shown in Table 5.

Additionally, for various spectral bands of the UAV HS sensor, the HS VIs minimized the spectral redundance, which was usually found in the hyperspectral data and also promoted the computational optimization [29,30]. Therefore, 17 prevalent formulas, composed of two, three, or four spectral bands using function of sum, difference, ratio, double difference, normalized difference, and hybrid, were regarded as the HS VIs and employed in the study [31]; the 17 formulas are shown in

Table 6. Band combinations of the UAV-derived hyperspectral images used in the study.

Number	Band Combination	Number	Band Combination
1	$b_1+b_2$	10	$\left(b_1-b_2\right)/b_3$
2	$b_1-b_2$	11	$\left(b_1+b_2\right)/b_3$
3	$b_1/b_2$	12	$\left(b_1+b_2\right)/\left(b_3+b_2\right)$
4	$b_1/\left(b_1+b_2\right)$	13	$\left(b_1-b_2\right)/\left(b_3+b_2\right)$
5	$b_1/\left(b_1-b_2\right)$	14	$\left(b_1+b_2\right)/\left(b_3-b_2\right)$
6	$\left(b_1-b_2\right)/\left(b_1+b_2\right)$	15	$\left(b_1-b_2\right)/\left(b_3-b_2\right)$
7	$b_1+b_2-b_3$	16	$\left(b_1-b_2\right)+\left(b_3-b_4\right)$
8	$b_1/\left(b_2+b_3\right)$	17	$\left(b_1-b_2\right)/\left(b_3-b_4\right)$
9	$b_1/\left(b_2-b_3\right)$

Note: ‘$b_1$’, ‘$b_2$’, ‘$b_3$’, and ‘$b_4$’: the average value of each band of the UAV-derived hyperspectral images, respectively.

. Band iteration was applied to all hyperspectral bands in each formulation.

Table 5. Multispectral vegetation indices used in the study.

Multispectral Vegetation Index	Formulation	Reference
Difference vegetation index	$\mathrm{DVI}=\mathrm{G}-\mathrm{B}$	[32]
Ratio vegetation index	$\mathrm{RVI}=\frac{\mathrm{NIR}}{\mathrm{R}}$	[33]
Green chlorophyll index	$\mathrm{GCI}=\frac{\mathrm{NIR}}{\mathrm{G}}-1$	[34]
Red-edge chlorophyll index	$\mathrm{RECI}=\frac{\mathrm{NIR}}{\mathrm{RE}}-1$	[34]
Normalized difference vegetation index	$\mathrm{NDVI}=\frac{\mathrm{NIR}-\mathrm{R}}{\mathrm{NIR}+\mathrm{R}}$	[35]
Green normalized difference vegetation index	$\mathrm{GNDVI}=\frac{\mathrm{NIR}-\mathrm{G}}{\mathrm{NIR}+\mathrm{G}}$	[36]
Green-red vegetation index	$\mathrm{GRVI}=\frac{\mathrm{G}-\mathrm{R}}{\mathrm{G}+\mathrm{R}}$	[33]
Green-blue vegetation index	$\mathrm{GBVI}=\frac{\mathrm{G}-\mathrm{B}}{\mathrm{G}+\mathrm{B}}$	[37]
Normalized difference red-edge	$\mathrm{NDRE}=\frac{\mathrm{NIR}-\mathrm{RE}}{\mathrm{NIR}+\mathrm{RE}}$	[38]
Normalized difference re-edge index	$\mathrm{NDREI}=\frac{\mathrm{RE}-\mathrm{G}}{\mathrm{RE}+\mathrm{G}}$	[39]
Simplified canopy chlorophyll content index	$\mathrm{SCCCI}=\frac{\mathrm{NDRE}}{\mathrm{NDVI}}$	[40]
Enhanced vegetation index	$\mathrm{EVI}=2.5\times \frac{\mathrm{NIR}-\mathrm{R}}{1+\mathrm{NIR}-2.4\times \mathrm{R}}$	[41]
Two-band enhanced vegetation index	$\mathrm{EVI}2=2.5\times \frac{\mathrm{NIR}-\mathrm{R}}{\mathrm{NIR}+2.4\times \mathrm{R}+1}$	[42]
Optimized soil adjusted vegetation index	$\mathrm{OSAVI}=\frac{\mathrm{NIR}-\mathrm{R}}{\mathrm{NIR}-\mathrm{R}+\mathrm{L}} \left(\mathrm{L}=0.16\right)$	[43]
Modified chlorophyll absorption in reflectance index	$\mathrm{MCARI}=\left[\left(\mathrm{RE}-\mathrm{R}\right)-0.2\times \left(\mathrm{RE}-\mathrm{G}\right)\right]\times \frac{\mathrm{RE}}{\mathrm{R}}$	[44]
Transformed chlorophyll absorption in reflectance index	$\mathrm{TCARI}=3\times \left[\left(\mathrm{RE}-\mathrm{R}\right)-0.2\times \left(\mathrm{RE}-\mathrm{G}\right)\times \frac{\mathrm{RE}}{\mathrm{R}}\right]$	[45]
MCARI/OSAVI	MCARI/OSAVI	[44]
TACRI/OSAVI	TACRI/OSAVI	[45]
Wide dynamic range vegetation index	$\mathrm{WDRVI}=\frac{\mathrm{a}\times \mathrm{NIR}-\mathrm{R}}{\mathrm{a}\times \mathrm{NIR}+\mathrm{R}} \left(\mathrm{a}=0.12\right)$	[46]

Note: ‘R’, ‘G’, ‘B’, ‘RE’, and ‘NIR’: the average value of the red, green, blue, red-edge, and near-infrared bands of the UAV-derived multispectral images, respectively.

Table 5. Multispectral vegetation indices used in the study.

Multispectral Vegetation Index	Formulation	Reference
Difference vegetation index	$\mathrm{DVI}=\mathrm{G}-\mathrm{B}$	[32]
Ratio vegetation index	$\mathrm{RVI}=\frac{\mathrm{NIR}}{\mathrm{R}}$	[33]
Green chlorophyll index	$\mathrm{GCI}=\frac{\mathrm{NIR}}{\mathrm{G}}-1$	[34]
Red-edge chlorophyll index	$\mathrm{RECI}=\frac{\mathrm{NIR}}{\mathrm{RE}}-1$	[34]
Normalized difference vegetation index	$\mathrm{NDVI}=\frac{\mathrm{NIR}-\mathrm{R}}{\mathrm{NIR}+\mathrm{R}}$	[35]
Green normalized difference vegetation index	$\mathrm{GNDVI}=\frac{\mathrm{NIR}-\mathrm{G}}{\mathrm{NIR}+\mathrm{G}}$	[36]
Green-red vegetation index	$\mathrm{GRVI}=\frac{\mathrm{G}-\mathrm{R}}{\mathrm{G}+\mathrm{R}}$	[33]
Green-blue vegetation index	$\mathrm{GBVI}=\frac{\mathrm{G}-\mathrm{B}}{\mathrm{G}+\mathrm{B}}$	[37]
Normalized difference red-edge	$\mathrm{NDRE}=\frac{\mathrm{NIR}-\mathrm{RE}}{\mathrm{NIR}+\mathrm{RE}}$	[38]
Normalized difference re-edge index	$\mathrm{NDREI}=\frac{\mathrm{RE}-\mathrm{G}}{\mathrm{RE}+\mathrm{G}}$	[39]
Simplified canopy chlorophyll content index	$\mathrm{SCCCI}=\frac{\mathrm{NDRE}}{\mathrm{NDVI}}$	[40]
Enhanced vegetation index	$\mathrm{EVI}=2.5\times \frac{\mathrm{NIR}-\mathrm{R}}{1+\mathrm{NIR}-2.4\times \mathrm{R}}$	[41]
Two-band enhanced vegetation index	$\mathrm{EVI}2=2.5\times \frac{\mathrm{NIR}-\mathrm{R}}{\mathrm{NIR}+2.4\times \mathrm{R}+1}$	[42]
Optimized soil adjusted vegetation index	$\mathrm{OSAVI}=\frac{\mathrm{NIR}-\mathrm{R}}{\mathrm{NIR}-\mathrm{R}+\mathrm{L}} \left(\mathrm{L}=0.16\right)$	[43]
Modified chlorophyll absorption in reflectance index	$\mathrm{MCARI}=\left[\left(\mathrm{RE}-\mathrm{R}\right)-0.2\times \left(\mathrm{RE}-\mathrm{G}\right)\right]\times \frac{\mathrm{RE}}{\mathrm{R}}$	[44]
Transformed chlorophyll absorption in reflectance index	$\mathrm{TCARI}=3\times \left[\left(\mathrm{RE}-\mathrm{R}\right)-0.2\times \left(\mathrm{RE}-\mathrm{G}\right)\times \frac{\mathrm{RE}}{\mathrm{R}}\right]$	[45]
MCARI/OSAVI	MCARI/OSAVI	[44]
TACRI/OSAVI	TACRI/OSAVI	[45]
Wide dynamic range vegetation index	$\mathrm{WDRVI}=\frac{\mathrm{a}\times \mathrm{NIR}-\mathrm{R}}{\mathrm{a}\times \mathrm{NIR}+\mathrm{R}} \left(\mathrm{a}=0.12\right)$	[46]

Note: ‘R’, ‘G’, ‘B’, ‘RE’, and ‘NIR’: the average value of the red, green, blue, red-edge, and near-infrared bands of the UAV-derived multispectral images, respectively.

2.4. Yield Estimation Model

Two repeats (Repeat 1 and Repeat 2) were employed for all cultivars to construct the winter wheat yield estimation model from each vegetation index. Pearson’s correlation coefficient (r) was used to demonstrate the yield estimation accuracy by comparing the yield measured in the field and estimated yield from the regression models below using Student’s t-test at a 95% confidence level. After excluding the outliers described in Section 2.2.1, 144 samples were used to establish the yield estimation model in the study.

The most commonly used remote-sensing based approaches for crop yield estimation involve the use of empirical statistical models, which demonstrate the relationship between yield and canopy spectrum characteristics in an intuitive way. Therefore, the simple linear regression function (SLR) was used to analyze the relationship between winter wheat yield and vegetation indices in individual growth stages (Equation (1)).

$$y=a\times x+b_1$$

(1)

where $y$ represents the winter wheat yield; $x$ represents the vegetation index values at the booting, heading, flowering, filling, or maturation stages; and $a$ and $b_1$ are parameters calculated from a least-squares fitting method.

In addition, the multiple linear regression function (MLR), which was used to combine the vegetation indices among the multiple growth stages to estimate crop yield, is presented in Equation (2).

$$y=a_1\times x_1+a_2\times x_2+a_3\times x_3+a_4\times x_4+a_5\times x_5+b_2$$

(2)

In this equation, $y$ represents the winter wheat yield; $x_1$, $x_2$, $x_3$, $x_4$, and $x_5$ represent the vegetation index values at the booting, heading, flowering, filling, and maturation stages, respectively; and $a_1$, $a_2$, $a_3$, $a_4$, $a_5$, and $b_2$ are parameters calculated from a least-squares fitting method.

2.5. Validation of the Crop Yield Estimation Model

After eliminating outliers, 70 samples of Repeat 3 for all cultivars were used to validate the regression model using root mean square error (RMSE) and mean absolute percentage error (MAPE). The RMSE and MAPE equations are presented in Equations (3) and (4).

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

(3)

$$MAPE=\frac{100\%}{n}{\displaystyle \sum }_{i=1}^n\left|\frac{{\hat{y}}_i-y_i}{y_i}\right|$$

(4)

In these equations, $y_i$ is the crop yield measured in the field of sample i, ${\hat{y}}_i$ is the crop yield predicted by the estimation models of sample i, and n is the number of valid samples.

3. Results

3.1. Winter Wheat Yield Estimation using Vegetation Indices from Individual Growth Stages

The relationships between winter wheat yield and the 19 multispectral vegetation indices at five different growth stages are displayed based on correlation coefficients in

Table 7. Correlation coefficient (r) and RMSE between winter wheat yield and each multispectral vegetation index at the five different growth stages (n = 144).

Multispectral Vegetation Index	Growth Stages
	Booting		Heading		Flowering		Filling		Maturation
	r	RMSE (t/ha)	r	RMSE (t/ha)	r	RMSE (t/ha)	r	RMSE (t/ha)	r	RMSE (t/ha)
DVI	0.26	1.21	0.31	1.18	0.12	1.24	0.48	1.10	0.43	1.13
RVI	0.36	1.16	0.44	1.12	0.56	1.03	0.61	0.98	0.44	1.12
GCI	0.07	1.24	0.25	1.21	0.63	0.97	0.56	1.03	0.34	1.17
RECI	0.01	1.25	0.35	1.17	0.67	0.93	0.66	0.93	0.33	1.18
NDVI	0.37	1.16	0.46	1.11	0.57	1.02	0.64	0.96	0.49	1.09
GNDVI	0.10	1.24	0.28	1.20	0.64	0.96	0.61	0.98	0.36	1.16
GRVI	0.14	1.24	0.21	1.22	0.24	1.21	0.49	1.09	0.49	1.09
GBVI	0.21	1.22	0.22	1.22	0.26	1.20	0.39	1.15	0.39	1.15
NDRE	0.00	1.25	0.36	1.16	0.67	0.92	0.67	0.93	0.33	1.18
NDREI	0.16	1.23	0.14	1.24	0.48	1.10	0.52	1.07	0.28	1.20
SCCCI	0.09	1.24	0.26	1.20	0.48	1.09	0.37	1.16	0.16	1.23
EVI	0.38	1.15	0.47	1.10	0.16	1.23	0.12	1.24	0.04	1.25
EVI2	0.37	1.16	0.45	1.11	0.57	1.02	0.63	0.96	0.47	1.10
OSAVI	0.16	1.24	0.33	1.20	0.44	1.13	0.54	1.06	0.16	1.23
MCARI	0.21	1.22	0.14	1.23	0.20	1.22	0.48	1.10	0.47	1.10
TCARI	0.03	1.25	0.09	1.24	0.08	1.24	0.19	1.22	0.53	1.06
MCARI/OSAVI	0.21	1.22	0.14	1.23	0.20	1.22	0.48	1.10	0.47	1.10
TACRI/OSAVI	0.03	1.25	0.09	1.24	0.08	1.24	0.19	1.22	0.53	1.06
WDRVI	0.37	1.16	0.45	1.12	0.57	1.03	0.63	0.97	0.45	1.11

Note: The best results in terms of the r and RMSE values derived from the use of the simple linear regression functions for each growth stage are in bold typeface.

. According to Table 7, although the vegetation indices showed various relationships with yield at different growth stages, there was no multispectral vegetation index that significantly correlated with winter wheat yield at a single stage (r > 0.70). Most multispectral vegetation indices showed lower correlations with winter wheat yield at the booting, heading, and maturation stages than at the flowering and filling stages. The ‘NDRE’ vegetation index at the flowering stage had the best correlation with winter wheat yield (r = 0.67; RMSE = 0.92 t/ha) (Figure 4a).

The highest correlation with winter wheat yield at each individual growth stage for each band combination type of the UAV hyperspectral images are shown in

Table 8. The highest correlation coefficient between winter wheat yield and each band combination type of hyperspectral sensor images at five different growth stages (n = 144).

Band Combination	Growth Stages
	Booting		Heading		Flowering		Filling		Maturation
	r	RMSE (t/ha)	r	RMSE (t/ha)	r	RMSE (t/ha)	r	RMSE (t/ha)	r	RMSE (t/ha)
$b_1+b_2$	0.45	1.11	0.31	1.18	0.67	0.93	0.60	1.00	0.44	1.12
$b_1-b_2$	0.57	1.02	0.61	0.99	0.73 **	0.86	0.75 **	0.83	0.66	0.93
$b_1/b_2$	0.57	1.03	0.61	0.99	0.75 **	0.82	0.75 **	0.83	0.65	0.95
$b_1/\left(b_1+b_2\right)$	0.57	1.03	0.61	0.99	0.75 **	0.82	0.74 **	0.83	0.65	0.95
$b_1/\left(b_1-b_2\right)$	0.56	1.04	0.56	1.03	0.74 **	0.84	0.70	0.90	0.68	0.91
$\left(b_1-b_2\right)/\left(b_1+b_2\right)$	0.57	1.03	0.61	0.99	0.75 **	0.82	0.74 **	0.83	0.65	0.95
$b_1+b_2-b_3$	0.60	1.00	0.65	0.95	0.76 **	0.81	0.70	0.89	0.70	0.88
$b_1/\left(b_2+b_3\right)$	0.61	0.99	0.67	0.93	0.79 **	0.77	0.76 **	0.82	0.72 **	0.86
${\mathrm{b}}_1/\left({\mathrm{b}}_2-{\mathrm{b}}_3\right)$	0.60	0.99	0.65	0.95	0.76 **	0.80	0.74 **	0.84	0.73 **	0.85
$\left(b_1-b_2\right)/b_3$	0.60	1.00	0.65	0.95	0.78 **	0.78	0.75 **	0.83	0.73 **	0.86
$\left(b_1+b_2\right)/b_3$	0.61	0.99	0.67	0.93	0.78 **	0.77	0.76 **	0.82	0.73 **	0.86
$\left(b_1+b_2\right)/\left(b_3+b_2\right)$	0.60	1.00	0.65	0.95	0.78 **	0.79	0.75 **	0.83	0.72 **	0.86
$\left(b_1-b_2\right)/\left(b_3+b_2\right)$	0.60	1.00	0.65	0.95	0.78 **	0.79	0.75 **	0.83	0.72 **	0.86
$\left(b_1+b_2\right)/\left(b_3-b_2\right)$	0.60	0.99	0.65	0.95	0.77 **	0.80	0.73 **	0.85	0.73 **	0.85
$\left(b_1-b_2\right)/\left(b_3-b_2\right)$	0.57	1.02	0.61	0.99	0.76 **	0.80	0.77 **	0.79	0.68	0.91
$\left(b_1-b_2\right)+\left(b_3-b_4\right)$	0.64	0.96	0.66	0.94	0.79 **	0.77	0.76 **	0.81	0.73 **	0.85
$\left(b_1-b_2\right)/\left(b_3-b_4\right)$	0.65	0.95	0.68	0.92	0.80 **	0.75	0.78 **	0.78	0.75 **	0.83

Note: ‘$b_1$’, ‘$b_2$’, ‘$b_3$’, and ‘$b_4$’: the average value of each band of the UAV hyperspectral images, respectively. ** indicates significance at p-value < 0.01. The best results in terms of the r and RMSE values derived from the use of the multiple linear regression functions for each growth stage are in bold typeface.

. Similar to the multispectral vegetation indices, the correlation coefficients between winter wheat yield and the hyperspectral vegetation indices were higher at the flowering and filling stages than at the booting, heading, and maturation stages. No VIs showed a significant correlation with grain yield at the early growth stages of booting and heading. Unlike the multispectral vegetation indices, the band combinations of the hyperspectral images at the flowering stage exhibited more significant correlations with yield than that at the filling stage. Moreover, all VIs, excluding ‘$b_1+b_2$’, were significantly correlated with yield at the flowering stage (r > 0.70). The band combination type of ‘$\left(b_1-b_2\right)/\left(b_3-b_4\right)$’ showed the best correlation with grain yield throughout the whole growth period among all VIs. Hence, the best performance in terms of yield estimation with hyperspectral-based VIs, according to the simple linear regression model used in the present study, was achieved by the VI of ‘$\left(b_1-b_2\right)/\left(b_3-b_4\right)$’ at the flowering stage, with an r value of 0.80 and RMSE of 0.75 t/ha (Figure 4b), and the corresponding central wavelengths of ‘$b_1$’, ‘$b_2$’, ‘$b_3$’, and ‘$b_4$’ were 838 nm, 722 nm, 710 nm, and 970 nm.

3.2. Winter Wheat Yield Estimation with Vegetation Indices Combining Multiple Growth Stages

The multiple linear regression models (Equation (2)) were employed to further investigate the relationship between winter wheat yield and the multi-temporal spectral vegetation indices.

The correlation coefficients between winter wheat yield and the multispectral vegetation indices at the five growth stages were determined via the use of the multiple linear regression models, and the results are presented in Figure 5. The results show that the vegetation indices combining multiple growth stages presented a higher correlation than the indices did at the individual growth stages. According to Figure 5, the MLR models combining multispectral vegetation indices of ‘DVI’, ‘RECI’, ‘GBVI’, ‘NDRE’, ‘MCARI’, ‘TCARI’, ‘MCARI/OSAVI’, and ‘TCARI/OSAVI’ across the five growth stages showed a significant correlation with winter wheat yield (r > 0.70). The optimal fit of the MLR model for winter wheat yield estimation based on the multispectral vegetation indices in this study was found at the index of ‘NDRE’ (r = 0.75; RMSE = 0.83 t/ha) (Figure 6a).

The best correlation coefficients of each band combination type of the UAV hyperspectral images including the five growth stages with winter wheat yield are shown in Figure 7. Comparing Table 8 and Figure 7, all band combination types demonstrated better regression accuracy of grain yield estimation based on multiple growth stages compared to when the growth stages were only considered singularly. According to Figure 7, all hyperspectral band combinations except ‘$b_1+b_2$’ presented a significant correlation with winter wheat yield when combining the five growth stages. The optimal MLR model for winter wheat yield estimation was based on the combination type of ‘$\left(b_1-b_2\right)/\left(b_3-b_4\right)$’, with r of 0.84 and RMSE of 0.69 t/ha (Figure 6b), which was significantly higher than the optimal MLR model of the multispectral vegetation index above. The central wavelengths of the best hyperspectral band combination, ‘$\left(b_1-b_2\right)/\left(b_3-b_4\right)$’, were 782 nm, 874 nm, 762 nm, and 890 nm, respectively.

3.3. Validation of the Regression Models for Winter Wheat Yield Estimation

The validation of the regression models for winter wheat yield estimation was conducted by using the independent dataset of ‘Repeat 3’ (n = 70).

3.3.1. Validation of the Simple Linear Regression Model at a Single Growth Stage

The RMSE and MAPE of the simple linear regression models at individual growth stages based on the multispectral vegetation indices and hyperspectral band combination types are shown in

Table 9. RMSE and MAPE for winter wheat yield estimation with different multispectral vegetation indices at each individual growth stage (n = 70).

Multispectral Vegetation Index	Growth Stages
	Booting		Heading		Flowering		Filling		Maturation
	RMSE (t/ha)	MAPE	RMSE (t/ha)	MAPE	RMSE (t/ha)	MAPE	RMSE (t/ha)	MAPE	RMSE (t/ha)	MAPE
DVI	1.35	14.09%	1.39	14.62%	1.31	13.44%	1.15	11.15%	1.36	13.94%
RVI	1.24	13.28%	1.10	11.43%	1.00	9.93%	0.97	9.99%	1.17	12.12%
GCI	1.34	13.87%	1.22	12.71%	0.96	9.34%	1.03	10.52%	1.22	12.44%
RECI	1.35	13.94%	1.21	12.45%	0.91	8.58%	0.87	8.69%	1.19	12.05%
NDVI	1.24	13.22%	1.08	11.32%	0.99	9.91%	0.95	9.74%	1.15	12.01%
GNDVI	1.33	13.83%	1.20	12.48%	0.95	9.23%	0.96	9.78%	1.21	12.36%
GRVI	1.36	14.08%	1.37	14.37%	1.24	12.91%	1.09	11.33%	1.23	12.93%
GBVI	1.36	14.10%	1.40	14.58%	1.25	12.72%	1.19	11.82%	1.36	13.92%
NDRE	1.35	13.93%	1.20	12.38%	0.84	8.38%	0.85	8.47%	1.18	12.00%
NDREI	1.34	13.88%	1.28	13.30%	1.17	11.97%	1.12	11.52%	1.31	13.57%
SCCCI	1.38	14.11%	1.28	13.13%	1.13	11.22%	1.13	11.54%	1.35	13.96%
EVI	1.23	13.11%	1.10	11.37%	1.37	13.94%	1.35	13.93%	1.35	13.93%
EVI2	1.24	13.24%	1.09	11.33%	0.99	9.90%	0.95	9.79%	1.16	12.03%
OSAVI	1.34	14.35%	1.22	12.63%	1.12	11.34%	1.03	10.55%	1.35	13.95%
MCARI	1.39	14.46%	1.33	13.93%	1.30	13.60%	1.16	11.98%	1.20	12.65%
TCARI	1.35	13.89%	1.33	13.67%	1.35	13.92%	1.29	13.55%	1.19	12.90%
MCARI/OSAVI	1.39	14.46%	1.33	13.93%	1.30	13.60%	1.16	11.98%	1.20	12.65%
TACRI/OSAVI	1.35	13.89%	1.33	13.67%	1.35	13.92%	1.29	13.55%	1.19	12.90%
WDRVI	1.24	13.26%	1.09	11.37%	0.99	9.90%	0.96	9.88%	1.17	12.07%

and

Table 10. RMSE and MAPE for winter wheat yield estimation with different band combination types of the UAV-derived hyperspectral sensor images at each individual growth stage (n = 70).

Band Combination	Growth Stages
	Booting		Heading		Flowering		Filling		Maturation
	RMSE (t/ha)	MAPE	RMSE (t/ha)	MAPE	RMSE (t/ha)	MAPE	RMSE (t/ha)	MAPE	RMSE (t/ha)	MAPE
$b_1+b_2$	1.21	12.69%	1.24	13.24%	0.94	9.50%	1.09	11.24%	1.24	12.61%
$b_1-b_2$	1.15	11.43%	1.21	11.94%	0.94	9.20%	0.97	9.80%	1.00	10.34%
$b_1/b_2$	1.27	12.58%	1.10	10.95%	0.90	8.84%	0.97	10.05%	0.95	9.71%
$b_1/\left(b_1+b_2\right)$	1.27	12.59%	1.10	10.94%	0.89	8.77%	0.97	10.04%	0.96	9.65%
$b_1/\left(b_1-b_2\right)$	2.21	17.83%	1.16	11.81%	0.87	8.79%	1.01	10.03%	0.91	9.26%
$\left(b_1-b_2\right)/\left(b_1+b_2\right)$	1.27	12.59%	1.10	10.94%	0.89	8.77%	0.97	10.04%	0.96	9.65%
$b_1+b_2-b_3$	1.09	11.07%	1.10	10.90%	0.83	7.81%	0.97	9.82%	0.93	9.39%
$b_1/\left(b_2+b_3\right)$	1.15	11.55%	1.02	10.29%	0.79	7.38%	0.98	10.11%	0.89	8.79%
${\mathrm{b}}_1/\left({\mathrm{b}}_2-{\mathrm{b}}_3\right)$	1.07	10.95%	1.09	10.81%	0.84	8.36%	0.93	9.13%	0.91	9.18%
$\left(b_1-b_2\right)/b_3$	1.09	11.27%	1.09	10.86%	0.80	7.64%	0.96	9.94%	0.94	9.51%
$\left(b_1+b_2\right)/b_3$	1.15	11.55%	1.02	10.30%	0.79	7.68%	0.98	10.11%	0.88	8.69%
$\left(b_1+b_2\right)/\left(b_3+b_2\right)$	1.08	11.25%	1.12	10.99%	0.84	7.87%	0.97	9.99%	0.93	9.35%
$\left(b_1-b_2\right)/\left(b_3+b_2\right)$	1.08	11.25%	1.12	10.99%	0.84	7.87%	0.97	9.99%	0.93	9.35%
$\left(b_1+b_2\right)/\left(b_3-b_2\right)$	1.08	11.12%	1.11	10.92%	0.84	8.40%	0.90	9.24%	0.91	9.03%
$\left(b_1-b_2\right)/\left(b_3-b_2\right)$	1.27	12.55%	1.10	10.92%	0.91	9.00%	0.83	8.08%	0.98	9.88%
$\left(b_1-b_2\right)+\left(b_3-b_4\right)$	1.10	11.30%	1.05	10.48%	0.79	7.40%	0.88	8.94%	0.95	10.02%
$\left(b_1-b_2\right)/\left(b_3-b_4\right)$	0.97	9.93%	1.02	10.10%	0.78	7.24%	0.84	8.47%	0.89	8.92%

, respectively. The VI of ‘NDRE’ derived from the UAV-derived multispectral images at the flowering stage, which achieved the best correlation with grain yield, showed the lowest RMSE of 0.84 t/ha and a MAPE of 8.38%. For the VIs calculated based on the UAV-derived hyperspectral images, the best RMSE (0.78 t/ha) and MAPE (7.24%) were achieved by ‘$\left(b_1-b_2\right)/\left(b_3-b_4\right)$’ at the flowering stage. Moreover, these two indices achieved the best correlation with crop yield for the UAV-derived multispectral and hyperspectral images at all individual growth stages considered in this study, respectively. The validation of the estimated winter wheat yield (estimated via the use of two simple linear regression models) based on these two indices at the flowering stage is illustrated in Figure 8.

3.3.2. Validation of the Multiple Linear Regression Models for Winter Wheat Yield Estimation Combining Five Different Growth Stages

The RMSE of the multiple linear regression model for winter wheat yield estimation based on UAV-calculated multispectral and hyperspectral vegetation indices combining five growth stages and the corresponding MAPE are shown in Figure 9. According to Figure 9a, the lowest MAPE for the multispectral VIs was 8.44% (found at ‘NDRE’; RMSE of 0.90 t/ha). For the hyperspectral VIs, the lowest MAPE was 6.56%, which was achieved by ‘$\left(b_1-b_2\right)/\left(b_3-b_4\right)$’, with an RMSE of 0.70 t/ha (Figure 9b). The validation of the estimated winter wheat yield (derived from the use of two multiple linear regression models) based on these two indices is illustrated in Figure 10. Similar to the validation of the simple linear regression models at the individual growth stages, the RMSE and MAPE were lower for the multiple linear regression model based on the hyperspectral vegetation index than that based on the multispectral vegetation index. Moreover, comparing the regression models constructed only by the vegetation indices of an individual growth stage, the regression models which combined multi-temporal vegetation indices information demonstrated their robustness in winter yield estimation. The multiple linear regression model based on the band combination of ‘$\left(b_1-b_2\right)/\left(b_3-b_4\right)$’ from the hyperspectral sensor performed reasonably well (achieving the highest correlation coefficient and lowest RMSE and MAPE for winter wheat yield estimation) and could be regarded as the best grain yield estimation model independent of the photography conditions in this study.

4. Discussion

Currently, NDVI is the most widely used vegetation index for crop yield estimation. However, because of the characteristics of NDVI, it has significant saturation under a high vegetation coverage level, thereby affecting estimation accuracy [47,48,49,50]. According to the spectral reflectance characteristics of the plant, the absorption of chlorophyll on the red-edge waveband is weaker than that of red band, with the red-edge region having stronger transmission ability with the crop canopy [51]. The use of a red-edge band rather than a red band in the vegetation index of NDVI can reduce the saturation phenomenon, improving crop yield estimation accuracy [51,52]. Additionally, the spectrum of each irrigation group for the winter wheat cultivator of ‘C9’ at the flowering stage is shown in Figure 11. According to the figure, there were significant differences among different irrigation groups at the near-infrared band. Therefore, in this study, vegetation indices composed of the red-edge and near-infrared bands, both for MS imagery (NDRE) and HS imagery (‘$\left(b_1-b_2\right)/\left(b_3-b_4\right)$’), demonstrated reasonable robustness in terms of their winter wheat yield estimation.

Remote crop yield estimation methods are commonly based on the high correlation between the crop yield and the vegetation index taken at a specific crop growth stage [50]. The success of the development of a vegetation index is dependent on the use of bands with different sensitivities to the key parameter that is to be monitored [51]. The potential of multi-spectrum and hyper-spectrum for winter wheat yield estimation was systematically compared in this study. The lower yield estimation accuracy based on multispectral vegetation indices is mostly due to the obvious shortcoming of the limited fixed bands with a wide resolution. The hyperspectral sensor captured much richer information and is more sensitive to crop canopy characteristics with the continuous acquisition of reflectance at narrow wavelengths [45,52]. Additionally, for all vegetation indices calculated from the UAV-hyperspectral imagery used in this study, with the number of bands including in the vegetation index, the yield estimation accuracy of winter wheat assumed a rising tendency. According to the study of Thenkabail et al. [53], four sensitive band combination-based optimum multiple narrow band reflectance models could explain up to 92% of the crop biophysical parameter variability. Hence, the four-band combination type based on hyperspectral imagery achieved a more significant correlation with winter wheat yield than the vegetation indices derived from multispectral imagery throughout the whole growth period of winter wheat in this study.

According to Qader et al., yield estimation models that use VIs from the crop’s critical growth stage can obtain a higher accuracy across remote sensing data [54]. Our results regarding winter wheat yield estimation based on singular growth stages confirmed this and strongly indicated that the flowering stage was a critical period for winter wheat yield estimation. Some studies have shown that the accumulative vegetation index can improve the stability of yield estimation and that adding one or more growth stages (except the critical growth stage) could improve estimation accuracy [28,55,56]. In this study, the correlation coefficient between crop yield and the hyperspectral vegetation index of ‘$\left(b_1-b_2\right)/\left(b_3-b_4\right)$’ increased from 0.80 to 0.84, and RMSE decreased from 0.75 t/ha to 0.69 t/ha when the booting, heading, filling, and maturation stages were added with the MLR model in the flowering stage of the simple linear model.

Although the MLR model combining multiple temporal hyperspectral vegetation indices calculated according to hyperspectral imagery acquired good yield estimation accuracy for winter wheat, machine learning algorithms have demonstrated the potential to retrieve crop characteristics using multispectral satellite data, aerial multispectral data, aerial hyperspectral data, and so on [57,58,59]. Therefore, future research should aim to explore machine learning regression models to strengthen crop estimation ability via multiple temporal UAV-derived hyperspectral datasets. Additionally, data fusion approaches which could integrate UAV data and satellite data-based vegetation index time-series curves together and improve temporal resolutions will also be investigated in the context of estimating crop yield in a future study.

5. Conclusions

This study assessed the accuracy of winter wheat yield estimation values based on UAV-derived multispectral and hyperspectral images from single and multiple growth stages. The results suggested that the proposed multiple linear regression model, constructed by the vegetation index of ‘$\left(b_1-b_2\right)/\left(b_3-b_4\right)$’ with central wavelengths of 782 nm, 874 nm, 762 nm, 890 nm, which were calculated from UAV-based hyperspectral images using the growth stages from booting to maturation, can be used as a fast and reliable method for winter wheat yield estimation to contribute to the breeding of drought-resistant varieties of winter wheat in a field plot scale over a short amount of time. Moreover, the red-edge and near-infrared bands are recommended for use in the context of crop yield estimation.

From a longer-term perspective, more in-depth investigations into crop yield estimation (including rice, maize, soybean, and other crops) based on UAV-mounted hyperspectral datasets (via not only linear regression models but also machine learning algorithms, data assimilation, and so on) are expected.