Evaluation of Soybean Drought Tolerance Using Multimodal Data from an Unmanned Aerial Vehicle and Machine Learning

Abstract

Drought stress is a significant factor affecting soybean growth and yield. A lack of suitable high-throughput phenotyping techniques hinders the drought tolerance evaluation of multi-genotype samples. A method for evaluating drought tolerance in soybeans is proposed based on multimodal remote sensing data from an unmanned aerial vehicle (UAV) and machine learning. Hundreds of soybean genotypes were repeatedly planted under well water (WW) and drought stress (DS) in different years and locations (Jiyang and Yazhou, Sanya, China), and UAV multimodal data were obtained in multiple fertility stages. Notably, data from Yazhou were repeatedly obtained during five significant fertility stages, which were selected based on days after sowing. The geometric mean productivity (GMP) index was selected to evaluate the drought tolerance of soybeans. Compared with the results of manual measurement after harvesting, support vector regression (SVR) provided better results (N = 356, R² = 0.75, RMSE = 29.84 g/m²). The model was also migrated to the Jiyang dataset (N = 427, R² = 0.68, RMSE = 15.36 g/m²). Soybean varieties were categorized into five Drought Injury Scores (DISs) based on the manually measured GMP. Compared with the results of the manual DIS, the accuracy of the predicted DIS gradually increased with the soybean growth period, reaching a maximum of 77.12% at maturity. This study proposes a UAV-based method for the rapid high-throughput evaluation of drought tolerance in multi-genotype soybean at multiple fertility stages, which provides a new method for the early judgment of drought tolerance in individual varieties, improving the efficiency of soybean breeding, and has the potential to be extended to other crops.

1. Introduction

Drought is the most complex global abiotic stress [1]. Reduced plant water content following drought stress severely affects plant growth and development, leading to a significant reduction in yield [2,3]. The average yield loss of cereals due to drought is up to 17%, and according to statistics, China loses 7–8 billion kilograms of grain production annually due to drought [4]. Plant root systems have difficulties absorbing enough water from the soil under drought conditions, causing various changes in the metabolism of the tissues and cells and leading to differences in the phenotypic traits of plants to a certain extent [5]. Water-sensitive soybean is more susceptible to drought stress.

Soybeans are an important legume grown worldwide [6]. Soybeans are rich in nutritional value, and soy protein is considered an essential source of protein for the human body [7]. Soybeans are also an important cash crop with significant pharmaceutical, food, and feed-processing value [8]. China is one of the world’s leading producers of soybeans. However, its domestic production is unable to meet the growing demand, and China remains the world’s largest soybeans consumer [9,10]. Soybeans have an underdeveloped root system and require adequate water throughout their growth, making them highly susceptible to increasing natural droughts [8,11]. In traditional soybean breeding, many soybean genotypes are evaluated under drought-stress environments, and experienced breeders often combine their experience and post-harvest yields at maturity to acquire drought-tolerant and high-yielding soybean varieties, which is subjective, time-consuming, and financially expensive [12]. There is an urgent need for a method to evaluate drought tolerance in multi-genotyped soybeans early, non-destructively, accurately, economically, and in a high-throughput manner.

Recently, plant phenomics technologies based on various sensors have been developed [13,14,15]. Researchers can use different sensors to acquire various data form plants and then combine them with machine learning to reflect multiple physiological parameters and phenotypic traits of plants [16,17]. Among these, unmanned aerial vehicle (UAV)-based aerial phenotyping platforms have been widely used for the precise phenotyping of field crops owing to their nondestructive detection, low cost, and high throughput [18,19], including the accurate and rapid estimation of plant height [20], leaf area index [21], biomass [22], plant density [23], relative maturity [24], flowering time [25], and other traits. Crops under different stress environments produce several morphological responses (e.g., reduced plant height, leaf area, and biomass) and physiological responses (e.g., reduced chlorophyll and relative water content) at various fertility stages [26]. These responses can potentially be detected and quantified by UAV remote sensing methods. Li et al. [27] identified previously unexpressed loci through UAV time-series data for a soybean canopy. Qiao et al. [28] improved the accuracy of leaf area index dynamic estimation by using morphological and spectral information over multiple fertility stages in maize. Zhou et al. [29] estimated flood-induced damage by monitoring soybeans after flooding stress with a UAV system. Jiang et al. [30] quantified the dynamic drought response of rice populations in a rice field environment by analyzing manual leaf roll scoring and plant content via a consumer-grade UAV platform. Liu et al. [31] differentiated the drought severity of maize by classifying predictions of maize drought stress through the visual assessment of drought stress from UAV system visible-light imagery. Zhou, Zhou, Ye, Ali, Nguyen, and Chen [29] demonstrated that UAV imagery can evaluate fast- and slow-wilt traits in soybean genotypes for drought tolerance. Soybeans under drought stress also produce several stress responses at various fertility stages, which makes it possible to accurately extract changes in multiple genotypes of soybeans under drought stress and analyze the differences between genotypes using UAV remote sensing methods.

