Fast Detection of Heavy Metal Content in Fritillaria thunbergii by Laser-Induced Breakdown Spectroscopy with PSO-BP and SSA-BP Analysis

Fast detection of heavy metals is important to ensure the quality and safety of herbal medicines. In this study, laser-induced breakdown spectroscopy (LIBS) was applied to detect the heavy metal content (Cd, Cu, and Pb) in Fritillaria thunbergii. Quantitative prediction models were established using a back-propagation neural network (BPNN) optimized using the particle swarm optimization (PSO) algorithm and sparrow search algorithm (SSA), called PSO-BP and SSA-BP, respectively. The results revealed that the BPNN models optimized by PSO and SSA had better accuracy than the BPNN model without optimization. The performance evaluation metrics of the PSO-BP and SSA-BP models were similar. However, the SSA-BP model had two advantages: it was faster and had higher prediction accuracy at low concentrations. For the three heavy metals Cd, Cu and Pb, the prediction correlation coefficient (Rp2) values for the SSA-BP model were 0.972, 0.991 and 0.956; the prediction root mean square error (RMSEP) values were 5.553, 7.810 and 12.906 mg/kg; and the prediction relative percent deviation (RPD) values were 6.04, 10.34 and 4.94, respectively. Therefore, LIBS could be considered a constructive tool for the quantification of Cd, Cu and Pb contents in Fritillaria thunbergii.


Introduction
Fritillaria thunbergii (Zhebeimu) is a perennial herb in the lily family [1] that is widely recognized for its medicinal properties; its main production areas are the Zhejiang, Jiangsu and Anhui provinces in China [2]. Chinese herbal medicines have been used in China for thousands of years and have spread worldwide [3]. Although herbal medicines are often considered natural and harmless, there are safety concerns, such as the potential of heavy metal overdose and toxicity [4]. Owing to the growth environment and characteristics of the herb species itself, herbs commonly contain a certain amount of harmful heavy metals, such as Pb, Cd, Hg, As, Cu and Cr. The dry bulbs of Fritillaria thunbergii are widely used as expectorants and antitussives [5], and their safety has become a matter of considerable concern. The unique biological properties of Fritillaria lead to the accumulation of heavy metals, especially Cd, when grown in acidic soils. Rapid determination of Cd, Cu and Pb contents is important for evaluating the safety and quality of Fritillaria thunbergii.
Laser-induced breakdown spectroscopy (LIBS) is an atomic emission spectroscopic technique that uses a laser as the excitation source [6]. Compared with traditional analytical techniques, such as inductively coupled plasma mass spectrometry (ICP-MS) [7], atomic absorption spectrometry (AAS) [8], and inductively coupled plasma atomic emission

Determination of the Best Preprocessing Algorithm
The spectral profile is shown in Figure 1, and the raw spectral data contain 1024 variables from 210 to 231 nm. Based on the National Institute of Standards and Technology (NIST) Atomic Spectroscopy Database (ASD), the characteristic spectral lines of heavy metals (Cd, Cu and Pb) were identified in this wavelength range as Cd II (214.44 nm, 226.50 nm), Cd I (228.80 nm), Cu II (217.94 nm, 222.89 nm, 224.26 nm), Pb I (217.00 nm) and Pb II (220. 35 nm). It was observed that the elements of the Fritillaria thunbergia sample were relatively complex with a heterogeneous distribution of elemental emission spectra. Therefore, LIBS data need to be further analyzed using chemometric methods.
To determine the optimal preprocessing algorithm, different prediction methods were applied to raw LIBS spectral data. PLSR models were developed using different preprocessed spectral data as input variables. The results of the PLSR models based on the raw and preprocessed data are shown in Table 1. For Cd prediction, the optimal performance with an R cv 2 of 0.950 and RMSECV of 7.294 mg/kg was obtained using the spectra preprocessed with MSC. Similarly, for Cu prediction, the optimal performance with an R cv 2 of 0.978 and RMSECV of 10.951 mg/kg was obtained using the MSC preprocessed spectra. For Pb prediction, the optimal performance with an R cv 2 of 0.883 and RMSECV of 20.370 mg/kg was obtained using the spectra preprocessed with SNV. The optimal input variables were employed for further calculations. Rcv 2 of 0.978 and RMSECV of 10.951 mg/kg was obtained using the MSC prepro spectra. For Pb prediction, the optimal performance with an Rcv 2 of 0.883 and RMSE 20.370 mg/kg was obtained using the spectra preprocessed with SNV. The optima variables were employed for further calculations.

