Samples were collected from the tomato plants after 65 days of nitrogen stress treatment. During leaf collection, the healthy 7-leaf pinnate compound leaves of the tomato that best reflect the growth state were selected and cut off. They were immediately placed in a sealed bag to maintain freshness, which was placed in a portable refrigerated incubator (Ni et al., 2021) to prevent the external environment from affecting it. Twenty leaf samples were selected for each nitrogen stress gradient, and a total of 80 samples were collected from four gradients. The samples were then placed in the THz-TDS measurement system for sample scanning to obtain the spectral information. Before the experiment, a dehumidifier was turned on, and the relative humidity in the sample detection box was reduced to below 5% to eliminate the interference of water vapor on the THz spectrum. Ten sampling points were scanned for each sample to obtain the spectral information, and the average value was taken as the data collected for the sample.

The nitrogen content of the collected test samples was determined by the Nessler reagent colorimetric method, and the colorimetry was performed with a spectrophotometer at 420 nm. The actual measurement value of the nitrogen content of the sample was calculated by drawing the standard curve of the nitrogen content concentration and the photometric values, as presented in Figure 3.

The Savitzky-Golay (S-G) smoothing algorithm is a commonly used algorithm in data preprocessing due to its simple, fast, and easy-to-use process. The principle is to first take a window with an odd number of points in width, and to use the least-squares method for fitting via the translation of the window. The original value is then replaced with the fitted value of the midpoint of the window to achieve the smoothing of the data (Zhao et al., 2018). In this study, the S-G smoothing algorithm was used to preprocess the data, and the window width was 7 points/time. Taking the power spectrum data as an example, the comparison of the effect before and after the smoothing of the THz power spectrum is shown in Figure 4. The results demonstrate that the algorithm can effectively reduce the interference signal and improve the modeling efficiency and model accuracy.

The division of the sample set is the key to the applicability of the model. If the selected calibration set has good representativeness, the predictive ability of the model can be enhanced.

The random sampling (RS) algorithm is a simple method of randomly extracting samples regularly or irregularly from the entire sample set. A portion of the samples can be randomly selected as the prediction set until the sampling is full, and the remaining portion is used as the calibration set. The Kennard-Stone (KS) method was jointly proposed by Kennard and Stone (He et al., 2018), and its sample screening process of the calibration set is as follows. By calculating the Euclidean distance between two samples, the two samples with the largest distance enter the calibration set. The distance between the two selected samples is selected, the shortest distance among them is selected, and the sample corresponding to the longest distance among these shortest distances is entered into the calibration set. This method can ensure that the calibration set samples are evenly distributed according to the spatial distance, and the Euclidean distance between the two vectors is calculated as follows:

where d_x (p,q) is the Euclidean distance of the spectral reflectance between samples p and q, is the reflectance of sample p at the j^th wavelength point, x_q (j) is the reflectance of q at the j^th wavelength the reflectivity of each wavelength point, and J is the number of wavelength points.

Table 1 reports the results of the power spectrum and absorbance after sample division by the RS and KS algorithms.

The data in the calibration set obtained by sample division by the KS algorithm had a higher determination coefficient and a lower root-mean-square error (RMSE) than the data obtained by the RS algorithm. The subsequent data analysis and processing were therefore carried out on the basis of the division of the sample data set by the KS algorithm.

Uninformative variable elimination (UVE) (Jiang et al., 2019) can filter out spectral variables that contribute less to the modeling, which allows for the selection of representative spectral variables. The filtered variables are called uninformative variables, and by filtering them out, the complexity of the model and the number of variables required for subsequent modeling can be reduced. The UVE algorithm is based on the partial least-squares (PLS) algorithm. During variable screening, artificial noise variables equal to the original variables are added to the PLS model. The variables are randomly numbered, and one is left by crossover. This method is used to obtain regression coefficients of variables including artificial noise. In the analysis, the random change generated by artificial noise is used as a reference, and the reliability of each variable is measured by the threshold and the stability value. When the absolute value of the stability is less than the threshold, the variables in this part are regarded as non-informative variables. The stability value is defined as follows:

where i is the sample variable number, S_i is the stability value of the variable numbered i, b_i is the regression coefficient of the variable numbered i, std(b_i ) is the regression coefficient of b_i , mean(b_i ) is the mean value of b_i , and m is the total number of variables.

