UAV-Based Multi-Temporal Thermal Imaging to Evaluate Wheat Drought Resistance in Different Deficit Irrigation Regimes

Abstract

Deficit irrigation is a common approach in water-scarce regions to balance productivity and water use, whereas drought stress still occurs to various extents, leading to reduced physiological performance and a decrease in yield. Therefore, seeking a rapid and reliable method to identify wheat varieties with drought resistance can help reduce yield loss under water deficit. In this study, we compared ten wheat varieties under three deficit irrigation systems (W0, no irrigation during the growing season; W1, irrigation at jointing; W2, irrigation at jointing and anthesis). UAV thermal imagery, plant physiological traits [leaf area index (LAI), SPAD, photosynthesis (Pn), transpiration (Tr), stomatal conductance (Cn)], biomass and yield were acquired at different growth stages. Wheat drought resistance performance was evaluated through using the canopy temperature extracted from UAV thermal imagery (CT-UAV), in combination with hierarchical cluster analysis (HCA). The CT-UAV of W0 and W1 treatments was significantly higher than in the W2 treatment, with the ranges of 24.8–33.3 °C, 24.3–31.6 °C, and 24.1–28.9 °C in W0, W1 and W2, respectively. We found negative correlations between CT-UAV and LAI, SPAD, Pn, Tr, Cn and biomass under the W0 (R² = 0.41–0.79) and W1 treatments (R² = 0.22–0.72), but little relevance for W2 treatment. Under the deficit irrigation treatments (W0 and W1), UAV thermal imagery was less effective before the grain-filling stage in evaluating drought resistance. This study demonstrates the potential of ensuring yield and saving irrigation water by identifying suitable wheat varieties for different water-scarce irrigation scenarios.

1. Introduction

Wheat is the main food crop grown in North China, and nearly half of the country’s wheat production comes from the North China Plain [1]. Although wheat plays a critical role in ensuring food security in China, wheat production consumes a lot of water resources, increasing the burden of water consumption [2]. In addition, the current severe groundwater overexploitation in the North China Plain has seriously damaged the local ecological environment [3]. In water-scarce regions, ‘deficit irrigation’ is a commonly used water-saving measure [4,5,6]. Deficit irrigation uses less water than is required for potential evapotranspiration and maximum yield, thereby saving limited irrigation water [7]. Deficit irrigation can increase water use efficiency because it reduces evapotranspiration while maintaining yield [8]. However, while crops are being irrigated in deficit, they will face the threat of drought stress, resulting in a decrease in yield [9]. Since different crop varieties have different sensitivities to water deficit, they will behave different drought resistance performances under different water deficit conditions [10]. Choosing wheat varieties with good drought resistance performance can minimize yield loss while carrying out different degrees of deficit irrigation, which is very important for water shortages such as drought or semi-arid [11].

In the past, the research on the drought resistance performance of crops was based on the investigation of physiological traits and yield, which requires destructive sampling and is time-consuming and labor-intensive [12,13,14]. However, Canopy temperature (CT) can be used to identify changes in physiological traits [15,16], enabling rapid screening of varieties with drought resistance [17]. For decades, canopy temperature has been commonly measured by a handheld infrared analyzer in most research [18,19,20]. The handheld infrared analyzer can quickly obtain the canopy temperature and then evaluate the physiological traits of crops, and it also showed that crop varieties with lower canopy temperature may have better physiological performance under drought stress due to their deep root or osmotic adjustment capacity [21,22,23]. However, using a handheld thermal infrared analyzer to obtain canopy temperature may have two problems; one of them is that as the number of samples collected increases, it remains time-consuming and laborious. Secondly, a long time span of the canopy temperature measurement will easily lead to a decrease in the accuracy of the canopy temperature acquisition and affect the final analysis result. This is because the canopy temperature of plants is easily affected by environmental factors such as ambient temperature [24].

Unmanned aerial vehicles (UAVs) equipped with different types of sensors become a promising tool in various fields of agricultural research, e.g., biomass and grain yield estimation, planting density assessment, pest and disease detection, and crop stress monitoring [25,26,27,28,29,30]. Among them, the thermal infrared cameras onboard unmanned aerial vehicles (UAV) have been increasingly used for crop drought stress detections [28,31,32]. UAV-based thermal images can provide sub-meter spatial resolution, making it possible to accurately obtain canopy temperature over large fields in a short time [33,34,35,36,37,38,39]. Accordingly, crop transpiration can also be predicted from crop canopy temperature to guide irrigation [40]. In addition, UAV thermal infrared images can also be used to predict some physiological traits of crops, such as water use efficiency and biomass [41,42], which brings new insights into the assessment of crop drought resistance under water-deficit and drought conditions. For instance, UAV thermal infrared imagery has been used for evaluating crop stress resistance in drought-stressed conditions [43,44,45], e.g., the canopy temperature extracted from UAV thermal infrared images was used to rapidly identify the drought resistance of different soybean varieties [46]. Moreover, UAV thermal infrared images enabled the evaluation of the stomatal conductance, plant moisture content and biomass of different wheat genotypes under salt stress, as well as to grade them for salt tolerance [47]. Despite that UAV thermal imagery has been applied to various stress monitoring or crop moisture monitoring, little is known about the evaluation of the drought resistance performance of wheat varieties in deficit-irrigation conditions, such as in the North China Plain.

