Vigna unguiculata (cowpea) and Chenopodium quinoa (quinoa) were cultivated in a Plant Growth Facility https://www.cranfield.ac.uk/facilities/plant-growth-facility) at Cranfield University (Cranfield, UK). The plants were grown in pots filled with compost with reconstituted nutrients. The compost was first washed with deionized water following the works of Masters-Clark et al. (2020) with some modifications. This was done by flooding one-part compost with five-part deionized water. A double-headed 0.5 µm sieve was used to drain off the soluble solution. The steps were repeated five times, ensuring the plant nutrients were depleted, and the washed compost was dried in an oven at 70°C for 72 hours. Based on the Letcombe nutrition solution (Masters-Clark et al., 2020), the required treatments were prepared as: high nitrogen high phosphorus (HNHP), high nitrogen low phosphorus (HNLP), low nitrogen high phosphorus (LNHP) and low nitrogen low phosphorus (LNLP), with five replications each for the treatments. The treatments had concentrations of 49.1 mM and 14.6 mM N for HN and LN and 13.4 mM and 3.3 mM P for HP and LP, respectively. By visual inspection, the different treatments resulted in changes in the leaf pigments, resulting in variations in the spectral characteristics during plant development. For both species, the leaves were green when the nitrogen and phosphorus levels were adequate, whilst deficiency in nitrogen was characterized by chlorosis progressing from light green to yellow to brown. Phosphorus deficiency inhibited shoot growth and decolorized leaves from a blue-green color to pale green/yellow in severely affected regions (McCauley et al., 2011).

Spectral images were collected 3.0 m above the canopy level using a Lemnatec Scanalyzer installed in the glasshouse. The Scanalyzer has a push broom hyperspectral camera (hyperspec^® inspectorTM Headwall Photonic) that operates within the 390 – 2500 nm region covering the visible–near-infrared (VNIR) and short-wave infrared (SWIR) regions. The sensor collects data with a 0.7 nm step (at the VNIR region) and 6 nm step (within the SWIR region) and an FWHM (full width at half maximum) image slit of 2.5 nm. This results in a hypercube of 1015 spectral bands for the VNIR data and 275 spectral bands for the SWIR data, with a dynamic range of 16 bits. Data collection was performed twice per week throughout the full development of the plants. Five different growth stages of cowpea and quinoa, hereafter referred as growth stages I, II, III, IV, and V, were considered for this study. They represent the 19, 29, 36, 50, and 69 DAT (days after transplanting) and 13, 30, 51, 60, and 72 DAT for the cowpea and quinoa plants, respectively. Supplementary Table 1 is the description of the selected growth stages based on the BBCH system for coding the phenological growth stages of plants (Meier et al., 2009).

Prior to any analysis, the raw data were pre-processed to normalize the spectral data from ambient illumination and to reduce noise and other artefacts that were produced during scanning. The pre-processing steps comprised: (i) radiometric calibration to remove illumination system effect, (ii) spectral down-sampling to remove redundant wavelengths, and (iii) noise reduction.

Radiometric calibration: Hyperspectral image acquisition suffers from radiometric errors caused by illumination from varying light exposure. Radiometric correction is essential to reduce the variable illumination effect and the influence of dark current on the spectral data. During scanning, a white panel (Zenith Lite™ Ultralight Targets 95%R, Sphereoptics^®) was imaged as the white reference data. Dark reference data with approximately 0% reflectance were collected in the night without any light source by completely covering the camera lens with an opaque cap (Zhu et al., 2013). The calibrated data were calculated using equation (1):

Where Ic is the calibrated data, I_raw is the raw spectral data, R_dark is the dark reference, and R_white is the white reference.

Down-sampling: Although the hyperspectral data contain essential plant information, the huge dataset poses computational challenges. Down-sampling helps to reduce these computational complexities and the noise generated during scanning (Sadeghi-Tehran et al., 2021). Moreover, down-sampling is conducted to map the spectral resolution to its reference target. In this study, an averaging window with a 2 nm spectral width was used to down-sample the spectral data.

