UAV- and Machine Learning-Based Retrieval of Wheat SPAD Values at the Overwintering Stage for Variety Screening

: In recent years, the delay in sowing has become a major obstacle to high wheat yield in Jiangsu Province, one of the major wheat producing areas in China; hence, it is necessary to screen wheat varieties are resilient for late sowing. This study aimed to provide an effective, fast, and non-destructive monitoring method of soil plant analysis development (SPAD) values, which can represent leaf chlorophyll contents, for late-sown winter wheat variety screening. This study acquired multispectral images using an unmanned aerial vehicle (UAV) at the overwintering stage of winter wheat growth, and further processed these images to extract reﬂectance of ﬁve single spectral bands and calculated 26 spectral vegetation indices. Based on these 31 variables, this study combined three variable selection methods (i.e., recursive feature elimination (RFE), random forest (RF), and Pearson correlation coefﬁcient (r)) with four machine learning algorithms (i.e., random forest regression (RFR), linear kernel-based support vector regression (SVR), radial basis function (RBF) kernel-based SVR, and sigmoid kernel-based SVR), resulted in seven SVR models (i.e., RFE-SVR_linear, RF-SVR_linear, RF-SVR_RBF, RF-SVR_sigmoid, r-SVR_linear, r-SVR_RBF, and r-SVR_sigmoid) and three RFR models (i.e., RFE-RFR, RF-RFR, and r-RFR). The performances of the 10 machine learning models were evaluated and compared with each other according to the achieved coefﬁcient of determination (R 2 ), residual prediction deviation (RPD), root mean square error (RMSE), and relative RMSE (RRMSE) in SPAD estimation. Of the 10 models, the best one was the RF-SVR_sigmoid model, which was the combination of the RF variable selection method and the sigmoid kernel-based SVR algorithm. It achieved high accuracy in estimating SPAD values of the wheat canopy (R 2 = 0.754, RPD = 2.017, RMSE = 1.716 and RRMSE = 4.504%). The newly developed UAV- and machine learning-based model provided a promising and real time method to monitor chlorophyll contents at the overwintering stage, which can beneﬁt late-sown winter wheat variety screening.