Sensitive Variable Selection
The LIBS spectrum collected in this study contained 1024 variables. These va have problems with covariance and high dimensionality, which can lead to redu information and a slow computing speed. To reduce these problems, PCA, SPA and were used for feature selection to achieve dimensionality reduction.
The results of the dimensionality reduction using the PCA algorithm are sho Figure 2. The predictions for Cd and Cu had the same dimensionality reduction b they have the same preprocessing method, and PCA involves unsupervised learni shown in Figure 2a, the first 6 principal components could explain 95.34% of th spectral data. For the Pb prediction, 9 principal components accounted for 95.33% cumulative explained variance, as shown in Figure 2b. In other words, these variab retain more raw data than the others. Therefore, we selected the sensitive variab further modeling based on the first six principal components for Cd and Cu and first nine for Pb.
The feature variables were extracted using SPA and CARS as shown in Fig  Based on the CARS method, 134, 118, and 91 variables were selected for the predic Cd, Cu, and Pb, which represented 13.09%, 11.52%, and 8.89% of the total wavele respectively. The sensitive bands selected by SPA for the predictions of Cd, Cu a were 4, 5 and 3, which were only 0.39%, 0.49% and 0.29% of the full variables, respec

Sensitive Variable Selection
The LIBS spectrum collected in this study contained 1024 variables. These variables have problems with covariance and high dimensionality, which can lead to redundant information and a slow computing speed. To reduce these problems, PCA, SPA and CARS were used for feature selection to achieve dimensionality reduction.
The results of the dimensionality reduction using the PCA algorithm are shown in Figure 2. The predictions for Cd and Cu had the same dimensionality reduction because they have the same preprocessing method, and PCA involves unsupervised learning. As shown in Figure 2a, the first 6 principal components could explain 95.34% of the total spectral data. For the Pb prediction, 9 principal components accounted for 95.33% of the cumulative explained variance, as shown in Figure 2b. In other words, these variables can retain more raw data than the others. Therefore, we selected the sensitive variables for further modeling based on the first six principal components for Cd and Cu and on the first nine for Pb.
The feature variables were extracted using SPA and CARS as shown in Figure 3. Based on the CARS method, 134, 118, and 91 variables were selected for the prediction of Cd, Cu, and Pb, which represented 13.09%, 11.52%, and 8.89% of the total wavelengths, respectively. The sensitive bands selected by SPA for the predictions of Cd, Cu and Pb were 4, 5 and 3, which were only 0.39%, 0.49% and 0.29% of the full variables, respectively.

BPNN Model Prediction Results
To select the best prediction model, BPNN models for predicting the Cd, Cu and Pb concentrations were developed using the characteristic variables selected by the PCA, CARS and SPA algorithms as input variables. The results of each model operation are shown in Table 2.
For the prediction of Cd, SPA and PCA selected significantly fewer variables than CARS. However, the SPA-BP model performed better with Rp 2 = 0.968, RMSEP = 5.951 mg/kg and RPD = 5.74. The PCA-BP exhibited overfitting; therefore, it required optimization. For Cu concentration prediction, the BPNN models based on SPA and PCA selected variables that were only 0.49% and 0.59% of the full spectrum, respectively. Their performance with Rp 2 = 0.943, 0.931, RMSEP = 19.367, 21.311 mg/kg and RPD = 4.27, 3.84, respectively. However, they yielded better results than the variables selected by CARS that accounted for 13.09% of the full spectrum. For the prediction of Pb, although the CARS-BP model achieved superior performance (Rp 2 = 0.950, RMSEP = 13.835 mg/kg and RPD = 4.55), the number of variables filtered out was 91, which was much larger than the 3 and 9 variables selected for the SPA and PCA algorithms, respectively, resulting in a long modeling duration.

