Predicting Equivalent Water Thickness in Wheat Using UAV Mounted Multispectral Sensor through Deep Learning Techniques

Abstract

The equivalent water thickness (EWT) is an important biophysical indicator of water status in crops. The effective monitoring of EWT in wheat under different nitrogen and water treatments is important for irrigation management in precision agriculture. This study aimed to investigate the performances of machine learning (ML) algorithms in retrieving wheat EWT. For this purpose, a rain shelter experiment (Exp. 1) with four irrigation quantities (0, 120, 240, 360 mm) and two nitrogen levels (75 and 255 kg N/ha), and field experiments (Exps. 2–3) with the same irrigation and rainfall water levels (360 mm) but different nitrogen levels (varying from 75 to 255 kg N/ha) were conducted in the North China Plain. The canopy reflectance was measured for all plots at 30 m using an unmanned aerial vehicle (UAV)-mounted multispectral camera. Destructive sampling was conducted immediately after the UAV flights to measure total fresh and dry weight. Deep Neural Network (DNN) is a special type of neural network, which has shown performance in regression analysis is compared with other machine learning (ML) models. A feature selection (FS) algorithm named the decision tree (DT) was used as the automatic relevance determination method to obtain the relative relevance of 5 out of 67 vegetation indices (Vis), which were used for estimating EWT. The selected VIs were used to estimate EWT using multiple linear regression (MLR), deep neural network multilayer perceptron (DNN-MLP), artificial neural networks multilayer perceptron (ANN-MLP), boosted tree regression (BRT), and support vector machines (SVMs). The results show that the DNN-MLP with R² = 0.934, NSE = 0.933, RMSE = 0.028 g/cm², and MAE of 0.017 g/cm² outperformed other ML algorithms (ANN-MPL, BRT, and SVM- Polynomial) owing to its high capacity for estimating EWT as compared to other ML methods. Our findings support the conclusion that ML can potentially be applied in combination with VIs for retrieving EWT. Despite the complexity of the ML models, the EWT map should help farmers by improving the real-time irrigation efficiency of wheat by quantifying field water content and addressing variability.

1. Introduction

Wheat production accounts for nearly 50% of China’s National Agricultural output. Increasing wheat consumption requires effective decision-making during the wheat growth period. Classically, field data collection has been employed to diagnose plant biophysical parameters, including equivalent water thickness (EWT) [1]. Improving crop water management requires the accurate and timely monitoring of water in the plant [2]. EWT has previously been used to derive the expected grain yield [3] by managing irrigation scheduling [4]. Although the classical method is accurate, it is time-consuming, laborious, and does not indicate the spatial variability of EWT in the entire field [5,6,7]. Airborne remote sensing has been used to detect plant water variability. Recently, UAV-mounted multispectral, hyperspectral, and thermal sensors have been used to acquire very high spectral and spatial resolutions [8]. Many studies have been carried out to explore accuracy in the retrieval of crop biophysical parameters using the UAV remote sensing platform [8,9,10,11,12,13]. The use of the UAV platform not only allows non-destructive and timely data acquisition, but it also captures variability at the field scale as well. Therefore, it is also possible to use UAV-based remote sensing techniques to increase the accuracy of EWT retrieval by reducing the assessment time.

In addition, EWT have been assessed using the broad waveband [14], narrow waveband in the near-infrared (NIR), and shortwave infrared (SWIR) regions [15,16,17,18]. Ceccato et al. [15,17] confirmed that the wavelength absorption bands of plant and water molecules are in the region of 900 to 2500 nm. The reflectance spectrum of green vegetation in this region is reported to be affected by strong liquid water absorption. Tucker [19] suggested that the 1555–1750 nm region was best-suited for the remote sensing of plant canopy water status using satellite platforms. In the past, data analysis techniques such as laboratory spectral reflectance and remote sensing platforms using hyperspectral sensors have been used to estimate EWT [1,20,21,22,23]. Physical-methods have shown that the simple ratio between 1600 and 820 nm is heavily influenced by EWT [15]. Recently, multispectral sensors with higher spatial and spectral resolutions have been used for site-specific crop-management, which provides a reliable and effective remotely sensed source of information for agriculture. For instances, multispectral data have been used to assess chlorophyll content [24,25] and plant nitrogen content [26,27] and their wavelengths are located in the 500 to 800 nm region. However, the most commonly used multispectral vegetation indices (VIs) have limitations when used for predicting plant water content via a linear or multilinear model. In a previous study undertaken by Thornton et al. [28], the change in plant nitrogen concentration was explained by changes in shoot water content, suggesting that changes in chlorophyll or nitrogen content affect plant water content. Poblete et al. [13] have suggested that several indices, when applied to information between 500 and 800 nm, can help to indirectly estimate water status using nonlinear methods.

Estimating EWT at the canopy level using optical remotely sensed data can be achieved via several methodologies, including the empirical method based on hyperspectral VIs [29]. However, the indices generated to estimate EWT in this context yield poor performance. Physical models have also been used to assess EWT from remotely sensed data, based on physical laws such as radiative transfer models that describe the transfer and interaction of radiation within the atmosphere and plant canopy [30,31]. The fact that the site-specific information requirement for a proper model parameterization, which is not always available, can be the biggest disadvantage of radiative transfer modeling. The VIs-based and physical models are critically limited in their ability to accurately estimate EWT. Consequently, alternative methods for the retrieval of EWT using multispectral VIs are needed. The machine learning (ML) method is an advanced computational method that has shown good performance in retrieving biophysical parameters from crops using multispectral and hyperspectral reflectance data. The ML models have been used to produce accurate and robust models in engineering [32], agriculture [33,34,35], hydrology [36], and forestry [37]. ML is a step forward to implement artificial intelligence without the need for explicit programming while deep learning, a subset of machine learning algorithm, makes intuitive and intelligent decisions using a neural network stacked layer-wise. ML algorithms such us neural networks (NN) are known to learn the underlying relationship between the input and the output data [38]. NN is one of the powerful learning algorithms for standing exemplary performance in regression problems [39,40,41]. The success of NN captured the attention of researchers by exploring the possibility to train NN with several hidden layers. However, the success of the NN was limited to one or two hidden layers. For instance, the difficulties of training NN such as artificial neural network (ANN) with several hidden layers lie in vanishing gradients or exploding gradients with the increment of the hidden layer number (depth of network). The introduction of a deep neural network (DNN) smoothed the issue of vanishing or exploding gradients, thereby speeding up training [42]. The attractive feature of DNN is its ability to exploit new activation functions and learning algorithms. In addition, deep architectures often have an advantage over ANN architectures when dealing with complex learning problems [38].

However, ML techniques have not yet been fully employed to bridge the knowledge gap between VIs from the visible infrared region and EWT. Consequently, it is important to define an ML method to accurately assess vegetation EWT using only VIs. Despite the advantages of using ML, they present some major drawbacks. ML models can include any superfluous input variable, increasing the complexity and providing erroneous information about the variables that actually affect the model’s performance [43,44]. To solve these problems, dimensionality reduction approaches, such as feature selection (FS) methods, are applied to reduce the number of input variables to only the relevant ones. As such, FS is used to clean up the redundant, irrelevant, and noisy data. As a result, the performance is boosted. In the FS, a subset of features is selected from the original set based on their redundancy and relevance [44]. In recent studies [45,46,47], FS methods such as filter, wrapper, and embedded methods have been used based on their interaction with the learning model. The filter and wrapper have been widely used in various remote sensing applications; nevertheless, their large computational time is a major drawback. Embedded methods, such as decision tree (DT) algorithms, use ensemble learning and hybrid learning methods for FS [48]. The DT offers excellent sparseness performance compared to the filter and wrapper methods. It is less computationally intensive than the other two methods [47].