Introduction
Wheat is one of the three major food crops in China. Jiangsu Province, located in the lower reaches of the Yangtze River, is one of the major wheat production areas in China. In recent years, the delay in rice maturity in this region has led to a significant delay in subsequent winter wheat sowing. For example, the percentages of late-sown (seven or more days later than local normal seed sowing date) wheat area in Jiangsu Province were 48.6%, 51.2%, and 59.6% in 2015, 2016, and 2017, respectively [1]. The delay in sowing has become a major obstacle to high wheat yield; therefore, it is necessary to screen wheat varieties suitable for late sowing. For example, optimal wheat varieties should be able to maintain a certain amount of growth even at the overwintering stage [2]. In addition, wheat in the lower reach of the Yangtze River often suffers from low-temperature frost damage at the overwintering stage, which severely affects wheat growth and development [3]; hence, good wheat varieties should have strong resistance to low temperatures. Therefore, accurately monitoring wheat growth status at the overwintering stage is critical for latesown winter wheat variety screening.
As an important pigment for photosynthesis, chlorophyll has a critical impact on a plant's ability to exchange material and energy with the external environment. Chlorophyll content can indicate crop growth status, primary productivity, and nitrogen use efficiency [4]. Exposure to various kinds of stresses may reduce crop chlorophyll content [5]; so, chlorophyll content can provide information on a wheat variety's tolerance capability to endure stresses due to late sowing, low temperature, insufficient nitrogen application and so on.
The traditional methods for measuring chlorophyll content are usually time-consuming, laborious, and destructive to crop leaves [4]. Although the Soil and Plant Analysis Development (SPAD) method is non-destructive, it can work only at limited measuring points, and cannot provide the spatially continuous distribution of SPAD [6].
Late-sown wheat variety screening requires a large number of experimental plots. Better monitoring wheat growth status may benefit from rapid, accurate, and non-destructive estimation of wheat chlorophyll content in each plot. Remote sensing provides a great potential for chlorophyll estimation over large regions [7].
Spectral vegetation indices (VIs) have been widely used to estimate vegetation chlorophyll content from spectral data [8]. Good vegetation indices are able to maximize sensitivity to the vegetation characteristics, while reducing the spectral effects due to atmosphere, soil background, topography, and sensor view angle [9,10].
The VI of modified chlorophyll absorption in reflectance index (MCARI) that was proposed by Daughtry et al. [11] is based on the reflectance of green, red, and red edge bands, and this VI is sensitive to leaf chlorophyll variations. Using the reflectance of green, red edge, and near-infrared (NIR) bands, Cao et al. [12] modified MCARI to propose a new VI of modified chlorophyll absorption reflectance index 1 (MCARI1). They indicated that the MCARI1 displayed quite significant correlations with rice's above-ground biomass and plant nitrogen uptake at each growth stage of rice.
The Medium Resolution Imaging Spectrometer (MERIS) terrestrial chlorophyll index (MTCI) that was designed by Dash and Curran [13] is suitable for the estimation of chlorophyll content from the MERIS data. MTCI has become a frequently used VI to monitor spatial variability of crop chlorophyll [14].
A 2-band VI of green chlorophyll vegetation indices (GCVI) that was proposed by Gitelson et al. [15], using green and NIR bands, had good correlations with chlorophyll content in maize and soybean canopy. The developed models could provide accurate estimation of canopy chlorophyll contents, although calibration coefficients were different for maize and soybean.
The chlorophyll vegetation index (CVI) that was developed by Vincini et al. [16] uses the reflectance of green, red, and NIR broad bands. Vincini et al. [16] indicated that CVI was specifically sensitive to leaf chlorophyll content at the canopy scale of sugar beet.
Since hyperspectral data are composed of a large number of continuous and narrow bands, a number of spectral indices have been proposed to estimate the chlorophyll content of vegetation [17]. Main et al. [18] developed two red edge derivative based indices (i.e., red edge position via linear extrapolation index and the modified red edge inflection point index), and found that the two indices were consistent and robust in chlorophyll content estimation in three crop species and a variety of savanna tree species. Jin et al. [19] developed a spectral index of double-peak canopy nitrogen index I (DCNI I), and they indicated that DCNI I produced accurate estimation of chlorophyll content in cotton. However, hyperspectral images are difficult to obtain. These hyperspectral VIs cannot be calculated from broadband multispectral data that are more readily available.
The spectral indices that were proposed for SPAD estimation use different spectral bands and have very different equations. Their performances vary among different studies, and none of them performed the best in all studies. More importantly, while most previous studies in remote sensing of SPAD only involved single variety or a very small number of varieties of a crop, very little research focused on the relationships between SPAD and spectral indices that were impacted by a large number of varieties. Accurate estimation of SPAD of various varieties is critical for late-sown wheat variety screening.
Traditional algorithms such as simple or multiple linear regressions have often been used in remote sensing for crop biophysical and biochemical variable retrieval. In recent years, machine learning algorithms (MLAs) have been employed increasingly. Different from traditional algorithms, MLAs are data-driven, and they are able to autonomously cope with linear correlations as well as solve strong nonlinear problems possessed by agricultural and remote sensing variables [20].
Among a number of powerful MLAs, support vector regression (SVR), random forest regression (RFR), and Artificial Neural Network (ANN) are most frequently used for agricultural remote sensing [21]. Yang et al. [22] estimated green leaf chlorophyll density of rice from hyperspectral reflectance measured over two experimental rice fields containing two cultivars treated with three levels of nitrogen application. They found that SVR, the regression version of support vector machine (SVM), largely improved the estimation accuracy in comparison with the stepwise multiple regression (SMR). They indicated that SVR deals better than traditional regression algorithms with non-linear processes that exist in the relation between green leaf chlorophyll density and spectral data.
Similarly, an SVR-based model developed using the canopy spectral reflectance of maize, measured by a handheld spectrometer, was demonstrated to be able to estimate the chlorophyll content of maize canopy in the field non-destructively and rapidly [23]. Cavallo et al. [24] used RFR to predict total chlorophyll content of fresh-cut rocket leaves from spectral data measured using a spectrophotometer. The developed RFR model could provide accurate estimation of total chlorophyll content.
Although ANN is also widely used to estimate agricultural variables, it is not as practical as SVR and RFR because it requires complex and time consuming procedures such as selecting the number and size of hidden layers, setting the learning rate, obtaining a large training dataset, and dealing with the problem of overfitting [20].
Prior to modeling, it is important to determine which variables should be included in the model, using a suitable variable selection method. A good variable selection method can select the smallest number and most efficient subset of variables from the original set, which can improve the estimation power of the model and speed up the model execution time [25].
Random forest (RF) is not only a regression algorithm but also a frequently used variable selection method. This is because the importance of each variable can be calculated using RF and ranked during RF model development [26]. Recursive feature elimination (RFE) is another widely used method for variable selection. The method uses a base model to perform multiple rounds of training, during which the weakest features are eliminated until a specified number of features is reached [27,28]. The base model could be SVM [27], RF [29], or other regression models. Nevertheless, there have been few previous studies that used variable selection in remote sensing of vegetation chlorophyll contents, although variable selection is critical. To fill this gap, the current study was designed to develop a good combination of variable selection methods with machine learning algorithms that can accurately estimate SPAD of wheat canopy at the overwintering stage, by comparative analysis. This study was based on the unmanned aerial vehicle (UAV) technology because the wheat variety screening experiments required high spatial resolution imagery.

