Mid-Infrared Laser Spectroscopy Detection and Quantiﬁcation of Explosives in Soils Using Multivariate Analysis and Artiﬁcial Intelligence

: A tunable quantum cascade laser (QCL) spectrometer was used to develop methods for detecting and quantifying high explosives (HE) in soil based on multivariate analysis (MVA) and artiﬁcial intelligence (AI). For quantiﬁcation, mixes of 2,4-dinitrotoluene (DNT) of concentrations from 0% to 20% w / w with soil samples were investigated. Three types of soils, bentonite, synthetic soil, and natural soil, were used. A partial least squares (PLS) regression model was generated for predicting DNT concentrations. To increase the selectivity, the model was trained and evaluated using additional analytes as interferences, including other HEs such as pentaerythritol tetranitrate (PETN), trinitrotoluene (TNT), cyclotrimethylenetrinitramine (RDX), and non-explosives such as benzoic acid and ibuprofen. For the detection experiments, mixes of di ﬀ erent explosives with soils were used to implement two AI strategies. In the ﬁrst strategy, the spectra of the samples were compared with spectra of soils stored in a database to identify the most similar soils based on QCL spectroscopy. Next, a preprocessing based on classical least squares (Pre-CLS) was applied to the spectra of soils selected from the database. The parameter obtained based on the sum of the weights of Pre-CLS was used to generate a simple binary discrimination model for distinguishing between contaminated and uncontaminated soils, achieving an accuracy of 0.877. In the second AI strategy, the same parameter was added to a principal component matrix obtained from spectral data of samples and used to generate multi-classiﬁcation models based on di ﬀ erent machine learning algorithms. A random forest model worked best with 0.996 accuracy and allowing to distinguish between soils contaminated with DNT, TNT, or RDX and uncontaminated soils.


Introduction
The intensive use of high explosives (HEs) in military operations and mining excavations has contributed to soil contamination. Providing that HEs and their decomposition products that are highly persistent, mutagenic, and classified as Group C human carcinogens threaten human health [1,2], research on the timely detection of HEs has continued to receive considerable attention over the past few years. The methods currently used to detect HEs include gas chromatography-mass spectroscopy (GC-MS), gas chromatography-chemiluminescence (GC-CL), ion mobility spectrometry (IMS) [3], immunosensors [4], electrophoresis [5], fluorescence [6], high-pressure liquid chromatography (HPLC) [7,8], HPLC/mass spectrometry [7], and photo-assisted electrochemical detection [9]. However, none of these methods provides the required speed or accuracy for in situ detection of HEs in the presence of solid interfering materials.
Soil is considered to be a challenging matrix of organic compounds that interfere with HEs, making the detection of HEs in soils a difficult task [10][11][12]. While remote sensing has been applied to soils, the proposed system is complex [13]. Other studies conducted by the same research group involved the characterization [14][15][16], interaction [17][18][19], and detection [20,21] of HEs using Raman and Fourier-transform infrared (FT-IR) spectroscopy [22,23]. The detection of HEs was only marginally possible in all these cases providing that crystalline samples of explosives must be found within the solid matrix by microscopy to achieve the detection. The in situ detection of HEs in real soils with high selectivity and sensitivity has not been reported yet. In this study, a small, portable, and easy to handle system employing mid-infrared (MIR) quantum cascade laser (QCL) source [24] was used for analyzing HEs in soils at a distance of 15 cm. Furthermore, multivariate analysis (MVA) and artificial intelligence (AI) techniques were employed to quantify and detect the analytes of interest on the soil samples studied dosed with complex matrices of organic compounds [25].
External cavity QCLs provide ample wavelength tunability, justifying studies on threat chemicals [26][27][28][29][30] and biological compounds and microorganisms [31] with broad absorption features, including solids. QCLs were first demonstrated in 1994 [32]. These devices offer multiple benefits such as room-temperature operation, small sizes, long lifetimes, low energy consumption, long-term power stability, and fine-tuning of the output frequency [33]. Moreover, they provide an opportunity to devise portable systems for remote testing with excellent sensitivity, particularly on samples with low reflectivity or high absorptivity in the MIR range, such as soils. Compared to the existing sources in the MIR range, QCL sources have a higher output power that allows remote sensing and analysis. The rapid detection of HEs on highly reflective substrates using this method was demonstrated and reported by other investigators [34][35][36][37][38]. This spectroscopic analysis can be achieved without contact and in a non-destructive way. A study conducted at Pacific Northwest National Lab [39] demonstrated that the possibility of remote detection of many important HEs, including trinitrotoluene (TNT), cyclotrimethylenetrinitramine (RDX), and Tetryl, occurs in the spectral fingerprint window from 800 to 1400 cm −1 . This spectral window corresponds to the fingerprint region of the HE samples deposited on painted metal car doors. Other studies have demonstrated diverse QCL approaches to the rapid identification and characterization of HEs, given the routine requirements of security checks [40,41].
QCLs can provide significant benefits for public safety, particularly in locations such as airports, railways, bus stations, sports stadiums, and marathon sites. Therefore, this study is focused on the quantification and detection of HEs in soils [25] using a QCL source to enable remote sensing [42][43][44][45][46][47][48][49][50][51][52][53]. The limit of detection (LOD) was calculated for QCL spectroscopy based on the calibration curves of one HE in solid mixtures. The statistical figure of merit (FM) was achieved using the partial least squares (PLS) multivariate algorithm. Other FMs were also calculated.
Furthermore, AI methods were employed to achieve HE detection. In the future, two main applications can be developed based on the outcomes of this research. First, the proposed methods can be used to detect explosives on ordinary, non-ideal, low-reflective, real-world solid matrices. Second, the methods can be employed to detect HEs on natural solid matrices; for example, they can be used to locate landmines in sites affected by war and conflicts throughout the world, both for military purposes and implementing humanitarian demining applications.

