Detection of protein, starch, oil, and moisture content of corn kernels using one-dimensional convolutional autoencoder and near-infrared spectroscopy

Dataset description

In this study, the corn dataset, which is commonly used in the literature, was used to test the validity of the proposed method. The corn dataset contains the spectra of 80 corn kernels measured on three different devices. These devices are called M5 (FOSS NIRSystems 5000), MP5 (FOSS NIRSystems 5000), and MP6 (FOSS NIRSystems 6000). The wavelength range covered in the dataset is 1,100–2,498 nm at 2 nm intervals. In addition, the reference values for moisture, oil, protein, and starch targets of each corn kernel are also included in the dataset. The mean spectra for each device in the corn dataset are given in Fig. 1. The corn dataset can be accessed from (https://eigenvector.com/resources/data-sets/).

Autoencoders

Autoencoders are generative and unsupervised neural network algorithms. In this learning algorithm, the main goal is to get output values equal to input values. An autoencoder framework consists of two main blocks: encoder and decoder. The encoder block compresses input data into a low-dimensional representation called latent variables, which contains valuable input data information. The decoder block takes these latent variables as input and attempts to obtain the original data. An autoencoder framework also includes a waypoint, called bottleneck, between the encoder and the decoder (Zhang, Liu & Jin, 2020a). A simple diagram of an autoencoder is shown in Fig. 2.

For a given input data x_i, i = 1, 2, …, N, latent variables, h_i, can be obtained as: (1)${\displaystyle h_i=\psi \left(w_i\ast x_i+b\right)}$

where w_i denotes coefficients, b denotes biases and $\psi \left(x\right)$ denotes the activation function of the encoder layer. After the encoding process, the decoding process starts with obtained h_i using (2): (2)${\displaystyle y_i=\psi \left({\tilde{w}}_i\ast h_i+\tilde{b}\right)}$

here, ${\tilde{w}}_i$ and $\tilde{b}$ denote coefficients and biases of the decoder layer. In an ideal autoencoder, x_i and y_i are expected to be equal. During the training phase, the created network model tries to minimize the loss function, $J\left(\theta \right)$, by searching for optimal values for the weight and bias parameters. (3)${\displaystyle J\left(\theta \right)=\frac{1}{N}\sum _{i=1}^N{\left(x_i-y_i\right)}^2.}$

1D convolutional autoencoders

Convolutional neural network (CNN) is a special form of neural networks that uses convolution operation in layers. Previous works show that convolutional layers are more successful than fully connected layers in retrieving high-level features (Kiranyaz et al., 2021). Because of this, CNN has taken the field of artificial intelligence to an advanced level by adding a different depth. Similarly, a convolution operation can be applied to layers of the autoencoder network. Thus, a convolutional autoencoder can extract high-level features that can be used for classification or regression. CNNs are mainly used for high-dimensional data such as images; however, they can also be applied to low-dimensional data such as signals or time series.

The mathematical convolution operation of two discrete signals in one dimension can be defined as follows: (4)${\displaystyle y\left[n\right]=x\left[n\right]\ast h\left[n\right]=\sum _{k=-\infty }^{\infty }x\left[k\right]h\left[n-k\right].}$

Where $x\left(i\right)$, $y\left(i\right)$ and $h\left(i\right)$ are input, output, and filter vectors, respectively.

For a given vector, x, and 1D filter, w, whose length is m, the convolution formula can be reorganized as: (5)${\displaystyle \mathrm{conv}{\left(x,w\right)}_k=\sum _{i=1}^m{w}_i.x_{k+i-1}.}$

In the forward propagation step of the 1D CNN network, we can generalize formula (5) for each neuron in each CNN layer. (6)${\displaystyle x_k^l=\psi \left(\mathrm{conv}\left(x_k^{l-1},w_k^{l-1}\right)+b_k^l\right)=\psi \left(\sum _{i=1}^m{w}_k^{l-1}.x_{k+i-1}^{l-1}+b_k^l\right).}$

Here, $b_k^l$ is the bias of the k th neuron at layer l, $x_k^{l-1}$ is the output of layer l − 1, $w_k^{l-1}$ is filter coefficients of the k th neuron at layer l − 1, ψ(x) is the activation function of the current layer. An activation function is used to ensure the nonlinearity of the system. The selection of the activation function is one of the essential phases in creating a network model. Sigmoid, hyperbolic tangent (tanh), and rectified linear activation function (ReLU) are the most popular activation functions. Of these, the hyperbolic tangent is preferred when the input and output are constrained to values between −1 and 1 (Patil & Kumar, 2021). The formula of the hyperbolic tangent activation function is given in Eq. (7). (7)${\displaystyle \psi \left(x\right)=\frac{e^x-e^{-x}}{e^x+e^{-x}}.}$

A convolution layer is generally followed by a pooling layer. The pooling layer reduces features without changing the number of channels. Pooling layers do not contain any parameters, so no learning occurs in this layer. Max pooling and average pooling are the most popular pooling algorithms. An example of the 1D-Max pooling used in this study is given in Fig. 3.

A loss function needs to be utilized to evaluate the learning progress of the network. Mean squared error is the most preferred loss function for regression tasks, and its formula is given in Eq. (8). (8)${\displaystyle J=\frac{1}{N}\sum _{i=1}^N{\left({\hat{y}}_i-y_i\right)}^2.}$

Here, ${\hat{y}}_i$ is the predicted output, and y_i is the target output of the network, which is also equal to the input, x_i. For backward propagation, gradients must be calculated and propagated from the output layer to the input layer using the chain rule. (9)${\displaystyle \frac{\partial J}{\partial w_k}=\sum _{i=1}^{N-m+1}\frac{\partial J}{\partial x_i^l}.\frac{\partial x_i^l}{\partial w_k}.}$

In Eq. (9), $\frac{\partial x_i^l}{\partial w_k}$ term can be calculated as in Eqs. (10) and (11)

(10)${\displaystyle \frac{\partial x_i^l}{\partial w_k^l}=\frac{\partial \left(\psi \left(\sum _{i=1}^m{w}_k^{l-1}.x_{k+i-1}^{l-1}+b_k^l\right)\right)}{\partial w_k^l}}$ (11)${\displaystyle \frac{\partial x_i^l}{\partial w_k^l}={\psi }^{{}^{\prime }}\left(x_{k+i-1}^{l-1}\right).}$

Substituting Eqs. (11) into (9), the gradients needed to update the weights are obtained as in Eq. (12). (12)${\displaystyle \frac{\partial J}{\partial w_k}=\sum _{i=1}^{N-m+1}\frac{\partial J}{\partial x_i^l}.{\psi }^{{}^{\prime }}\left(x_{k+i-1}^{l-1}\right).}$

Various optimization algorithms in the literature have been proposed to update weights and biases, such as Stochastic Gradient Descent (SGD), RMSProp, Adam, and Adadelta. Among these, the Adam optimizer was used in our study. The Adam optimizer is based on adaptive moment estimation and combines Momentum and RMSProp (Kingma & Ba, 2014). To apply the Adam optimizer, firstly, the moving averages of the gradients and the moving averages of the square gradients, m^t and v^t, needs to be calculated using formulas Eqs. (13) and (14).

(13)${\displaystyle m^t={\beta }_1{m}^{t-1}+\left(1-{\beta }_1\right)\frac{\partial J}{\partial w_k}}$ (14)${\displaystyle v^t={\beta }_2{v}^{t-1}+\left(1-{\beta }_2\right)\frac{\partial J^2}{\partial w_k^2}.}$

Where, β₁ and β₂ are decay rate parameters. Using Eqs. (13) and (14), we can calculate bias corrected m^t and v^t.

(15)${\displaystyle {\hat{m}}^t=\frac{m^t}{\left(1-{\beta }_1^t\right)}}$ (16)${\displaystyle {\hat{v}}^t=\frac{v^t}{\left(1-{\beta }_2^t\right)}.}$

Using Eq. (17), we can update weights and biases. (17)${\displaystyle w_k^{\left(t\right)}=w_k^{\left(t-1\right)}-\eta \frac{{\hat{m}}^t}{\sqrt{{\hat{v}}^t}+\in }.}$

Here, η is named as the learning rate, another critical hyperparameter affecting the learning speed of the network.

Proposed model

In Fig. 4, the proposed 1D-CAE model is shown. This model consists of two convolutional layers, two pooling layers, and two dense (fully connected) layers in the encoder sub-model and three convolutional layers, two upsampling layers, and one dense layer in the decoder sub-model. The main reason for choosing two stacked convolutional layers in our model is that some previous works show that two or three layers are sufficient for CNN-based NIRS applications (Zhang et al., 2020b). Using random search, the optimal number of filters for each convolutional layer in the encoder and decoder model was determined and given in Tables 1 and 2. The hyperbolic tangent function was chosen to provide nonlinearity. The filter weights were randomly initialized. Training of the autoencoder model was done using randomly chosen samples. The reference values in the dataset were not used in this process which is unsupervised learning. The backpropagation algorithm was used to update the convolution filter weights. Although two, four, eight, 16, and 32 neurons were tried as the number of latent variables, no remarkable success was achieved in models containing fewer than 32 neurons, forcing us to select 32 neurons in our model. The optimization of unsupervised learning was done using the ADAM optimizer (Kingma & Ba, 2014). Selected values for learning rate, β₁ and β₂ were 0.001, 0.9, and 0.999, respectively. After the unsupervised learning training process, 32 latent variables for each corn kernel were exported for further analysis.

In the second part of the proposed model, multiple linear regression was employed to establish linear relations between latent variables and reference outputs. For each reference output, moisture, oil, protein, and starch, different MLR models were developed using the same latent variables. In addition, all the processes mentioned above were performed for three devices in the corn dataset to confirm that our results are device independent.

Dataset processing

In order to create a reliable model and to make accurate comparisons with known methods, the spectral data in the dataset were divided into two sets: one for calibration and one for prediction by the random division method. This way, the same samples from different devices were used for calibration and prediction. Of 80 samples, 60 were labeled as the calibration set, and the remaining 20 were labeled as the prediction set. While the calibration set was used to train the model, the prediction set was used to evaluate the model’s performance. The reason for splitting the dataset into two subsets is that for small datasets, the additional splitting can result in a smaller training set which can be subject to overfitting (Ashtiani et al., 2021; Féré et al., 2020). Statistics of the split dataset are given in Table 3.

Hyperparameters

In machine learning, tuning the hyperparameters of a model is an essential step that determines the performance of the model. In this work, the optimization of hyperparameters was carried out using random search with a lookup table. This table includes kernel size, the number of latent variables, and batch size. The lookup table is given in Table S1. During training, the maximum number of epochs was set to 20. Validation loss was tracked along the training of the model. When validation loss increased for two consecutive epochs, the training of the network was stopped to avoid overfitting.

Test environment

In this article, autoencoder and regression models were implemented in Python (version 3.7.13) using Keras (version 2.9.0), which is a high-level neural networks library (Chollet, 2022) and scikit-learn (version 1.0.2) which provides regression models and model evaluation tools (Lemaitre, 2021). All training and testing processes were performed using a computer with Intel i7 10870H CPU, 16GB Ram, and Nvidia RTX 2070 GPU.

Performance evaluation criteria

The coefficient of determination (R²), root mean squared error (RMSE), and root mean squared percentage error (RMSPE) indicators were used to test the evaluation of our proposed model. R² and RMSE indicators were used for overall model evaluation, while RMSPE was used to determine the performance of the autoencoder model in reconstructing the input spectrum. The formulas of R², RMSE, and RMSPE are given in Eqs. (18), (19) and (20) (Ashtiani, Salarikia & Golzarian, 2017; Chen et al., 2008; Miles, 2005).

(18)${\displaystyle R^2=1-\frac{\sum _{i=1}^N{\left(y_i-{\hat{y}}_i\right)}^2}{\sum _{i=1}^N{\left(y_i-{\overline{y}}_i\right)}^2}}$ (19)${\displaystyle \mathrm{RMSE}=\sqrt{\frac{1}{N}\sum _{i=1}^N{\left(y_i-{\hat{y}}_i\right)}^2}}$ (20)${\displaystyle \mathrm{RMSPE}=100\ast \sqrt{\frac{1}{N}\sum _{i=1}^N{\left(\frac{y_i-{\hat{y}}_i}{y_i}\right)}^2}.}$

Here N is the sample size, y_i, ${\hat{y}}_i$ and ${\overline{y}}_i$ are the actual output, the predicted output, and the mean value of actual outputs, respectively. As one can understand from Eq. (18), R² indicator is the proportion of the dependent variable variation explained by the independent variables and takes values between 0 and 1. RMSE, another indicator often used in regression tasks, is equal to the standard deviation of the residuals. Similarly, RMSPE gives the ratio of the error to the input spectrum in percent. Values closer to 0 are preferable for RMSE and RMSPE. Although most studies include the ratio of performance to deviation (RPD) metric to show model quality, some articles argue that RPD is not different from R² (Minasny & McBratney, 2013). For this reason, we did not find it necessary to include both metrics.

The 1D-CAE model was created separately for M5, MP5, and MP6 datasets in the first experiment. In this stage, the main objective was to obtain a reliable model that reconstructs the spectrum like the input spectrum and to obtain meaningful latent variables. The RMSPE indicator was utilized to show the reconstruction performance of the model, and the obtained results are given in Table 4. Besides, sample input and reconstructed spectra for each dataset are shown in Fig. 5. The 1D-CAE model successfully reconstructed the spectrum and obtained a mean RMSPE value of 1.90% for calibration and 2.27% for prediction.

The most common regression methods used in NIR systems, PLSR and PCR, were used as comparison methods. Another popular regression model, MLR, was not used because the sample number is lower than the feature number, which is necessary for MLR models. The latent variable and principal component parameters of the PLSR and PCR methods were selected as the optimal value between 1 and 10. RMSE and 5-fold cross-validation were used to find the optimal value for the latent variable and principal component parameters. Together with the original spectrum, four different preprocessing methods were applied to spectral data to increase the accuracy of these methods. Besides, the proposed method, 1D-CAE+MLR, was also applied to spectrum data. R² and RMSE values were calculated separately for each combination. The block diagram summarizing the whole process is given in Fig. 6.

The M5 dataset is the first used for the test, and the obtained results are given in Table 5. A careful examination of these results indicates that the proposed method exhibits satisfactory performance compared to other methods, as evidenced by the minimum R² values of 0.9560 and 0.9012 in the calibration and prediction sets, respectively. Furthermore, the proposed method demonstrated superior performance when predicting oil and starch content, as it yielded higher R² values and lower RMSE values. However, when analyzing the prediction of moisture and protein content, it was found that the PLSR method yielded higher R² values and lower RMSE values.

The MP5 dataset was employed in another experiment, and the results are presented in Table 6. The proposed method demonstrates superior performance, as evidenced by the higher R² values for all targets in the calibration set. Similar to the M5 dataset, the proposed method outperforms conventional methods when predicting oil and starch content, as it yields higher R² values and lower RMSE values. However, when assessing the prediction of moisture and protein content, it is observed that the PLSR method yields higher R² values and lower RMSE values. A notable difference is observed, particularly in the oil parameter, with an increase of 20.9% in the R² metric.

The MP6 dataset was utilized in the final experiment, and the test procedure was applied in the same manner as in previous experiments. The obtained results are presented in Table 7. Although the R² values were lower than those obtained in the M5 and MP5 datasets, the proposed method showed improved performance on all targets in the MP6 dataset, as evidenced by the higher R² values and lower RMSE values. Additionally, when analyzing the prediction of the oil and starch parameters, it was found that conventional methods were unable to establish a viable model, as the R² value was below 0.7.