Spectra smoothing and denoising: This is a common pre-processing practice that involves some numerical operations on the raw spectral data to reduce spectral noise levels. This eliminates spikes and smooths spectral curves whiles isolating essential features, which may be obscured by the presence of noise across the different wavelengths. During spectral smoothing and denoising, the original shape and features of the spectra are normally preserved. In our study, we applied the Savitzky–Golay filter (a commonly used low-pass filter), an effective and computationally fast filter, to smooth the spectra. For each data point in the spectrum, the filter selects an odd-sized window of spectral points and fits a least-square regression with a polynomial of higher order. During this operation, the data points in question are eventually replaced with the corresponding values of the fitted polynomial. We used a window of size 13- and second-degree polynomial as the optimal parameters. It should be noted that a small window size leads to the emergence of large artefacts in the smoothed spectra, whilst the larger the window size, the smaller the distinction between full and moving window processing (Rinnan et al., 2009).

Segmentation: After pre-processing, the HSI data were segmented using a selected spectral ratio combined with Otsu thresholding. To segment the VNIR dataset, 705 nm and 750 nm wavelengths were extracted to create a red-edge normalized difference vegetation index (RENDVI) image, which was then combined with automatic Otsu thresholding for segmentation. The RENDVI generally differentiates plant vegetation from non-vegetation regions. Hence, the combination with Otsu thresholding creates a binary image where the plant images are labeled as one and all other items labeled as zeros. The SWIR dataset segmentation used a similar approach, where the NDVI image was derived from wavelengths 1375 nm and 1141 nm, as proposed by Williams et al. (2017).

Although HSI contains important information that relates objects to its absorptance and reflectance, the measured spectra are subject to spurious signals such as multiplicative and additive effects. Hence, further pre-processing steps are applied to minimize the effect of undesirable occurrences on spectral measurement such as artifacts, particle size effects, and light scattering (Amigo, 2010). In this study, the effects of different pre-processing steps on the performance of classification models were studied by transforming the HSI data using standard normal variate (SNV), second derivative (SDeri), and linear discriminant analysis (LDA). In addition to the raw reflectance spectra, three datasets were generated from the SNV, SDeri, and LDA, respectively, to train the models.

Hyperspectral data contain a large amount of information that could be highly redundant and multi-collinear within adjacent wavelengths, resulting in computational complexities during processing and application (Hennessy et al., 2020). Feature mining methods that extract the most relevant and sensitive informative wavelengths from the spectral data are important to reduce multi-collinearity in the data and enhance model robustness. In this work, a correlation-based feature selection (CFS) method was applied to the best group of spectral data. CFS is one of the most popular data engineering methods for selecting a sensitive set of features to build a discriminative model for a specific purpose. It works on the principle that good and informative wavelengths are those that are highly correlated with a particular class but uncorrelated with each other. To use the CFS algorithm, a heuristic search algorithm was applied along with a correlation function (Pearson’s correlation) to assign high scores and select the best subset of features that had high predictive ability of the class labels with poor correlation with each other. A greedy stepwise (GS) search strategy was applied during the application of the CFS algorithms to select the space attributes of the variable subsets.

To evaluate the performance of the proposed model against the classification models, the overall accuracy (OA), the F- Score, and the kappa coefficient (Kappa) evaluation metrics were considered. OA represents the number of correctly classified samples in the overall test samples, which is given by equation 5. The F-score evaluates the accuracy of the model on the entire dataset and is useful when there is an uneven class distribution. It is given by the expression shown in equation 6. The kappa coefficient is a statistical measure that describes the mutual information between the ground truth map and the predicted classified map. Kappa values ranges between 0 and 1, such that a coefficient of 1 means perfect agreement with predicted class and ground truth and vice versa when the coefficient is 0.

where OA is the overall accuracy, N is the number of all samples, x is the number of class labels present in the dataset, and dii is obtained from the diagonal element of the confusion matrix.