Until now, most of the reference standards in drought research based on UAV remote sensing are the results of visual evaluations at maturity [29,30,31]. UAV data are analyzed and visual evaluations are conducted using characteristics based on canopy growth, which can lead to idealized results. However, yield is considered the primary criterion for selecting high-yielding genotypes under different environmental conditions in practical research [32]. Yield performance of all the genotypes under different environments was inconsistent, and the performance of the same genotype under drought stress and full irrigation was uncertain [33]. Therefore, breeding experts have proposed the drought tolerance index (DTI), which is based on the mathematical relationship between yield under normal and stress conditions, for a more accurate and intuitive evaluation of drought tolerance in crops [33]. Common DTIs include mean productivity (MPI) [34], geometric mean productivity (GMP) [35], stress tolerance (STI) [35], stress sensitivity (SSI) [36], yield (YI) [37], tolerance (TOL) [34], and yield stability indices (YSIs) [38]. Some researchers have used UAV remote sensing methods to predict yield [39,40,41]. Still, the results showed that they are highly influenced by genotype, environment, geographic location, irrigation conditions, and mostly at maturity. This approach limits the timeliness and robustness of drought tolerance evaluation. However, the relationship between UAV traits and DTIs has not been analyzed in soybean drought tolerance studies, and there have been no attempts to evaluate drought tolerance at multiple fertility stages. This study aimed to explore the potential of UAV phenotypes for drought tolerance evaluation in multi-genotyped soybeans at multiple fertility stages. The specific objectives of this study were as follows:

(1) Evaluating the differences between the UAV-based canopy trait data of different soybean genotypes under well water (WW) and drought stress (DS).

(2) Evaluating the correlation between DTIs and UAV-based canopy trait data during maturity and yield.

(3) Evaluating the ability of UAV-based canopy trait data and machine learning to regress the GMP and categorize the DISs.

2. Materials and Methods

2.1. Field Experiments and Experimental Equipment

As shown in Figure 1 the experiment was conducted at the Yazhou and Jiyang trial sites in Sanya City, Hainan Province, China, using uniform soil throughout the field. Hainan, which receives low rainfall in winter, sowed 215 and 510 soybean varieties in Jiyang and Yazhou, respectively, in the winters of 2021 and 2022. Owing to a force majeure, 427 and 509 soybean data points were eventually retained, respectively. Four blocks were set up in each field, including two WW and two DS replicates. In 2021, planting was performed in double rows with plant spacing of 10 cm and ridge spacing of 60 cm, and in 2022, double rows with plant spacing of 15 cm and ridge spacing of 80 cm. Each cultivar was sown using 16 seedlings. All blocks were irrigated before seedling emergence. WW continued to maintain irrigation after seedling emergence for moisture in the soil layer from 0 to 50 cm to reach 85 ± 5% of the field water holding capacity to ensure that the soybean was in a moisture-suitable condition during its whole life cycle. The DS group was not irrigated after seedling emergence, and further management was consistent with that of WW. Harvesting began 75 days after sowing (DAS) and was completed at DAS105. Breeders selected three consecutive plants of each soybean variety to determine the yields through a digital seed examiner.

2.2. UAV Image Data Collection