The steps of using the UVE algorithm to filter the feature frequencies are as follows.

Stability competitive adaptive reweighted sampling (SCARS) (Liu et al., 2014) is a common algorithm used to screen the optimal feature combination. During calculation, the measured THz spectral data can be set as a matrix X_N×P . The number of samples is N, P is the number of variables, and the specific operation steps of SCARS are as follows.

(1) The stability value c_j of each frequency-band variable is calculated as follows:

where c_j is the stability value of the j^th variable during Monte Carlo (M) sampling, b_j is the value of the j^th variable during M sampling, and s(b_j) is the standard deviation of the j^th variable during M sampling.

(2) The adaptive reweighted sampling method is combined with forced frequency band selection to screen groups with large stability values. They are combined into a subset of variables, and the ratio of variables to the whole frequency band is determined by the exponential decay function (EDF).

(3) These two steps are repeated in turn to obtain variable subset K. A PLS regression (PLSR) model is obtained, and the obtained variable subset is then evaluated through tenfold cross-validation. The K value is the number of operations in the SCARS algorithm, and the RMSE of cross-validation (RMSECV) value can be used as the judgment basis for whether the variable subset is a feature variable subset. The feature variable subset can be obtained at the smallest value of the RMSECV.

Interval PLS (iPLS) (De Lima et al., 2012) is a commonly used interval filtering algorithm for characteristic variables that was proposed by Norgaard at the beginning of the 20^th century. Based on the PLSR model, the algorithm divides the overall intervals into equal intervals to be filtered. The PLSR model of each equally spaced interval n is respectively established. By comparing the model accuracy of each subinterval, the subinterval with the best accuracy is selected as the modeling candidate frequency interval.

Radial basis function neural networks (RBFNNs) are developed based on the multi-dimensional spatial interpolation of radial basis functions. They are feed-forward neural networks with good performance and can be understood as function approximation or curve fitting in high-dimensional space. The proposal of the RBFNN provided new ideas and methods for the application and research of neural networks in various fields. In practical applications, this neural network has the advantages of a fast learning speed, no local minimum problem, and the ability to establish a corresponding network topology according to different types of data (Gao et al., 2014).

The RBFNN is mainly composed of an input layer, a hidden layer, and an output layer. When using this neural network to establish the model, P can be set as the sample input matrix of the correction set, and T is the sample output matrix of the correction set. The calculation formula of hidden-layer neurons can be obtained as follows:

where Q is the number of calibration set samples, and C is the center of the radial basis function. The connection weight between the hidden layer and input layer is set as W. Moreover, b₂ is the threshold value of N output layer neurons, and the following formula can be obtained.

By solving Eq. (5), the threshold value b₂ and the connection weight value W between the output layer and the hidden layer can be obtained as follows.

Backpropagation neural networks (BPNNs) (Li et al., 2020; Liu et al., 2021) can achieve the goal of not solving the relevant mapping equation in advance by learning the relationship between the input and output. Because of its high fault tolerance and parallel processing ability, the BPNN has a wide range of applications in target classification, recognition, and optimization.

The BPNN algorithm is mainly composed of forward calculation and error-reverse calculation. The topology includes an input layer, a hidden layer, and an output layer. During forward calculation, when the signal received by the output layer is an unexpected output value, a reverse propagation error signal will be generated, and the propagation path is the initial connection path. Compared with the RBFNN, the BPNN has the following characteristics.

Figure 6 shows the results of the uninformative variable screening.

In Figure 6, the left side is the THz spectrum variable, the right side is the artificially generated noise variable, and the abscissa and the ordinate are respectively the serial number and stability index corresponding to the variable. The larger the absolute value of the stability index corresponding to the variable, the greater the correlation with the model. The dotted line in the figure is the screening threshold set based on artificial noise. Based on this threshold, a total of 17 characteristic frequency bands were screened, as shown in Figures 7, 8.

In this study, when using the SCARS algorithm for filtering, the number of cycle sampling instances was 50 and THz power spectrum data were taken as the object. After 50 rounds of sampling had been reached, each index value had reached a stable state. The processing result of the SCARS algorithm is shown in Figure 9.