BPNN Model Prediction Results
To select the best prediction model, BPNN models for predicting the Cd, Cu and Pb concentrations were developed using the characteristic variables selected by the PCA, CARS and SPA algorithms as input variables. The results of each model operation are shown in Table 2.
For the prediction of Cd, SPA and PCA selected significantly fewer variables than CARS. However, the SPA-BP model performed better with Rp 2 = 0.968, RMSEP = 5.951 mg/kg and RPD = 5.74. The PCA-BP exhibited overfitting; therefore, it required optimization. For Cu concentration prediction, the BPNN models based on SPA and PCA selected variables that were only 0.49% and 0.59% of the full spectrum, respectively. Their performance with Rp 2 = 0.943, 0.931, RMSEP = 19.367, 21.311 mg/kg and RPD = 4.27, 3.84, respectively. However, they yielded better results than the variables selected by CARS that accounted for 13.09% of the full spectrum. For the prediction of Pb, although the CARS-BP model achieved superior performance (Rp 2 = 0.950, RMSEP = 13.835 mg/kg and RPD = 4.55), the number of variables filtered out was 91, which was much larger than the 3 and 9 variables selected for the SPA and PCA algorithms, respectively, resulting in a long modeling duration.

BPNN Model Prediction Results
To select the best prediction model, BPNN models for predicting the Cd, Cu and Pb concentrations were developed using the characteristic variables selected by the PCA, CARS and SPA algorithms as input variables. The results of each model operation are shown in Table 2.
For the prediction of Cd, SPA and PCA selected significantly fewer variables than CARS. However, the SPA-BP model performed better with R p 2 = 0.968, RMSEP = 5.951 mg/kg and RPD = 5.74. The PCA-BP exhibited overfitting; therefore, it required optimization. For Cu concentration prediction, the BPNN models based on SPA and PCA selected variables that were only 0.49% and 0.59% of the full spectrum, respectively. Their performance with R p 2 = 0.943, 0.931, RMSEP = 19.367, 21.311 mg/kg and RPD = 4.27, 3.84, respectively. However, they yielded better results than the variables selected by CARS that accounted for 13.09% of the full spectrum. For the prediction of Pb, although the CARS-BP model achieved superior performance (R p 2 = 0.950, RMSEP = 13.835 mg/kg and RPD = 4.55), the number of variables filtered out was 91, which was much larger than the 3 and 9 variables selected for the SPA and PCA algorithms, respectively, resulting in a long modeling duration. For the three heavy metal predictions, the BPNN models based on the SPA and PCA algorithms were more suitable for the next step of the calculation. This is because these two algorithms greatly simplified the model operations and had a short modeling duration and high accuracy. Because the BPNN can easily fall into local optimal solutions, overfitting and other shortcomings, it can lead to unsatisfactory results, such as the low prediction accuracy of PCA-BP for Cd and Pb concentrations and the relatively poor prediction accuracy of SPA-BP for Pb concentrations. Therefore, further optimization of the BPNN model is required to improve the prediction accuracy.

PSO-BP and SSA-BP Models' Prediction Results
The PSO-BP and SSA-BP models were established based on variables selected by the PCA and SPA methods, respectively, to further explore their ability to predict heavy metal concentrations in Fritillaria thunbergii. The performance of each model is listed in Table 3. It can be seen from the results that the PSO-BP and SSA-BP models had better prediction performance than the raw BPNN model (Table 2). Therefore, the two optimization algorithms enhanced the performance, and the accuracy of the models based on different characteristic variables and predictions of different heavy metals was improved in both cases. The performance evaluation metrics of the PCA-PSO-BP and PCA-SSA-BP models were better than those of SPA-PSO-BP and SPA-SSA-BP (Table 3), indicating that PCA is more suitable for joint modeling with PSO-BP and SSA-BP. In addition, the optimal models for all three heavy metal predictions were based on the variables selected using PCA. For the optimal PCA-PSO-BP models of Cd and Pb, R p 2 = 0.982 and 0.964, RMSEP = 4.461 and 11.791 mg/kg, and RPD = 7.72 and 5.47. For the prediction of Cu, the best result was obtained by the PCA-SSA-BP model, and the optimal model had an R p 2 of 0.991, RMSEP of 7.810 mg/kg and RPD of 10.34. Scatter plots of the models with the best results in the prediction set for the prediction of Cd, Cu and Pb are shown in Figure 4.

