Use of HSI Coupled with Machine Learning for the Identification of Water Stress in Pear Seedling Leaves

1. Material

QingzhenD3’ pear root-stocks with around 12 leaves per plant were grown in the National Agricultural and Forestry Science and Technology incubator seedling base in Zhucheng, China. Pear seedlings of similar growth were chosen. Pear seedlings were transported to greenhouses based at Qingdao Agricultural University. Two weeks after transplanting, all pear seedlings were treated with the complete nutrient solution and supplied with all the essential nutrients. A known amount of nutrient solution was provided to each plant using a trickle nozzle.

The temperature of the greenhouse was 25℃. The humidity of the greenhouse was around 95%. Each pot contained one seedling. The substrate was composed of peat and pastoral soil (1:1). The pastoral soil was procured from Qingdao Agricultural University and carefully sieved. The experiment took place between July and September 2021.

2. Experiment design

Treatment groups were established as the prediction models and simulate the differing water stress situations. There were three treatment groups: the excessive water treatment group, the drought group, and the control group-normal watering group. For each treatment, 30 pear seedlings were used for hyperspectral imaging collection. Specifically, there were 30 repetitions for each treatment group. A pre-experiment was conducted to conclude that the amount of water a pot of pear seedling needs was 10mL. The night before the experiment, the pear seedlings of the three treatment groups were watered thoroughly. Once the pear seedlings of the drought group were watered, no further watering was required during the experiment. The excessive water treatment group maintained water at a depth of 1.5cm, and the water was uniformly replenished at 6 pm. In addition, 10mL of water was added to pear seedlings in the regular watering group every day at 6 pm.

One leaf in the middle of each plant (usually the fourth leaf from the bottom) was employed for hyperspectral data collection (Figure 1). According to the research method requirement, the middle leaf is relatively stable and representative of the whole plant. The acquisition time was 0 days, one day, three days, five days, and seven days after water treatment before apparent symptoms such as wilting and yellowing appeared on pear leaves. Images were taken every other day to ensure the internal environment of the plant was recorded consistently and accurately. On day 0, the starting point of water treatment, only 30 hyperspectral images were collected (10 for each treatment). From day 0, day 1, day 3, day 5, day 7, 30 images were gathered for each treatment. The sample collection period was from 9:00 am to 11:00 am. Table 1 indicates the number of pictures collected. Three areas of interest were extracted from each hyperspectral cube, therefore, we have 765 areas of interest in total.

A workflow consisting of three phases was proposed, from data analysis to model development. In order to clearly demonstrate the experiment process, the workflow of this experiment is shown in Figure 1.

Figure 1. Workflow of this research.

Figure 1. Workflow of this research.

Phase1: region of interest procured from the acquired hyperspectral image (with leaf as an example); Phase 2: data cube acquired after performing dimension reduction (SPA key wavelengths selection); Phase 3: three types of information input into the neural network; Phase 4: machine learning analysis of the input

The proposed workflow has four phases: region-of-interest (ROI) extraction, feature selection, three types of information input into the neural network (whole data cube, main spectral, and CNN feature extraction of key wavelength images), and machine-learning analysis. Region of interest (ROI) selection is usually used in the study to remove the irrelevant target regions. Researchers usually use the application ENVI to extract pixel points in ROI or average spectra of hyperspectral images as pixel-level or object-level data for analysis and processing (Kang et al.,2022). In phase 1, the region-of-interest leaf area was chosen using ENVI. In phase 2, feature wavelengths were selected using a successive projections algorithm (SPA). In the third phase, we procured three types of information, whole data cubes, key wavelengths obtained using SPA, and CNN features of images of key wavelengths. In the fourth phase, three machine-learning models were employed and compared to model plant water stress detection. The final output provided results in the early detection of dry, normal and overwatering plants from the feature wavelengths selected.

2.1. Hyperspectral image acquisition