A UAV (Matrice 300 RTK, DJI, Shenzhen, Guangdong, China) equipped with visible and multispectral sensors was used to simultaneously collect both types of remote sensing images. The visual sensor (P1, DJI, Shenzhen, Guangdong, China) had a resolution of 8192 × 5460 pixels. The multispectral sensor (Rededge-MX; MicaSense, Seattle, WA, USA) consisted of five bands with a wavelength range of 400–900 nm and resolution of 1280 × 960. The flight path was planned using a DJI plot; the altitude of each flight was 30 m, the speed was 1 m/s, and the imaging interval was one frame/s. The heading and side overlap rates were 90% and 80%, respectively, for the visible images and 80% and 70% for the multispectral images. Two reflector plates with reflectance values of 5% and 15% were used to calibrate the multispectral camera, and 12 GCPs evenly distributed around and in the middle of the field were used for calibration during the UAV image reconstruction. The GCP geographic coordinates were acquired using hemispherical real-time kinematic (RTK) differential global navigation satellite system (GNSS) measurements. Jiyang UAV data were acquired at 71 days after sowing (DAS) when more than 90% of the varieties were at maturity. Yazhou UAV data were first acquired at 21DAS during the nutrient growth stage, and then at 10-day intervals after entering fertility growth. To ensure that the weather conditions at the time of acquisition fulfilled the requirements, image data were acquired on five separate dates: 35DAS, 46DAS, 56DAS, and 64DAS. The final data collection was conducted when more than 90% of the species were in the maturity stage. The data were collected between 11:00 and 14:00, and the weather was sunny, breezy, and cloudless.

2.3. Image Processing and Canopy Trait Extraction

The RGB and multispectral (MS) images captured by the UAV were reconstructed by Structure in Motion Recovery (SFM) using Agisoft photoscan 1.4.3 (Agisoft LLC, St. Petersburg, Russia) software to obtain dense point clouds and visible orthophotographs from the RGB images and multispectral orthophotographs from the MS images [42]. The sensor size, focal length, and pixel size were used to improve the accuracy of stitching. After importing the UAV images, the photographs were first aligned based on the feature points between the images, and the positions of the 12 GCPs were manually adjusted to improve the horizontal and vertical accuracies of the output orthophotos and dense point clouds. The image pixel values of the MS orthophotos were converted into reflectance values using a calibration board. In Figure 2a,b shows the process of extracting canopy length, width, and cover, Figure 2c shows the process of extracting plant height and canopy volume, and Figure 2d shows the process of extracting vegetation index.

2.3.1. Canopy Coverage, Length, and Width

A fixed-size rectangular box was used to intercept the RGB orthophotos based on the species to acquire a soybean plot, and each plot was divided into three bands: blue, green, and red. Exg is often used to segment green plants from shadows and soil [43]. The acquired Exg image was first subjected to Otsu segmentation and then binarized to acquire a mask, from which the soybean region of interest (ROI) was extracted [44]. Three canopy traits were extracted from the processed RGB images: canopy coverage, length, and width. Canopy coverage was defined as the ratio of soybean pixels to total pixels in the ROI, and the numbers of rows and columns of non-0-pixel rows and non-0-pixel columns of the ROI after removing the background were used as the canopy length and canopy width [45]. To better segment soybeans from backgrounds, such as shadows and soil, Exgh was defined based on Exg. In testing, it was found that the segmentation of soybeans was best when λ = 3, μ = 2.1, and ω = 1.

Exg = 2 × G − R − B

(1)

$$\mathrm{Exgh}=\mathsf{\lambda }\times \mathrm{G}-\mathsf{\mu }\times \mathrm{R}-\mathsf{\omega }\times \mathrm{B}$$

(2)

In the above equations, G, R, and B are the reflectance values of the green, red, and blue spectral bands, respectively.

2.3.2. Plant Height and Canopy Volume

A 3D point-cloud model of plants can directly and effectively reflect their morphological characteristics. The random sample consensus (RANSAC) algorithm has been used to segment plant point clouds and has shown better results [46]. The RANSAC algorithm has a wide range of applications in plant point-cloud segmentation. In this experiment, the RANSAC algorithm was used to separate the soybean point cloud from the soil, and all operations were performed using the open-source 3D data-processing library Open3D. First, the planting area of the soybean was cut out from the dense point cloud acquired by reconstruction, and then, the point cloud was filtered to remove noise and reduce the point-cloud density. The filtering used uniform downsampling filtering, uniform downsampling reduced the number of points under the premise of retaining the point-cloud traits, and Open3d provided the function of uniform downsampling to the point cloud. The RANSAC algorithm was applied to the filtered point-cloud data to segment the soybean into two parts, the above-plane and the below-plane parts, through a plane. The above-plane part was the soybean, and finally, it was segmented into a single ROI region. Two canopy traits were acquired from the processed soybean point-cloud data: plant height and canopy volume. Plant height was defined as the average height between the maximum and 70th percentiles of the point-cloud heights in the ROI plots, and the canopy volume was generated by calculating the surface difference score of the volume between the highest and lowest points of the canopy [47].

2.3.3. Vegetation Index