Comparison with a Conventional Method (PLSR)
We compared the PSO-BP and SSA-BP methods with the commonly used spectral data modeling method, PLSR. Each method used the same data segmentation, preprocessing, and sensitive variable selection methods for modeling, and the results are shown in Table 4. The PLSR method achieved the best results for Cu estimation, as indicated by the following statistical results: Rp 2 = 0.951, RMSEP = 17.996 mg/kg and RPD = 4.50. However, compared with the results in Table 3, the prediction accuracy of the PLSR model for Cu was not as good as that of the PSO-BP and SSA-BP models. For Cd and Pb concentration predictions, the prediction results of PLSR were far inferior to those of the two optimization algorithms with the BPNN. The prediction results of the three models were compared to better visualize the effectiveness of the SSA-BP and PSO-BP prediction models. Figure 5 shows the relative errors between the predicted and true values. For the prediction of Cd and Cu concentrations, the relative errors of both SSA-BP and PSO-BP were lower than those of PLSR at high concentrations, with a better prediction performance of SSA-BP at low concentrations and relatively larger errors of both PSO-BP and PLSR. Regarding the prediction of Pb concentrations, SSA-BP and PSO-BP provided better predictions than PLSR at both low and high concentrations, but the effect of PSO-BP was not obvious. In summary, the two optimization models, PSO-BP and SSA-BP, were found to be more accurate than the commonly used linear PLSR model for the prediction of heavy metals in Fritillaria thunbergii.

Comparison with a Conventional Method (PLSR)
We compared the PSO-BP and SSA-BP methods with the commonly used spectral data modeling method, PLSR. Each method used the same data segmentation, preprocessing, and sensitive variable selection methods for modeling, and the results are shown in Table 4. The PLSR method achieved the best results for Cu estimation, as indicated by the following statistical results: R p 2 = 0.951, RMSEP = 17.996 mg/kg and RPD = 4.50. However, compared with the results in Table 3, the prediction accuracy of the PLSR model for Cu was not as good as that of the PSO-BP and SSA-BP models. For Cd and Pb concentration predictions, the prediction results of PLSR were far inferior to those of the two optimization algorithms with the BPNN. The prediction results of the three models were compared to better visualize the effectiveness of the SSA-BP and PSO-BP prediction models. Figure 5 shows the relative errors between the predicted and true values. For the prediction of Cd and Cu concentrations, the relative errors of both SSA-BP and PSO-BP were lower than those of PLSR at high concentrations, with a better prediction performance of SSA-BP at low concentrations and relatively larger errors of both PSO-BP and PLSR. Regarding the prediction of Pb concentrations, SSA-BP and PSO-BP provided better predictions than PLSR at both low and high concentrations, but the effect of PSO-BP was not obvious. In summary, the two optimization models, PSO-BP and SSA-BP, were found to be more accurate than the commonly used linear PLSR model for the prediction of heavy metals in Fritillaria thunbergii.

Comparison of PSO-BP and SSA-BP
As seen in Table 3, the performances of both the PSO-BP and SSA-BP models are good for the prediction of the 3 heavy metals--the values of Rp 2 and RMSEP are very similar, and RPD values are greater than 3, indicating that both models are valid and stable. However, the modeling duration of the models differed (the average value of 10 operation times was used as the final reference basis)--the modeling duration of the PSO-BP model was approximately 3 times longer than that of the SSA-BP model; thus, SSA-BP significantly outperformed PSO-BP in terms of computing speed. Figures 6 and 7 show the relative errors between the first two sets of data predicted by the BPNN, PSO-BP and SSA-BP models and the true values. For the prediction of low concentrations, the relative errors of the SSA-BP model are all lower than those of the BPNN model, whereas the PSO-BP model is not stable for the prediction of low concentrations and is sometimes less accurate than the BPNN model. The relative errors of the SSA-BP model are lower in most cases compared with the PSO-BP model; therefore, the prediction results of the SSA-BP model are more informative in practical use.
In summary, the SSA-BP model was used to predict the concentrations of three heavy metals (Cd, Cu and Pb) in Fritillaria thunbergii, and a stable, accurate and practical reference result was obtained in a short period of time.

