Application of Computational Intelligence Methods for the Automated Identification of Paper-Ink Samples Based on LIBS

Laser-induced breakdown spectroscopy (LIBS) is an important analysis technique with applications in many industrial branches and fields of scientific research. Nowadays, the advantages of LIBS are impaired by the main drawback in the interpretation of obtained spectra and identification of observed spectral lines. This procedure is highly time-consuming since it is essentially based on the comparison of lines present in the spectrum with the literature database. This paper proposes the use of various computational intelligence methods to develop a reliable and fast classification of quasi-destructively acquired LIBS spectra into a set of predefined classes. We focus on a specific problem of classification of paper-ink samples into 30 separate, predefined classes. For each of 30 classes (10 pens of each of 5 ink types combined with 10 sheets of 5 paper types plus empty pages), 100 LIBS spectra are collected. Four variants of preprocessing, seven classifiers (decision trees, random forest, k-nearest neighbor, support vector machine, probabilistic neural network, multi-layer perceptron, and generalized regression neural network), 5-fold stratified cross-validation, and a test on an independent set (for methods evaluation) scenarios are employed. Our developed system yielded an accuracy of 99.08%, obtained using the random forest classifier. Our results clearly demonstrates that machine learning methods can be used to identify the paper-ink samples based on LIBS reliably at a faster rate.