A grayscale map of Exgh was acquired from the red, green, and blue bands of the MS orthophoto. A mask of the entire soybean region was acquired after the Otsu threshold segmentation and binarization. Soybean images in the red, green, blue, near-infrared, and red-edge bands were extracted separately using a mask, and each soybean image was segmented into several ROIs based on the species. The grayscale mean values were calculated for each ROI and converted into reflectance values using a reflectance plate. Three vegetation indices, the normalized difference vegetation (NDVI), atmospherically resistant vegetation (ARVI), and red-edge chlorophyll (RECI) indices, were calculated from five multispectral bands [48,49,50]. These indices have been widely used and reported to be effective indicators for assessing crop drought stress [51,52].

NDVI = (NIR − RED)/(NIR + RED)

(3)

ARVI = (NIR − (2 × R − B))/(NIR + (2 × R − B))

(4)

RECI = (NIR/R) − 1

(5)

In the above equations, R, B, NIR, and RED are the reflectance values of the red, blue, near-infrared, and red-edge spectral bands, respectively.

2.4. Data Analysis

All statistical analyses were performed using Python 3.8. The differences between canopy traits under WW and DS and their correlations with yield were analyzed, and the dynamics between canopy traits at different fertility stages of soybeans under WW and DS and their relationship with yield were evaluated. Using the extracted canopy traits, two algorithms of ensemble learning, extreme gradient boosting (XGBoost) and random forest (RF), as well as the classical linear-kernel support vector machine (linear-SVM) algorithm, were used to estimate the DTI and DIS [53,54,55]. Ensemble learning combines multiple weak learners to acquire better-performing learners. XGBoost is an iterative algorithm based on boosting that implements a gradient-boosting decision-tree architecture, which leads to higher accuracy in regression and classification. The core innovations of XGBoost include optimized gradient-boosting algorithm, support for multiple loss functions, feature importance assessment, support for parallelization, and the prevention of overfitting. The optimized gradient-boosting algorithm improves the efficiency and accuracy of the gradient-boosting algorithm through a series of optimisation strategies, including the second-order Taylor expansion of the objective function, the introduction of regularization terms, and the use of approximate greedy algorithms. RF is an ensemble learning algorithm based on bagging, which improves overfitting tolerance after averaging the decision trees. Different training sets are constructed by randomly putting back multiple samples from the original dataset, and then, multiple decision trees are trained on each of these training sets. The random selection of features is also introduced when training each decision tree, and this randomness increases the variability between each tree and reduces the possibility of overfitting. SVM is a classical supervised machine learning algorithm. SVM has better results in both linear and nonlinear problems, in dealing with nonlinear problems by choosing a suitable kernel function, and in mapping the data into a higher dimensional space, thus achieving linear differentiability, SVM has strong robustness and generalization ability in high-dimensional space, and it performs well for small sample datasets. In this study, a linear kernel was chosen in SVR because it was found to be the most efficient among many kernel functions, while the Gaussian kernel was the most efficient in SVC [56].

When modeling through machine learning, real data are still required for training. To minimize the time spent manually checking the seeds, a model with excellent generalization and accuracy should ideally be able to excel when trained with a small number of datasets and generalize to a large number of datasets. To verify the accuracy of the model on different data, we kept the number of data from Yazhou and Jiyang as close as possible, which gives a stronger indication that our model has a better ability to migrate to data from different times and locations. In the regression training process of DTIs, the data from Yazhou 64DAS were divided into training and test sets at a ratio of 3:7, and the model was then migrated to the data from Jiyang. The model hyperparameters were optimized using a nested 5-fold cross-validated grid search. In the classification training process of DISs, the data of Yazhou 64DAS were divided into training and test sets in a ratio of 7:3, and then, all the data of 21DAS, 35DAS, 46DAS, 56DAS, and 64DAS were used to test the performance of the model. The classification accuracy of the test set was calculated to evaluate the performance of the proposed model. All modeling processes in this experiment were performed on a PC equipped with an NVIDIA RTX 3060 GPU and a 3.60 GHz Intel Core i7-11700k CPU with 32 GB of RAM. The R², RMSE, and relative RMSE (RRMSE) between the measured and predicted GMP were calculated, and these three metrics were used to evaluate the model’s accuracy and robustness.

$${{\displaystyle \mathrm{R}}}^2=1-{\displaystyle {\sum }_{\mathrm{i}}^{\mathrm{n}}{{\displaystyle ({{\displaystyle \mathrm{y}}}_{\mathrm{i}}-{{\displaystyle \widehat{\mathrm{y}}}}_{\mathrm{i}})}}^2/{\displaystyle {\sum }_{\mathrm{i}}^{\mathrm{n}}{{\displaystyle ({{\displaystyle \mathrm{y}}}_i-\overline{\mathrm{y}})}}^2}}$$