Sample Preparation
Samples containing from 0% to 20% DNT w/w and a total mixture mass of approximately 0.20 g were prepared for the quantification analysis. The initial set of 10 samples and replicates entailed mixtures of several matrices with DNT. The matrices consisted of KBr, BC, SYN-S, and NAT-S. The SYN-S comprised 37% bentonite, 27% alundum™ cement, 16% montmorillonite, 10% silica gel, 8% washed sea sand, and 2% activated charcoal. The samples were prepared by grinding the HE into a fine powder using a mortar and pestle, followed by mixing in a mini vortex mixer for approximately 10 s at 3000 rpm. The mixed samples were ground again and mixed in the mini vortex mixer for a second time. Mixtures of KBr and NAT-S with lower concentrations of DNT, approximately from 0% to 3% by weight, were made and named KBr Low and NAT-S Low, respectively. Only these matrices were prepared and tested for comparing the variation of the results with the change of the concentration range. The preparation of these samples was similar to that of the others, except for using different DNT concentrations.
Mixtures containing other highly intrusive compounds such as explosives (TNT, RDX, and PETN) and non-explosives (BA and IBP) were prepared using the same preparation method with only one concentration per interference analytes, namely, 10% w/w. These samples were tested using the PLS model of DNT mixed with NAT-S. For the pattern recognition analysis based on AI, samples of mixtures containing from 0% to 20% of DNT, TNT, and RDX by weight in three soil types (BC and two NAT-S samples) were prepared.

Soil Characterization
Soil samples from Mayagüez, PR (USA) were characterized using thermogravimetric analysis (TGA) to measure the water content, sand and clay percentages, and total organic matter (TOM) [54]. TOM values were determined by oxidation with hydrogen peroxide [55,56] and calcination [55]. The total dissolved solids levels were measured by dissolving soil in water, followed by its filtering and drying. The percent distributions of components in the Mayagüez-PR NAT-S samples are illustrated in the Supplementary Materials: Table S1.

Data Acquisition and QCL System
The background spectrum of a KBr substrate was obtained before measuring the QCL spectra of the samples. Providing the lack of MIR signals from the employed solid matrix, this background spectrum provided an excellent and smooth reference trace. The first samples were then placed in the wells of metal holders (1.1 cm in diameter and 3 mm deep). Duplicate spectra were collected at 10 different locations on the sample surfaces, resulting in a total of 20 spectra per sample. Sixteen ranges were used for calibration and internal validation or cross-validation (CV), while the remaining four spectra were used for testing or external validation. This process was repeated for each concentration. The spectra were obtained in the reflectance mode at a distance of approximately 15 cm using a LaserScan™ MIR pre-dispersive spectrometer (Block Engineering, Marlborough, MA, USA) equipped with three tunable MIR lasers diodes having a tuning range from 990 to 1111 cm −1 , 1111 to 1178 cm −1 , and 1178 to 1600 cm −1 The spectral linewidth was < 2 cm −1 and the scan time was approximately 1.5 s for each of the diodes. The average power typically varied between 0.5 and 10 mW across the entire tuning range of ≈ 600 cm −1 with 100:1 Transverse Electromagnetic Mode (TEM 00 ) polarization and a beam divergence of < 2.5 mrad in the x-axis and < 5 mrad in the y-axis. The spectrometer had a 7.6-cm diameter ZnSe lens, which was used to focus the MIR beam to collect the reflected light and focus the light onto a thermoelectrically cooled mercury-cadmium-telluride (MCT) detector. The wavelength accuracy and precision were 0.5 and 0.2 cm −1 , respectively. The spectroscopic system worked best at a distance to the target of 15 ± 3 cm, with each laser producing an elliptical spot with diameters of 4 and 2 mm in the same space at a distance of 15 cm due to the difference of beam divergence in the axes (Galan-Freyle et al. [30]).