Resonon imager PIKA L was used for image acquisition (Figure 2). Resonon imaging spectrometers were line-scan imagers, meaning they collected data one line at a time. PIKA L covered the Visible and Near-Infrared (VNIR) spectral range of 400-1000nm. It had 281 spectral channels, with a spectral bandwidth of 2.1nm. The spectral resolution-FWHM was 3.3 nm, while the spatial pixels were 900. Its max frame rate (fps) was 249. The application SpectrononPro can collect data cubes with a Benchtop System. In addition, because the water condition of the live plants can be influenced by the outside environment, such as light, heat, etc., the water content of leaves is prone to change, which may affect the differentiation of leaves under different treatments. Therefore, the hyperspectral camera must be carefully calibrated as frequently as needed to obtain a valid data cube.

The lifting platform was placed on the moving stage. The vertical and horizontal directions were parallel to the corresponding sides of the moving stage. In the meantime, the height of the stage was adjusted so that the distance from the bottom edge of the light source to the whiteboard was about 35cm. The position of the stage was changed so that the inner side of the stage coincided with the longitudinal axis of the lens. And the horizontal axis of the lens was located at 1/3 of the whiteboard on the right. Next, the pear seedling was placed on the moving stage by adjusting the height and spatial orientation of the pear seedling. The target leaf was placed in the longitudinal centre of the whiteboard along the horizontal axis of the lens. It should be sure that the leaf was on the same horizontal plane as the whiteboard. During the entire experiment, the position of the lifting platform was kept fixed, to make sure of the clarity of the image.

The pear seedling leaf was kept on the black background of stiff paper and imaged to acquire the image. Image capture parameters used were 256 frames and 45 steps to get the best resolution image between 400 and 100nm. The exposure time used was 24ms^-1. The parameters were kept consistent for all the photos. The resulting hyperspectral images was a particular block of 900×500×300 reflectance image, representing a 3-D image with X-axis and Y-axis coordinate information and the other representing the spectral information at 300 different wavelengths. This information was stored for subsequent analysis.

2.2. Image acquisition and image correction

The image acquisition was conducted in a dark room to avoid undesired light and at a controlled temperature and humidity of 20℃ and 65%, respectively. Hyperspectral datacubes were acquired in the 400-1000nm spectral range with 5nm intervals between contiguous bands. The acquired raw images were corrected with two reference images using the following Equation (1).

$$\mathrm{R}=\frac{R_0-D}{W-D}$$

(1)

Where: R was the relative reflectance image of the sample, $R_0$ was the raw image of the sample, W was the white reference image acquired from a uniform, stable, and high reflectance ceramic tile (reflectance), and D was the current dark image acquired by completely covering the camera lens with its non-reflective opaque black cap (Wang et al., 2018).

2.3. Feature wavelength selection

Feature selection methods in HSI aim to reduce dimensionality while preserving relevant information for later classification (Wang et al., 2004). The selection of wavelength variables is integral to establishing a nondestructive testing model based on hyperspectral imaging. The model can be simplified by screening the characteristic wavelengths or wavelength ranges. Various feature selection techniques have been used to select important variables for spectroscopic data (Wang et al., 2018).