Therefore, the main goal of this study was to test the feasibility of using UAV thermal infrared imaging to identify drought-resistant varieties suitable for different deficit-irrigation regimes in water-scarce regions in China. The specific objectives are (1) to understand the relations of canopy temperature from UAV-thermal imaging with physiological traits and (2) to evaluate the drought resistance performance of ten wheat varieties by combining physiological trait measurements, yield and UAV-based canopy temperature.

2. Materials and Methods

2.1. Study Site and Experimental Design

This experiment was conducted in the Experimental Station of China Agricultural University (37°41′N, 116°37′E) in Yaozhuang Village, Wuqiao County, Hebei Province, from 2019 to 2020 (Figure 1). The site is located in The Heilong Port Basin in the central part of Hebei Province, with an altitude of about 20 m above mean sea level. It has a continental monsoon climate in summer, with an average annual temperature of 12.5 °C, an average annual rainfall of 550 mm and an average rainfall of 124 mm in wheat season. Soil nutrients in the surface layer (0–20 cm) were as follows: organic matter content of 10.4 g·kg⁻¹, the total nitrogen content of 0.98 g·kg⁻¹, alkali hydrolyzed nitrogen of 60.2 mg·kg⁻¹, available phosphorus of 39.5 mg·kg⁻¹, and available potassium of 92 mg·kg⁻¹.

The experiment used a split-plot design, with three irrigation times for the primary plots and fifteen varieties for the subplots. Three irrigation rates (times during the growing season) were applied, including the W0 treatment (no irrigation), W1 treatment (one irrigation at jointing) and W2 treatment (two irrigations at jointing and anthesis). The same amount of water (750 m³/ha) was used for each irrigation. The varieties included the Gaoyou 2018 (GY2018), Hengmai 4399 (H4399), Guanmai 35 (G35), Gaoyou 5766 (GY5766), Shinong 086 (SN086), Shimai 22 (SM22), Jimai 418 (JM418), Zhongmai 1062 (ZM1062), Cangmai 6005 (C6005), Nongda 3636 (ND3636), and five new untested genotypes. The experiment consisted of two replications and in total 90 plots. The plot size was 10 m × 1.5 m, with a row spacing of 15 cm. The experiment was sown on 16 October 2019, and the harvest of each treatment was conducted on the 7th, 9th, and 12th of June 2020, respectively, for W0, W1 and W2. The five untested genotypes (30 plots) did not show unstable physiological characteristics, and thus they were excluded from further data analysis (60 plots). Protective plots were planted between irrigation treatments to avoid the influence of soil water translocation.

2.2. UAV Systems and Thermal Infrared Camera

In this study, a UAV thermal remote sensing system (Figure 2a) was developed with a Matrice M600 Pro (DJI Inc., Shenzhen, China). The total weight of the UAV, including the batteries, is 9.5 kg, and its maximum takeoff weight is 15.5 kg. The UAV can use the GNSS system to locate coordinates in real-time and give information of position data. The UAV was controlled using the DJI Matrice 600 series remote controller and a Xiaomi tablet (Xiaomi Inc., Beijing, China) with the DJI Pilot app (DJI Inc., Shenzhen, China). The UAV requires six charged 4500 mAh batteries for operation. It has a maximum flight time of about 20 min within its payload.

The thermal infrared camera is Zenmuse XT2 (FLIR Systems Inc., Wilsonville, OR, USA). The camera consists of two lenses (Figure 2b): a thermal infrared lens and an RGB lens. So, it can obtain thermal infrared images and RGB images at the same time. The thermal infrared lens has a field of view (FOV) of 25° × 20° and a resolution of 640 × 512 pixels. The camera’s sensor is an uncooled Vanadium Oxide (VOx) microbolometer detector with a detector pitch of 17 μm measuring in the spectral range of 7.5–13 μm. The specified temperature range of the measurement objects is −40 °C to +550 °C and it has a claimed accuracy of ±5 °C and thermal sensitivity of 0.05 °C. In addition, the thermal infrared camera also has a temperature self-calibration function, which can ensure accurate access to the ground temperature information. The RGB lens has a field of view (FOV) of 57° × 42° and a resolution of 3840 × 2160 pixels, supporting both image optimization and digital image enhancement to avoid image distortion.

2.3. Acquisition of UAV Thermal Images and RGB Images

The UAV flew and acquired the thermal infrared images that were acquired from 14:00 to 15:00 at noon local time on 5 May, 11 May, 20 May, 25 May and 2 June 2020. The thermal infrared camera lens was nadir angle to the test plots. The flight altitude was 40 m, and the side overlap and front overlap were set at 90% and 80%, respectively. Under this condition, the ground resolution of thermal infrared images and RGB images is 7.5 cm/pixel and 1.2 cm/pixel, respectively. DJI Pilot app (DJI Inc., v. 1.7.2) was used for UAV route planning, and the flight speed is set at 1.2 m/s, and the time of a single flight mission is about 10 min. The camera was preheated for 15 min before each take-off to reduce systematic error. All missions were conducted under sunny, cloudy-free and windless conditions to minimize the impact of the environment on temperature acquisition.

2.4. UAV Thermal Images Processing