It can be seen from Figure 9 that with the increase of the number of iterations of the SCARS algorithm, the number of frequency bands retained generally decreased. However, the speed of reduction slowed down, indicating that the SCARS algorithm used coarse screening in the early stages of screening the characteristic frequency bands and fine screening in the later stages. It can be seen that with the gradual increase of the number of runs, the RMSECV value presented a decreasing trend. The RMSECV value was the minimum when the number of runs reached 43, which indicates that the frequency band with less correlation with the sample nitrogen content had been removed before. After 43 runs, the value had rebounded, and, combined with the change in the regression coefficient path, this indicates that the characteristic frequency bands related to the sample nitrogen content may have been removed by mistake. When the RMSECV value was the lowest, the subset of selected characteristic frequency bands was the best. Five characteristic frequencies were obtained, namely 0.574, 0.624, 0.642, 0.704, and 0.817 THz, and they were used as alternative characteristic frequency bands. The selected THz spectral bands are shown in Figures 10, 11.

When using the iPLS algorithm to filter the optimal frequency band interval, the number of subintervals to be divided must first be determined. Too many or too few subintervals will affect the subsequent selection of the optimal frequency band interval. In this study, the pretreated THz spectral intervals were divided into 10-45 equidistant intervals. As shown in Figure 12, when the RMSECV value was the minimum, the overall THz spectral interval was divided into 22 equidistant subintervals.

According to the number of subintervals shown in Figure 12, the THz spectral interval in the range of 0-1.4 THz was divided into 22 equal parts, and the PLS model index values in each subinterval were calculated. Table 2 reports the THz spectral frequency division range of each subinterval and the corresponding RMSECV value. Figures 13, 14 show the operation results of the iPLS algorithm.

It can be seen from Figure 13 that the index values of the PLS model established in each subinterval presented a trend of first decreasing and then increasing, and were closely related to the power spectrum curves of each stress gradient. When the power spectrum curve density was large, the corresponding modeling effect was poor; on the contrary, the modeling effect was better. Finally, the subinterval with the interval number of 12 was selected as the characteristic frequency subinterval, and there were 31 characteristic frequency points in total.

An RBFNN with strict logic was established by using the function package in MATLAB software. The expansion speed is the key to determining the quality of the RBFNN, and the setting of this parameter depends on the target object. Figure 15 exhibits the change of the RMSE of correction (RMSEC) of the RBFNN, iPLS, UVE, and SCARS models with an expansion speed of 0.1-1 and an interval of 0.1.

It can be seen from Figure 15 that when the diffusion speed was 0.5 and 0.6, the RBFNN model had the best RMSEC value. The established neural network model was tested to verify the prediction set, and the prediction results are reported in Table 3.

According to the results reported in Table 3, when the RBFNN algorithm was combined with the characteristic bands screened by SCARS, the detection model had the best effect. The RMSEC and the RMSE of prediction (RMSEP) were respectively 0.1322% and 0.1855%, and the determination coefficients of the correction set (R_c²) and prediction set (R_p²) were respectively 0.8714 and 0.8463. The scatter plot of the optimal results of the model based on the RBFNN algorithm is presented in Figure 16.

The sample data were trained using the scaled conjugate gradient function in the BPNN tool in MATLAB software. The specific training process included using the newff function to construct the basic structure of the BPNN. Moreover, the architecture based on the BPNN trained the model through the Trainlm function, and the Sim function output model was used to train the process. Taking the combination of characteristic bands screened by the iPLS algorithm as an example, the training process of its BPNN is shown in Figure 17.

It can be seen from Figure 17 that the number of network iterations epochs was 22. Moreover, the learning rate was set to 0.01, and the number of hidden-layer neurons was set to 6. The sample data filtered by other algorithms were input into the network after normalization, and the training results are reported in Table 4. By comparing the results in the table, it is evident that the modeling effect obtained by combining the SCARS filtering algorithm was the best.

According to the results reported in Table 4, when the BPNN algorithm was combined with the characteristic bands screened by SCARS, the RMSE values of the correction and prediction sets were respectively 0.1721% and 0.1843%, and the R² values of the correction and prediction sets were respectively 0.8447 and 0.8375. The scatter plot of the optimal results of the BPNN model is shown in Figure 18.

The analysis of the nitrogen detection model based on the THz spectrum shows that the models built by the two algorithms had good prediction effects. In terms of accuracy, the RBFNN model was slightly better than the BPNN model.