(6)

$$\mathrm{RMSE}=\sqrt{{\displaystyle {\sum }_{\mathrm{i}}^{\mathrm{n}}{{\displaystyle ({{\displaystyle \mathrm{y}}}_i-{{\displaystyle \widehat{\mathrm{y}}}}_{\mathrm{i}})}}^2/\mathrm{n}}}$$

(7)

$$\mathrm{RRMSE}=\mathrm{RMSE}/\overline{\mathrm{y}}\times 100$$

(8)

$$\mathrm{Accuracy}=\frac{\mathrm{Number} \mathrm{of} \mathrm{samples} \mathrm{classified} \mathrm{correctly} \mathrm{in} \mathrm{a} \mathrm{test} \mathrm{set}}{\mathrm{Total} \mathrm{Number} \mathrm{of} \mathrm{samples} \mathrm{in} \mathrm{a} \mathrm{test} \mathrm{set}}\times 100\%$$

(9)

In the above equations, y_i and ${\widehat{\mathrm{y}}}_i$ are the actual and predicted values for the ith sample, respectively, and $\overline{\mathrm{y}}$ is the mean value of the total number of samples.

In

Table 1. Descriptions of DTIs based on grain yield used in this study.

DTI	Formula	Reference
YI	${\mathrm{Y}}_{\mathrm{DS}}{/\overline{\mathrm{Y}}}_{\mathrm{DS}}$	[38]
SSI	$(1-{\mathrm{Y}}_{\mathrm{WW}}{/\mathrm{Y}}_{\mathrm{DS}})/(1-{\overline{\mathrm{Y}}}_{\mathrm{WW}}{/\overline{\mathrm{Y}}}_{\mathrm{DS}})$	[37]
STI	${(\mathrm{Y}}_{\mathrm{WW}}\times {\mathrm{Y}}_{\mathrm{DS}})/{({\overline{\mathrm{Y}}}_{WW})}^2$	[36]
TOL	YWW − YDS	[35]
MPI	(YWW + YDS)/2	[35]
GMP	(YWW × YDS)0.5	[36]
YSI	YDS/YWW	[39]

, Y_WW and Y_DS are the grain yields of soybean germplasm grown under WW and DS, respectively. ${\overline{\mathrm{Y}}}_{\mathrm{WW}}$ and ${\overline{\mathrm{Y}}}_{\mathrm{DS}}$ are the mean grain yields of all soybean germplasms under WW and DS, respectively.

Table 2. Classification of DISs based on GMP.

GMP	DIS
[13, 80]	1
(80, 147]	2
(147, 214]	3
(214, 256]	4
(256, ∞)	5

shows the process of calculating DIS based on GMP.

3. Results and Discussion

3.1. Differences between Soybean Canopy Traits under WW and DS

3.1.1. Differences between Canopy Traits at Maturity in Soybeans under WW and DS

As shown in Figure 3, significant differences in canopy coverage, plant height, length, width, canopy volume, NDVI, ARVI, and RECI at maturity between soybeans grown under WW and DS were observed. For most varieties, the absolute values of each soybean trait were significantly higher in the WW than in the DS group [57]. To reflect the differences in the same canopy trait of soybeans under the two moisture conditions, the mean absolute percentage error (MAPE) was calculated for the same canopy trait under DS relative to WW. The MAPE values for canopy coverage and plant height were 39.49% and 48.73%, respectively, and the trends were similar. Generally, a plant grows, expands upward, and is synchronized in all directions. Length showed no significant difference between WW and DS, with a MAPE of only 12.32%. This is related to the planting method and its calculations. The MAPE value of 26.86% for width showed a more significant difference than that for length. As it was noted that canopy volume characterizes plant morphology in the horizontal and vertical directions, it makes the differences in traits under WW and DS even greater [47]. Thus, the MAPE value for volume was the largest of all the traits related to biomass (78.71%). The RECI showed more significant differences among the three vegetation indices, but the differences in NDVI were minor under WW and DS. The trends of the three vegetation indices tended to be similar because they were all derived from multispectral bands.