Experimental Site and Experimental Design
The experiment was conducted at the Shatou Agricultural Experimental Farm located in Yangzhou, Jiangsu, China (32 • 18 39.96 N, 119 • 32 45.34 E) during the 2020/2021 wheat growing season (Figure 1). In total, 24 winter wheat varieties were sown under four rates of nitrogen fertilization through top-dressing (i.e., N0: 0 kg ha −1 , N210: 210 kg ha −1 , N270: 270 kg ha −1 , and N330: 330 kg ha −1 pure nitrogen) on 15 November 2020, which was 10 days later than the local normal date for seed sowing. There were 96 plots in this experiment, and the size of each plot was 3 m long and 3 m wide. Each plot received the same irrigation and field managements.
The 24 winter wheat varieties are listed in Table 1. They had different optimal growth habits [30].  The imagery of all plots in the experimental site was collected by the Phantom 4 Multispectral UAV (SZ DJI Technology Co.; Shenzhen, China) under clear weather conditions on 13 January 2021, at the overwintering stage of the winter wheat. The UAV was designed for precision agriculture equipped with a real-time kinematic (RTK) module capable of increasing the accuracy of Global Positioning System (GPS) location data.
Utilizing a gimbal with a 6-camera array, the UAV carried a RGB camera as well as a five-single-sensor camera that captured 2MP images in blue (450 nm ± 16 nm), green (560 nm ± 16 nm), red (650 nm ± 16 nm), red edge (730 nm ± 16 nm), and NIR (840 nm ± 26 nm) spectral regions. This study only used the images acquired by the five single sensors.
A gray gradient calibration panel, having 10 calibration plates with different gray gradients, was placed on the ground while the UAV was flying, for reflectivity correction [31]. The reflectance of these plates were spectrally measured in situ using an Analytical Spectral Device (ASD) FieldSpec ® 3 Full-Range spectroradiometer (Analytical Spectral Devices, Inc.; Boulder, CO, USA).
The flight routes were designed using the software of the ground station (DJI GS PRO). Both the forward and side overlaps were set at 80%. The flight was conducted at 10:00 AM local time. The flight speed was fixed at 3 m s −1 , and the flight altitude above ground level (AGL) was 15 m, which returned a spatial image resolution (i.e., ground sample distance) of 0.59 cm.