Feature extraction can reduce the dimension of spectral data and improve regression model performance. The successive projections algorithm (SPA) was used to obtain the optimum wavelength. The successive projection algorithm (SPA), a variable-selection technique, was proposed for constructing multivariate calibration models (Woldgiorgis et al., 2021). Later, it was extended to address classification problems (Kandpal et al., 2016). SPA is a one-way selection algorithm that mainly uses collinear minimisation to select the optimal variables (Arau ́jo et al., 2001). The advantage of this algorithm is that the variable group with the least redundant information can be chosen from more spectral variables. In addition, collinearity between variables within a variable group is minimised (Pontes et al., 2005). The SPA is executed in Matlab 2023 (The Math Works, Natick, USA) using a SPA toolbox (available at http://www.ele.ita.br/~kawakami/spa/). The test_size is 0.4, m_min is 2, m_max is 28.

According to Figure 3 and Figure 7 wavelengths were chosen, substituting the 300 wavelengths of the original sample. The optimal wavelengths selected were listed in Table 2. Next, the eigenvectors of the image corresponding to the characteristic wavelength were extracted and put into neural networks for analysis. For these seven wavelengths, we obtain seven images. After applying CNN, for each image, we have 4096 features. After obtaining all features relevant to all the data points, we put into neural networks for analysis.

Different data input is shown in Table 3. In the data input part, the height of each hyperspectral imaging data is 300, the total treatment samples is 30+90×4=390, after SPA major wavelength selection, 7 spectra were selected to represent the original 300 wavelengths. Since we use AlexNet for feature selection, we procured 4096 features for each HSI data cube. Different parameters of neural networks are shown in Table 4.

According to Table 4, for the RR-MLR, the regularization parameter (gam) is 10, for the Bayesian Linear Regression, the covariance (sigma_squared) is 0.01, for the Elman Neural Network, the maximum number of iterations (epoch) is 2000, the target error (goal) is 1e-5, and the learning rate (lr) is 0.01.

3. Machine learning methods

The acquired hyperspectral data was analysed using machine learning methods. The data set was randomly divided into training dataset (80% of the hyperspectral data) and testing dataset (20% of the hyperspectral data).

Machine learning algorithms are effective for expressing complex relationships (Chang, 2007). Convolutional Neural Network (CNN) is used for feature extraction. Models such as Data Classification for Multiple Linear Regression Based on Ridge Regression, Bayesian Linear Regression, and Elman Neural Network are introduced in our research. These machine learning models are classic algorithms known for their relative high efficiency and accuracy, which can be used for image classification. In our research, we designed these three models under the same frame.

Multiple Linear Regression Based on Ridge Regression is a valuable technique for improving the performance and stability of Multiple Linear Regression models, particularly when dealing with multicollinearity and the risk of overfitting (Luo & Liu, 2017). It strikes a balance between bias and variance, making it a useful tool in many real-world regression problems (Matdoan et al., 2021). Regression models for predicting rice yield and protein content using unmanned aerial vehicle-based multispectral imagery was employed by Kang et al. (2021).

Bayesian Linear Regression offers several advantages over traditional frequentist linear regression methods: Bayesian Linear Regression allows to incorporate prior information or beliefs about the model parameters (Baldwin & Larson, 2017). Bayesian Linear Regression provides a full posterior distribution over the model parameters, not just point estimates (Kong et al., 2020). Bayesian methods can be effective when dealing with small datasets, as they naturally incorporate prior information, which becomes more influential when data is limited (Barbier et al., 2021). Bayesian model averaging was used to improve the yield prediction in wheat breeding trials (Fei et al., 2023).

Elman networks can be designed with varying degrees of complexity by adjusting the number of hidden units in the recurrent layer (Wang et al., 2021). This flexibility allows you to control the model’s capacity, making it suitable for both simple and complex tasks (Thilagaraj et al., 2021). Elman networks can handle input sequences of varying lengths (Sriram et al., 2018). This is advantageous when dealing with data where the number of time steps may change from one example to another. Elman neural network is used as a rapid prediction method of moisture content for green tea fixation (Lan et al., 2022).

3.1. Convolutional neural network (CNN) feature extraction

Convolutional neural networks, as a deep feedforward network, are commonly used to process multiple arrays of data, such as time series, images, and audio spectrograms (Chen et al.,2023). This application showed that CNN is capable of learning the features of the hyperspectral data well when properly trained. A typical convolutional neural network consists of a convolutional layer, a pooling layer, and a fully connected layer, each with a different function. The convolution layer convolves the input vector through the convolution kernel to generate a feature vector, and all units in the same feature value share the same filter. The pooling layer is located behind the convolutional layer and is divided into maximum pooling and average pooling, which compresses the input feature information and simplifies the computational complexity of the network (Gu et al.,2018).

The spectral information is calculated by the convolutional layer and the maximum pooling layer (Li et al.,2021). The network model effectively extracts the local and global features that cannot be directly obtained from the original spectral data and convert the extracted feature inputs into vectors by the global average pooling layer (Yamashita et al.,2018). The fully connected layer is used to convert the inputs into vectors and to realise the classification function (O’Shea & Nash, 2015).

The structural system of the convolutional neural network is mainly composed of two parts, namely, the feature extractor and the classifier (Albawi et al.,2017). The feature extractor usually consists of a stack of several convolutional layers and a maximum pooling layer, and the classifier is usually a fully connected softmax layer (Yang & Li,2017).

At a convolution layer, the previous layers’ feature maps are convolved with learnable kernels and put through the activation function to form the output feature map. Each output map may combine convolutions with multiple input maps. In general, we have that

$$X_j^{\mathcal{l}}=\mathrm{f}({\sum }_{i\in M_j}{X}_i^{\mathcal{l}-1} \ast K_{ij}^{\mathcal{l}}+b_j^{\mathcal{l}})$$

(2)

Where $M_j$ represents a selection of input data, and the convolution is of the “valid” border handling type when implemented in MATLAB (Bouvrie,2006). Each output map is given an additive bias b, however for a particular output map, the input maps will be convolved with distrinct kernels. That is to say, if output map j and map k both sum over input map i, then the kernels applied to map i are different for output maps j and k.

3.2. Baseline model

3.2.1. Data Classification for Multiple Linear Regression Based on Ridge Regression

In Hoerl & Kennard (1970) it was shown that the optimal value of k, i.e., the value that minimizes the mean square error (MSE), is the following:

$$k_i=\frac{{\sigma }^2}{{\alpha }_j^2}$$

(3)

Based on the results obtained from the optimal value of k in that same article the authors sugested the following estimator of the ridge parameter:

$$\mathrm{K}1={\widehat{k}}_{HK}=\frac{{\widehat{\sigma }}^2}{{\widehat{\alpha }}_{max}^2}$$

(4)

Where ${\widehat{\sigma }}^2$=${\sum }_{i=1}^n{u}_i^2/(n-p-1)$ and $u_i$ are the residuals obtained from the OLS regression and ${\widehat{\alpha }}_{max}$ is the maximum element of $\widehat{\alpha }$. Hence, for this estimator, ${\sigma }^2$ and ${\alpha }_j^2$ are simly replaced by their unbiased estimators. Further developments were then made by Kibria (2003) where the following estimators were proposed:

$$\mathrm{K}2={\widehat{k}}_{GM}=\frac{{\widehat{\sigma }}^2}{{({\prod }_{j=1}^p{\widehat{\alpha }}_j^2}^{(\frac{1}{p})}}$$

(5)

and

$$\mathrm{K}3={\widehat{k}}_{MED}=\mathrm{Median}\left\{m_j^2\right\}$$

(6)

Where $m_j=\sqrt{\frac{{\widehat{\sigma }}^2}{{\widehat{\alpha }}_j^2}}$.

Khalaf and Shukur (2005) suggested a new method to estimate the ridge parameter k, as a modification of $k_{HK1}$:

$$\mathrm{K}4={\widehat{k}}_{KS}=\frac{t_{max}{\widehat{\sigma }}^2}{\left(n-p\right){\widehat{\sigma }}^2+t_{max}{\widehat{\alpha }}_{max}^2}\mathrm{}$$

(7)

Where $t_{max}$ is the maximum eigenvalue of $X^{\prime }\mathrm{X}$. Using the same idea as in Kibria (2003), Khalaf & Shukur (2005), and Alkhamisi et al. (2006) we have the following estimators for k:

$$\mathrm{K}5=k_{s3}={\widehat{k}}_{max}^{KS}=\max (s_i)$$

(8)

and

$$\mathrm{K}6={\widehat{k}}_{md}^{KS}=\mathrm{median}(s_j)$$

(9)

Where $s_j$ = $\frac{t_j{\widehat{\sigma }}^2}{\left(n-p\right){\widehat{\sigma }}^2+t_j{\widehat{\alpha }}_j^2}$.

In this article, we propose a modification of all of these estimators by multiplying them by the amount $w_i=\frac{t_{max}}{{\sum }_{j=1}^p{\left|\widehat{\alpha }\right|}_j}$ (10). The eigenvalues of the matrix of cross-products equal to one when the regressors are independent.

To investigate the performance of the RR and OLS, we calculate the MSE using the following equation:

$$\mathrm{MSE}=\frac{{\sum }_{i=1}^R{(\widehat{\beta }-\beta )}_i^{\prime }{(\widehat{\beta }-\beta )}_i}{R}$$

(11)

Where $\widehat{\beta }$ is the estimator of β obtained from OLS or RR, and R equals 2000 which corresponds to the number of replicates used in the Monte Carlo simulation (Khalaf et al., 2013).

3.2.2. Bayesian linear regression

Provided a dataset ${\left\{x_i,t_i\right\}}_{i=1}^N,x\in R^d$ is the input variable, $t_i\in \mathrm{R}$ is the corresponding target value, N is the number of data samples (Kong et al., 2020).

Regression aims at providing a specific predictive value $\mathrm{y}\left(x_i;w\right)$ given the input variable $x_i$. Form of linear regression mainly contains two types: standard linear model and kernelized model, which are given by Equations (12), (13), respectively (Bishop, 2006). The standard linear model is a linear combination of the elements of the input variable while the kernelized model is a linear combination of a set of nonlinear functions of the input variable.

$$\mathrm{y}\left(x_t;w\right)={\sum }_{i=1}^d{w}_i{x}_t^i+w_0$$

(12)

$$\mathrm{y}\left(x_t;w\right)={\sum }_{i=1}^N{w}_{i}K(x_t,x_i)\mathrm{}+w_0$$

(13)

where $w_i$ is the ith element of the weight vector w, $x_t^i$ is the ith element of the input variable $x_t$, $K(x_t,x_i)$ is the kernel function (Tipping, 2003).

The above two regression models can be uniformly expressed by:

$$\mathrm{y}={\mathsf{\Phi }}_w$$

(14)

where $\mathrm{y}={(y_1,y_2,\dots ,y_N)}^T$, $\mathsf{\Phi }={[\varphi \left(x_1\right),\varphi \left(x_2\right),\dots ,\varphi \left(x_N\right)]}^T$ is the $\mathrm{N}\times (\mathrm{M}+1)$ ‘design’ matrix, $\mathrm{w}={(w_0,w_1,\dots ,w_M)}^T·y_i$=y($x_i;\mathrm{w}),\mathsf{\varphi }(x_i$)is referred to as ‘basis function’.

$\mathrm{y}={\mathsf{\Phi }}_w$ is the standard linear model under the conditions of M=d and $\mathsf{\varphi }\left(x_i\right)=[1;x_i]$.

$\mathrm{y}={\mathsf{\Phi }}_w$ is the kernelized model under the conditions of M=N and $\mathsf{\varphi }\left(x_i\right)=\left[1,K\left(x_i,x_1\right),K\left(x_i,x_2\right),\dots ,K\left(x_i,x_N\right)\right].$

3.2.3. Elman Neural Network

The network is a single recursive network that has a context layer as an inside self-reference layer, see Figure 6. During operation, both current input from the input layer and previous state of the hidden layer saved in the context layer activate the hidden layer. Note that there exists an energy function associated with the hidden layer, context layer, and input layer (Liou, 2006). With successive training, the connection weights can load the temporal relations in the training word sequences.

Figure 6. the Elman network (Liou & Lin, 2006).

Figure 6. the Elman network (Liou & Lin, 2006).

The context layer carries the memory. The hidden layer activates the output layer and refreshes the context layer with the current state of the hidden layer. The back-propagation learning algorithm (Rumelhart & McClelland,1986) is commonly employed to train the weights in order to reduce the difference between the output of the output layer and its desired output. In this study, the threshold value of every neuron in the network is set to zero. Let $L_0,L_h,L_c,$ and $L_i$ be the number of neurons in the output layer, the hidden layer, the context layer, and the input layer, respectively. In the Elamn network, $L_h$ is equal to $L_C{,that is, L}_h=L_c$. In this study, the number of neurons in the input layer is equal to that in the output layer and is also equal to the number of total features, that is, R=${ L}_0=L_i$ (15)

Let $\left\{w_n,n=1-N\right\}$ be the code set of different words in a corpus. The corpus, D, contains a collection of all given sentences. During training, a sentence is randomly selected from the corpus and fed to the network sequentially, word by word, starting from the first word of the sentence. Let |D| be the total length of all the sentences in the corpus, D. |D| is the total number of words in D. Usually, |D| is several times the number of different words in the corpus, so |D|>N. Initially, t=0, and all weights are set to small random numbers. Let w(t) be the current word in a selected sentence at time t, i.e.,

$$\mathrm{w}\left(\mathrm{t}\right)\in \mathrm{D},\mathrm{}\mathrm{w}\left(\mathrm{t}\right)\mathsf{ϵ}\left\{w_n,n=1-N\right\},t=1-T$$

(16)

where $\mathrm{w}\left(\mathrm{T}\right)$ is the last word of a training epoch. In this study, we set T=4|D| in one epoch. This means that in each epoch, we use all the sentences in the corpus to train the Elman network four times. Let the three weight matrices between layers be $U_{oh},U_{hc},$ and $U_{hi}$, where $U_{oh}$ is an $L_h$ by $L_0$ matrix, $U_{hc}$ is an $L_c$ by $L_h$ matrix, and $U_{hi}$ is an $L_i$ by $L_h$ matrix, as shown in Figure 6. The output vector of the hidden layer is denoted as $\mathrm{H}(\mathrm{w}\left(\mathrm{t}\right))$ when $\mathrm{w}(\mathrm{t})$ is fed to the input layer. $\mathrm{H}(\mathrm{w}\left(\mathrm{t}\right))$ is an $L_h$ by 1 column vector with $L_h$ elements. Let $\mathrm{E}(\mathrm{w}\left(\mathrm{t}+1\right))$ be the output vector of the output layer when $\mathrm{w}(\mathrm{t})$ is fed to the input layer. $\mathrm{E}(\mathrm{w}\left(\mathrm{t}+1\right))$ is an $L_0$ by 1 column vector.

The function of the network is

$$\mathrm{H}\left(\mathrm{w}\left(\mathrm{t}\right)\right)=\mathsf{\phi }(U_{hi}w\left(t\right)+U_{hc}H(w\left(t-1\right))$$

(17)

Where $\mathsf{\phi }$ is a sigmoid activation function that operates on each element of a vector (Rumelhart & McClelland,1986). We use the sigmoid fucntion $\mathsf{\phi }\left(\mathrm{x}\right)=1.7159\ast \tanh (\mathrm{x}\ast \frac{2}{3})$ for all neurons in the network. This function gives a value roughly between $+$1.7159 and $-$1.7159.

4. Performance metrics

4.1. Confusion matrix

Performance measurement is vital in defining the effectiveness of a program. Confusion matrixes are a common evaluation tool used in machine learning (An, 2020). Generally, they consist of a n $\times$ n table plotting actual class against predicted class (n denoting the number of classes, so a binary classifier would utilize a 2$\times$ 2 table), to which the true and false (determined by the actual classes) positives and negatives (determined by the predicted classes) fit within (Visa et al., 2011).

4.2. Accuracy

In practical applications, we should take the accuracy of the classifier into consideration. Because scientists and farmers are more concerned with the situations where the classifier sorts the drought pear seedling leaves as sound ones if the classifier makes the wrong decision, which will hinder the timely watering of the plants, leading to more significant potential economic losses than discarding the plants.

$$\mathrm{Accuracy}=\frac{correctly classified samples}{total number of samples}$$

(18)

Accuracy (total correct divided by the total number of assessments), however, does not consider the significance of misidentified class (Halimu et al., 2019) and tends to be an overaly optimistic performance indicator.

4.3. evaluation metrics for classification algorithms

Commonly used evaluation metrics for classification algorithms include recall, precision, F1 score and mse-loss. Recall measures the ability to identify positive samples, precision measures the accuracy of positive sample predictions, and F1 score is a metric that combines recall and precision. MSE-loss is a criterion that measures the mean squared error between each element in the input and target.These metrics can be selected and weighted based on specific requirements. The formulas for these metrics are provided accordingly.

$$\mathrm{Recall}=\frac{TP}{TP+FN}$$

(19)

$$\mathrm{Precision}=\frac{TP}{TP+FP}$$

(20)

$$\mathrm{F}1 \mathrm{score}=\frac{2\times P\times R}{P+R}$$

(21)

Where TP is true positive, TN is true negative, FP is false positive, FN is false negative, P is precision, and R is recall (Shu et al., 2023).

All hyperspectral data processing was conducted in ENVI and statistical analysis were executed in MATLAB 2023. All experiments were performed under a Windows 10 OS on a machine with CPU Intel Core i7-7820HK @ 2.90 GHz, GPU NVIDIA GeForce 1080 with Max-Q Design, and 8GB of RAM.