Data Quantification Analysis
The OPUS™ data acquisition and analysis software (v. 4.2 and v. 6.0, Bruker Optics, Billerica, MA, USA) were used to perform the multivariate data analysis. Calibration curves were generated from the PLS models of DNT in the studied solid mixtures, and the uncertainties and FM for this model were estimated. The accuracies of the multivariate calibration curves were evaluated using the root-mean-square error of estimation (RMSEE), the root-mean-square error of CV (RMSECV), and rootmean-square error of prediction (RMSEP) for external validation. These parameters were used as criteria for evaluating the quality of the proposed method.
The linearity of the calibration curves was evaluated using the values of the coefficient of determination (R 2 ), which indicates the percentage of variance present in the true component values reproduced by the PLS regression model. In contrast, the sensitivity (SEN) of multivariate methods can be estimated as the net analyte signal (NAS) [57] at a unit concentration as follows [58]: where b denotes the regression vector of the PLS model. The analytical sensitivity (γ) can be defined similarly to univariate calibrations [59] as where ε denotes instrumental noise. In the absence of interferences, the NAS would be equal to the intensity of the total analyte signal. The noise level was measured by collecting 20 spectra of a blank (KBr) and calculating the average of the standard deviations for all wavenumbers. The LOD was estimated using Equation (3) [60,61]: where h o denotes the distance of the predicted sample to the mean of the calibration set at zero concentration and ∆(α,β) is a statistical parameter correlated with the α and β probabilities of falsely stating the presence or absence of the analyte. ∆(α,β) = 3.3 was used to compute the LOD values providing that the value for the degrees of freedom was >25. The residual prediction deviation or relative predictive determinant (RPD) was used to represent how the calibration model predicts a specific set. This value can also be used to evaluate the performance of a model in absolute terms. The RPD values can be calculated as where SD denotes the standard deviation, SEP is the standard error of prediction, N stands for the number of samples, and Diff indicates the difference in concentration values of the analytes between the predicted and reference sets.

Pattern Recognition Analysis
An algorithm for comparing different machine learning (ML) methods for classification was developed in Python 3 using the sklearn 3.2 library [62]. Ten ML methods for classification were employed: K-neighbors classifier, support vector machine for classification (SVM), decision tree classifier, random forest classifier (RFC), AdaBoost classifier (ABC), gradient boosting classifier, linear discriminant analysis, and quadratic discriminant analysis. A basic description of each ML method used in this study is included in Supplementary Materials: Table S2. Each classifier was trained using 53% of data, while the remaining 47% of data were reserved for testing. The method achieving the highest accuracy value and lowest log-loss value on the test data was considered to be the most efficient.

Artificial Intelligence Scheme
The spectra were preprocessed before they were evaluated by the ML methods. The preprocessing scheme employed in this study is shown in Figure 1. First, all spectra were normalized using vector normalization (VN) preprocessing. Next, the spectral data were reduced using principal component analysis (PCA) from 3057 points to 20 PCs. Separately, the preprocessed spectra were compared with those in a database of soils. Subsequently, a classical least squares preprocessing (Pre-CLS) [36] was applied to the original spectra to extract a parameter that was proportional to the percentage of soil from the database spectra. Finally, PCA and the extracted parameter were used to generate the ML models.
Appl. Sci. 2020, 10, x 19 of 19 The criteria used to evaluate the performance of the classification models were recall, log-loss, precision, f1-score, weighted average, support, and accuracy. For binary classification, the recall of the positive class is also known as sensitivity, while the recall of the negative class is called specificity. The log-loss function is used in (multinomial) logistic regression and its extensions, such as neural networks. It is defined as the negative log-likelihood of the true labels given the predictions of a probabilistic classifier. The log-loss is defined only for two or more labels. For a single sample with the true label yt and estimated probability yp that yt = 1, the log-loss is Precision is often called the positive predicted value and calculated as the ratio TP/(TP + FP),

Vector
Normalization preprocessing Pre-CLS is based on a linear model of classical least squares represented as where f ϕ j , β j , ω i denotes the normalized calculated spectrum (CS) of soil derived from a mixture of several normalized spectra of soils ϕ(ω i ) j recorded in the database of soils spectra (DBS); ω i denotes the wavenumber and β j is a parameter indicating the fraction or proportion of ϕ(ω i ) j of a particular component in the CS. The model assumes that there are no binding interactions among the components in the mixture, which implies that the intensity contributions are additive. The β j parameters can be calculated by finding the minimum of the square difference of the normalized intensity between the real spectrum and the CS. The minimum value of the sum of the squares of d i (E) with respect to β j can be found by equating the first-order partial derivatives with respect to β j to zero and finding the β j values. Since the model contains n parameters, n partial derivative equations are generated as follows: The value of ϕ(ω i ) j can be determined using a simple forward selection algorithm (FSA). According to this algorithm, an empty model is generated first, and then variables (soils spectra from the database) are added one by one. At each next step, the E value shows the improvement of the model. The process is stopped when E cannot be further improved.
The criteria used to evaluate the performance of the classification models were recall, log-loss, precision, f1-score, weighted average, support, and accuracy. For binary classification, the recall of the positive class is also known as sensitivity, while the recall of the negative class is called specificity.
The log-loss function is used in (multinomial) logistic regression and its extensions, such as neural networks. It is defined as the negative log-likelihood of the true labels given the predictions of a probabilistic classifier. The log-loss is defined only for two or more labels. For a single sample with the true label yt and estimated probability yp that yt = 1, the log-loss is Precision is often called the positive predicted value and calculated as the ratio TP/(TP + FP), where TP denotes the number of true positives. FP denotes the number of false positives. Intuitively, precision is the ability of a classifier to not label a negative sample as positive. The f1-score is also known as the balanced f-score or the f-measure. The f1-score can be defined as a weighted average of precision and recall, where the f1-score best value is 1, and the worst is 0. The relative contributions of precision and recall to the f1-score are equal. The f1-score in the multi-class and multi-label cases is the average of the f1-score for each class with weighting that depends on the average between precision and recall. The f1-score can be calculated as The support is the number of used records, i.e., the number of spectra for each class. The classification accuracy score in multi-label classification computes the accuracy subset: the set of labels predicted for a sample should exactly match the corresponding set of actual labels in an ideal case.