The EWT varies with vegetation growth, and thus estimating EWT across the field at different growth stages using the ML method and VIs could offer the farmer critical time- and location-specific information for monitoring their crops and managing farming activities, thus enabling them to increase yields. Therefore, the objective of the present study was to investigate the feasibility of different ML algorithms, such as DNN, ANN, boosted regression tree (BRT), and support vector machines (SVMs) models, to predict wheat EWT using multispectral single-band images and VIs. The key idea behind using several ML methods is to take advantage of each method’s capability to predict EWT. We thus investigated an FS algorithm-based DT to determine best VIs for use as input parameters in the different machine learning models.

2. Materials and Methods

2.1. Description of The Study Area and Experimental Design

Three experiments with varied water and nitrogen rate were conducted from September 2019 to June 2020 at the Agricultural Station of the Chinese Academy of Agricultural Sciences (CAAS), located in Qiliying county of the North China Plain (Figure 1). For all experiments, UAV and destructive sampling were performed from March to May in 2020. This area is categorized by a subtropical warm climate with rainfall varying between 900 and 1200 mm, about 65 to 70% of which occurs between June and September. The average solar radiation is about 4900 MJ m⁻² yr⁻¹. The mean annual temperature is about 14.5 °C. Exp. 1 was conducted in a rain-out shelter facility using a split-plot design with four treatments of different irrigation quantities (0, 120, 240, 360 mm) and two nitrogen (N) levels (75 and 255 kg N ha⁻¹). In Exp.1, twenty-four plot experiments were conducted in the irrigation treatment involving one wheat cultivar, Zhoumai27. All of the plot experiments adopted a randomized complete block design. The size of each plot was 3.5 m × 1.9 m and this did not change during the study period. N fertilizer was applied via a two-plot split design with three replications, and the N was only applied at the top dressing. In addition to the plot experiments, two field experiments, Exps. 2–3 were conducted under field conditions to precisely compare the wheat water statuses. The plot size for each treatment varied from 14 × 10 m to 15 m × 10 m. The wheat cultivar was Zhoumai22 and each treatment was performed once. All the plots received a basic N application of 75 kg N ha⁻¹, and at the top dressing the application rate varied from 0 to 150 kg N ha⁻¹. Previous field plot experiments have had various objectives, but this study took advantage of the varying N status to evaluate different UAV remote sensing-based water status estimation methods. Areas of 0.36 m² were selected as sampling areas in each experiment. However, the samples in Exps. 2–3 were taken at three different locations, and the average value was taken.

The description of all the abbreviated parameters and description of parameters and variables used have been summarized in

Table 1. Description of the abbreviations.

Name	Description
ANN-MLP	artificial neural networks–multilayer perceptron
BRT	boosted tree regression
DT	decision tree
DNN-MLP	deep neural network–multilayer perceptron
FS	feature selection
ML	machine learning
NIR	near-infrared region
NN	neural network
SWIR	short-wavelength infrared region
SVMs	support vector machines
UAV	unmanned aerial vehicle
VIs	vegetation indices

and

Table 2. Description of parameters and variables.

Parameters/Variables	Description	Unit
DW	dry weight	t·ha−1
EWT	equivalent water thickness	g·cm−2 or cm
FW	fresh weight	t·ha−1
MAE	mean absolute error	g·cm−2 or cm
MLR	multiple linear regression	g·cm−2 or cm
NSE	Nash–Sutcliffe efficiency
RMSE	root means square error
R2	determination coefficient
$X_n$	input variable
${\beta }_n$	regression coefficients associated $X_n$ input variable
$X_{norm}$	normalized value of the input variable
$X_i$	real value of the input variable
$X_{min}$	minimum input variable
$X_{max}$	maximum input variable
$Y_{norm}$	denormalized value of the output variable	g·cm−2 or cm
$Y_i$	real value of the output variable	g·cm−2 or cm
$Y_{min}$	minimum output variable	g·cm−2 or cm
$Y_{max}$	maximum output variable	g·cm−2 or cm
$f\left(ne{t}_i\right)$	transfer function of neural network
xi	input from
$w_{ij}$	weight of the connection between unit i and unit j
bi	bias
$e_i$	Backpropagation (error)
δ	summation index that enforces j > i
e	products errors
ε	injected errors
$\widehat{y_i}$	predicted values	g·cm−2 or cm
$\overline{\widehat{y_i}}$	actual values	g·cm−2 or cm
$\overline{y_i}$	mean of the observed values	g·cm−2 or cm
$y_i$	mean of predicted values	g·cm−2 or cm
n	number of data points

respectively.

2.2. Sampling and Measurements

Destructive sampling was undertaken immediately after UAV remote sensing data collection. Wheat plants were harvested from 0.36 m² quadrats in all experiments. The flights and the ground samplings occurred on 7 March, 4 April, and 28 May 2020, reflecting the three stages of growth: the early stem elongation, late stem elongation, and anthesis growth stages. The leaves and stems of harvested wheat plants were separated and weighed, followed by oven drying for 24 h at 105 °C before the dry weight was determined. The canopy equivalent water thickness ($EW{T}_{Canopy}$_y) corresponded to the hypothetical thickness of a single layer of water averaged over the whole ground area (Ag). Wocher et al. [21] defined $EW{T}_{Canopy}$ as the sum of the EWT of the leaf, stem, and fruit over one square meter of ground. To derive an unbiased estimation of EWT, we considered the weight of the fresh and dry masses of the leaf, stem, and ear for the calculation of $EW{T}_{Canopy}$ over one square meter of ground.

$$TotalEW{T}_{Canopy}=\sum \left(F{W}_{leaf+stem+ears}-D{W}_{leaf+stem+ears}\right)\ast A{g}^{-1}\left[{\mathrm{g}\mathrm{cm}}^2\right]or\left[\mathrm{cm}\right]$$

(1)

where Ag denotes the ground area, FW is the fresh sample weight, and DW is the oven dry weight. The $EW{T}_{leaf}$, $EW{T}_{stem}$, $EW{T}_{ear}$, and $EW{T}_{canopy}$ per cm² were calculated from specific water content per ground area.

2.3. UAV Data Collection and Processing