After obtaining the thermal infrared images and RGB images, Pix4D software (Pix4D Inc., Lausanne, Switzerland) was used for mosaic processing. During the photo stitching process, there are ground control points to calibrate the geo-reference of thermal and RGB orthophoto, and the coordinates of these ground control points (Figure 2d) were measured through D-RTK2 mobile GNSS station (DJI Inc., China, Figure 2c). Next, we used a method that combined RGB images and thermal infrared images to avoid the influence of a soil background on canopy temperature extraction [39]. The method had two key steps (Figure 3): First, we used UAV RGB orthophoto to separate the wheat canopy pixels from the soil background. Second, we binarized the segmented wheat cover images for wheat binary maps. In order to better extract canopy pixels, we combined the excess green vegetation index (ExG) with Otsu’s method to set a separate threshold for each plot; and second, we resampled the wheat cover binary maps from 1.2 cm/pixel to 7.5 cm/pixel to match the size of thermal orthophoto through using the nearest-neighbor interpolation algorithm. Finally, the thermal infrared orthophoto was multiplied by the wheat binary maps to produce a new thermal orthophoto in which the soil background has been removed. All of the image processing was completed by MATLAB (v. 2018b). Then, the canopy temperature was extracted from each plot using the open-source GIS software QGIS (v. 3.20.3-Odense). In order to reduce the influence of shadow and brightness, we normalized three bands—R, G and B [48,49,50]. The ExG was calculated through using Equations (1)–(4):

$$\mathrm{Rn}=\frac{\mathrm{R}}{\mathrm{R}+\mathrm{G}+\mathrm{B}}$$

(1)

$$\mathrm{Gn}=\frac{\mathrm{G}}{\mathrm{R}+\mathrm{G}+\mathrm{B}}$$

(2)

$$\mathrm{Bn}=\frac{\mathrm{B}}{\mathrm{R}+\mathrm{G}+\mathrm{B}}$$

(3)

$$\mathrm{ExG}=2\mathrm{Gn}-\mathrm{Rn}-\mathrm{Bn}$$

(4)

where R, G and B represent the raw digital numbers of red, green and blue bands, Rn, Gn and Bn represent the normalized digital numbers of red, green and blue bands.

2.5. Acquisition of CT by A Handheld Thermal Infrared Analyzer

Immediately after the UAV mission, the thermometry measurements were made with an AZ8895 infrared thermometer (Heng Xin Inc., Taiwan, China, Figure 2e). The infrared thermometer had the specified temperature range of the measurement objects is −40 °C to +816 °C, with a temperature resolution of 0.1 °C. To avoid the impact of soil, it scanned a certain area (D:S = 12:1, about 0.0125 m²) of the wheat canopy at an angle of 75° in the horizontal direction and 15 cm above the canopy with the emissivity of 0.95 [51]. Each plot was measured three times in different locations and then calculated its average canopy temperature. Then, the canopy temperature data will be used to calibrate the canopy temperature that was extracted from UAV thermal imagery. The relationship between the canopy temperature extracted from the UAV thermal imagery (CT-UAV) and the canopy temperature obtained by the handheld infrared analyzer (CT-handheld) and the calibration of the canopy temperature of the UAV are shown in Table S1.

2.6. Physiological Parameters Measurements

2.6.1. Photosynthetic Data

Li-6400 photosynthetic instrument was used to measure the net photosynthetic rate (Pn), transpiration rate (Tr), and stomatal conductance (Cn) of flag leaves of wheat on 11 May and 20 May from 8:00 to 11:00 in the morning, with a built-in LED light source with a light intensity of 1500 μmol·m⁻²·s⁻¹ and an airflow rate of 400 μmol·m⁻²·s⁻¹. Three wheat flag leaves with the same size and direction were measured repeatedly in each plot.

2.6.2. SPAD and Leaf Area Index

Chlorophyll content (SPAD) measurements were made using a SPAD-502 Minolta chlorophyll meter (Spectrum Technologies Inc., Aurora, IL, USA) from 11:00 to 12:00 afternoon on 5 May, 11 May, 20 May, 25 May, and 2 June. Five flag leaves of similar size, orientation and shape were selected for measurement in each plot, measured once on both sides of each leaf and again in the middle position, for a total of three times, and the average value was taken.

Leaf area index (LAI) measurements were taken from 16:00 to 17:00 afternoon on 5 May, 11 May, 20 May, 25 May, and 2 June using Li-COR 2200C (LI-COR Inc., Lincoln, NB, USA) canopy analyzer. Five measurements will be made within each plot, with the five sample points arranged in a Z-shape, and the results will be averaged.

2.6.3. Biomass and Yields

After the wheat ripened, the wheat was harvested on 7, 11, and 12 June 2020 under W0, W1 and W2 treatments, respectively. Two subplots of 1 m² (1 m × 1 m area) were collected from each plot, and all the wheat culms in the ground were harvested (excluding the roots). After they were all dried, the dry matter weight was taken as the biomass. After that, the wheat was threshed, and the grain water content was obtained, which was finally converted into 13% water content to be the final grain yield of each variety. Finally, take the average of the two numbers measured in each plot as the biomass and yield of this plot.

2.6.4. The Corresponding Days and Growth Stage after Sowing

The sampling and physiological trait measurement dates are indicated for the days after sowing (DAS) and corresponding growth stages in

Table 1. The sampling or measuring date corresponds to the days after sowing and the growth stages.