Quantification Analysis
All spectra were converted to the Kubelka-Munk function [63] because the reflectance values for all samples were very low. This transformation is appropriate for the diffuse reflectance spectra of powders [64] with R values < 60% [65]. The models were generated using this transformation. The graphs plotting the predicted versus true values for the samples used in CV and testing are shown in Figure  of powders [64] with R values <60% [65]. The models were generated using this transformation. The graphs plotting the predicted versus true values for the samples used in CV and testing are shown in  The data set employed in CV was generated using the leave-one-out CV (LOO-CV). It can be noticed from Figure 2 that the predicted values for the CV and test data sets closely follow the best performance line (y = x). This alignment can also be confirmed by the low values of the RMSEE, RMSECV, and RMSEP listed in Table 1.  The data set employed in CV was generated using the leave-one-out CV (LOO-CV). It can be noticed from Figure 2 that the predicted values for the CV and test data sets closely follow the best performance line (y = x). This alignment can also be confirmed by the low values of the RMSEE, RMSECV, and RMSEP listed in Table 1. All models were generated using the complete spectral window (1000-1600 cm −1 ). In this study, VN was used as a preprocessing step as it proved to be better than the other tested preprocessing steps except for KBr. Other preprocessing steps were tested, such as mean centering (MC), linear offset subtraction, straight-line subtraction, minimum and maximum normalization, multiplicative scatter correction, first derivative, second derivative, and no preprocessing step. When applying VN, the average intensity is calculated first, and then this value is subtracted from the spectrum. Next, the sum of the squared intensities is calculated, and the spectrum is divided by the square root of this sum. Models for DNT in KBr were generated to evaluate the detection in the absence of interferences: one at high concentrations (0-20%) and another at low concentrations (0-3%). The errors for these models are listed in Table 1. The most effective preprocessing method for the KBr models was MC. This result suggests that VN is a suitable preprocessing step only when interferences from the matrix are present. In a model free from the interfering matrix (KBr), applying VN for generating samples with low concentrations had the same signal intensity as samples with high concentrations. This is reflected in the lousy prediction and high uncertainty for samples with low concentrations (see Supplementary Materials: Figure S3).
The spectra of the neat matrices, matrices with 20% DNT, and the corresponding standard reference spectra of DNT in KBr for transmittance are shown in Figure 3a-c, respectively. Figure 3d shows the resulting spectrum for the matrix with 2% DNT in NAT-S minus the spectrum of the neat matrix and corresponding standard reference spectra of DNT in KBr for transmittance. For clarity, DNT reference spectra are represented with the red, dotted lines. From these graphs, it is evident that DNT signals are observed on the background matrix signals in all cases, including 2% DNT in NAT-S.
All models were generated using the complete spectral window (1000-1600 cm −1 ). In this study, VN was used as a preprocessing step as it proved to be better than the other tested preprocessing steps except for KBr. Other preprocessing steps were tested, such as mean centering (MC), linear offset subtraction, straight-line subtraction, minimum and maximum normalization, multiplicative scatter correction, first derivative, second derivative, and no preprocessing step. When applying VN, the average intensity is calculated first, and then this value is subtracted from the spectrum. Next, the sum of the squared intensities is calculated, and the spectrum is divided by the square root of this sum. Models for DNT in KBr were generated to evaluate the detection in the absence of interferences: one at high concentrations (0-20%) and another at low concentrations (0-3%). The errors for these models are listed in Table 1. The most effective preprocessing method for the KBr models was MC. This result suggests that VN is a suitable preprocessing step only when interferences from the matrix are present. In a model free from the interfering matrix (KBr), applying VN for generating samples with low concentrations had the same signal intensity as samples with high concentrations. This is reflected in the lousy prediction and high uncertainty for samples with low concentrations (see Supplementary Materials: Figure S3).
The spectra of the neat matrices, matrices with 20% DNT, and the corresponding standard reference spectra of DNT in KBr for transmittance are shown in Figure 3a-c, respectively. Figure 3d shows the resulting spectrum for the matrix with 2% DNT in NAT-S minus the spectrum of the neat matrix and corresponding standard reference spectra of DNT in KBr for transmittance. For clarity, DNT reference spectra are represented with the red, dotted lines. From these graphs, it is evident that DNT signals are observed on the background matrix signals in all cases, including 2% DNT in NAT-S.