where Recall = TP TP  +  FN and Precision = T P TP + FP and TP is true positive, FN is false negative, and FP is false positive.

where P_o is the probability of observed agreements and P_c is the expected agreement.

Figures 3 and 4 represent the average spectral curves for cowpea and quinoa, respectively, at growth stage III. For both species, although the spectral curves are similar in shape, there are visible differences between the four treatments. There is a characteristic peak in the visible spectra observed, especially for wavelengths between 550 nm and 560 nm, for both plant species. This indicates a high amount of chlorophyll absorption in this region (Siedliska et al., 2021), which shows that chlorophyll could be a responsive feature for the dynamics of N and P in both plants. The region between 600 nm and 700 nm (red band) shows a clear and continuous distinction between some variants of the experiment for both species. The low-nitrogen treatments exhibited higher reflectance in contrast to the high-nitrogen treatments, which had low reflectance in this region. HNLP had the lowest reflectance, while LNLP had the highest reflectance in cowpea. Similarly, LNLP and HNHP had the highest and lowest reflectance in the red-edge region for quinoa. This agrees with Clevers and Gitelson (2013), who demonstrated rapid change in the reflectance of the plant canopy in the red-edge region in response to high levels of nitrogen. The regions between 750 and 1100 nm (near infrared) exhibited the highest differences between the spectral curves of the two species. The HNHP had the highest reflectance in both species at above 0.62 nm and above 0.64 nm for cowpea and quinoa, respectively. This agrees with Zhai et al. (2013), who showed that reflectance of healthy plants (in this case, HNHP plants) in the NIR region has a strong correlation with the chemical and cellular architecture of plants and exhibits high reflectance. In the SWIR region, there was a characteristic drop in reflectance from 1250 to 1500 nm and a gradual rise and fall in reflectance across the spectra from 1500 to 2000 nm. Moreover, there was comparably low reflectance in the two - related treatments at 1450 nm, which agrees with Siedliska et al. (2021), who suggest that this region has a quantitative relationship between light reflectance and P treatments. From wavelengths 2000 nm to 2500 nm, there was an observable increase in reflectance for both species.

Figures 5 and 6 are box and whisker plots summarizing the accuracy distribution scores for the different models of the two species. Tables 3 , 4 summarize the performance of each model based on the evaluation metrics used in this study. From Figure 5 , the pre-processing methods had an effect on the performance of the different models. SDeri datasets had higher accuracies in all four models, achieving the highest OA (99.24%), kappa (98.64%), and F-score (99.19%) for cowpea and 99.18% OA, 98.85% kappa, and 98.76% F-score for quinoa. The high OA and F-score show the good generalization capacity of the models in classifying the nutrient status. In contrast, the RAW dataset produced the lowest performance for all the models, with the hybrid CNN model achieving the highest performance for cowpea (98.39%, 98.37%, and 97.81% OA, F-score, and kappa coefficient, respectively). In quinoa, the HybridSN model achieved the best performance when trained on the raw dataset with 98.57%, 98.43% and 97.76% OA, F-score, and kappa coefficient, respectively. On the contrary, the 2DCNN model (trained with raw dataset) had the lowest performance with 92.77% OA for cowpea. Furthermore, the 2DCNN model had the lowest performance (91.15% OA) when trained with the LDA dataset for quinoa. There was a 3.58% and 1.57% reduction in accuracy performances when the hybrid model was trained with the SNV and LDA datasets, respectively, for cowpea as compared to the performance when trained with the SDeri dataset. Similarly, the proposed model decreased in performance when trained with the SNV and LDA datasets, achieving 96.32% and 91.45% OA scores, respectively, for quinoa. This suggests that the SDeri dataset had important learnable wavelengths to discriminate the treatments for both species. Although the LDA dataset had reduced training samples, all the LDA-based models displayed lower performance compared to the SDeri dataset-trained model’s performance.