SPAD Measurements
SPAD values of wheat leaves were measured in situ in each plot using a non-destructive and portable SPAD-502plus handheld chlorophyll meter (Minolta Camera Co.; Osaka, Japan) immediately after UAV overpass. Within each plot, 10 leaves of the top layer were randomly selected. For each selected leaf, SPAD readings were recorded at one, three, and five sixths of leaf length of the leaf, and then the average was calculated and used for this leaf. The SPAD readings of the 10 selected leaves were then averaged to represent the plot.

Machine Learning Regression Algorithms Used
As mentioned above, RFR, SVR, and ANN are most frequently used for agricultural remote sensing [21], but ANN is not as practical as RFR and SVR [20]. Hence, this study employed RFR and SVR in wheat canopy SPAD estimation.

Random Forest Regression (RFR)
The random forest (RF) model is a classification or regression tree-based machine learning method proposed by Breiman [32]. The RF model uses bootstrap resampling methods to extract multiple bootstrap datasets from the original training dataset, develops a decision tree for each bootstrap dataset, and then combines these multiple parallel decision trees into a single RF model for prediction. The final prediction is obtained by applying a majority voting decision mechanism to the output of all individual decision trees. Stochasticity is introduced by the bootstrap strategy, which improves the resistance of the RF to overfitting [33]. When RF is used for regression, it is referred to as random forest regression (RFR). For RFR models, parameter tuning is required to optimize its two main parameters, namely the number of trees to grow in the forest (ntree) and the number of randomly selected predictor variables at each node (mtry) [34].

Support Vector Regression (SVR)
A support vector machine (SVM) or its regression version, support vector regression (SVR), is a supervised machine learning algorithm, which is powerful and flexible for classification and regression [35]. Using an optimal kernel, SVM maps the input data into different classes in a hyperplane in multidimensional space in an iterative manner until it finds a maximum marginal hyperplane, in which the differences among classes are maximized so as to minimize the error of classification [34].
Among the model parameters of SVR, gamma (kernel parameter) and C (penalty coefficient) have great influences on the estimation results, so they should be optimized through parameter tuning [17,27,36]. In addition, the kernels of linear, radial basis function (RBF) and sigmoid are commonly used; they were applied and tested in this study.

UAV Image Processing and Index Extraction
The acquired UAV images were imported into the DJI Terra software to generate a digital surface model (DSM) and an orthomosaic image of the experimental plots in GeoTIFF format with UTM projection. The DN values were extracted from the orthomosaic image using the Environment for Visualizing Images (ENVI) software (ITT Exelis; Boulder, CO, USA) for each calibration plate, and then regression relationships between the DN values and their reflectance values were constructed for each of the five bands. Based on these regression relationships, DNs were converted to reflectance using the simplified empirical line approach proposed by Wang and Myint [37].
Then, the reflectance values at the blue, green, red, red edge, and NIR bands were extracted for each pixel in the entire UAV image of the experimental site. These reflectance values were used to further calculate 26 spectral indices. These spectral indices were often used to estimate crop chlorophyll contents and monitor crop growth status in the literature. In this study, a total of 31 variables, including 26 spectral indices and reflectance of five single bands, were used to develop a machine learning-based model to estimate SPAD values of wheat canopy. These 31 spectral variables and related computation formulas are listed in Table 2.
Difference between green and blue (DifGB) In general, the pixels with normalized difference vegetation index (NDVI) values of 0.2-0.3 were mixed pixels of non-vegetation and vegetation [55][56][57]. A careful visual check of the image revealed that the NDVI value of 0.25 could better discriminate wheat from soil, so the pixels with NDVI values lower than 0.25 were masked out in this study. Thus, the 31 variables were extracted and averaged as the variable values of the plot, respectively.