Figures of Merit (FM)
The accuracies of the models were determined from the values of the RMSEE, RMSECV, and RMSEP. The results for each model (labeled by matrix) are listed in Table 1. The precision was inferred through the values of the relative standard deviation (RSD). The RSD values for the prediction of the same sample in the same site were calculated for different concentrations and their average to measure the repeatability (RSDr). The RSD values for the prediction of the same sample in various locations on the sample surfaces were calculated for different concentrations and their average to measure the mixture homogeneities (RSDh). The RSD values for different samples at the same concentration were calculated for different concentrations and their average to measure the reproducibility (RSDrd). From these values, which are listed in Table 2, it can be concluded that the technique has excellent repeatability, while the samples have excellent homogeneity. However, only good reproducibility was determined for the models. Another FM used for measuring precision was the RPD. As previously mentioned, the RPD is the ratio of the variation in the validation samples and the size of probable errors occurring during predictions. The RPD values are also listed in Table 2.
The best models were the ones with the lowest RMSEP values and the highest RPD values. The values obtained for the RPD were excellent. It has been suggested that models with RPD values greater than 5 can be considered suitable for quality control. On the other hand, models with RPD values greater than 6.5 can be used for process monitoring. Moreover, models with RPD values greater than 8 can be used in any application [66]. All the models tested in this study had RPD values higher than 8. This result indicates that the proposed technique can be used for the direct analysis of HE in soils. At the same time, the detection of explosives and their respective effects over different soil types have not been studied in detail so far. We plan to expand the present study to include reliable models over various types of soils as part of our future work.
In the context of determining HEs, sensitivity refers to the ability of a method or instrument to detect an analyte at a specified concentration. In contrast, the sensitivity of an analytical approach as defined in this study, is the capability of the method to discriminate between small differences in concentrations or masses of a target analyte. The sensitivity of each method studied was calculated according to Equations (1) and (2), and the results of these calculations are presented in Table 3. The SEN value calculated for a model with VN preprocessing should be interpreted differently from the one without preprocessing or employing other preprocessing steps. Thus, the sensitivity values were derived from two types of models: models with VN as the preprocessing step and models without preprocessing steps. When using VN preprocessing, the spectra were normalized such that the sum of squares of intensity is equal to one. In other words, the spectra are converted to unit vectors. The SEN parameters were calculated from the magnitude of the vector calibration functions b, which, in turn, were derived from the spectra and their respective concentrations of standards. VN directly affects the magnitude b and value of SEN as a consequence. However, a better parameter for sensitivity can be obtained by calculating γ because this parameter is only affected by instrumental noises. The noise level was measured by collecting 20 spectra of a blank (target) and calculating an average of the standard deviations for all wavenumbers. The resulting noise for the models was different when considering 20 normalized spectra with VN (see Table 3). Otherwise, the γ-values calculated from the two types of models were very similar, which is an indication that γ is not affected by VN preprocessing or any other types of preprocessing. Table 3. Sensitivity, analytical sensitivity, the limit of detection (LOD), and noise of the models, with and without vector normalization (VN) preprocessing. The inverse of γ (denoted as γ −1 ) provides an estimation of the minimum concentration difference (resultant) discernible by the model considering the instrumental noise as the only source of error. In the case of NAT-S Low, γ was 0.003% for the model with VN preprocessing and 0.002% for the model without preprocessing, with the difference being statistically insignificant. It is not possible to make a comparison of the sensitivities between the modes with matrix and without matrix because the magnitude of b depends on the number of signals and number of latent variables (LV) (see Supplementary Materials: Figure S5). For the models of spectra with many signals (BC, SYN-S, and NAT-S), the magnitude of b is higher than that for the models with low-intensity signals (KBr models). A better FM for this comparison is the LOD. For the models with the matrix (BC, SYN-S, and NAT-S), the LOD values are close to that of the model without the matrix (KBr models). Curves for the samples of low concentrations of DNT in NAT-S were generated to determine whether the employed concentration range influenced the LOD values of the curves. In these cases, the LOD decreased from 0.8% to 0.3%. DNT concentrations between 0.3% and 0.8% in the NAT-S Low model can be quantified with higher uncertainties and higher probabilities of false positives and missed detections because RSD should be between 10% and 33%.