The proposed model for both species exhibited good performance when trained with SNV- and SDeri-transformed datasets. Although the HybridSN and 3D CNN models had good performance when trained with SDeri datasets, their training time (caused by the high number of trainable parameters) made the implementation of this model a challenge. The proposed hybrid CNN model had reduced training parameters (< 200,000 trainable parameters), subsequently affecting the computational processing time. The SDeri-trained hybrid CNN model had a processing time reduced by over 1000 s as compared to the 3DCNN-SDeri model, which had a 2141.63 processing time for cowpea. For quinoa, the SDeri-hybrid CNN model had about 54.16% reduced training time when compared to the 3D CNN model (trained with a raw dataset, which had the highest training time of over 2000 s). For both species, the HybridSN-based models had a significantly high training time irrespective of the training dataset.

Although the proposed model with the SDeri dataset produced the best performance in classifying plant nutrient status, the high dimensionality of the hypercube posed computational challenges. Hence, effective wavelengths that accurately distinguish between the N and P levels were selected and used to retrain the models. From Figure 7 , the selected wavelengths for cowpea were localized in the blue, red-edge, near-infrared, and short-wave infrared regions (having the highest number of selected wavelengths). Sixteen (16) wavebands were selected as the most sensitive wavelengths, with three selected from the blue regions (411, 431, and 455 nm). Moreover, two were selected from the red-edge regions (683 and 691 nm), four from the near-infrared regions (791, 871, 919, and 923 nm), and seven from the SWIR regions (1438, 1450, 1461, 1485, 1508, 1602, and 1660 nm).

In quinoa, the wavelengths were selected within the blue, near-infrared, and SWIR regions. From Figure 8 , twelve (12) wavebands were selected. These include four from the blue region (571, 575, 583, and 599 nm), one from the red spectrum (607 nm), four the near-infrared regions (723, 731, 924, and 967 nm), and three from the SWIR regions (1680 nm and 17520 nm). The differences in the selected wavebands could be due to the variability in the plant structure and the changes in the nutrient composition of the plants.

After the waveband’s selection, the full plant (with full spatial resolution and selected spectral wavelength) in the second derivative format was used as a dataset to retrain the different models to predict the plants’ N and P statuses at different growth stages. Datasets were generated using the same approach as that in subsection 2.4.2. The training datasets were augmented to artificially enlarge the number of training images using rotation and flipping methods.

Tables 5 and 6 show the performance of the models when trained with the selected wavelengths for cowpea and quinoa, respectively. All the models performed well (with 79-100% accuracy on the test data) irrespective of the growth stage under consideration. This performance was better than that found in Siedliska et al. (2021), where a 40-100% classification accuracy rate was achieved for predicting different P levels of three plant species at five growth stages. Comparing the four models, the proposed hybrid model had the best performance, attaining above 94% and 97% test accuracy across the selected growth stages for cowpea and quinoa, respectively. The 2DCNN model had the highest misclassification, with 84.08% in cowpea at growth stage IV and 79.23% in quinoa at growth stage II in the test dataset. For each species, there were high classification errors at the early growth stages, especially in the LNHP and HNLP plants. However, as the plants developed, there were clear distinctions between the spectral signatures of the different treatments, which resulted in an increase in the classification accuracy.

Table 7 is the summary of the confusion matrix for the proposed model, which explains the accuracy and misclassification of the individual classes for cowpea and quinoa. For both species, the model accurately classified HNHP treatment across all the five growth stages. The model had difficulties in classifying some variants of the treatments, especially those with low phosphorus. In cowpea, HNLP (5% error) and LNHP (> 10% error) were misclassified at growth stage I. Moving to growth stage II, although a similar trend occurred, there was an improvement in the classification accuracy of LNHP (3.64% increment). In the subsequent growth stages, the model improved its performance, with > 95% at growth stages III, IV, and V. In quinoa, there was a 100% accuracy in the HNHP and HNLP classification. However, the model had lower performance in correctly classifying the LNHP treatment 92.82%. Moving across the growth spectrum, the model performance increased, achieving above 98% classification for HNHP and HNLP at growth stages IV and V.