Sampling or Measuring Date	The Days after Sowing (DAS)	Growth Stages
5 May 2020	202 days	Flowering
11 May 2020	208 days	6 days after flowering
20 May 2020	217 days	Early stage of grain-filling
25 May 2020	222 days	Middle stage of grain-filling
2 June 2020	230 days	Late stage of grain-filling
7 June 2020	235 days	Ripe stage, W0
9 June 2020	237 days	Ripe stage, W1
12 June 2020	240 days	Ripe stage, W2

.

2.7. Statistical Analysis

In this study, the canopy temperature and physiological traits measured in all different stages were used for principal component analysis to determine important variables and the relationship between canopy temperature and physiological traits.

Hierarchical cluster analysis (HCA) characterizes the similarity between samples by examining the distances between points representing all possible pairs of samples in a high-dimensional space [52]. When using the algorithm to execute, it does not require any information about the number of clusters in advance. Furthermore, it is easy to implement, less complex and tends to give the best results in some cases and can be used to comprehensively rank observations across multiple parameters [53]. In this study, we conducted the Ward’s method hierarchical cluster analysis by computing the Euclidean distance between the observations (rows) of the data matrix. The distance between two cluster centers was calculated through using Equation (5). Next, we clustered the drought resistance of different wheat varieties into groups (good/high, moderate, bad/low) based on the measured physiological traits (SPAD, LAI, Pn, Tr, Cn and biomass), grain yield and CT-UAV across growth stages, respectively, for physiological-trait-based clustering, yield-based clustering and CT-based clustering.

$$\mathrm{d}(\mathrm{x},\mathrm{y})=\sqrt{{{(\mathrm{x}}_1{-\mathrm{y}}_1)}^2{+(\mathrm{x}}_2{-\mathrm{y}}_2{)}^2{+\dots +(\mathrm{x}}_{\mathrm{k}}{-\mathrm{y}}_{\mathrm{k}}{)}^2}$$

(5)

where d(x,y) represents the Euclidean distance, x_k and y_k represent the observed data for two different subjects in k-th variable, and k represents the number of variables.

In addition, SPSS 19.0 (SPSS, Inc., Chicago, IL, USA) was also used to conduct Duncan’s multiple range test to detect physiological differences of wheat under different moisture conditions. In addition, in order to explore the differences between the canopy temperatures under different water deficit treatments, this study used the R language (v. 4.1.2) to perform a multiple-test method on the canopy temperature. In this study, the correlation analysis between canopy temperature and physiological traits is to use Pearson correlation analysis to determine the correlation coefficient between variables. Finally, the canopy temperatures of the different variety groups were compared for their drought resistance performance using the t-test (α = 0.05) and visualized using the R language.

3. Results

3.1. Physiological and Yield Performance of Wheat under Different Water Treatments

Physiological traits of wheat showed significant differences under different water treatments (

Table 2. Physiological traits at five growth stages under different water treatments.

Physiological Traits	Treatment	n	5 May	11 May	20 May	25 May	2 June
SPAD	W0	20	52.3 b	50.4 b	47.8 b	45.5 b	10.7 c
	W1	20	54.2 a	53.3 a	52.8 a	48.5 a	22.1 b
	W2	20	54.6 a	54.3 a	54.4 a	51.5 a	31.4 a
LAI	W0	20	4.7 b	4.5 b	3.6 c	3.1 c	2.7 c
	W1	20	6.5 a	5.5 a	4.5 b	3.8 b	3.2 b
	W2	20	7.1 a	6.2 a	5.0 a	4.2 a	3.7 a
Pn (μmol CO2·m−2·s−1)	W0	20	-	17.3 b	11.4 c	-	-
	W1	20	-	19.7 a	13.9 b	-	-
	W2	20	-	20.7 a	15.6 a	-	-
Tr (mmol H2O·m−2·s−1)	W0	20	-	2.73 b	2.21 c	-	-
	W1	20	-	2.82 ab	2.64 b	-	-
	W2	20	-	3.09 a	3.63 a	-	-
Cn (mmol H2O·m−2·s−1)	W0	20	-	0.153 b	0.099 c	-	-
	W1	20	-	0.160 b	0.156 b	-	-
	W2	20	-	0.227 a	0.208 a	-	-

Duncan’s multiple range test was used to compare the mean values of physiological traits under different water treatments. Within each trait, any two treatments showing different letters are significantly different in the mean values, at p < 0.05.

). First, the physiological performances of wheat such as SPAD, LAI, biomass, Pn, Tr and Cn were better with increasing irrigation, with the best performance under W2 treatment and the worst under W0 treatment. For the W0 treatment, compared with the W2 treatment, these physiological traits decreased significantly in all stages and the decline has been increasing over time, while for the W1 treatment, the decrease in physiological performance was smaller and did not reach a significant level before 11 May.

With the progress of the growth stages, SPAD, LAI, Pn, Tr and Cn further decreased, reaching a significant level compared with the W2 treatment after 20 May, suggesting that drought stress will cause the further deterioration of wheat physiological traits as the growth stage progresses. Figure 4. shows the physiological performance of each wheat variety in different water treatments on 20 May. The physiological performance of almost all wheat cultivars decreased continuously with the intensification of the water deficit. At the same time, there are also great differences in all measured physiological traits in different wheat varieties under the same water treatment, indicating the degree of drought impact of each wheat variety was different and different wheat varieties have different abilities to cope with drought stress.