Comparison of PSO-BP and SSA-BP
As seen in Table 3, the performances of both the PSO-BP and SSA-BP models are good for the prediction of the 3 heavy metals-the values of R p 2 and RMSEP are very similar, and RPD values are greater than 3, indicating that both models are valid and stable. However, the modeling duration of the models differed (the average value of 10 operation times was used as the final reference basis)-the modeling duration of the PSO-BP model was approximately 3 times longer than that of the SSA-BP model; thus, SSA-BP significantly outperformed PSO-BP in terms of computing speed. Figures 6 and 7 show the relative errors between the first two sets of data predicted by the BPNN, PSO-BP and SSA-BP models and the true values. For the prediction of low concentrations, the relative errors of the SSA-BP model are all lower than those of the BPNN model, whereas the PSO-BP model is not stable for the prediction of low concentrations and is sometimes less accurate than the BPNN model. The relative errors of the SSA-BP model are lower in most cases compared with the PSO-BP model; therefore, the prediction results of the SSA-BP model are more informative in practical use.
In summary, the SSA-BP model was used to predict the concentrations of three heavy metals (Cd, Cu and Pb) in Fritillaria thunbergii, and a stable, accurate and practical reference result was obtained in a short period of time.

Comparison of PSO-BP and SSA-BP
As seen in Table 3, the performances of both the PSO-BP and SSA-BP models are good for the prediction of the 3 heavy metals--the values of Rp 2 and RMSEP are very similar, and RPD values are greater than 3, indicating that both models are valid and stable. However, the modeling duration of the models differed (the average value of 10 operation times was used as the final reference basis)--the modeling duration of the PSO-BP model was approximately 3 times longer than that of the SSA-BP model; thus, SSA-BP significantly outperformed PSO-BP in terms of computing speed. Figures 6 and 7 show the relative errors between the first two sets of data predicted by the BPNN, PSO-BP and SSA-BP models and the true values. For the prediction of low concentrations, the relative errors of the SSA-BP model are all lower than those of the BPNN model, whereas the PSO-BP model is not stable for the prediction of low concentrations and is sometimes less accurate than the BPNN model. The relative errors of the SSA-BP model are lower in most cases compared with the PSO-BP model; therefore, the prediction results of the SSA-BP model are more informative in practical use.
In summary, the SSA-BP model was used to predict the concentrations of three heavy metals (Cd, Cu and Pb) in Fritillaria thunbergii, and a stable, accurate and practical reference result was obtained in a short period of time.