Variable Selection
The variable selection method of recursive feature elimination (RFE) starts working by searching a subset of variables from all variables in the training dataset, removes the least important variables, and then uses the remaining variables to refit the base model [27,58]. This process is repeated until a specific number of variables are retained, and the order in which variables are eliminated in this process is the ranking of the variables. This study tested using RFR and linear kernel SVR as a base model to conduct the RFE method, respectively. For RFE, its base model cannot be RBF kernel-based or sigmoid kernel-based SVR models [59].
RF is a fast and efficient machine learning algorithm, even when dealing with noisy variables [26]. RF is also a variable selection method [60]. It provides an easy way to assess variable importance. As a decision tree-based machine learning algorithm, RF can output a feature importance attribute; hence, all variables can be ranked in order of their own importance values calculated in the process of RF modeling [26].
The Pearson correlation coefficient (r) indicates a linear correlation between a predictor variable and a target variable, and its value ranges between −1 and 1. The predictor variables with higher absolute value of r have stronger linear correlation with the target variable, and hence they are selected prior to those predictor variables with lower absolute value of r in the variable selection method of Pearson correlation coefficient (r).
This study employed RFE, RF, and r to select variables. They are also most commonly used in previous studies [27].

Modeling, Cross Validation, and Performance Assessment
The 31 spectral variables (Table 2) were employed for feature selection after being standardized using the standardized method of StandardScaler. The standardized method was from a Python library named "sklearn.preprocessing" [61].
To obtain the most effective model for estimating the SPAD of wheat canopy through comparative analysis, this study combined three variable selection methods with four machine learning algorithms, resulting in 10 SPAD estimation models, which included seven SVR models (i.e., RFE-SVR_linear, RF-SVR_linear, RF-SVR_RBF, RF-SVR_sigmoid, r-SVR_linear, r-SVR_RBF, and r-SVR_sigmoid) and three RFR models (i.e., RFE-RFR, RF-RFR, and r-RFR).
Parameter tuning plays an important role in achieving the best performance of machine learning-based models [35], and this study combined grid search techniques with cross-validation to find the best parameter values. For each of the 10 models, the optimal number of variables (n_features_to_select) was determined by searching values from 1 to 31 in an interval of 1 and selecting the number of variables that achieves the best cross-validation accuracy.
To optimize the RFR models, the number of input variables was used as the mtry value. Nine values (i.e., 50, 100, 200, 700, 1000, 1100, 1200, 1300, and 2000 trees) were used as the ntree value, respectively, and the ntree value resulting in the best cross-validation accuracy was selected as the optimal ntree.
To optimize the SVR models, the penalty C was evaluated using values from 1 to 500 with an interval of 1, and the parameter gamma was searched with values from 0.0001 to 0.01, with an interval of 0.0002. The C and gamma achieving the best cross-validation accuracy were selected as the optimal parameter values. Figure 2 summarizes the methodology used. The entire process of variable selection, modeling, cross-validation, and performance evaluation was performed with a desktop computer using a Python program written for this study based on the Scikitlearn package [62], which is an open source Python module. The optimal values of the above mentioned key parameters such as n_features_to_select, mtry, ntree, C, and gamma were determined simultaneously according to the SPAD estimation accuracy, using crossvalidation and the grid search strategy. Correspondingly, the SPAD estimation model with the highest estimation accuracy was considered as the optimum SPAD estimation model of wheat canopy. Among various cross-validation techniques, this study selected the leave-one-out cross-validation (LOOCV). The LOOCV [35] selects one sample each time for validation and uses all the other samples to develop a model, and then uses the selected one sample to calculate the estimation errors until all samples are involved in the cross validation. The LOOCV can assess the generalization capability of the models [35] and eliminate overfitting [63], and it was used to evaluate the model reliability and stability.
The accuracy of models was evaluated by achieved R 2 , RMSE, and relative RMSE (RRMSE) in SPAD estimation. The LOOCV results of R 2 , RMSE, and RRMSE were calculated following Wang and Lu [64].
In general, a higher R 2 value and lower RMSE and RRMSE values indicate better model performance [65]. In addition, this study also calculated the ratio of percentage deviation (RPD) to assess the predictive capability of models. RPD is the ratio of the standard deviation of the measured values to the RMSE of the cross-validation [66]. According to Viscarra Rossel et al. [67], RPD < 1.4 indicates very poor or poor estimations, 1.4 ≤ RPD < 1.8 indicates fair estimations, 1.8 ≤ RPD < 2.0 indicates good estimations, 2.0 ≤ RPD < 2.5 indicates very good estimations, and RPD ≥ 2.5 indicates excellent estimations.  The 96 plots were divided to four groups (i.e., N0 plots, N210 plots, N270 plots, and N330 plots) according to their rates of nitrogen fertilization (i.e., 0, 210, 270, and 330 kg ha −1 pure nitrogen), respectively. Table 3 demonstrates that the mean SPAD values increased significantly from 33.80 to 38.31 from N0 plots to N210 plots, and then increased slightly to 39.73 of N270 plots and to 40.54 of N330 plots. In addition, the standard deviation of N0 plots was also smaller than those of N210 plots, N270 plots, and N330 plots.