VN
To test the NAT-S model, a map of three new types of samples was analyzed, and a map of 10 × 10 mm 2 (100 points) for each sample was generated. Two samples comprised the same soil as that used in the NAT-S models and were contaminated with 10% of DNT. The first sample consisted of a simple mixture (NAT-S-M). The second sample involved mixing and macerating the components (NAT-S-MM). The concentrations were predicted using the NAT-S model with three preprocessing methods: VN, MC, and no preprocessing. The map for the NAT-S model with VN is included in Supplementary Materials: Figures S6-S9, while the predictions for the 100 points for each map using different preprocessing methods are shown in Figure 4a,b. VN preprocessing applied to both samples (NAT-S-M and NAT-S-MM) resulted in better results, with the predicted values being close to the true value on average (10% DNT, see Figure 4c). MC preprocessing worked better in NAT-S-MM than NAT-S-M. This indicates that while the macerated process homogenized the size of the particles, MC was not able to compensate for the difference in the particle size. In contrast, the VN preprocessing method was able to compensate for these differences. The prediction of DNT in NAT-S-M shows peaks of high DNT concentrations (>10% DNT). This is because the particles of DNT were not homogenized. The models with no preprocessing provided bad predictions due to the difference in the baseline of the spectra. The above discussion demonstrates that VN preprocessing corrects the spectral variation due to changes in the particle size. To demonstrate the change in spectrum with the particle size, the tested soil was sieving for three particle sizes (d): d > 0.85 mm, 0.85 mm > d > 0.25 mm, and d < 0.25 mm. One hundred spectra were acquired at various locations on the sample surface for each value of d, and the average and standard deviation of the spectra were determined (see Supplementary Materials: Figures S10 and S11). The background offset spectral decreases with d, whereas the standard deviation increases with d; however, this pattern is not consistent throughout the spectra. It is higher in the 1000-1200-cm −1 and 1400-1600-cm −1 regions. This can be explained by the fact that MC is not able to correct the background offset spectral completely in contrast to VN; VN is better because it scales the spectrum to unit vectors, whereas MC only changes the baseline.
Appl. Sci. 2020, 10, x 19 of 19 the difference in the baseline of the spectra. The above discussion demonstrates that VN preprocessing corrects the spectral variation due to changes in the particle size. To demonstrate the change in spectrum with the particle size, the tested soil was sieving for three particle sizes (d): d > 0.85 mm, 0.85 mm > d > 0.25 mm, and d < 0.25 mm. One hundred spectra were acquired at various locations on the sample surface for each value of d, and the average and standard deviation of the spectra were determined (see Supplementary Materials: Figure S10-S11). The background offset spectral decreases with d, whereas the standard deviation increases with d; however, this pattern is not consistent throughout the spectra. It is higher in the 1000-1200-cm −1 and 1400-1600-cm −1 regions. This can be explained by the fact that MC is not able to correct the background offset spectral completely in contrast to VN; VN is better because it scales the spectrum to unit vectors, whereas MC only changes the baseline. In the third sample, another type of natural soil (NAT2-S) was used to evaluate the NAT-S model. Two samples of NAT2-S were contaminated with 10% of DNT and mixed. Mappings of 10 × 10 mm 2 (100 points) were generated with the %DNT predicted from the NAT-S model using VN preprocessing. The mappings are present in Supplementary Materials: Figure S12, while predictions for each point are illustrated in Figure 4d. It can be noticed from the figure that the predictions were lower than the true values. This indicates that the NAT-S model quantifies below the true value of 10% DNT; however, it is capable of predicting the existence of an explosive in a soil. This indicates that the technique should be used in known soil to have good quantification.
To challenge the NAT-S model with other interferences, mixtures of other analytes in soil were prepared with a concentration of 10% (median of the model). The analytes used as interferences In the third sample, another type of natural soil (NAT2-S) was used to evaluate the NAT-S model. Two samples of NAT2-S were contaminated with 10% of DNT and mixed. Mappings of 10 × 10 mm 2 (100 points) were generated with the %DNT predicted from the NAT-S model using VN preprocessing. The mappings are present in Supplementary Materials: Figure S12, while predictions for each point are illustrated in Figure 4d. It can be noticed from the figure that the predictions were lower than the true values. This indicates that the NAT-S model quantifies below the true value of 10% DNT; however, it is capable of predicting the existence of an explosive in a soil. This indicates that the technique should be used in known soil to have good quantification.
To challenge the NAT-S model with other interferences, mixtures of other analytes in soil were prepared with a concentration of 10% (median of the model). The analytes used as interferences were BA, IBP, PETN, RDX, and TNT. BA and IBP do not have nitro groups but have aromatic rings and thus many common signals with DNT in the range of 1000-1600 cm −1 . PETN and RDX are nitro aliphatic explosives. TNT is very similar to DNT but considered to be a more challenging interference. Predictions for these samples were generated using the calibration curves for DNT/NAT-S. While the predicted values should have been zero or close to zero because the samples did not contain DNT, the average predicted concentrations were 8.8% (BA), 3.1% (IBP), -8.0% (PETN), 2.3% (RDX), and 25.8% (TNT) (see Table 4). The objective was to measure the model's capability of discriminating against these interferences. To improve the model with respect to the recognition of interferences, an optimization procedure was applied. This procedure involved implementing the optimization of the most significant regions of the spectra after using various preprocessing steps. During optimization, the spectral region was divided into equal spectral sub-regions. Then, the optimum combination of sub-regions was determined by starting with 10 sub-regions and successively excluding one sub-region. This procedure continued until the values of the cross-validation errors did not improve further. The RMSECV values were calculated for each combination of the preprocessing steps, and the models with the lowest values of errors were selected. Twenty spectra of each interference were introduced together with the validation spectra set off with 0% of DNT true value. Then, optimization was performed by minimizing the RMSECV value. This was labeled as an optimization of the validation set (OPT-Val).
In OPT-Val, all interferences were stabilized with the correct rejections. The parameters for the NAT-S OPT-Val model were approximately similar to those of the NAT-S model except for the LV. In particular, nine LV were added to the NAT-S OPT-Val model to stabilize the interferences (Table 5). In the optimization process, the entire region was selected, and the best pretreatment was found to be VN. Therefore, the only difference between the two models was that the optimized model had more LVs. This procedure was optimum for the elimination of the interferences; however, it assumed that the interferences do not interact with the analyte or the matrix. As part of our future work, samples of DNT in soil with interferences interacting with DNT can be added to the model to remove these potential errors. For the NAT-S OPT-Val model, 19 LV were required to obtain a good RMSECV value. Figure 5 shows the dependence of the RMSECV values for the prediction of each interference as a function of the LV. The NAT-S OPT-Val and NAT-S models were generated to determine how the RMSECV values decrease with the LV and compare the two models. The NAT-S model showed the minimum RMSECV at L = 10. Adding more LV resulted in worse RMSECV values. For the NAT-S OPT-Val model, the minimum RMSECV was achieved with 19 LV. Nevertheless, at 10 LV, the RMSECV values for other interferences were also sufficiently close to the RMSECV value for the model. For that reason, TNT and IBP interferences did not allow the RMSECV value of the NAT-S OPT-Val model to approach the RMSECV value of the NAT-S model.
For the NAT-S OPT-Val model, 19 LV were required to obtain a good RMSECV value. Figure 5 shows the dependence of the RMSECV values for the prediction of each interference as a function of the LV. The NAT-S OPT-Val and NAT-S models were generated to determine how the RMSECV values decrease with the LV and compare the two models. The NAT-S model showed the minimum RMSECV at L = 10. Adding more LV resulted in worse RMSECV values. For the NAT-S OPT-Val model, the minimum RMSECV was achieved with 19 LV. Nevertheless, at 10 LV, the RMSECV values for other interferences were also sufficiently close to the RMSECV value for the model. For that reason, TNT and IBP interferences did not allow the RMSECV value of the NAT-S OPT-Val model to approach the RMSECV value of the NAT-S model.