Table 3. The yield of different varieties under different deficit irrigation systems.

Varieties	Treatment	Yield (kg/ha−1)	Varieties	Treatment	Yield (kg/ha−1)
GY2018	W0	6017.0 c	ND3636	W0	5871.9 c
	W1	7261.4 b		W1	7088.0 b
	W2	7788.4 a		W2	7568.4 a
G35	W0	6520.5 c	H4399	W0	6315.3 c
	W1	7352.0 b		W1	7877.2 b
	W2	7917.8 a		W2	8784.6 a
GY5766	W0	6825.8 c	SN086	W0	5973.6 c
	W1	7170.9 b		W1	7840.1 b
	W2	7840.1 a		W2	8481.8 a
SM22	W0	6629.3 c	JM418	W0	7103.8 c
	W1	7657.3 b		W1	8150.6 b
	W2	8573.0 a		W2	9030.4 a
ZM1062	W0	7045.0 c	C6005	W0	6779.3 c
	W1	7538.4 b		W1	7452.0 b
	W2	8211.0 a		W2	8573.0 a

Duncan’s multiple range test was used to compare the mean values of yield in different water treatments. The means followed by different letters are significantly different at p < 0.05.

shows the yield of different wheat varieties under different deficit irrigation systems. There were significant yield differences among wheat varieties, in W0, W1 and W2 conditions, yield ranges from 5871.9 kg/ha⁻¹ to 7103.8 kg/ha⁻¹, 7088.0 kg/ha⁻¹ to 8150.6 kg/ha⁻¹, and 7568.4 kg/ha⁻¹ to 9030.4 kg/ha⁻¹, respectively. With the decrease in irrigation times, the yield decreased significantly, and the JM418 variety had the highest yield, while the ND3636 variety had the lowest yield in the three deficit irrigation treatments. These results proved that the water deficit affects the physiological performance and ultimately crop yield.

3.2. The HCA Results Based on Physiological Traits and Yield

The ten wheat varieties were artificially divided into three categories (Figure 5) which were the first group (with good drought resistant performance), the second group (with moderate drought resistant performance) and the third group (with poor drought resistant performance). For the W0 treatments, C6005, ZM1062 and JM418 were grouped into the first cluster (green color group) due to their better and similar physiological performance, while the second cluster mainly included SM22 and GY5766, and these two clusters were re-aggregated into one group when the Euclidean distance was 6.7, and the third cluster was SN086, G35, GY2018, H4399 and ND3636. The third cluster and the other two groups were aggregated into one group when the Euclidean distance was 10.9. For the W1 treatment, the first cluster included JM418, SN086, H4399 and SM22, the second cluster included ZM1062, G35, C6005 and GY5766 and the wheat varieties of the final cluster were GY2018 and ND3636. Then, the first and second clusters were re-aggregated into one group when the Euclidean distance was 6.4 and all of the groups were clustered into one group as the Euclidean distance was 10.7. For the W2 treatment, the first cluster included JM418, H4399, SN086 and C6005, the second cluster mainly included ZM1062, GY2018 and SM22, and the third cluster was GY5766, G35 and ND3636. Among them, the third cluster and the second cluster were merged into one larger group when the Euclidean distance was 6.6, and all of the clusters were aggregated into one larger cluster when the Euclidean distance was 8.4.

Table 4. Yield-based hierarchical clustering in different deficit irrigation treatments and the average value of each group was also displayed.

Treatment	Group	The Average Yield (kg/ha−1)	Varieties
W0	High	6970.3 a	JM418, ZM1062, C6005, GY5766
	Moderate	6488.3 b	SM22, G35, H4399
	low	5954.2 c	GY2018, SN086, ND3636
W1	High	7956.0 a	JM418, H4399, SN086
	Moderate	7532.5 b	SM22, ZM1062, C6005
	low	7218.1 c	G35, GY2018, GY5766, ND3636
W2	High	8907.5 a	JM418, H4399
	Moderate	8421.9 b	SN086, ZM1062, C6005
	low	7822.0 c	SM22, G35, GY2018, GY5766, ND3636

Duncan’s multiple range test was used to compare the mean values of yield under different water treatments. The means followed by different letters are significantly different at p < 0.05.

shows the yield clustering results and the centroid value of each group, 10 wheat varieties were also divided into three different yield categories. For the W0 treatment, the high-yield wheat varieties include JM418, ZM1062, C6005 and GY5766, the moderate varieties were SM22, G35 and H4399 and the low-yield varieties were GY2018, SN086 and ND3636. Their centroid values of yield were 6970.3 kg/ha⁻¹, 6488.3 kg/ha⁻¹ and 5954.2 kg/ha⁻¹, respectively. As for the W1 treatment, JM418, H4399 and SN086 wheat varieties belong to the high yield group; the centroid value of yield was 7956.0 kg/ha⁻¹, and the moderate group includes SM22, ZM1062 and C6005, their centroid value was 7532.5 kg/ha⁻¹. The low-yield wheat varieties were G35, GY2018, GY5766 and ND3636; their centroid value was 7218.1 kg/ha⁻¹. For the W2 treatment, the high-yield group included JM418 and H4399, and the centroid value of yield was 8907.5 kg/ha⁻¹. The moderate group included SN086, ZM1062 and C6005, and their centroid value was 8421.9 kg/ha⁻¹. The low-yield group included SM22, G35, GY2018, GY5766 and ND3636, and the centroid value of this group was 7822.0 kg/ha⁻¹. The centroid yield values of the three groups differed significantly between the three deficit irrigation treatments.