The above data analysis showed that DS caused a smaller canopy, a shorter height, and a smaller vegetation index in soybean plants, reflecting the fact that drought stress leads to a reduction in the physiological activity of soybean plants, which, in turn, affects the morphological characteristics and physiological parameters of soybeans and ultimately leads to a decrease in yield. This has been confirmed in previous studies [26,58]. However, some varieties that showed opposite trends under WW and DS were observed; that is, DS was superior to WW in some traits, which may be related to the adaptive ability of varieties from different latitudes in the tropics or possible differences in material composition among separate compartments in the experimental area [59]. In addition, the level of each trait under WW conditions was not the same among varieties; therefore, acquiring accurate traits under WW and DS conditions and conducting correlation analyses to accurately characterize each variety’s drought tolerance are necessary. The analysis of soybean canopy traits at maturity under WW and DS conditions revealed the final performance of each trait. This provided the most direct basis for evaluating drought tolerance. The drought tolerance of soybean plants is also reflected in the dynamics of each stage, and assessing drought tolerance through early plant performance during the breeding process is essential.

3.1.2. Dynamic Changes in Canopy Traits at Different Fertility Stages of Soybean under WW and DS

Drought stress continuously affects soybeans at all fertility stages, and the drought tolerance of different soybean genotypes varies. This resulted in an additional growth state for each genotype at each fertility stage, which manifested in various distribution ranges and specific values of canopy traits. Breeders recorded the fertility stages of all experimental samples, and UAV data were acquired during several major stages of soybean fertility. Soybeans gradually increased in canopy size after seed germination and were clearly distinguishable from weeds (21DAS). Soybeans entered the flowering stage when flowers opened at any node of the main stem (35DAS), and then entered the podding stage when young pods gradually formed (46DAS). The flowering and podding stages represented the most vigorous soybean growth and development stages. Subsequently, the soybean entered the bulging stage, the seed began to expand gradually, and the bulging stage was the stage in which the dry matter of the soybean accumulated the most (56DAS). After entering the maturity, the soybean leaves began to turn yellow and fell off, and the seeds were dehydrated (64DAS) [60].

Figure 4 shows boxplots of the Yazhou canopy traits under WW and DS. As the soybeans grew and developed, biomass-related traits (canopy coverage, plant height, length, width, and volume) increased under WW and DS conditions, with little change at 56 and 64DAS, which is consistent with the result of Roth et al. [61]. The length and width exhibited a wider distribution range during the early stages. The distribution range of length slowly narrowed as the soybean grew. The width distribution range initially narrowed and then increased. The distribution ranges of canopy coverage, plant height, and volume were narrow in the early stages, and continued to widen as the soybeans grew. Similar results can be found in the study of Zhou, Zhou, Ye, Ali, Chen, and Nguyen [41], whereby the differences between the sample populations in the field consistently became more significant with increasing drought stress. Biomass-related traits under WW increased rapidly before 46DAS but changed insignificantly afterwards. The vegetation indices (NDVI, ARVI, and RECI) reached their highest values at 46DAS. The vegetation indices were slightly higher under DS than under WW in the early stage and were significantly lower after 46DAS. The general trends in soybean canopy traits were similar across the fertility stages. However, the distribution ranges of the population samples changed over time during different fertility stages and in different environments. This provides data support for subsequent DTI and DIS studies using canopy traits.

Changes in Pearson’s correlation between Yazhou soybean yield and canopy coverage, plant height, length, width, volume, NDVI, AVRI, and RECI during different fertility stages are shown in Figure 5. Under WW, the Pearson correlation of the biomass-related traits gradually increased, reaching a maximum at 56DAS and fluctuating with a minimal magnitude at 64DAS. The three vegetation indices showed some correlation in the early stages, and all three correlations with yield decreased at 35DAS. As growth and development proceeded, the correlation increased to a certain level beyond 0.6 and reached its highest value close to 0.8 at 64DAS. Under DS, both biomass-related traits and vegetation indices showed some correlation and gradually increased; however, biomass-related traits reached their highest values at 46DAS. The correlations between canopy coverage and plant height, width, and volume showed a pronounced decline after 46DAS. Consistent with the research of Zhou, Zhou, Ye, Ali, Chen, and Nguyen [41], the correlations for the vegetation index consistently increased and reached a maximum at the final time point.

When WW and DS were combined, the strongest correlations for biomass-related traits tended to occur at 64DAS under WW and at 46DAS under DS. The vegetation indices showed similar characteristics under WW and DS, both of which gradually increased, which is consistent with the observation of da Silva et al. [62]. Notably, at 21DAS and 35ADS, all trait and yield correlations were stronger under DS than under WW conditions. However, all canopy traits and yield correlations were significantly better under WW than under DS during the later stages. The changing relationships among traits during different stages provide theoretical evidence for the inversion of DTIs using UAV data and the evaluation of drought tolerance levels during different stages.