Pattern Recognition Analysis
All sample spectra were normalized using VN preprocessing. Then, they were evaluated using Pre-CLS. The soil database used for the evaluation initially contained QCL reflectance spectra of four different natural soils, sand, and clay standards. It was expanded using linear combinations of sand and clay by applying the concept of Self-Simulated Learning Artificial Intelligence (SSLAI). SSLAI consists of an Artificial Intelligence method that uses minimal information to develop the Machine Learning models. In this case, other types of soils are simulated using those that are available, and these new simulated soils are used for the analysis. FSA algorithm for Pre-CLS was applied for spectra of the samples, and one spectrum of the database was added to the model (see Equation (7)) for most spectra. Next, the sum of (SUM(β)) parameters was calculated without including the bias ( ), and a simple binary discrimination model was generated and evaluated (see Table 6) using this parameter to distinguish between contaminated and clean soils. Figure 6a shows the probability distribution of SUM(β) for the two classes of the binary model. In this figure was observed a broad distribution of SUM(β) for the samples contaminated by explosives (EXP), it is due to the range of concentrations. Additionally, when samples EXP has the highest concentration of HEs, the value SUM(β) approaches 0 because it has the lowest proportion of soil. On the other hand, in the NONE sample, the distribution is sharpest, and the SUM(β) is highest.

Pattern Recognition Analysis
All sample spectra were normalized using VN preprocessing. Then, they were evaluated using Pre-CLS. The soil database used for the evaluation initially contained QCL reflectance spectra of four different natural soils, sand, and clay standards. It was expanded using linear combinations of sand and clay by applying the concept of Self-Simulated Learning Artificial Intelligence (SSLAI). SSLAI consists of an Artificial Intelligence method that uses minimal information to develop the Machine Learning models. In this case, other types of soils are simulated using those that are available, and these new simulated soils are used for the analysis. FSA algorithm for Pre-CLS was applied for spectra of the samples, and one spectrum of the database was added to the model (see Equation (7)) for most spectra. Next, the sum of β j (SUM(β)) parameters was calculated without including the bias (β 0 ), and a simple binary discrimination model was generated and evaluated (see Table 6) using this parameter to distinguish between contaminated and clean soils.  Figure 6a shows the probability distribution of SUM (β) for the two classes of the binary model. In this figure was observed a broad distribution of SUM (β) for the samples contaminated by explosives (EXP), it is due to the range of concentrations. Additionally, when samples EXP has the highest concentration of HEs, the value SUM (β) approaches 0 because it has the lowest proportion of soil. On the other hand, in the NONE sample, the distribution is sharpest, and the SUM (β) is highest. The evaluation criteria (precision, recall, f1-score, and accuracy) and the confusion matrix for the model are shown in Table 6. This model has high precision in predicting clean soils (NONE), whereas it is moderate in predicting contaminated soils (EXP). This is because soils with low explosive concentrations are predicted as clean, increasing the number of false negatives. This happens at low concentrations (<2%), where the model does not work well. Following the scheme illustrated in Figure 1, PCA was applied to the normalized data providing 20 principal components accounting for 99.99% of the variance. These components were normalized using an auto-scale preprocessing, and the SUM (β) was added as an additional principal component. Then, various ML methods were used to generate multi-classification models to discriminate between soils contaminated with three types of explosives (DNT, TNT, and RDX) and clean soil (NONE). These models were evaluated using accuracy and log-loss calculated over the test data (Figure 6b). The best method was RFC. This method constructs a multitude of decision trees during the training process and outputs the class that is the mode of the classes of the individual trees. The evaluation criteria (precision, recall, f1-score, and accuracy) and the confusion matrix calculated for the RFC model over the test data are listed in Table 7. The test data were classified almost entirely. Other strategies such as the generation of models without SUM (β) were also evaluated (data not shown); however, these provided poor results with a high number of false positives, i.e., many samples of clean soil were classified as containing explosives. This was resolved by adding the SUM (β) parameter. The evaluation criteria (precision, recall, f1-score, and accuracy) and the confusion matrix for the model are shown in Table 6. This model has high precision in predicting clean soils (NONE), whereas it is moderate in predicting contaminated soils (EXP). This is because soils with low explosive concentrations are predicted as clean, increasing the number of false negatives. This happens at low concentrations (<2%), where the model does not work well.
Following the scheme illustrated in Figure 1, PCA was applied to the normalized data providing 20 principal components accounting for 99.99% of the variance. These components were normalized using an auto-scale preprocessing, and the SUM (β) was added as an additional principal component. Then, various ML methods were used to generate multi-classification models to discriminate between soils contaminated with three types of explosives (DNT, TNT, and RDX) and clean soil (NONE). These models were evaluated using accuracy and log-loss calculated over the test data (Figure 6b). The best method was RFC. This method constructs a multitude of decision trees during the training process and outputs the class that is the mode of the classes of the individual trees. The evaluation criteria (precision, recall, f1-score, and accuracy) and the confusion matrix calculated for the RFC model over the test data are listed in Table 7. The test data were classified almost entirely. Other strategies such as the generation of models without SUM (β) were also evaluated (data not shown); however, these provided poor results with a high number of false positives, i.e., many samples of clean soil were classified as containing explosives. This was resolved by adding the SUM (β) parameter.