Introduction
During the last decade, there have been important developments in laser induced breakdown spectroscopy (LIBS). This atomic emission spectroscopy technique, also known as laser induced plasma spectroscopy (LIPS), is used for qualitative and quantitative chemical analysis of samples in all states of matter [1,2]. In this technique, high-power, short-duration laser pulse causes an ablation of the analyzed material, which due to its high temperature (10,000 K) dissociates into excited ions and atoms. When plasma cloud cools down this excited species revert to lower energy states and emit optical radiation which can be recorded and analyzed, revealing information about the elemental composition of the sample [2]. LIBS spectra are generally very rich in emission lines coming from excited atoms and ions occurring in a high temperature plasma cloud. Physical and chemical phenomena behind LIBS are not fully understood yet, as they are very complex in nature [3,4]. Nevertheless, LIBS applications have recently been rapidly growing due to the number of advantages of this method. The most important ones are minimally destructive measurements with little or no sample preparation, no efficiency, or no possibility to analyze in real-time all elements in a single laser shot in the samples of all three states of matter [5]. For solids, both mapping (2D) and depth profiling (3D) can be obtained. LIBS can be used for quantitative chemical analysis, material identification, and discrimination [1]. This method can be applied in laboratories and industrial plants, and even at stand-off distances of tens of meters [6]. These properties predispose the LIBS technique for use in many fields: space exploration [7], remote analysis of hazardous materials [8], on-line quality control in various industries [9], cultural heritage studies [10], forensic chemistry [11,12], geology [13], weld quality assurance [14], robotics [15], and many others [16][17][18][19].
Although currently LIBS is already an established technology, spectrochemical LIBS analysis is not straightforward [20]. Identification of elemental constituents in the sample is usually based on the strongest lines present in the LIBS spectrum (so called persistent lines [1]), which are compared with a literature database collected for all the elements from the periodic table. This analysis is cumbersome and time-consuming because the emission spectrum is determined also by the properties of the plasma, not only by the composition of the examined sample [21].
At present, the great advantages of LIBS seem to be impaired to some extent by the main drawbacks: problems coming from often poor reproducibility, the impact of sample composition on spectra recorded for individual components (matrix effects), and the difficulty of performing an overall reliable data analysis (sometimes over 500 spectral lines need to be interpreted). Thus, a question is immediately raised as to whether the LIBS method can be supported by the modern achievements in the field of computational intelligence [22] in an attempt to overcome the limitations of the current LIBS spectrum analysis methodology. The analogous problem of classification of a one-dimensional series of data is well known in computer science. The methods developed for this problem have been successfully applied in many areas [23][24][25][26][27][28][29][30][31]. Therefore, the application of machine learning methods to LIBS spectra of samples may certainly give a strong impact for further developments and applications.
Motivated by the aforementioned arguments, we decided to apply computational intelligence methods to the classification of paper-ink samples for forensic purposes, based on LIBS spectra of the samples. The problem of discrimination of different paper-ink samples has already been addressed in previous studies based on the LIBS spectrum analysis. Trejos et al. showed [32] that the highest discrimination power (DP) (96.4%) was obtained when comparisons were done qualitatively by spectral overlap of the regions of interest (three different emission lines monitored per element) and quantitatively followed by pairwise comparisons (one emission line per element) using ANOVA (analysis of variance). Kula et al. [33] presented the LIBS method as a useful tool in qualitative elemental differentiation of ink samples. They obtained the discrimination power (DP) coefficient of 61%, 82%, and 83% for red, black, and blue inks, respectively. Elsherbiny and Nassef [34] studied the dependence of the obtained spectra of various black gel inks on the wavelength of laser excitation, reporting the DP in the range of 88-91%. In the next study [35], the elemental analysis with the use of LIBS was performed with inkjet inks and toners on office paper. The LIBS results supported by pairwise comparison analysis (ANOVA with Tukey's post hoc test) provided discrimination power of 98.4% (3 indistinguishable/190 compared pairs) for the toners and 100% for the inkjet inks. Moreover, these three undifferentiated toner pairs were discriminated using a Student's t-test at a 95% confidence limit. The authors claimed that LIBS as a tool for the determination of elemental composition of sample can be a part of a procedure for questioned document examination.
The problem of paper-ink classification, on the other hand, is less common in the literature. To the best of authors' knowledge there are only two articles dealing with LIBS spectra of paper-ink samples employed in classification/fitting problems by means of computational intelligence techniques. In a paper by Hoehses et al. [36], the benefit of applying several independent chemometric methods to LIBS data was demonstrated. A consecutive methodology of applying soft independent modeling of class analogy (SIMCA) and partial least-squares discriminant analysis (PLS-DA) enabled the step-wise classification of data and separation of inks that were not identified by principal component analysis (PCA). The support vector machine (SVM) yielded a correct classification rate of 87%, and cross-validation accuracy amounted to 81%. In the second paper [37], multiple methods such as three comparative functions (linear correlation, overlapping integral, and sum of squared deviations) and two advanced statistical methods (multivariate curve resolution alternating least squares (MCR-ALS) with classification tree and discriminant analysis (DA)) were applied to statistically evaluate LIBS spectra. The newly introduced MCR-ALS/DA methodology showed identification of the paper and printer type with an accuracy of 96.3% and 83.3%, respectively. Computational intelligence methods analyzed in this study have a larger hypothesis space than SIMCA, PLS-DA, MCR-ALS, or SVM; therefore, they can fit the data more accurately.
In the present paper, we show that the similar problem can be solved using computational intelligence methods. The classification problem discussed in this paper is very difficult due to similarities between LIBS spectra from different paper-ink samples and differences between such spectra from one single class. Additionally, we operate over spectra with a large number of lines with significant noise. To solve this classification problem, we tested many preprocessing ways and computational intelligence methods, but the paper presents only selected and best ones.

Materials
Fifty ballpoint pens (10 items of each of 5 models) produced by four different manufacturers from Germany and Poland were purchased in Poland: Bic Orange fine blue (denoted B), Rystor Kropka 0.5 1986 (denoted R), Staedtler Stick 430F (denoted S), Staedtler Ball 432F (denoted SB), and Toma Sunny fine 050 (denoted T). Fifty sheets (10 sheets of each of 5 types) of five Canadian certified reference papers were used (papers denoted A, D, L, N, and O from the set "Fillers in paper" supplied by A.S.O. Design Canada).
Each particular ink from each of 10 pens of 5 types was deposited (as straight lines) on each of 10 sheets of 5 types of papers, using standard hand pressure. All 2500 (50 sheets of papers with 50 deposits each) paper-ink samples were placed in plastic bags and stored in darkness at room temperature.
In the experiment, data for 30 classes (A, A + B, A + R, A + S, A + SB, A + T, D, D + B, . . . , O + SB, O + T) based on the combination of ink, paper, and empty papers were recorded. They are summarized in Table 1.

LIBS
The analysis of all paper-ink samples was carried out using a laser induced breakdown spectroscopy system LIBS-6 (Applied Photonics, Skipton North Yorkshire, UK). It consists of an integrated Q-switched Quantel Ultra Nd-YAG laser (Quantel, Les Ulis, France) working at λ = 1 064 nm emitting a maximum energy equal to 150 mJ due to one laser pulse (6 ns), and an Avaspec-2048-2-USB2 fibre optic Czerny-Turner spectrometer (6-channel) with a CCD detector (Avantes, NS Apeldoorn, The Netherlands). The system was also equipped with a camera enabled to observe the analyzed object and a movable sample table. The LIBS-6 system was operated by LIBSoft V6.0.1 software (Applied Photonics, Skipton North Yorkshire, UK). Under normal conditions, because of no moving elements inside, a wavelength calibration of the spectrometer was not required. Every measurement were conducted in air under atmospheric pressure. The spectrum of a new background was collected prior to the analysis of each new sample. Table 1. List of papers and ballpoint pens examined in this study.

Pens Papers
A Q-switch delay time of 165 µs, the integration delay time of 1.27 µs, and integration time of 1.2 ms were used. Samples were analyzed directly without any special preparation. They were placed at the sample table at the focal point of the focusing lens (at a distance determined by the nozzle about 70 mm from the optical head). The diameter of the ablation spot, ranging from 0.6 mm to 1 mm, was dependent on the analyzed material.
The results of LIBS analysis is an emission spectrum-a two-dimensional dependence between the intensity and the wavelength (nm) of electromagnetic radiation emitted by an atom or molecule making a transition from a high energy state to a lower energy state. There are many possible electron transitions for each atom, and each has a specific energy difference. These different transitions, leading to different radiated wavelengths, create a unique emission spectrum. The emission spectra were collected in the UV-vis range (185 nm to 904 nm, spectral resolution 0.1 nm).
In total, 30 classes of experimental data listed above were recorded. For each class, 100 LIBS spectra were acquired, resulting in total in 3000 LIBS emission spectra, each one constituting a single sample. There was one spectrum collected from each line on a piece of paper. A sample spectrum is shown in Figure 1. A detailed version of Figure 1 with descriptions of the most prominent lines is in supplementary materials. The spectrum of each sample is available on our website [38].

Spectral Line Identification
In the spectral line identification process (Table 2), the standard protocol was applied [1]. Different elemental compositions for different types of papers and inks were identified. In order to find which spectral line in a particular spectrum of a paper-ink sample originates from ink, a comparison with a spectrum of pure paper was made. The general observation was that the elemental profile of all papers is in accordance with the information (included in the certificate of analysis) about fillers used in analyzed standard papers. Thus, the identified elements were mainly corresponding to the following fillers: calcium carbonate (CaCO 3 ), kaolinite (Al 4 [Si 4 O 10 ](OH 8 )), titanium dioxide (TiO 2 ), and talc (Mg 3 (OH) 2 Si 4 O 10 ). In spectra obtained for papers containing TiO 2 (L and N) due to the ubiquity of Ti spectral lines, no lines from ink components were observable. The best searching of ink spectral lines of elements originating from inks was in the case of O paper in which kaolinite was the main filler. Concerning ballpoint pen inks, four of them have relatively complicated elemental composition (at least three elements). In one of them (ink T) only one element was found.

Data Analysis
We first performed a visual inspection of data after using PCA. In supplementary materials, we include 14 figures presenting data after PCA for all classes on truncated and normalized data. As can be seen, various types of paper (A, D, L, N, and O) separated well. In contrast, various types of inks (B, R, S, SB, and T) were not distinct from each other. This visual analysis confirms the difficulty of the issue-the automated identification of paper-ink samples.
The analysis of the LIBS spectrum consisted of the following steps: 1. independent preprocessing of LIBS spectra; 2. selection of data for the cross validation and testing sets; 3. data analysis based on computational intelligence methods; 4. evaluation of the results.
Steps of the experiment are presented in the flowchart in Figure 2.

Signal Preprocessing
The initial preprocessing of the raw data was performed to reduce the number of data points corresponding to individual LIBS spectra. This preprocessing step was applied because the spectral range of LIBS is very broad (11,746 data points in this case). From the viewpoint of computational intelligence applied to classification tasks, it contains irrelevant information. Moreover, the removed spectral ranges are also not chemically relevant-there is only noise (the intensities close to zero) and there are no spectral lines. Therefore, these ranges do not provide any relevant information about the elemental composition of the sample. It can be expected that the truncation of the datasets enhances extraction of characteristic features of the distinctive classes (types of material) [39]. Four preprocessing steps were taken under consideration: • data truncation by removing the data points from the beginning (first 3746 data points) and from the end (last 1000 data points) of the LIBS spectrum, which do not contain relevant information (7000 data points are left for further analysis); • normalization of intensity values to the interval [0, 1]; • standardization of intensity values (therefore, the mean value becomes equal to 0 and standard deviation becomes equal to 1).
Four preprocessing ways based on these steps were constructed and evaluated: data truncation and either standardization (also known as standard normal variate) or normalization. These steps are depicted at the bottom of Figure 2.
The Simple Intuitive Language for Experiment Modeling (SILEM) [40] software was used for data preprocessing tasks.

Cross-Validation
The set of 3000 LIBS emission spectra was divided into two subsets [41]: a training subset containing 90% (2700) of spectra and a test subset containing 10% (300) of spectra. Selection of spectra for both subsets was performed in a stratified way, that is the numbers of spectra from different classes had the same proportions in each subset. The cross-validation subset was used for parameter optimization, and the test subset was then used for final evaluation of different types of classifiers with optimized parameters.
The 5-fold stratified cross-validation method was applied to build the training and validation data sets. In each of the five analyzed combinations, a training set of 72 LIBS spectra from each of 30 classes (2160 in total) and a validation set of 18 LIBS spectra from each of 30 classes (540 in total) were used.
Evaluation of the classifiers with optimized parameters was performed separately for each training subset using the test subset. Results of evaluation based on each training subset were obtained by averaging the folds.
Each method is potentially dependent on a set of parameters that are either quantitative of qualitative. They have an influence on the overall performance of the method. These parameters were separately optimized for each machine learning algorithm to receive the lowest number of erroneous classes. The basic categorical parameters for each method were set to constants listed in Table 3. Methods based on an adaptive neuro-fuzzy inference system (ANFIS) [48] and Gaussian process [49] were also tested, but they are computationally intensive and take a long time to receive results even with a reduced data set.
The SILEM software uses the scikit-learn Python library [50] as a source of implementations of employed computational intelligence algorithms.

Evaluation Criteria
Evaluation of classification process was based on methodology described in [51,52]. The 5-fold cross-validation strategy was adopted in our experiment. Five performance parameters (figures of merit), namely accuracy (ACC), sensitivity (SEN), specificity (SPE), mean values, and Cohen's kappa (κ), were calculated. The mean values were calculated to estimate the overall performance of the computational intelligence methods used in this study separately for the task of recognition of each of the different classes of LIBS spectra.
To test whether there are classes which are classified with higher accuracy, specificity, or sensitivity, we calculated for each class S the number of true positives TP(S), false positives FP(S), true negatives TN(S), and false negatives FN(S) separately. Then, the accuracy ACC(S), sensitivity SEN(S), and specificity SPE(S) with respect to class S are defined as averages over five folds of cross-validation (K = 5 is the number of folds): where TP i (S), FP i (S), TN i (S), and FN i (S) are, respectively, the number of true positives, false positives, true negatives, and false negatives for the ith fold of cross-validation with respect to class S, i = 1, 2, . . . , K.
The overall values of ACC, SEN, and SPE for the classification system are the arithmetic means of ACC(S), SEN(S), and SPE(S) over all classes.
To evaluate the degree of discrimination between two different samples discrimination power (also known as discrimination accuracy) coefficient (denoted DP) was proposed [53,54]. DP is the ratio of the number of correctly identified pairs of test samples (identified as from different classes) to the number of all possible pairs of test samples. It can be calculated by the following equation: where • D-the number of differentiated pairs, that is the number of pairs correctly identified as belonging to the same class or correctly identified as belonging to different classes;

Results
An example of a raw spectrum and the output from successive stages of spectrum preprocessing (data truncation and normalization) are shown in Figure 3. The data truncation stage is optional and normalization can be replaced by standardization.   Visualization of data preprocessing for single sample spectrum from Class A. Figure (a) depicts an example of a raw spectrogram with marked truncation points, (b) shows normalized data.
Parameter selection is a key part in reaching the optimal overall performance of a classification system. There are many possible options available, so it is important to analyze the outcomes of experiments for different values of parameters to demonstrate the possibilities of machine learning methods. Particular classification methods depend on the various basic parameters set as listed in Table 3 and some other parameters that were optimized. Optimization of these parameters was performed in two steps. The first step was to find out the general range of values for each parameter for the fine-tuning procedure. Then, the detailed grid search of the selected ranges of parameters was performed.
The classification algorithms, the tuning parameters and range of these parameter values are described below.

•
Decision trees-the range of the number of features to consider when looking for the best split: from 100 to 7000 with a step equal to 100. • Random forest-two parameters were optimized during preprocessing. The number of trees in the forest was optimized in range from 10 to 200 with a step of 10 and from 200 to 1000 with a step equal to 50, and the number of features to consider when looking for the best split was optimized in the same range. MLP-the number of neurons was optimized in a range from 10 to 200 with a step equal to 10. The activation function was selected from "identity," "logistic," "tanh," "relu". The solver for weight optimization was chosen from "lbfgs" (an optimizer in the family of quasi-Newton methods), "sgd" (a stochastic gradient descent), or "adam" (a stochastic gradient-based optimizer proposed in [55]).
These optimal parameters were selected in a series of experiments to maximize the mean ACC, SEN, and SPE values described in Section 2.4.4.
The best results of the experiments with different preprocessing methods and machine learning algorithms based on 5-fold cross-validation (3000 spectra were divided into training set containing 2160 spectra, a validation set with 540 spectra, and test sets containing 300) are presented in Table 4. The preprocessing involves data truncation and standardization or normalization. Generally, the classification results were worse when the preprocessing included normalization instead of standardization. In the table, we show the values of ACC, SEN, SPE (and their mean value denoted MEAN), κ, and DP. Table 5 shows the values of ACC(S), SEN(S), and SPE(S) for all S classes and all classifiers with appropriate preprocessing methods. It can be seen from the table that the most easily distinguishable (highest ACC(S), SPE(S), and SEN(S)) spectra are in six classes: A, A + T, D, N, O, and O + T (ACC(S), SPE(S), and SEN(S) reached the value of 100%) using the random forest method. Additionally, spectra of samples from 2 classes-L and N + T-are well recognized by SVM and MLP methods as well. The least distinguishable classes are L + S and O + B with a value of SEN(S) less than 50% by the random forest method, but the SPE(S) value was over 99%. This means that the spectra of samples from these classes are often assigned to other classes.  Table 5. Results of classification per class.

Discussion
The results of experiments confirm that the computational intelligence methods can be used to analyze LIBS data and obtain accurate classification of paper-ink samples (please see Table 4). We have obtained an accuracy of 99.08% for the best classifier. The results of the present study cannot be directly compared with those obtained in [33,34] focusing on the problem of discrimination between two samples which were not assigned to a priori classes. However, the methods used in [33,34] used Plasus SpecLine 2.13 software for the identification of spectral lines and National Institute of Standards and Technology (NIST) spectral database, gave values of the DP in the range from 85 to 92%. The same standard methods with visual spectra was applied to the data collected in our study gave DP equal to 90.6%. After adopting our algorithms to the problem of discrimination between two different samples, we achieved the DP of 99.0% with SVM and more than 98% with decision trees, random forest, PNN, and MLP classifiers. The limitation of applying our solution to the problem of discrimination is that the classes must be predefined.
The detailed analysis of Table 4 shows that the best method for the class of problems addressed in the present study is the random forest classifier. Parameters for this classifier were optimized using the 5-fold cross-validation using standardized data. The classifier during test over independent test set was able to correctly recognize 1294 cases of LIBS spectra out of 1500 cases and had no erroneous recognitions in training sets (10,800 cases). The cases are summed across all five folds of cross-validation. This classifier achieved an accuracy of 99.08%, a sensitivity of 86.27%, a specificity of 99.53%, a kappa of 85.79% and a DP of 99.08% for the test set.
Taking into account the number of erroneous classifications under different settings, the worst classifier was GRNN.
Two preprocessing methods were applied on the raw data prior to the classification. Hence, we have shown that an appropriate selection of preprocessing procedures influences the performance of classification. However, although in most cases the actual selection of the pair preprocessing-classifier determines the overall performance of the classification method, for some specific selections of classifier the final classification performance depends only marginally on the preprocessing. For this reason, it is not possible to select a single best preprocessing method.
We have taken two other methods into consideration: resampling and wavelet decomposition. The former resamples the input data at given wavelengths of the original sequence. The latter performs one dimensional Haar wavelet decomposition [56]. The initial results were promising, but this method requires more detailed investigation in the future. This technique needs further research to determine which type of wavelet decomposition structure and which parameters will lead to the best results.
The results of the present study show that modern computational intelligence methods can be used for the classification of LIBS spectra into predefined classes, which solves a broad class of problems related to LIBS applications. The main limitation of this study is that with the chosen settings, another broad class of LIBS applications, that is analysis of elemental composition, cannot be solved. The problem of classification addressed in the present study is certainly simpler than the problem of detailed analysis of the composition of materials. However, a successful solution to this problem encourages further research in this area. We suppose that, because the problem of analysis of composition of materials based on a LIBS spectrum may be equivalent to the recognition of a presence of specific spectral lines, this problem can be brought to a series of classification subproblems.

Conclusions
The aim of this study was to solve the difficult problem of distinguishing paper-ink samples. Hence, we have designed a computational intelligence system to solve this problem using LIBS spectra. The difficulty of this problem is caused by the high variability of spectra within a single class of paper-ink samples and strong similarity of spectra from samples of different classes. A machine learning system is developed and based on results presented in Tables 4 and 5, we can conclude that the system performs well. The random forest classifier provided the most accurate classification with an accuracy of 99.08%, a sensitivity of 86.27%, a specificity of 99.53%, a kappa of 85.79% and a DP of 99.08% for the test set.
Our developed system can be used for other applications, e.g., the identification of metal alloys. The main limitation of this method is that we do not list particular elements of the tested material, but it may be considered as future work. Additional feature extraction and classification methods, such as principal component analysis [57], can be applied in the future. Hence, the sensitivity of the described method can also be further improved.