The drought resistance clustering based on physiological traits was similar to the clustering results based on yield, with similarities of 80% (two varieties with inconsistent classification), 70% (three varieties with inconsistent classification) and 60% (four varieties with inconsistent classification) in W0, W1 and W2 treatments, respectively. Results showed a strong link between drought resistance and yield.

3.3. Canopy Temperature from Different Treatments and Multiple Growth Stages

Figure 6 shows the distribution of canopy temperature among treatments in different growth stages. There were significant differences in canopy temperature among the three treatments. In the five growth stages, the canopy temperature ranged from 24.8 °C to 33.3 °C under the W0 treatment, from 24.3 °C to 31.6 °C for the W1 treatment, and from 24.1 °C and 28.9 °C for the W2 treatment. The canopy temperature of all treatments increased continuously with the growth process (except on 11 May, there was precipitation before), and showed the maximum on 2 June, which were 33.3 °C, 31.6 °C and 28.9 °C, for W0, W1 and W2, respectively. W0 treatment had the highest canopy temperature and was 4.7 °C higher than W2 on 20 May, indicating that W0 treatment was most seriously affected by water deficit. W1 treatment differed greatly from W2 and showed the largest difference of 2.6 °C on 2 June.

Figure 7 shows the canopy temperature difference between varieties on 20 May. For the W0 treatment, there were apparent differences in canopy temperature among wheat varieties, in which JM418 and ZM1062 had lower canopy temperatures of 31.4 °C and 31.8 °C, respectively, while SN086 and ND3636 had higher canopy temperatures, 33.5 °C and 33.6 °C. For the W1 treatment, the canopy temperature differences among wheat varieties were still large, in which JM418, H4399 and SN086 all had lower canopy temperatures. GY5766 and ND3636 had higher canopy temperatures, 29.7 °C and 30.2 °C. For the W2 treatment, JM418 varieties had a lower canopy temperature that was 27.4 °C and GY5766, GY2018 and SM22 varieties had a higher canopy temperature of 28.5 °C. Moreover, in W0, W1 and W2 treatments, the maximum canopy temperature difference among wheat varieties was 2.2 °C, 1.9 °C and 1.1 °C, respectively.

3.4. The PCA of Physiological Traits and CT-UAV in Different Deficit Irrigation Treatments

According to the scree plot (Supplementary Figure S1), the first two principal components were selected for analysis for both W0 and W1 treatments, because the first two principal components in them could explain 87.4% and 83.3% of the data variance, respectively. For the W0 treatment, the first principal component explains 79.1% of the variance, while the second principal component explains 8.3% of the variance (Figure 8a). For the W1 treatment, the first and second principal components can contribute 71.6% and 11.7% of data variance (Figure 8b). Among them, for the canopy temperature part, CT-UAV has more contributions in the first principal component and the total contribution of CT-UAV was still relatively high at any growth stage.

In addition, there was a negative correlation between canopy temperature and physiological traits such as LAI, SPAD, Pn, Tr, Cn, biomass and yield by PCA under both water deficit treatments. At the same time, the first principal component can explain more than 70% of the data variation in all water stress treatments. In contrast, for the W2 treatment, the first two principal components could only explain 64.4% of the data variance, and there were no strong correlations between canopy temperature and physiological parameters.

3.5. The Evaluation of CT-UAV on Physiological Traits at Different Growth Stages

Through the linear regression model, the negative relationship between canopy temperature and physiological traits in each growth stage was determined, as shown in

Table 5. The linear relationship between CT-UAV and physiological traits in different deficit irrigation treatments.

Physiological Traits	Treatment	5 May	11 May	20 May	25 May	2 June
	W0	0.41 **	0.53 ***	0.50 ***	0.51 ***	0.79 ***
SPAD	W1	0.28 *	0.46 **	0.49 ***	0.56 ***	0.72 ***
	W2	0.01	0.05	0.00	0.00	0.05
	W0	0.41 **	0.60 ***	0.64 ***	0.50 ***	0.72 ***
LAI	W1	0.22 *	0.28 *	0.65 ***	0.54 ***	0.68 ***
	W2	0.13	0.01	0.03	0.00	0.07
	W0	-	0.48 ***	0.68 ***	-	-
Pn (μmol CO2·m−2·s−1)	W1	-	0.35 **	0.58 ***	-	-
	W2	-	0.03	0.18	-	-
	W0	-	0.44 **	0.66 ***	-	-
Tr (mmol H2O·m−2·s−1)	W1	-	0.40 **	0.58 ***	-	-
	W2	-	0.03	0.25 *	-	-
	W0	-	0.41 **	0.71 ***	-	-
Cn (mmol H2O·m−2·s−1)	W1	-	0.28 *	0.60 ***	-	-
	W2	-	0.01	0.24 *	-	-
	W0	0.42 ***	0.43 ***	0.47 ***	0.60 ***	0.62 ***
Biomass (ton/ha−1)	W1	0.43 **	0.37 **	0.45 **	0.32 **	0.42 **
	W2	0.19	0.00	0.01	0.08	0.65 ***