Spectral Characteristics of Wheat Canopy SPAD Values
The correlation coefficients between 31 UAV-derived variables and SPAD values of wheat canopy were calculated ( To compare spectral curves for varying SPAD values, this study averaged the SPAD and spectral reflectance values for the four groups (i.e., N0 plots, N210 plots, N270 plots, and N330 plots). As illustrated in Figure 4, the spectral reflectance of the four groups had similar variations with increasing wavelength, although they had different mean SPAD values. Moreover, in general, reflectance increases with increasing SPAD, especially in the red edge and near-infrared regions.  (Table 2) and SPAD values of wheat canopy. The orange colour represents the variables that included only the visible spectral bands, the green colour represents the variables that included the near-infrared spectral band but not the red edge spectral band, and the blue colour represents the variables that included the red edge spectral band.

SPAD Inversion Models
For each model, the best model parameters were retrieved using the grid search technique, according to the resultant lowest RMSE by the LOOCV. These optimal parameters are listed in Table 4. The selected optimal variables of each model are listed in Table 5. The optimal variables were quite different among the 10 models. MCARI1 was the only variable that was commonly selected by all 10 models. The three RFR models had three common selected variables (i.e., MCARI1, MTCI, and CIRE). The seven SVR models had five common selected variables (i.e., MCARI1, VARI, GI, RGRI, and MCARI).  Table 6 displays the accuracies of the 10 models for estimating the SPAD. RF-SVR_sigmoid achieved the highest estimation accuracy (Optimal number of variables = 27, R 2 = 0.754, RPD = 2.017, RMSE = 1.716, and RRMSE = 4.504%), followed by r-SVR_sigmoid (Optimal number of variables = 31, R 2 = 0.743, RPD = 1.973, RMSE = 1.754, and RRMSE = 4.605%), and RFE-SVR_linear (Optimal number of variables = 7, R 2 = 0.730, RPD = 1.926, RMSE = 1.797, and RRMSE = 4.718%). Table 6. Optimal number of variables and estimation accuracies of the developed 10 SPAD estimation models from the LOOCV.

Model
Optimal To further analyze the modeling accuracy, Figure 5 presents the scatterplots of the measured SPAD against the estimated SPAD by the LOOCV, when the three models that achieved the best estimation performance were applied. For all three models, the data points are well distributed along the 1:1 relationship, demonstrating good agreements between the measured and estimated SPAD values.
This study compared the accuracies of the selected best model, RF-SVR_sigmoid, in estimating wheat canopy SPAD for plots under four different rates of nitrogen fertilization (Table 7). RMSE ranged from 1.500 (N0 plots) to 2.017 (N210 plots), and RRMSE ranged from 4.054% (N270 plots) to 5.265% (N210 plots). Table 7. Comparison of SPAD estimation accuracy between plots under four different rates of nitrogen fertilization. The selected best model, RF-SVR_sigmoid, was used to estimate SPAD by the LOOCV.

Plots
Number This study also compared the SPAD estimation accuracies of the RF-SVR_sigmoid model for wheat varieties that have three different optimal growth habits ( Table 8). RMSE ranged from 1.486 to 2.088, and RRMSE ranged from 3.866% to 5.176%, with the lowest in the N270 plots and the highest in the N210 plots. The semi-winterness varieties had the lowest RMSE and RRMSE, followed by the springiness varieties and the weak springiness varieties.

Optimal SPAD Estimation Model for Late-Sown Winter Wheat Variety Screening
This study involved up to 24 wheat varieties, four rates of nitrogen fertilization, and 96 plots in a variety screening experiment, which resulted in a very complex relationship between wheat canopy SPAD values and spectral variables. Therefore, this study investigated the feasibility of using machine learning-based models rather than conventional ones to estimate SPAD of wheat canopy from UAV data at the overwintering stage. Seven SVR models and three RFR models were developed, cross-validated, and compared. Of these 10 models, the RF-SVR_sigmoid model, which was constructed by combining the RF variable selection method and the sigmoid kernel-based SVR algorithm, was identified as the best model. It achieved the highest accuracy in estimating SPAD values by the LOOCV (R 2 = 0.754, RPD = 2.017, RMSE = 1.716, and RRMSE = 4.504%). According to Viscarra Rossel et al. [67], the model of RF-SVR_sigmoid (RPD = 2.017) achieved very good SPAD estimation of wheat canopy. In contrast, r-SVR_sigmoid (RPD = 1.973) and RFE-SVR_linear (RPD = 1.926) only produced good estimation, but their RPD values were very close to 2.
Despite the fact that the study involved up to 24 different wheat varieties, the newly developed RF-SVR_sigmoid model was able to achieve high accuracy in SPAD estimation. Moreover, the model worked well for both plots under different rates of nitrogen fertilization and plots with springness, weak springness, and semi-winterness varieties. It is particularly important for wheat variety screening that the developed model can be used for all the wheat varieties involved, considering that the best varieties are usually screened from a large number of varieties. In contrast, most previous studies of SPAD remote sensing involved only a single variety or a very small number of varieties, which is not suitable for variety screening [68][69][70].
Given the encouraging performances of the SVR, however, it is too early to say that the SVR model always outperforms the RFR model in estimating vegetation chlorophyll content. For the estimation of other agricultural variables, some studies had reported that the RFR model outperformed the SVR model. For example, Osco et al. [71] noted that, using UAV multispectral images, the RFR model was able to predict leaf nitrogen content (LNC) of maize more accurately than the SVR model. Similarly, Zha et al. [65] reported that the RFR algorithm performed better than the SVR and ANN models in estimating the nitrogen nutrient index (NNI) of rice from drone data.
Yang et al. [22] indicated that appropriate kernels could better avoid overfitting. The RBF kernel was considered to be the most frequently used kernel function [72]. Some previous studies (e.g., Ahmad et al. [73]; Chen and Hay [74]) reported that the RBF kernel performed better relative to linear and sigmoid kernels. However, this study found that among the seven SVR models, the two sigmoid kernel-based SVR models performed the best, followed by the RFE-SVR_linear model, while the two RBF kernel-based SVR models produced lower accuracy. Hence, the sigmoid kernel seemed more appropriate for SVR in terms of SPAD estimation.
In addition, although the RFE-SVR_linear model was not as good as the two SVR models with sigmoid kernel, its SPAD estimation was also good. Moreover, the linear kernel-based SVR model runs much faster than the RBF kernel-based SVR model [75].

Influence of Different Variable Selection Methods on Model Estimation Performance
Little previous research has investigated employing variable selection in machine learning-based remote sensing of vegetation chlorophyll contents. This research found that the optimal combination of variable selection methods and machine learning algorithms could produce more accurate SPAD estimation of wheat canopy. Among various combinations of three variable selection methods and four machine learning algorithms, the combination of RF and SVR_sigmoid demonstrates the best capability of SPAD estimation.
In this study, the RFR and SVR_linear models using the RFE variable selection method provided overall higher R 2 and more robust results than the RFR and SVR_linear models using RF or r variable selection methods. In addition, the number of optimal variables selected using RFE was much smaller than those selected using RF or r. Results from this study disagree with the superiority of RF over RFE as a variable selection method, as reported by Chen et al. [27]. More research should be conducted to further evaluate these variable selection methods in agricultural remote sensing.
Using RF to select the optimal variables resulted in more accurate SPAD estimates compared with using the r variable selection method (Table 6). However, the differences caused by using different variable selection methods are not as large as those caused by using different machine learning algorithms.
Three models produced the lowest RMSE and RRMSE (i.e., RFE-SVR_linear, RF-SVR_sigmoid, and r-SVR_sigmoid), as displayed in Table 6. They had seven common optimal variables (i.e., green, VARI, GI, RGRI, RI, MCARI1, and MCARI) ( Table 5). This indicates that the seven variables are particularly important for accurate estimation of wheat canopy SPAD. The seven variables were selected from the 31 variables (Table 2), that included 10 variables (i.e., NPCI, GCVI, CVI, MCARI1, MCARI, MTCI, NCARI, TCI, TCARI, and CIRE) proposed and widely used for monitoring chlorophyll or SPAD in previous studies [16,40,54]. It is a bit surprising that, among the 10 variables, only two (i.e., MCARI1 and MCARI) are commonly selected by the three best SPAD estimation models in this study.

Limitations and Future Research
Although this preliminary study proposed a promising method for SPAD monitoring, there were still some limitations. This study found that the NDVI value of 0.25 could discriminate wheat from soil. When the new method is applied on other dates or on other farms, the threshold should be determined through a very careful visual check of the entire image.
Besides chlorophyll content, canopy reflectance is also sensitive to other influencing factors, such as canopy structure and leaf area index [9,15]. The error in the SPAD estimation of the RF-SVR_sigmoid model could be partially attributed to the omission of these influencing factors. Considering these influencing factors may further improve the SPAD estimation accuracies in future research.
This study applied 96 plots, and future research will employ a larger data base. In addition, this study involved 18 springness wheat varieties, and in contrast it involved only two weak springness varieties and four semi-winterness varieties. Hence, the wheat variety shares should be balanced in future research.
This study used only data from one growing season (i.e., the overwintering growth stage). Future research should collect data at more growth stages. In addition, this study involved only a single year, but variety screening often needs multiple-year experiments. Hence, the conclusions drawn in this study should be further evaluated in future studies with data from multiple growth stages and more years on various experimental farms.

Conclusions
This study demonstrates that the combination of different variable selection methods with different machine learning algorithms can impact the estimation accuracy quite significantly. The optimal combination of variable selection methods and machine learning algorithms could produce more accurate SPAD estimation. The newly developed RF-SVR_sigmoid model, which is the combination of the RF variable selection method and the sigmoid kernel-based SVR algorithm, can accurately estimate SPAD of wheat canopy at the overwintering stage. The results of LOOCV prove that the estimation is very good, and the newly developed model can be used for all 24 wheat varieties involved.
The spectral data used in this study were acquired using a multispectral UAV. This study demonstrates the potential of UAV remote sensing for late-sown winter wheat variety screening. Use of UAVs is more efficient, time-saving, and convenient than traditional in situ measurements of SPAD. In comparison to satellite remote sensing, use of low-altitude UAVs is more flexible, and it provides imagery with much higher spatial resolutions, which is necessary considering the size of plots in variety screening experiments. The widespread use of UAV technology nowadays allows a fast and accurate monitoring of wheat canopy SPAD at the overwintering stage using the newly developed RF-SVR_sigmoid model, and this is critical for late-sown winter wheat variety screening.