3.2. Correlation between Canopy Traits at Maturity and DTIs of Soybean Yield

Various DTIs have been proposed for agriculture based on the application requirements in different scenarios. Pearson’s correlations between canopy traits, yield, and DTIs were calculated under WW and DS conditions based on data collected within two years for each trait. As illustrated in Figure 6, the canopy trait Pearson’s correlations under YI and DS ranged from 0.35 to 0.56 with high correlations, and the correlations under WW were significantly lower at only 0.13 to 0.23. The opposite was true for SSI, YSI, and TOL, with all three DTIs showing high correlations under WW and low correlations under DS. STI, MPI, and GMP showed high correlations under both conditions. The strong correlation between YI and yield under DS implies that YI is more suitable for screening high-yielding varieties under adverse environmental conditions. Similarly, the higher correlation between SSI, TOL, YSI, and yield under WW indicates that they are more suitable for screening high-yielding varieties in a typical environment. STI, MPI, and GMP, with higher correlations under both DS and WW conditions, can be used for screening high-yielding tolerant varieties in both adverse and typical environments, with GMP showing a better balance. This is consistent with the results of drought tolerance research on chickpeas and wheat [32,63]. The drought tolerance judgment of multi-genotype soybeans over multiple fertility stages requires the comparison of data under WW and DS; therefore, a more balanced metric, GMP, would be the optimal choice.

3.3. GMP Regression of Soybean at Maturity Based on Machine Learning

UAVs can quickly acquire crop traits over a large area and differentially reflect growth changes in each area. However, when modeling using machine learning, actual measurement data are required for training. To minimize the time spent on manual measurements, a model with excellent generalization ability and accuracy should ideally acquire extraordinary results in small-dataset training and be generalized to many datasets. From the previous analysis, it was found that GMP performs well in screening drought-tolerant soybean varieties under WW and DS, and the accurate prediction of GMP is vital for selecting drought-tolerant and high-yielding soybean varieties. Figure 7 shows the GMP regression prediction results using the canopy traits of Yazhou 64DAS achieved by the ensemble learning models, XGBoost, RF, and support vector regression (SVR). Thirty percent of the Yazhou data was used to train the model, and the remaining seventy percent was used to test the model’s generalization over large areas. All data from Jiyang 71DAS were used to test the model’s generalizability over different years and geographic regions. SVR exhibited the best performance on both the Yazhou and Jiyang datasets. For R², the Yazhou dataset showed a higher R² than the Jiyang dataset, with values of 0.74–0.75 and 0.63–0.68, respectively. For the RMSE and RRMSE, the Jiyang dataset was better than the Yazhou dataset, and the average values of the RMSE and RRMSE for the three models decreased by 14.01 g/m² and 10.30%, respectively. The main reason for this result is that the total amount of data in Jiyang was small, and the distribution was narrower than that of Yazhou. Another reason was that the intercropping pattern in Jiyang was slightly more intensive, which led to the accuracy of the extracted data being slightly lower than that of the Yazhou data. The prediction results at the ripening stage showed that machine learning and UAV methods have high accuracy for evaluating the drought tolerance of soybeans. However, maturity is reached close to harvest time, and breeders like to roughly evaluate drought tolerance in soybeans at an early stage.

3.4. Drought Injury Score (DIS) Classification of Soybean at Multiple Fertility Stages Based on Machine Learning

In practice, drought tolerance grades are more intuitive and widely used than DTIs. Previous studies have been conducted on UAV-based waterlogging grades for soybeans and drought tolerance grades for rice, while there are still fewer studies on drought tolerance grades for soybeans [30,45]. A DIS with five grades was assigned to soybean varieties according to the GMP. Most of the samples in the five grades were distributed in grades “2” and “3”, accounting for 71.9% of the total number of samples. The canopy traits of Yazhou 64DAS were used to classify the DIS using XGBoost, RF, and support vector classification (SVC), and the results are shown in Figure 8a–c. The model results on 153 test sets showed that XGBoost had the highest accuracy of 77.12%, which was better than those of RF (75.16%) and SVC (72.55%). The DIS was predicted using the trained XGBoost model for Yazhou at 64DAS, and the results are shown in Figure 8d. The accuracy of the classification was 80.55%; only eight samples were misclassified as non-adjacent grades, and more than 98% of the samples were correctly classified or classified as adjacent grades.