Significance analysis is based on Pearson correlation, 0.01 < p < 0.05, is “*”, 0.001 < p < 0.01, is “**”, p < 0.001, is “***”, “-” is no value available.

. CT-UAV was closely correlated with SPAD and LAI in W0 and W1 treatments, with R² values ranging from 0.28 to 0.79 and 0.22 to 0.72, respectively. The correlations were strong in the late growth stages, with the maximum R² values appearing on 2 June. CT-UAV also correlated closely with Pn (R² = 0.35 to 0.68), Tr (R² = 0.40 to 0.66) and Cn (R² = 0.28 to 0.71). Biomass was always correlated with CT-UAV under W0 and W1 on each of the five sampling dates, with the exception that biomass in W2 was correlated with CT-UAV only on the latest sampling date (2 June). Overall, CT-UAV was correlated with biomass and physiological traits mainly under the W0 and W1 treatments, and the correlations were strong in the late growth stages.

3.6. Evaluation of Drought Resistance Performance Based on Multi-stages CT-UAV

Figure 9 shows the variety clustering groups under the three irrigation treatments based on CT-UAV. The 10 wheat varieties were clustered for the first group (good: with lower CT), the second group (moderate: with moderate CT) and the third group (poor: with higher CT). For the W0 treatment, JM418, ZM1062 and C6005 were clustered into the good group, GY5766, SM22, G35 and H4399 were classified as the moderate group, and SN086, ND3636 and GY2018 were considered as the poor group. For the W1 treatment, the wheat varieties of the good group include JM418, H4399, SN086 and SM22, the moderate group included ZM1062, C6005, G35 and GY2018 and the poor group was GY5766 and ND3636. As for the W2 treatment, JM418, H4399, SN086 and C6005 were clustered into the good group. The moderate group mainly included ZM1062, GY2018 and SM22, and GY5766, G35 and ND3636 were considered as the poor group.

The CT-UAV-based clustering results of drought resistance were similar to those based on physiological traits, with similarities of 90% (only one variety with inconsistent classification), 80% (two varieties with inconsistent classification) and 50% (half of the varieties with inconsistent classification) in the W0, W1 and W2 treatments, respectively. Moreover, there were also 90%, 70% and 50% similarities between the clustering based on CT-UAV and yield under W0, W1 and W2 treatments, respectively.

In order to determine the effectiveness of CT acquisition time, we further compared the differences between the CT-UAV clustered groups. As shown in Figure 10, for the W0 treatment, according to the canopy temperature on 5 May and 11 May, varieties in the good performance group were significantly different from the moderate performance group and poor performance group. However, the canopy temperatures of the moderate and the poor groups were not significant on these two days. Following 20 May, canopy temperature showed a significant difference between any two of the three clustering groups.

In the W1 treatment, on 5 May, CT-UAV differed significantly only between the good and the poor groups. Until 11 May, CT-UAV showed significant differences between any two of the three clustering groups. Similar to the W0 treatment, CT-UAV differed significantly between all the groups from 20 May onwards. Overall, CT-UAV differed largely in the late stages compared to the early stages.

4. Discussion

4.1. The Different Clustering of Drought Resistance Performance in Different Water Deficits

Under the three water deficit regimes, we obtained different drought resistance clustering groups of wheat varieties based on physiological performance (Figure 5). It emphasizes the complexity of drought resistance at different drought stress levels.

In mild or moderate water deficit, drought stress can make crop physiological performance deteriorate. However, for some high-yielding varieties, such a level of water stress may have a limited impact on the physiological traits, such as green leaf area, chlorophyll content and photosynthetic rate [54,55,56,57]. These varieties may still maintain physiological performance and high yield in mild or moderate water deficit conditions [58,59]. In addition, in this study, under W1 and W2 treatments, the high-yielding varieties JM418, H4399 and SN086 varieties had better drought resistance performance and grain yield. In contrast, with the further aggravation of water stress such as W0, the physiological performance of some high-yielding varieties will decline rapidly, and the yield will also shrink [60]. Some drought-resistant varieties can maintain relatedly stable productivity and physiological performance by their regulation or deep root advantage [22,61,62]. For example, the ZM1062 and C6005 varieties, while they do not always have outstanding physiological and yield performance under regular water supply or mild water stress conditions.

In addition, the wheat crop has specific water-sensitive stages, such as the jointing stage, booting stage, and flowering and grain-filling stage [63], and different wheat cultivars have different sensitivities to drought stress during these stages. Simultaneously, applying drought treatment to wheat at different growth stages has different effects on physiological traits and yield [58,64,65]. For instance, drought stress during the flowering causes a grain yield reduction of 11% to 39% and the water stress at the heading reduces grain yield by 57% [66,67]. Similarly, a yield reduction of 58% to 92% has been measured when severe drought stress was imposed at grain-filling [68]. Collectively, drought resistance performance of the same varieties may differ in different growth stages and deficit irrigation conditions. Therefore, in different deficit irrigation systems, appropriate varieties should be used according to the degree of water stress to ensure the final yield, so as to improve water use efficiency as much as possible on the premise of saving irrigation water.

4.2. Canopy Temperature Variation under Different Deficit Irrigation Treatments