Fritillaria thunbergii Material and Sampling
Eight different brands of Fritillaria thunbergii commonly found on the market were purchased and used in this experiment (Table 5). Samples with different Cd, Cu and Pb contents were prepared in the laboratory by treatment with cadmium nitrate, copper nitrate and lead nitrate solutions.
Cadmium nitrate, copper nitrate and lead nitrate powders were used to prepare a solution of 0.01 mol/L each. The 8 types of Fritillaria thunbergii were crushed using a highspeed crusher and each type of Fritillaria thunbergii was equally divided into 8 groups and coded. A concentration gradient (for Cd: 0, 5, 10, 20, 25, 50, 80 and 100 mg/kg; for Cu: 0, 20, 40, 60, 80, 100, 200 and 300 mg/kg; and for Pb: 0, 5, 20, 40, 60, 80, 100 and 200 mg/kg) was predetermined between each group. The purpose was to produce samples whose true values could also be re-examined by 1-8 ICP-MS. Each group contained 5 g of Fritillaria thunbergii powder; 1 group was the control group, and the rest were the treatment groups. Different amounts of the solution were added to each group (for Cd cadmium nitrate solution: 22 . All samples were dried in an oven at 60 °C and crushed separately using a pulverizer (AQ-180E, Nail, Ningbo, China) to obtain a homogeneous mixture of each group of samples and to ensure consistency in the content of each heavy metal. The ground samples were formed into 1.5 cm diameter pellets using a tablet compressor. For each group, three parallel samples were prepared, and a total of 192 samples were obtained.

Fritillaria thunbergii Material and Sampling
Eight different brands of Fritillaria thunbergii commonly found on the market were purchased and used in this experiment (Table 5). Samples with different Cd, Cu and Pb contents were prepared in the laboratory by treatment with cadmium nitrate, copper nitrate and lead nitrate solutions. Cadmium nitrate, copper nitrate and lead nitrate powders were used to prepare a solution of 0.01 mol/L each. The 8 types of Fritillaria thunbergii were crushed using a high-speed crusher and each type of Fritillaria thunbergii was equally divided into 8 groups and coded. A concentration gradient (for Cd: 0, 5, 10, 20, 25, 50, 80 and 100 mg/kg; for Cu: 0, 20, 40, 60, 80, 100, 200 and 300 mg/kg; and for Pb: 0, 5, 20, 40, 60, 80, 100 and 200 mg/kg) was predetermined between each group. The purpose was to produce samples whose true values could also be re-examined by 1-8 ICP-MS. Each group contained 5 g of Fritillaria thunbergii powder; 1 group was the control group, and the rest were the treatment groups.  482.63 µL). All samples were dried in an oven at 60 • C and crushed separately using a pulverizer (AQ-180E, Nail, Ningbo, China) to obtain a homogeneous mixture of each group of samples and to ensure consistency in the content of each heavy metal. The ground samples were formed into 1.5 cm diameter pellets using a tablet compressor. For each group, three parallel samples were prepared, and a total of 192 samples were obtained.
In our experiment, the Nd: YAG laser was used as an excitation source. The laser beam was focused 2 mm below the samples with a 100 mm focal length optical lens, forming an ablation spot. To increase the signal-to-noise ratio, the following parameters were optimized: λ = 532 nm, energy = 40 mJ, delay = 2 µs, and gate width = 10 µs. The repetition rate was 1 Hz. Each spectrum was cumulatively acquired, and three laser pulses were taken. In total, 7 LIBS spectra were acquired in 7 different positions, yielding 21 (3 × 7) spectra. To reduce the fluctuation between the laser points, the first of the 7 points was removed, and the average of 18 (3 × 6) spectra was recorded as the LIBS data of the sample.

Reference Method for Heavy Metal Content Determination
The concentrations of Cd, Cu and Pb in all Fritillaria thunbergii samples, as determined by inductively coupled plasma mass spectrometry (ICP-MS) analysis (ELAN DRC-e, PerkinElmer, Waltham, MA, USA), were used as reference values. Ginseng (GBW10020) and French beans (GBW10021) were used to assure the accuracy of the results. The results of the Cd, Cu and Pb content determination are shown in Table 6. Table 6. Inductively coupled plasma mass spectrometry (ICP-MS) measurement of heavy metal (Cd, Cu, and Pb) content in Fritillaria thunbergii samples with eight different treatments (range: the range of the heavy metal content, mean ± SD: the mean and standard deviation of the heavy metal content).

Cd
Cu Pb

Sample Division and Spectral Preprocessing
The 8 concentration gradients of Fritillaria thunbergii samples in this experiment were numbered 1-8, and each sample contained 3 sets of experimental data, for a total of 192 spectral sets of data. We used 2 randomly selected sets of data from each sample as the training set for model calculation, and the remaining set of data was used as the test set for a total of 128 sets of training data and 64 sets of test data.
Before the calibration process, it was necessary to employ preprocessing to reduce irrelevant information and remove high-frequency noise and baseline variations. Four preprocessing methods, including multiplicative scatter correction (MSC), standard normal variate (SNV), wavelet transform (WT) and Savitzky-Golay smoothing (SG), were applied separately and compared. The coefficient of determination for cross-validation (R cv 2 ) and the root mean squared error cross-validation (RMSECV) were used as evaluation indicators to select the optimum preprocessing strategies.