The DIS was predicted for 21DAS, 35DAS, 46DAS, and 56DAS data by the established XGBoost model and the results are in

Table 3. DIS classification accuracy from 21DAS, 35DAS, 46DAS, and 56DAS.

Comparison	Accuracy	21DAS	35DAS	46DAS	56DAS
Manual result	Correct DIS	27.90%	42.44%	63.26%	62.48%
	Correct and adjacent DIS	73.48%	89.00%	94.89%	95.48%
Predicted result at 64DAS	Correct DIS	26.72%	47.15%	68.57%	69.94%
	Correct and adjacent DIS	76.62%	90.37%	96.27%	97.05%

. The prediction results are shown in Figure 8e–h. The prediction accuracy significantly increased with a change in the fertility stage. At 21DAS, the accuracy was only 27.90%, and most plants were assigned incorrect grades, with a large number of grade “2” plants even assigned to grade “1”. Soybeans were in the seedling stage at DAS21; the impact of drought stress was limited, and a slight difference was observed between soybean varieties under WW and DS, which could only classify categories with minor differences. After the seedling stage, 35DAS and 46DAS were at the flowering and podding stages, when plant growth and development were the most rapid, and classification accuracy appeared to have significantly improved (42.44% and 63.26%, respectively). Especially at 46DAS, more than 70% of grade “2” plants were correctly classified, while 94.89% of the samples were correctly classified or classified as an adjacent grade. At 56DAS, most grade “2” plants were correctly classified, with a further reduction in the number of grades classified as non-adjacent compared to 46DAS, where 95.48% of the samples were correctly classified or classified as an adjacent grade.

The predicted DIS was compared to the 64DAS predicted DIS for the 21DAS, 35DAS, 46DAS, and 56DAS data by the established XGBoost model. The comparison results are shown in Figure 8i–l. The overall trend changes are similar to those in Figure 8e–h; however, a more noticeable improvement in the prediction accuracy was observed. At 21DAS, the accuracy was only 26.72%, and more than 76% of the grade “2” plants were classified into grade “1”, and none of the grade “4” and grade “5” plants were correctly classified. Faster growth and development of 35DAS and 46DAS soybeans resulted in a rapid increase in prediction accuracy to 47.15% and 68.57%, respectively. It is noteworthy that at 35DAS, some grade “5” plants could be correctly classified. At 46 DAS, some grade “4” plants could be correctly classified, and the accuracy of being classified correctly or classified as an adjacent grade reached 96.27%. The accuracy further improved at 56 DAS, reaching 69.94%.

4. Conclusions

This study investigated the potential of using multimodal UAV data and machine-learning models to assess the drought tolerance of multi-genotype soybeans over multiple fertility stages. Hundreds of soybean canopy trait data were acquired in Jiyang (N = 427) and Yazhou (N = 509) under WW and DS conditions, respectively. Differences in the canopy traits of soybeans at maturity were analyzed under WW and DS. All canopy traits showed some correlation with yield. For biomass-related traits, the highest correlation under WW tended to occur at the maturity stage (64DAS), and under DS, at the podding stage (46 DAS). The vegetation index showed similar trends under WW and DS, both of which gradually increased.

The GMP was calculated to evaluate drought tolerance in soybeans and select drought-tolerant varieties. Thirty percent of the Yazhou data were selected for modeling to predict the GMP, and SVR showed the best performance on both the Jiyang and Yazhou test sets. The soybean DISs were classified based on GMP. Seventy percent of the Yazhou data were used for modeling to predict the DIS, and XGBoost had the highest accuracy. The predictions of the DIS and 64DAS predicted DIS were compared for the 21DAS, 35DAS, 46DAS, and 56DAS data, and the two datasets showed similar variations. It is worth mentioning the validation of the model’s advantages for extending the model to large areas and sites in different years, with a more visual representation of the drought tolerance of each soybean variety based on the grade classification. The studied UAV system will be an essential tool for screening drought-tolerant plants; it eliminated the need for the destructive sampling of the biomass or yield of all plants in the experimental plots and allows for the determination of whether plants are drought-tolerant and high-yielding early in life, which significantly reduces costs and saves time. In the future, we will conduct genetic analyses such as genome-wide association studies (GWAS) to identify functional gene loci for drought tolerance in soybeans and better explain the physiological mechanisms of drought in soybeans. In addition, we will continue to conduct UAV testing on various crops under various biotic and abiotic stresses to extend this research to more application scenarios.