Canopy temperature consistently increased with the decrease in irrigation (Figure 6), whereas a declined transpiration was observed (Table 2). Under water stress conditions, the water obtained by plants from the soil decreases, which reduces the transpiration of crop. Then, the heat exchange will decrease between the air and canopy or leaf surface and finally increases the canopy temperature [69,70]. It is worth noting that the canopy temperature difference between W0 and W2 was already large (4.1 °C) in the early stage, whereas the canopy temperature difference between W1 and W2 was more visible in the late stage (Figure 6). This also indicates that the W0 treatment was subjected to water stress and occurred earlier than the W1 treatment. Irrigation during the jointing stage (W1) can effectively reduce the impact of drought stress due to insufficient rainwater supply [71], and it has a great impact on the development of physiological traits. This irrigation enables crops to have a luxuriant plant canopy [72], which reduces bare soil and decreases soil evaporation to preserve water for transpiration and maintain lower canopy temperatures.

Our results also showed that the maximum temperature differences between varieties were 2.2 °C, 1.9 °C and 1.1 °C under W0, W1 and W2 treatments, respectively (Figure 7). The large extent of differences between varieties, especially in W0 and W1, provides the opportunity for UAV thermal infrared imagery to capture canopy temperature differences due to drought stress.

4.3. CT-UAV in Relation to Ground-based Physiological Traits in Different Deficit Irrigations

The close correlations between canopy temperature and SPAD value, LAI, Pn, Tr, Cn and biomass, particularly under W0 and W1, made it possible to evaluate the drought-resistant performance of wheat varieties using UAV thermal infrared images. The physiological performance of plants will deteriorate in drought stress conditions, while the canopy temperature is also sensitive to water deficit. Thus, there is also a certain relationship between canopy temperature and some physiological traits, e.g., chlorophyll, leaf area index, transpiration rate and stomatal conductance [39,47,73,74,75]. However, there was a poor correlation between canopy temperature and the above physiological traits in the W2 treatment. Compared with W2, W0 and W1 had more deficient irrigation, leading to more severe water stress responses and gradient changes in physiological traits [76,77,78]. Therefore, close correlations between canopy temperature and the above physiological traits were observed under the W0 and W1 treatments.

The effects of drought on crops would increase with the growth progress and organ senescence at later stages [14,79,80]. As a result, the differences in physiological traits such as SPAD, LAI, Pn, Tr and Cn among wheat cultivars were gradually significant after the anthesis stage, and the canopy temperature also increased due to the decrease in irrigation amount [81,82,83]. In contrast, the transpiration water of the crops was sufficient in W2 due to the irrigation during the flowering stage, and thus canopy temperature differences among varieties were small [75,84]. Recent studies also highlighted that no significant or weak correlation between canopy temperature and physiological traits, e.g., leaf chlorophyll content, LAI and Pn, under normal or mild water stress [85,86,87]. In this study, these results suggest that canopy temperature alone could not reflect the physiological changes of crops well under mild water stress.

Overall, the canopy temperature-based evaluation and the physiology-based evaluation (Figure 5 and Figure 9) produced consistent drought resistance clustering results for the W0 and W1 treatments, highlighting the potential of UAV thermal infrared image to identify drought resistant performance of wheat varieties in deficit irrigation conditions.

4.4. Canopy Temperature Acquisition Stage for Drought Resistance Assessment

Canopy temperature differences among the three drought-resistant performance groups became larger in the grain-filling stage (Figure 10), indicating that the canopy temperature obtained after this stage was more representative and can be more accurate to distinguish drought resistant performance among different wheat varieties. In deficit irrigation conditions, a large amount of soil moisture was lost due to evaporation in the early stage along with plant transpiration, resulting in running out of soil moisture during the grain-filling stage and further water stress [88]. This, in turn, leads to greater covariations in crop physiological parameters and canopy temperature among different wheat cultivars [14]. In this study, the correlations between canopy temperature and physiological traits (SPAD, LAI, Pn, Tr and Cn) were stronger under the W0 and W1 treatments and after 20 May. Therefore, the drought resistance can be well reflected by the canopy temperature, especially in the late stage of crop growth (grain-filling stage and later). Similarly, canopy or leaf temperatures at grain-filling stages and later were found to be better related to the physiological traits of the varieties, such as leaf area index, photosynthetic parameters, chlorophyll content and yield in previous studies [73,89,90,91]. Therefore, this study allows the recommendation for obtaining UAV thermal infrared images to evaluate drought resistant performance of wheat varieties during the grain-filling stage.

5. Conclusions

We obtained canopy temperature from UAV thermal infrared images and measured physiological traits (LAI, SPAD, Pn, Tr, Cn and biomass) of 10 wheat varieties under three deficit irrigation regimes, across multiple growth stages to evaluate their drought resistance. This study showed that irrigation deficit induced declines in the measured physiological traits but an increase in the canopy temperature. Moreover, the results showed closely negative correlations between canopy temperature and physiological traits under the W0 and W1 treatments, whereas almost no correlation was observed under the W2 treatment. The observed correlations between canopy temperature and physiological traits were relatively high following the grain-filling stage.

The UAV thermal imagery combined with the hierarchical clustering method successfully identified the drought-resistant performance (good, moderate, and poor) of 10 wheat varieties under W0 and W1 treatments, which was in line with the drought resistance evaluation results based on physiological traits. In contrast, the results were poor under the W2 treatment. The grain-filling stage is found to be more suitable than the early stages for drought resistance evaluation based on UAV thermal imaging, which should be further validated.