Sensitive Variable Selection
The raw LIBS spectrum of a Fritillaria thunbergii sample contained large amounts of data, including multiple element spectral lines and unnecessary information, such as redundancy, collinearity, and background information. In some cases, suitable methods can select effective wavelengths containing useful information and reduce the collinearity and high-dimensionality problems of full-spectrum data to reduce the input, develop simpler models, and increase speed. In this study, competitive adaptive reweighted sampling (CARS), successive projection algorithm (SPA), and principal component analysis (PCA) were used to select characteristic wavelengths.
3.6. Discriminant Analysis Method 3.6.1. BPNN A BPNN is self-organizing, self-learning, and self-adaptive, and its principles are simple and easy to implement. BPNNs have been used in several applications. A BPNN is a multilayer feed-forward neural network that includes an input layer, implicit layer, and output layer [27]. The main feature of the BPNN algorithm is that the signal is transmitted forward, and the weights and thresholds of the entire network are continuously adjusted by back propagation to reduce errors in the network output [28].
BPNN can achieve arbitrarily complex nonlinear mapping [29] with a certain degree of robustness and fault tolerance [30]. However, the BPNN algorithm has inherent shortcomings, including unstable training results that tend to fall into local optimum solutions, slow convergence, and data overfitting. Therefore, the initial weights and thresholds of the BPNN algorithm were optimized by introducing the PSO and SSA algorithms with the aim of iterative convergence to the optimal global solution.

PSO-BP
The PSO algorithm is used instead of gradient descent to determine the optimal weights and thresholds of the BPNN [31], which can improve the generalization performance. The PSO algorithm that is combined with the BPNN is called PSO-BP.
PSO, which is a swarm intelligence optimization initially presented by Eberhart and Kennedy (1995) [32], is derived from the predatory behavior of birds in nature. This algorithm is the simplest way for each bird to follow the bird closest to the food and search its surrounding area to obtain food. The basic principle is that each problem solution is considered to be a bird in the search space, called a "particle" whose position, velocity, and fitness values represent the particle characteristics. Assuming that in a d-dimensional space, there is an initialized population of potentially optimal particles, the particles adjust their velocity and position according to their individual and global extremes and update their fitness values to achieve the goal of finding an optimal solution.

SSA-BP
The SSA is a novel algorithm with strong optimization capability and fast convergence. The SSA is used to optimize the initial weights of a BPNN, establishing the SSA-BP model, which not only makes full use of the mapping capability of the BPNN, but also has the fast global convergence and learning capability of the SSA. The SSA algorithm that is combined with the BPNN algorithm is called SSA-BP.
The SSA is a novel smart swarm optimization algorithm based on the foraging and antipredatory behavior of sparrows [33]. Sparrows are classified as discoverers, joiners, and scouts in the foraging process [34]. The basic process of SSA is to initialize the sparrow population; calculate individual fitness values and determine the best and worst suited individuals; update the positions of discoverers, joiners, and scouts in turn; and update through successive iterations until the termination conditions are met.

Model Evaluation and Software
The performance of the models was evaluated using seven parameters: the correlation coefficient for calibration (R c 2 ), cross-validation (R cv 2 ) and prediction (R p 2 ), the root mean square error of calibration (RMSEC), the root mean square error of cross-validation (RMSECV) and the root mean square error of prediction (RMSEP), and relative percent deviation (RPD). A good model should have high R c 2 , R cv 2 and R p 2 values and low RMSEC, RMSECV and RMSEP values. An RPD greater than 1.5 is considered a good prediction; an RPD between 2.0 and 2.5 indicates that it is a satisfactory prediction model; and an RPD greater than 3.0 is considered a valid prediction model. Spectral data preprocessing, extraction, and calculation of the SPA, CARS, PCA and PLSR algorithms were performed using Python3.7.3. The BPNN, PSO-BP, and SSA-BP were constructed using MATLAB R2019b (Math Works, Natick, MA, USA).

Conclusions
This study demonstrated that LIBS can be used to predict the Cd, Cu and Pb contents in Fritillaria thunbergii. For the estimation of the content of heavy metals, optimal spectral estimation models were established using BPNN, PSO-BP, SSA-BP and PLSR models. The models established by the PSO-BP and SSA-BP models were better and had higher precision than the other models. The SSA-BP model is superior to the PSO-BP model because it runs faster, has higher prediction accuracy at low concentrations, and is more practical to apply. Data Availability Statement: The data, models, and code generated and/or employed during the study are available from the corresponding author upon request.