The multirotor UAV Spreading Wings S900 (DJI-Innovations Inc., Shenzhen, China) with six rotors, GPS, and flight control stabilizers was used in this study. To ensure maximum overlap in the fields, a flight path was set before the flights. The multispectral camera, Micasense Red Edge-MX (MicaSense, Seattle, WA, USA) (https://micasense.com/rededge-mx/, accessed on 3 November 2021) mounted on the UAV has five bands in the VIS-NIR spectral range (Red, Green, Blue, Nir, and Red Edge). The details of the multispectral sensor are shown in

Table 3. Specifications of sensors used in the present study.

Band	Bandwidth	Wavelength	Picture Resolution
Blue	20	475	1280 × 960
Green	20	560	1280 × 960
NIR	40	840	1280 × 960
Red	10	668	1280 × 960
Red Edge	717	10	1280 × 960

. Images were acquired across the entire field at the nadir viewing angle. The data were collected at 11 am–2 pm to minimize the shading effect of the canopy, at a speed of 2 ms⁻¹ and an altitude of 30 m (

Table 4. Flight details for the automated unmanned aerial vehicle imagery system during wheat growing season 2020.

Period	Flight Altitude (m)	Speed (ms−1)	Snapshot Interval (s)	Growth Stages
7 March 2020	30 m	2.5	2.5	Early stem elongation
4 April 2020	30 m	2.5	2.5	Late stem elongation
28 May 2020	30 m	2.5	2.5	Anthesis

). Radiometric calibration was performed using images of a spectral white panel before flight. The calibration image was used to calibrate each band during image processing in the Pix4D mapper (PIX4d, Lausanne, Switzerland) (https://www.pix4d.com/, accessed on 3 November 2021). The following steps were used for UAV image processing after UAV data collection. The pix4D mapper was used to process all the images into one large image using calibration images. The individual images were combined to form a large ortho-mosaic. The ortho-mosaic images were then radiometrically corrected and converted into reflectance. The resulting large ortho-mosaic images were Blue, Green, Near infrared (NIR), Red, and Red Edge (RE). Large ortho-mosaic images or their converted reflectance values were used to calculate the different VIs used in this study. The UAV data acquisition, processing methodology, image segmentation, and VIs calculation are described in Figure 2.

An automatic image segmentation algorithm called the Otsu segmentation algorithm [49] was implemented to minimize the effects of soil background and other non-leaf materials on canopy information pixels. The normalized difference vegetation index (NDVI) was first calculated and then used with the Otsu algorithm to separate soil from the canopy. The NDVI was employed for separating vegetation from other materials (soil included) [50]. The result of this was used in each plot in the form of a threshold value representing the border between soil and canopy. Soil or other background pixels were defined as 0 and wheat pixels as 1 during the extraction of the NDVI map. The integration of NDVI with the Otsu algorithm yielded efficient results in separating canopy from soil. The entire process of the NDVI–Otsu method was implemented using MATLAB and ArcMap.

2.4. $EW{T}_{canopy}$ Regression Model Development

2.4.1. Machine Learning for Regression

In this subsection, we present the ML regression method used to predict $EW{T}_{canopy}$ from the VIs. Figure 3 illustrates a typical workflow of a supervised ML algorithm for regression in MATLAB. In the first step, the FS method based on DT is used to reduce the input data to only five relevant input parameters. The next step is to normalize the data using the minimum–maximum normalization techniques. Seventy (70) percent of the normalized data were used as a training dataset to model EWT, while 15% were used for validation, and the remaining 15% were used for testing all the machine learning models used in this study. Next, the VIs and measured $EW{T}_{canopy}$ values pass into the training phase, where machine learning algorithms are used to identify a good model that can map the inputs to desired outputs. The validation and testing phases provide feedback to the ML phase so as to improve model accuracy. The training process is repeated until the desired accuracy level achieved. Once a model is constructed, it is used to predict $EW{T}_{canopy}$ from the new VIs data. FW and DW are estimated separately using vegetation indices and SVMs, BRT, ANN-MLP, and DNN-MLP algorithms. The resulting predicted FW and DW are used to estimate $EW{T}_{canopy}$. MATLAB MathWorks was used to simulate all the machine learning models used in this study. Constructing the $EW{T}_{canopy}$ maps involved two steps. The first was to convert the calculated VIs map (G, MTVI2, RE, OSAVI, NIR) to an “n × m” matrix with n rows and m columns. In order to be fed into the ML model, the matrix has to be presented as “n × 1”, with n rows and one column. In our study, we used the reshape function in Matlab to reshape the matrix. The output data are the predicted $EW{T}_{canopy}$ values, given in the form of an “n × 1” matrix. The second step is to reshape the “n × 1” to an “n × m” matrix and convert the matrix to a MATLAB image using the inwrite function.

2.4.2. Multiple Linear Regression

Multiple regression is generally used to explain the relationships between multiple independent/inputs variables and one dependent/target variable. The general form of the multiple regression equation is:

$$y=a+\beta { }_{1}X{ }_1+\beta { }_{2}X{ }_2+\cdots +\beta { }_{n}X{ }_n+\epsilon$$

(2)

where y is the dependent variable; a is the intercept; and the parameters $\beta { }_1$, $\beta { }_2$,…, $\beta { }_n$ are the regression coefficients associated with $X{ }_1$, $X{ }_2$,…, $X{ }_n$, respectively, while ε is the regression residual reflecting the difference between the observed and fitted linear relationships. The independent variables are VIs, represented as $X{ }_1$, $X{ }_2$,…, $X{ }_n$.

2.4.3. Support Vector Machine

The support vector machine (SVM), introduced by Boser et al. [51], is one of the most commonly applied supervised learning methods for regression as well as classification problems. When used for regression problems, the SVM model is known as support vector regression (SVR), which is used to predict a target using input variables. SVR is a powerful algorithm with the flexibility to be tolerant of an error margin (ϵ), and through tuning the tolerance, it can be made to fall outside the acceptable error rate. The kernel function is important for SVM analysis. Fan et al. [52] used SVM to estimate daily maize transpiration, and they found that SVM, alongside decision tree and deep learning models, can successfully estimate daily maize transpiration. Durbha et al. [53] retrieved the leaf area index using a multiangle imaging spectroradiometer.

In the present study, the SVM analyses were performed using different kernel functions, such as the Gaussian radial basis function (RBF) and polynomial kernel functions.

2.4.4. Boosted Regression Tree

Regression tree methods are commonly used to construct a model that can predict and explain target data from input data. The regression tree can capture the effect of each input variable on the target variable. However, the gradient boosted method combines several simple models to improve the prediction performance of one single model [54]. For instance Zhang et al. [55] modeled upland rice yield responses to climate factors using boosted regression tree in Sahel. A boosted regression tree (BRT) was employed to model the $EW{T}_{canopy}$ using spectral indices in this study. BRT incorporates the strengths of regression and boosted algorithms [55]. BRT is commonly used because it does not require data transformation or outlier elimination.

2.4.5. Artificial Neural Network Regression Model

The artificial neural network model is inspired by the structure of neural networks in the brain. Zha et al. [39] evaluated ANN other ML methods for estimating rice (Oryza sativa L.) aboveground biomass, plant N uptake, and N nutrition index. The ANN-MLP structure used in this study is a three-layer learning network consisting of an input layer, a hidden layer, and an output layer. The model minimizes the error based on the mean square error minimum value during the training process, using a tangent and sigmoid transfer function [55]. A maximum epoch of 10,000 iterations was set during the training, and we employed gradient descent with momentum and a learning rate of 0.2. A Bayesian regularized neural network model with a Levenberg–Marquart (LM) backpropagation algorithm was used during the training process to improve the generalization of the model. Sigmoid and logistic sigmoid functions were used for the activation function in each neuron, while a linear transfer function was used to calculate the network output. Five spectral indices were used as the inputs in this study, while $EW{T}_{canopy}$ was regarded as a target. The data were divided into three subsets: training, validating, and testing. The “dividerand” function of Matlab was used to randomly divide the data.

$$\mathrm{Forward}\mathrm{equation}\left(\mathrm{Transfer}\right){\mathrm{y}}_{\mathrm{i}}=f\left(ne{t}_i\right)=f\left({\displaystyle \sum }_j{w}_{ij}{x}_j+b_i\right)$$

(3)

where $f\left(ne{t}_i\right)$ is the transfer function, with a transfer threshold defined by [0, 1] for the sigmoid logistic and [−1, 1] for the sigmoid tangent transfer; xi is the input from i; $w_{ij}$ is the weight of the connection between unit i and unit j; bi is the bias.

$$\mathrm{Backpropagation}\mathrm{equation}\left(\mathrm{error}\right)e_i=-{\epsilon }_i+{\displaystyle \sum }_{j>1}{w}_{ij}{\delta }_j$$

(4)

where δ is the summation index that enforces j > i, and e and ε are the products and injected errors. The error is propagated from the output layer to the input layer in order to update the weights of connections using a gradient equation.

2.4.6. DNN-MLP Model Deployment

Deep learning algorithms were used for image classification [56] and regression analyses [34,52,57]. The addition of more data during deep learning improves the performance of the model, and this makes this technique superior to other learning techniques, such as ANNs, SVM, and RF, which reach a plateau in performance after a certain quantity of data is fed into the model. This study employed a DNN-MLP model with a Relu transfer function to stimulate $EW{T}_{canopy}$ using spectral indices by setting a maximum epoch at 10000 interactions. Adaptive moment estimation (Adam) was used as an optimization algorithm in the DNN-MLP model. The DNN-MLP structure is a three-layered feed-forward neural network consisting of an input layer, a hidden layer, and an output layer (Figure 4). The “dividerand” function in MATLAB was also used to randomly divide the data.

2.5. Data Pre-Processing Techniques

2.5.1. Data Normalization

The input variables were scaled to the same range (0, 1) using the minimum–maximum normalization techniques. Data normalization is carried out as part of data preparation for ML. The goal of normalization is to change the values of the inputs and output variables to a common scale [58]. $X_{min}$ and $X_{max}$ are the minimum and maximum values of the $ith$ attribute involved in the normalization process. Each input variable (VIs, $EW{T}_{canopy}$, FW, and DW) was normalized using the following Equation (5):

$$X_{norm}=(X_i-X_{min})/(X_{max}-X_{min})$$

(5)

where $X_{norm}$, $X_i$,$X_{min}$, and $X_{max}$ represent the normalized value, the real value of the input variable, the minimum input variable, and the maximum input variable, respectively. The real predicted value was denormalized after training, validation, and testing, according to the following Equation (6):

$$Y_i=Y_{min}+Y_{norm}\left(Y_{max}-Y_{min}\right)$$

(6)

where $Y_{norm}$, $Y_i$,$Y_{min}$, and $Y_{max}$ represent the normalized value, the real value of the output variable, the minimum output variable, and the maximum output value, respectively.

2.5.2. Feature Selection

Feature selection is a pre-processing technique used in ML to reduce the under- and over-fitting problems [46]. It is used to remove the irrelevant input variables and thus improve the learning accuracy of ML algorithms [58]. Several feature selection methods, such as functional discriminate analysis, principal component analysis (PCA), and sensitivity analysis, have been implemented. Feature selection based on a DT was used in this study to score the importance of each input variable in the model. It is easy to determine the contribution of each feature to the regression and its relative significance based on whether a leaf node, i.e., the output (DW, FW, and $EW{T}_{canopy}$), is higher or lower in the tree using DT. The five best VIs were used to model DW, FW, and $EW{T}_{canopy}$.

2.5.3. Model Performance

The determination coefficient (R²), Nash–Sutcliffe efficiency (NSE) [59], root mean square error (RMSE), and mean absolute error (MAE) were used to minimize the error and to assess the predictive accuracy of the regression models used in this study. The statistical indices are presented as follows:

$${\mathrm{R}}^2=\frac{{\left[{{\displaystyle \sum }}_{i=1}^n\left(y_i-\overline{y_i}\right)\left(\widehat{y_i}-\overline{\widehat{y_i}}\right)\right]}^2}{\left[{{\displaystyle \sum }}_{i=1}^n{\left(y_i-\overline{y_i}\right)}^2{{\displaystyle \sum }}_{i=1}^n{\left(\widehat{y_i}-\overline{\widehat{y_i}}\right)}^2\right]}$$

(7)

$$\mathrm{NSE}=1-\left[{{\displaystyle \sum }}_{i=1}^n{\left(y_i-\widehat{y_i}\right)}^2/{{\displaystyle \sum }}_{i=1}^n{\left(y_i-\overline{y_i}\right)}^2\right]$$

(8)

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

(9)

$$\mathrm{MAE}=\frac{1}{n}{\displaystyle \sum }_{i=1}^n\left|\widehat{y_i}-y_i\right|$$

(10)

where $\widehat{y_i} \mathrm{and} y_i$ are the predicted and actual values; $\overline{y_i} \mathrm{and} \overline{\widehat{y_i}}$ are the mean of the observed and predicted values; and n is the number of data points. The larger the values of R² and NSE, and the smaller the values of RMSE and MAE, the greater the precision and accuracy of the model in predicting $EW{T}_{canopy}$.

3. Results

3.1. Dynamic Changes of FW, DW and $EW{T}_{canopy}$ during the Growth Stage

Fresh weight (FW) and dry weight (DW) are amongst the most important crop growth indices. The accurate estimation of wheat FW and DW at different crop growth stages leads to effective agricultural field management. The dynamic changes in FW, DW, and $EW{T}_{canopy}$ observed at different growth stages are shown in Figure 5.

Table 5. Statistics (range, mean, standard deviation) for in situ measured $EW{T}_{leaf}$, $EW{T}_{stem}$, $EW{T}_{ear}$ and $EW{T}_{canopy}$.

Exps.	Exp. (1)	Exp. (2)	Exp. (3)
$EW{T}_{leaf}$ Range [g cm−2]	[0.003–0.296]	[0.014–0.186]	[0.001–0.158]
Mean (std) (g cm−2)	0.084 (0.072)	0.077 (0.041)	0.056 (0.05)
$EW{T}_{stem}$ Range [g cm−2]	[0.013–0.175]	[0.041–0.382]	[0.004–0.259]
Mean (std) (g cm−2)	0.081 (0.037)	0.166 (0.094)	0.091 (0.062)
$EW{T}_{ear}$ Range [g cm−2]	[0.01–0.112]	[0.039–0.108]	[0.006–0.082]
Mean (std) (g cm−2)	0.055 (0.016)	0.077 (0.019)	0.039 (0.022)
$EW{T}_{canopy}$ Range [g cm−2]	[0.03–0.442]	[0.09–0.567]	[0.011–0.459]
Mean (std) (g cm−2)	0.184 (0.089)	0.268 (0.126)	0.173 (0.113)

also summarizes the statistics (range, mean, standard deviation) used for in situ measured $EW{T}_{leaf}$, $EW{T}_{stem}$, $EW{T}_{ear}$, and $EW{T}_{canopy}$. From the early stem elongation stage to the late stem elongation stage, there was an increase in FW, DW, and $EW{T}_{canopy}$. However, from the late stem elongation to anthesis growth stages, there was a decrease in $EW{T}_{canopy}$ that can be explained by an increase in DW in favor of FW. The mean value of FW did not change significantly from the late elongation growth stage to the anthesis growth stage. With regard to FW and DW, we observed a gradual increase from the late stem elongation to the anthesis growth stage. Our results show that wheat growth was more vigorous in the late stem elongation growth stage compared to other stages. The mean values of FW at the early stem elongation, late stem elongation, and anthesis growth stages were 14.89, 28.71, and 29.05 t ha⁻¹, respectively, while the mean DW at the same growth stages was 2.67, 6.49, and 10.40 t ha⁻¹, respectively. The FW and DW increased by 92% and 143% from the early to late stem elongation growth stages, respectively (Figure 5). Conversely, FW and DW increased by 1.2% and 60% from the late stem elongation to the anthesis growth stage, respectively.

3.2. Variable Inputs’ Effects on FW, DW, and $EW{T}_{canopy}$ Estimation

Decision tree (DT) methods were implemented in this study to reduce the VIs to the most useful input variables only. The selected input variables were used in different ML models to predict the target variable ($EW{T}_{canopy}$, FW, and DW). It is very important to reduce the number of input variables in order to improve the performance of the model. The importance of VIs for estimating $EW{T}_{canopy}$, FW, and DW is shown in Figure 6. In this study, the sixty-seven individual calculated VIs were fed into the DT algorithm to select the five most relevant input variables on the basis of their feature weights (Appendix A 

Table A1. The vegetation indices evaluated in this study, B, G, R, RE, and NIR, indicate blue, green, red, red edge, and near-infrared band reflectance.

Vegetations (VIs)	Formulas	References
Blue Normalized Difference Vegetation Index (BNDV)	(NIR − B)/(NIR + B)	[72]
Green Chlorophyll Index (CIg)	NIR/G − 1	[73]
Red Edge Chlorophyll Index (CIre)	NIR/RE − 1	[25]
DATT Index (DATT)	(NIR − RE)/(NIR + R)	[24]
Excess Blue Vegetation index (ExB)	(1.4 × B − G)/(G + R + B)	[74]
Excess Green minus Excess Red (EXGR)	ExR − ExG	[75]
Excess Green index (ExG)	(2 × G − R − B)	[76]
Excess Red Vegetation index (ExR)	(1.4 × R − G)/(G + R + B)	[77]
Green Difference Vegetation Index (GDVI)	NIR − G	[78]
Green Leaf Index (GLI)	(2×G–R–B)/(– R − B)	[79]
Green Normalized Difference Vegetation Index (GNDVI)	(NIR − G)/(NIR + G)	[80]
Green Optimal Soil Adjusted Vegetation Index (GOSAVI)	(1 + 0.16) (NIR − G)/(NIR + G + 0.16)	[81]
Green Re–normalized Different Vegetation Index (GRDVI)	(NIR − G)/SQRT (NIR + G)	[26]
Green Ratio Vegetation Index (GRVI)	(G − R)/(G + R)	[82]
Green Red Vegetation Index (GRVI_Ratio)	NIR/G	[82]
Green Soil Adjusted Vegetation Index (GSAVI)	1.5× ((NIR − G)/(NIR + G + 0.5))	[81]
Green Wide Dynamic Range Vegetation Index (GWDRVI)	(0.12×NIR − G)/(0.12 × NIR + G)	[83]
Kawashima Index (IKAW)	(R − B)/(R + B)	[84]
Modified Chlorophyll Absorption in Reflectance Index 1 (MCARI1)	((NIR − RE) − 0.2 × (NIR − G)) × (NIR/RE)	[85]
Modified Chlorophyll Absorption in Reflectance Index 2 (MCARI2)	1.5 × (2.5 × (NIR − RE) –1.3× (NIR − G))/SQRT (SQ (2 × NIR + 1)) − (6 × NIR − 5 × SQRT(RE) − 0.5)	[86]
Modified Chlorophyll Absorption in Reflectance Index 3 (MCARI3)	((NIR − RE) − 0.2 × (NIR − R))/(NIR/RE)	[26]
Modified Chlorophyll Absorption in Reflectance Index 4 (MCARI4)	1.5 × (2.5 × (NIR − G) –1.3 × (NIR − RE))/SQRT (SQ (2 × NIR + 1)) − (6 × NIR − 5 × SQRT(G) − 0.5)	[26]
Modified Double Difference Index Green (MDD)	(NIR − RE) − (RE − G)	[87]
Modified Double Difference Index Red (MDD)	(NIR − RE) − (RE − R)	[27]
Modified Green Red Vegetation Index (MGRVI)	(SQ(G) − SQ (R))/(SQ(G) + SQ (R))	[88]
Modified Nonlinear Index (MNLI)	1.5 × (SQ(NIR) − R)/(SQ(NIR) + R + 0.5)	[89]
Modified Red Edge Difference Vegetation Index (MREDVI)	RE − R	[26]
Modified Red Edge Soil Adjusted Vegetation Index (MRESAVI)	0.5 × (2 × NIR + 1 − SQRT (SQ (2 × NIR + 1)) − 8 × (NIR − RE))	[90]
Modified Red Edge Transformed Vegetation Index (MRETVI)	1.2 × (1.2 × (NIR − R) − 2.5 × (RE − R))	[86]
Modified Soil Adjusted Vegetation Index (MSAVI)	0.5× (2×NIR + 1 − SQRT (SQ (2×NIR + 1) –8× (NIR − G)))	[90]
Modified Simple Ratio (MSR)	(NIR/R − 1)/SQRT (NIR/R + 1)	[91]
Modified Green Simple Ratio (MSR_G)	(NIR/G − 1)/SQRT (NIR/G + 1)	[91]
Modified Red Edge Simple Ratio (MSR_RE)	((NIR/RE) − 1)/SQRT ((NIR/RE) − 1)	[91]
Modified Transformed Chlorophyll Absorption in Reflectance Index (MTCARI)	3× ((NIR − RE) − 0.2 × (NIR − R) × (NIR/RE))	[92]
Modified Red Edge Soil Adjusted Vegetation Index (MRESAVI)	0.5 × (2 × NIR + 1 − SQRT (SQ (2 × NIR + 1)) − 8 × (NIR − RE))	[90]
Modified Red Edge Transformed Vegetation Index (MRETVI)	1.2 × (1.2 × (NIR − R) − 2.5 × (RE − R))	[86]
Modified Soil Adjusted Vegetation Index (MSAVI)	0.5× (2×NIR + 1 − SQRT (SQ (2×NIR + 1) –8× (NIR − G)))	[90]
Modified Simple Ratio (MSR)	(NIR/R − 1)/SQRT (NIR/R + 1)	[91]
Modified Green Simple Ratio (MSR_G)	(NIR/G − 1)/SQRT (NIR/G + 1)	[91]
Modified Red Edge Simple Ratio (MSR_RE)	((NIR/RE) − 1)/SQRT ((NIR/RE) − 1)	[91]
Modified Transformed Chlorophyll Absorption in Reflectance Index (MTCARI)	3× ((NIR − RE) − 0.2 × (NIR − R) × (NIR/RE))	[92]
Modified Triangular Vegetation Index (MTVI2)	1.5 × (1.2 × (NIR − G) − (2.5 × R–G))/SQRT (SQ (2 × NIR + 1) − (6 × NIR − 5 × SQRT(R)) − 0.5)	[86]
Normalized Difference Red Edge (NDRE)	(NIR − RE)/(NIR + RE)	[93]
Normalized Difference Vegetation Index (NDVI)	(NIR − R)/(NIR + R)	[94]
Normalized Green Index (NGI)	G/(NIR + RE + G)	[81]
Nonlinear Index (NLI)	(SQ(NIR) − R)/(SQ(NIR) + R)	[95]
Normalized NIR Index (NNIR)	NIR/(NIR + RE + G)	[81]
Normalized Near Infrared Index (NNIRI)	NIR/(NIR + RE + R)	[27]
Normalized Red Edge Index (NREI)	RE/(NIR + RE + G)	[81]
Normalized Red Edge Index (NREI)	RE/(NIR + RE + R)	[27]
Normalized Red Index (NRI)	R/(NIR + RE + R)	[27]
Optimized SAVI (OSAVI)	(1 + 0.16) × (NIR − R)/(NIR + R + 0.16)	[96]
Renormalized Difference Vegetation Index (RDVI)	(NIR − R)/SQRT (NIR + R)	[97]
Red Edge Difference Vegetation Index (REDVI)	NIR − RE	[26]
Red Edge Normalized Difference Vegetation Index (RENDVI)	(RE − R)/(RE + R)	[98]
Red Edge Optimal Soil Adjusted Vegetation Index (REOSAVI)	(1 + 0.16) × (NIR − RE)/(NIR + RE + 0.16)	[96]
Red Edge Renormalized Different Vegetation Index (RERDVI)	(NIR − RE)/SQRT (NIR + RE)	[26]
Red Edge Ratio Vegetation Index (RERVI)	NIR/RE	[99]
Red Edge Soil Adjusted Vegetation Index (RESAVI)	1.5× ((NIR − RE)/(NIR + RE + 0.5))	[81]
Red Edge Simple Ratio (RESR)	RE/R	[100]
Red Edge Transformed Vegetation Index (RETVI)	0.5× (120× (NIR − R) – 200 × (RE − R))	[91]
Optimized Red Edge Vegetation Index (REVIopt)	100 × (Ln (NIR) − Ln (RE))	[101]
Red Edge Wide Dynamic Range Vegetation Index (REWDRVI)	(0.12×NIR − RE)/(0.12 × NIR + RE)	[83]
Red Green Blue Vegetation Index (RGBVI)	(SQ(G) − (B × R))/(SQ(G) + (B × R))	[88]
Ratio Vegetation Index (RVI)	NIR/R	[102]
Soil–Adjusted Vegetation Index (SAVI)	1.5× (NIR − R)/(NIR + R + 0.5)	[103]
Transformed Normalized Vegetation Index (TNDVI)	SQRT ((NIR − R)/(NIR + R) + 0.5)	[104]
Optimal Vegetation Index (VIopt)	1.45×(SQ(NIR) + 1)/(R + 0.45)	[101]
Wide Dynamic Range Vegetation Index (WDRVI)	(0.12×NIR − R)/(0.12 × NIR + R)	[83]

and

Table A2. Mean, standard deviation, minimum, and maximum values of vegetation indices used in the study.

Vegetation Indices/Bands	Mean ± SD	Min.	Max.	Vegetation Indices/Bands	Mean ± SD	Min.	Max.
B	0.03 ± 0.02	0.08	0.01	MSR_RE	1.25 ± 0.55	2.12	0.51
BNDV	0.82 ± 0.11	0.95	0.62	MTCARI	0.04 ± 0.16	0.23	–0.39
Cig	6.48 ± 4.33	14.99	1.59	MTVI2	0.34 ± 0.35	0.84	–0.17
CIre	1.88 ± 1.41	4.51	0.26	NDRE	0.41 ± 0.2	0.69	0.11
DATT	0.48 ± 0.25	0.8	0.12	NDVI	0.64 ± 0.29	0.95	0.2
ExB	–0.13 ± 0.07	–0.05	–0.26	NGI	0.11 ± 0.04	0.18	0.05
ExG	0.01 ± 0.04	0.06	–0.07	NIR	0.38 ± 0.08	0.53	0.23
ExGR	0.09 ± 0.31	0.55	–0.33	NLI	0.31 ± 0.5	0.91	–0.46
ExR	0.1 ± 0.28	0.49	–0.29	NNIR	0.63 ± 0.12	0.8	0.46
G	0.07 ± 0.03	0.14	0.03	NNIRI	0.62 ± 0.15	0.83	0.41
GDVI	0.31 ± 0.1	0.48	0.14	NREI_G	0.26 ± 0.08	0.37	0.15
GLI	–0.57 ± 0.62	0.21	–1.68	NREI_R	0.25 ± 0.06	0.32	0.15
GNDVI	0.69 ± 0.15	0.88	0.44	NRI	0.13 ± 0.1	0.27	0.02
GOSAVI	0.46 ± 0.12	0.65	0.25	OSAVI	0.55 ± 0.25	0.85	0.18
GRDVI	0.59 ± 0.14	0.79	0.34	R	0.08 ± 0.07	0.27	0.01
GRVI	0.07 ± 0.28	0.47	–0.32	RDVI	0.43 ± 0.2	0.69	0.14
GRVI_Ratio	7.48 ± 4.33	15.99	2.59	RE	0.16 ± 0.05	0.33	0.09
GSAVI	0.49 ± 0.13	0.69	0.26	REDVI	0.23 ± 0.12	0.42	0.06
GWDRVI	–0.14 ± 0.28	0.31	–0.53	RENDVI	0.44 ± 0.26	0.78	0.09
IKAW	0.25 ± 0.21	0.54	–0.06	REOSAVI	0.36 ± 0.19	0.63	0.1
MCARI1	0.61 ± 0.55	1.77	0.03	RERDVI	0.41 ± 0.32	1.36	0.09
MCARI2	0.29 ± 0.6	1.16	–0.75	RERVI	2.77 ± 1.57	5.51	0.09
MCARI3	0.05 ± 0.01	0.08	0.03	RESAVI	0.42 ± 0.29	1.23	0.09
MCARI4	–0.1 ± 0.64	0.89	–1.21	RESR	3.63 ± 2.27	8.1	1.19
MDD_G	0.13 ± 0.14	0.36	–0.09	RETVI	10.35 ± 6.55	21.96	2.01
MDD_R	0.15 ± 0.1	0.34	0.02	REVIopt	92.05 ± 51.11	170.3	22.94
MGRVI	0.12 ± 0.51	0.77	–0.58	REWDRVI	–0.51 ± 0.18	–0.21	–0.74
MNLI	0.14 ± 0.24	0.49	–0.2	RGBVI	0.36 ± 0.3	0.75	–0.06
MREDVI	0.07 ± 0.02	0.12	0.03	RVI	13.69 ± 13.08	44.02	1.5
MRESAVI	–0.89 ± 0.48	–0.22	–1.69	SAVI	0.46 ± 0.21	0.73	0.15
MRETVI	0.2 ± 0.15	0.48	0.03	TNDVI	1.06 ± 0.14	1.21	0.84
MSAVI	0.49 ± 0.16	0.76	0.23	VIopt	3.21 ± 0.53	3.98	2.39
MSR	2.65 ± 2.01	6.38	0.32	WDRVI	–0.06 ± 0.5	0.67	–0.7
MSR_G	2.04 ± 0.92	3.63	0.84

). The FS method was separately applied to the input variables, VIs, and the target variables $EW{T}_{canopy}$, FW, and DW. The five VIs chosen as input variables by the algorithm, according to their scores, were G, MTVI2, RE, OSAVI, and NIR, and their relative importance to the $EW{T}_{canopy}$ varied from 2.4 to 78%. The top five VIs selected by the DT algorithm while modeling FW were G, RE, IKAW, R, and RESR, and their relative importance ranged between 1.48 and 58.54%. MTVI2, RE, RESAVI, R, and NNIR were selected while modeling DW, and the relative importance varied from 2.3 to 80.34%. Red edge (RE) was consistently selected among the top five VIs for estimating $EW{T}_{canopy}$, FW, and DW. MTVI2 was among the top five VIs for $EW{T}_{canopy}$ and DW. The R band was also among the top five VIs used for FW and DW.

3.3. $EW{T}_{canopy}$ Responses from the Multiple Regression Model

Multiple linear regression (MLR) has been extensively used to determine the relationships between target and input variables. VIs and EWT from the three experiments were used to develop a linear model in this study. Figure 7a shows the scatter plot between $EW{T}_{canopy}$ and the input VIs using the MLR model. The relationship between $EW{T}_{canopy}$ and the VIs shows acceptable prediction performance. As such, VIs could be used to predict $EW{T}_{canopy}$ by MLR. The R², NSE, RMSE, and MAE values of the MLR were 0.843, 0.843, 0.0433 g/cm², and 0.0313 g/cm², respectively. Nevertheless, MLR demonstrated poor performance for predicting $EW{T}_{canopy}$ as compared to the other ML models (

Table 6. Summary statistics of multiple linear regression models.

Statistics	Values
Multiple R	0.9185
R Square	0.8436
Adjusted R Square	0.8406
Standard Error	0.0439
Intercept	0.4312
Beta:
$\mathrm{G}(\beta { }_1$)	−0.5454
$\mathrm{MTVI}2(\beta { }_2$)	1.4150
$\mathrm{RE}(\beta { }_3$)	0.8637
$\mathrm{OSAVI}(\beta { }_4$)	−1.3980
$\mathrm{NIR}(\beta { }_5$)	−0.3555

and

Table 7. Statistical performance indices of different machine learning (ML) algorithm during the training, validation, and testing stages for $EW{T}_{canopy}$, FW, and DW.

Model	Variables	Training				Cross-Validation				Testing
		R2	ENS	RMSE	MAE	R2	ENS	RMSE	MAE	R2	ENS	RMSE	MAE
ANN-MLP	$EW{T}_{canopy}$(%)	0.916	0.916	0.032	0.021	0.922	0.915	0.034	0.025	0.905	0.894	0.0334	0.019
	FW (t/ha)	0.905	0.904	3.889	2.583	0.909	0.907	3.944	2.395	0.954	0.952	2.629	1.976
	DW (t/ha)	0.868	0.868	0.680	0.508	0.875	0.864	0.574	0.449	0.924	0.922	0.512	0.391
DNN-MLP	$EW{T}_{canopy}$(%)	0.938	0.937	0.027	0.015	0.933	0.930	0.030	0.021	0.913	0.909	0.034	0.022
	FW (t/ha)	0.934	0.934	3.215	1.953	0.914	0.897	4.085	2.413	0.903	0.902	3.943	2.659
	DW (t/ha)	0.900	0.900	0.571	0.413	0.894	0.893	0.508	0.421	0.882	0.881	0.701	0.531
BRT	$EW{T}_{canopy}$(%)	0.948	0.947	0.026	0.017	0.893	0.868	0.035	0.025	0.872	0.868	0.039	0.027
	FW (t/ha)	0.928	0.928	3.377	2.175	0.885	0.883	4.204	2.922	0.902	0.900	3.972	2.738
	DW (t/ha)	0.917	0.917	0.537	0.423	0.778	0.758	0.854	0.690	0.814	0.803	0.759	0.531
SVM-Gaussian	$EW{T}_{canopy}$(%)	0.955	0.950	0.024	0.015	0.908	0.904	0.032	0.026	0.915	0.907	0.035	0.025
	FW (t/ha)	0.937	0.936	3.192	1.873	0.880	0.878	4.299	3.041	0.925	0.898	4.007	2.937
	DW (t/ha)	0.922	0.920	0.516	0.352	0.864	0.864	0.626	0.502	0.924	0.922	0.519	0.434
SVM-Polynomial	$EW{T}_{canopy}$(%)	0.900	0.899	0.035	0.023	0.846	0.843	0.043	0.027	0.902	0.899	0.036	0.023
	FW (t/ha)	0.892	0.892	4.097	2.769	0.857	0.857	4.783	2.978	0.852	0.850	4.976	2.868
	DM (t/ha)	0.861	0.860	0.684	0.514	0.821	0.812	0.735	0.571	0.894	0.888	0.623	0.505

). The multiple regression model for $EW{T}_{canopy}$ (y) was expressed as follows:

y = 0.4311 − 0.5454 × G + 1.4150 × MTVI2 + 0.8637 × RE − 1.3980 × OSAVI − 0.3555

(11)

3.4. Modeling $EW{T}_{canopy}$ Using DNN-MLP, ANN-MLP, BRT, and SVM

To predict $EW{T}_{canopy}$, we employed regression analysis using the DNN-MLP model with the ReLu transfer function and Adam optimizers algorithm; the ANN-MLP model with sigmoid and sigmoid tangent functions; and BRT and SVM using Gaussian and polynomial kernel functions. In this study, 267 data samples were considered to train the model, of which 15% were used for validation and 15% were used for testing. The number of epochs was set to 10000. Here, we compare the results of DNN-MLP with the MLR, ANN-MLP, BRT, and SVMs models. Table 7 shows the statistical indices comparing the models used to calibrate the FW, DW, and $EW{T}_{canopy}$ to VIs. It is clear that the performances of all the ML methods (Table 7) were better than the MLR model (Table 6). The SVM-Gaussian regression method had an R² of 0.941, an NSE of 0.937, an RMSE of 0.0274 g/cm², and an MAE of 0.0181 g/cm². DNN-MLP had an R² of 0.934, an NSE of 0.933, an RMSE of 0.0283 g/cm², and an MAE of 0.0165 g/cm². The ANN-MLP-MPL had an R² of 0.914, NSE of 0.914, RMSE of 0.0321 g/cm², and MAE of 0.0211 g/cm². The BRT had an R² of 0.926, NSE of 0.926, RMSE of 0.0298 g/cm², and MAE of 0.0194 g/cm². The SVM-Polynomial had an R² of 0.892, NSE of 0.891, RMSE of 0.0362 g/cm², and MAE of 0.0231 g/cm².

3.5. Relationship between Measured and Predicted $EW{T}_{canopy}$

The scatter plots of the $EW{T}_{canopy}$ predicted by the models versus the observed $EW{T}_{canopy}$ are given in Figure 7. However, for the ML models, the ability to accurately predict $EW{T}_{canopy}$ depends on the typology of the ML algorithm. SVM-Gaussian (R² = 0.941) was found to be the best for predicting $EW{T}_{canopy}$ using VIs. The performance of DNN-MLP (R² = 0.934) was also better than BRT (R² = 0.926), ANN-MLP (R² = 0.914), and SVM-Polynomial (R² = 0.892). Of the ML networks, the SVM-Polynomial produced the weakest model. A cross-comparison of all the models employed showed that, in terms of performance, they ranked in the order SVN-Gaussian, DNN-MLP, BRT, ANN-MLP, SVM-Polynomial, and multiple regression (

Table 8. Comparative performance statistics of the machine learning (ML) models employed in the study.

Performance Rank	Model	R2	NSE	RMSE	MAE
1	SVM-Gaussian	0.942	0.937	0.027	0.018
2	DNN-MLP	0.934	0.933	0.028	0.017
3	BRT	0.926	0.926	0.030	0.019
4	ANN-MLP	0.914	0.914	0.032	0.021
5	SVM-Polynomial	0.892	0.891	0.036	0.023

). The scatter plots of the $EW{T}_{canopy}$ predicted by the models versus the observed $EW{T}_{canopy}$ are given in Figure 7. The performances of the SVM-Gaussian and DNN-MLP models were judged as satisfactory, as the slope of the regression line is close to 1. Figure 8 shows a curve between the predicted models and the observed $EW{T}_{canopy}$. From this figure, we can conclude that the results of the SVM-Polynomial and DNN-MLP are the closest to the observed $EW{T}_{canopy}$.

3.6. Calculating $EW{T}_{canopy}$ using Predicted DW and FW

The DW and FW predicted using the DNN-MLP and VIs were used to calculate $EW{T}_{canopy}$. The calibration, validation, and testing processes were carried out using the DW and FW collected during the three experiments. The $EW{T}_{canopy}$ calculation results show that the direct estimation of $EW{T}_{canopy}$ performs well compared to the indirect estimation using predicted DW and FW. The results of the indirect method show a decrease in the performance of R², RMSE, and MAE (Figure 9). However, the planned comparisons reveal that the assessment of $EW{T}_{canopy}$ using MLR and SVM-Polynomial models and FW and DW predicted from VIs performed slightly better than the $EW{T}_{canopy}$ estimated directly from the VIs. The results clearly demonstrate that VI_S are suitable for $EW{T}_{canopy}$ assessment. They proved to be sensitive to water stress, which is the major factor influencing $EW{T}_{canopy}$ retrieval. Additionally, VIs are suitable for assessing $EW{T}_{canopy}$ instead of using the predicted FW and DW.

3.7. Model Visualization

The results of the ML models and multi-linear regression were used to generate the $EW{T}_{canopy}$ map in Exp. (2). Although SVM-Gaussian achieved the best prediction for $EW{T}_{canopy}$, the $EW{T}_{canopy}$ map of Figure 10 shows that SVM-Gaussian cannot accurately assign $EW{T}_{canopy}$ because of the overfitting problem. Therefore, in terms of practically applying the ML in estimating $EW{T}_{canopy}$, DNN-MLP gave results close to the observed $EW{T}_{canopy}$ in regard to the study areas, and achieved the highest level of accuracy in terms of visual comparison. The results of this study (Table 8) suggest that DNN-MLP achieved the highest level of accuracy and is better used in the assessment of $EW{T}_{canopy}$ using VIs. We therefore found it possible to model $EW{T}_{canopy}$ using the DNN-MLP methodology, finding this to be a powerful computational tool to model $EW{T}_{canopy}$ using multispectral data. The results indicate that the model accurately identified the low and high regions of $EW{T}_{canopy}$. The results will allow us to identify water deficiencies in plants, and to take the appropriate action for irrigation. The $EW{T}_{canopy}$ map will allow agricultural decision-makers to remotely quantify plant water content ($EW{T}_{canopy}$) and address the variability in this value so as to improve input efficiency and irrigation management. The ANN-MPL, DNN-MLP, BRT, and SVM-Polynomial models were successfully used to quantify $EW{T}_{canopy}$ in the field, but the SVM-Gaussian model failed to estimate $EW{T}_{canopy}$ because of the overfitting problem. The map indicates that all the models except SVM-Gaussian will assist farmers in managing irrigation.

4. Discussion

4.1. Dynamic Changes in FW, FW and Ewtcanopy during the Growth Stage

In this study, we evaluated dynamic changes within wheat during the early and late stem elongation and anthesis growth stages. $EW{T}_{canopy}$ is dependent on FW and DW, and so has a higher value during the late stem elongation stage. These results reflect those of Jin et al. [34], who also found that the DW increased steadily with crop development. The changes in FW affect the value of $EW{T}_{canopy}$. At the anthesis stage, the wheat dries and holds a small amount of water. The increases in DW, FW, and EWT during late stem elongation are explained by the fact that wheat growth is vigorous during this period. However, in this study, during the anthesis growth stage, there was a decrease in EWT and an increase in DW, while the FW did not change significantly. These results agree with the study of Wocher et al. [21], showing that EWT decreased in the late growth stage.

4.2. Performance of the Machine Learning Models

The performances of the SVM-Gaussian and DNN-MLP models were remarkable in the training, cross-validation, and testing periods. The calibration, cross-validation, and testing processes were carried out using the observed $EW{T}_{canopy}$, FW, and DW values. Other studies employed the DNN- and ANN-MLP models to predict crops’ biophysical parameters, such as crop biomass and yield, finding that DNN was the best model for these purposes [34,45,60]. In addition to MLR, four different ML algorithms were applied to predict $EW{T}_{canopy}$ status indicators in this study. The ML model performed significantly better than the models based on MLR. Our results are consistent with previous studies [55], which used BRT and ANN to forecast upland rice yield under climate change conditions in the Sahel. The MLR models can only model a linear combination of predictors, while the ML models can also model nonlinear relationships. DNN-MLP is a regression with a great capacity for supporting many hidden layers [61,62]. This study’s results corroborated those of Jin et al. [34], who indicated that the DNN-MLP algorithm could be used to accurately estimate plant biomass from VIs. The SVM algorithm is based on statistical learning, and provided accurate results [63]. ANN-MLP regression is a nonparametric, nonlinear model that creates a neural network between inputs and target data. Fan et al. [52] confirmed the utility of SVM, extreme gradient boosting (XGBoost), ANN-MLP, and DNN-MLP models for estimating the daily temperature of maize in Northwest China. $EW{T}_{canopy}$ estimation using FW and DW was not as successful as that performed with the direct estimation approach. This may be because the indirect estimation approach led to the accumulation of errors.

4.3. Feature Selection Methods

In this study, 67 VIs were used in the decision algorithm, and only the 5 best-scoring (1.48 to 8.34) were selected for the estimation of $EW{T}_{canopy}$, FW, and DW. The RE band was one of the five most relatively important when applying the DT. MTVI2 was among the five top VIs when using the feature selection method for FW and DW, while the R band was among the five best VIs for $EW{T}_{canopy}$ and DW. Other VIs were also important for $EW{T}_{canopy}$ estimation. This is consistent with other studies that showed that different feature selection methods increase the performance of the model and decrease the computational time [48,58,64,65]. Haq et al. [64] stated that feature selection methods significantly increase model accuracy. In the feature selection method based on DT, the process of constructing the DT is the process of feature selection. The main advantages of the DT algorithm are its high classification and regression accuracy and its strong robustness [46,48,65]. Consistent with previous works [46,66], this study also found that feature selection improves the performance of the model by reducing the number of input variables. These results further support the idea that the feature vector plays a significant role in the performance of the ML model by including or excluding input variables. This is essential to the processes of training, testing, and validation.

4.4. Advantages and Limitations of Machine Learning

In this study, the multispectral data and VIs were obtained using a multi-rotor UAV remote sensing platform, while the $EW{T}_{canopy}$ distribution maps were created based on DNN-MLP, ANN-MLP, BRT, and SVM model prediction. The $EW{T}_{canopy}$ map can be used to guide farmers in the application of irrigation water. The use of an UAV remote sensing platform combined with the ML method could help overcome the limitations of airborne and satellite remote sensing platforms, and provide a reliable data source for $EW{T}_{canopy}$ assessment during the growth stage. In addition to the commonly used red, green, blue, red edge, and NIR bands, other spectral regions should be used to diagnose $EW{T}_{canopy}$, such as shortwave infrared (SWIR)-based indices or hyperspectral cameras. Other studies found that a combination of multispectral and thermal images via the ML method refines the estimation of plant chlorophyll concentration [67], improves vegetation monitoring [68], and can help predict water stress. In the future, thermal and SWIR- hyperspectral remote sensing data should be used to improve the performances of ML models.

Although the use of UAVs combined with ML enabled the prediction of $EW{T}_{canopy}$, there are several limitations that prevent their wider use. The multispectral camera was relatively expensive. However, these costs could be reduced by using low-cost RGB imaging instead of a multispectral camera. Sánchez-Sastre et al. [69] successfully used RGB VIs to estimate chlorophyll content in sugar beet leaves. The same results may be achieved by using RGB VIs to assess $EW{T}_{canopy}$. Hence, further work is needed to assess the applicability of RGB imaging for assessing $EW{T}_{canopy}$ using RGB VIs and ML models. In relation to the main limitations of this ML method, it should be noted that an average farmer may require training on how to operate the UAV platform and process the data using the ML methods, which may be costly. This fact may prohibit the adoption of UAV technologies for individual farmers with only small agricultural fields. This may affect the adoption of the UAV remote sensing technology reported in the literature [70]. We propose that scientists and experts should assist farmers in the field. Finally, crop biophysical parameters are known to be influenced by UAV flight height, as reported by Oniga et al. [71]. Further research should be undertaken to determine the influence of height in the retrieval of $EW{T}_{canopy}$. Another drawback of commercial UAVs technology is the short flight time, which ranges from 20 min to 1 h and can thus only cover a very restricted area with every flight. In addition, effective UAVs cannot be used on a very windy or rainy day, meaning flights must be postponed. In addition, feature selection is essential for optimizing the accuracy of the model and for enhancing model interpretability. However, in the FS method, all 67 VIs were blindly fed into the FS algorithm to select the 5 best VIs for $EW{T}_{canopy}$ modeling. Other FS methods need to be investigated in order to assess the effects of other VIs on $EW{T}_{canopy}$ assessment. RS scientist and/or progressive farmer with larger farms size are able to generate these valuable maps on demand using the methods prescribed in this study despite the complex nature of ML and DL models. Furthermore, the newly developed models might also assist the policy makers and agriculture extensionist for making the recommendations that can guide smallholder farmers to accomplish water management. However, further investigations with different regions and wheat cultivars are recommended to test the applicability of newly developed models for crop water diagnosis.

5. Conclusions

In our study, the ML algorithms DNN-MLP, ANN-MLP, BRT, SVM-Gaussian, SVM-Polynomial, and MLR were used to improve the estimation accuracy of $EW{T}_{canopy}$. The results show that only five VIs could be used to accurately estimate $EW{T}_{canopy}$ when using a feature selection algorithm. This study has also demonstrated the power of the DT algorithm used for automatic relevance determination when evaluating the more relevant input parameters in modeling $EW{T}_{canopy}$ using ML models. SVM-Gaussian gave the best performance, followed by DNN-MLP, BRT, ANN-MLP, and SVM-Polynomial. In terms of the $EW{T}_{canopy}$ map, the SVM-Gaussian model performed poorly when compared to the other ML models, making DNN-MLP the most suitable ML to assist farmers in irrigation. The linear model performed poorly in $EW{T}_{canopy}$ estimation. These findings contribute to the accurate estimation of $EW{T}_{canopy}$, and thus improve irrigation efficiency and grain yield. However, other advanced ML and thermal imagery models may be used to model $EW{T}_{canopy}$. Besides this, other feature selection methods may be implemented to improve the retrieval accuracy of $EW{T}_{canopy}$ using ML. More studies are needed to further improve these ML-based models by using thermal, SWIR, and hyperspectral images for irrigation and crop management.