Conclusions
This paper presented methods for quantifying DNT and detecting HEs in natural and synthetic soil matrices. The remote detection of HEs at a distance of 15 cm was achieved using multiple benefits of QCL spectroscopy, providing the evidence that it can be used in direct field applications. The obtained LOD values were adequate for the intended application. This study demonstrated that it is possible to get proper detection in the MIR range using diffuse reflection in the back-reflection mode for matrices containing components with low reflectivity such as soils. The remote detection of HEs in soils is now much more viable due to the high-power QCLs that are commercially available.
In the presented experiments, DNT signals could be observed in the spectra of solid mixtures with complex matrices, enabling the detection of explosives and quantification of DNT. For quantifying DNT in soil matrices, PLS models were useful for generating calibration curves. The best-performing models were obtained using VN preprocessing for all spectra. This preprocessing procedure proved to be suitable for samples, in which signals from analytes spectra overlapped with signals from matrices (natural soil, synthetic soil, and bentonite) spectra and relative intensity differences between spectra were small. In models that involved mixtures of explosives with KBr, VN preprocessing is not suitable because the KBr matrix does not have interfering signals in the spectroscopic window of interest. Therefore, this preprocessing step was not required. Hence, it can be concluded that VN is a useful preprocessing step for soil samples that have highly interfering matrices such as other explosives (TNT, RDX, and PETN) and non-explosives (BA and IBP). The strategy of adding samples of other analytes as possible interferences produced more selective models, making them more robust. The interferences used were organic compounds with many signals in the studied region. The interferences included explosives such as TNT, PETN, and RDX, as well as organic compounds such as BA and IBP. The NAT-S model was able to evade all interferences when it was trained with samples contaminated with these compounds. For the model to be accurate for other soils and possible interferences, it should be re-trained with corresponding samples.
Furthermore, an AI-based on an ML pipeline scheme (see Figure 1) was applied to obtain an RFC model for HEs detection in soil matrices. This scheme included Pre-CLS assisted by SSLAI and ML methods.
The expansion of the soil database generated by the linear combinations of QCL reflectance spectra of four different soils (natural soils, sand, and clay standards) allowed the application of Self-Simulated Learning Artificial Intelligence (SSLAI) concept. At the same time, it was possible to couple a simple forward selection algorithm (FSA) that facilitated the compression of the model when using simple multivariate analyses; CLS as additional preprocessing and PCA to reduce the number of spectral variables. Therefore, the combination of PCA, Pre-CLS, and SSLAI is an innovative ML pipeline that allows enriching the prediction of the model. Almost all tested samples of the studied soils and explosives (DNT, RDX, and TNT) were correctly classified by the model (see Table 7), providing an accuracy score of 0.996.

Supplementary Materials:
The following supplemental materials are available online at http://www.mdpi.com/ 2076-3417/10/12/4178/s1: Table S1: Composition percent of soil (NAT-S) samples; Figure S1: KM Spectra of 20%and 3%-DNT from the NAT-S sample; Figure S2: Vector Normalization of KM Spectra of 20%-and 3%-DNT from the NAT-S sample; Figure S3: PLS models for DNT in KBr using VN prepossessing; Figure S4: (a) b regression vector for models without preprocessing and with one Latent variable (LV), (b) b regression vector for models with vector normalization preprocessing and with one LV, (c) b regression vector for models without preprocessing and with eight LV, (d) b regression vector for models with vector normalization preprocessing and with eight LV; Figure S5: Plot of # Latent variable (LV) vs. analytical sensitivity for PLS models for DNT in KBr and Soil with VN preprocessing and without preprocessing; Figure S6: Map of % of DNT from the NAT-S model (VN) of the NAT-S-M sample; Figure S7: Map of % of DNT from the NAT-S model (VN) of the NAT-S-MM sample; Figure S8: Map of % of DNT from the NAT-S model (VN) of the NAT2-S sample; Figure S9: Map of % of DNT from the NAT-S model (VN) of the NAT2-S sample; Figure S10: QCL average spectra for different particle sizes of soil samples; Figure S11: QCL standard deviation spectra for different sizes of soil samples; Figure S12: QCL average spectrum for each spectrum in differently sized soil samples; Table S2: Basic description of ML methods for classification.