A Back Propagation Neural Network Model Optimized by Mind Evolutionary Algorithm for Estimating Cd, Cr, and Pb Concentrations in Soils Using Vis-NIR Diffuse Reflectance Spectroscopy

Abstract

Visible and near infrared spectroscopy is an effective method for monitoring the content of heavy metals in soil. However, due to the difference between polluted soil with phytoremediation and without phytoremediation, the common estimation model cannot meet accuracy requirements. To solve this problem, combined with an ecological restoration experiment for soil contamination using the plant Neyraudia reynaudiana, this study explored the feasibility of using a hyperspectral technology to estimate the heavy metal content (Cd, Cr, and Pb) of soil under phytoremediation. A total of 108 surface soil samples (from depths of 0–20 cm) were collected. Inversion models were established using partial least squares regression (PLSR) and the back propagation neural network optimized by a mind evolutionary algorithm (MEA-BPNN). The results revealed that: (1) modeling with derivative-transformed spectra can effectively enhance the correlation between soil spectral reflectance and heavy metal content. (2) Compared with the BP neural network model, the estimation accuracy (R²) was improved from 0.728, 0.737, and 0.675 to 0.873, 0.884, and 0.857 using the MEA-BP neural network model. The residual prediction deviation (RPD) values for the three heavy metals Cd, Cr, and Pb using the MEA-BPNN model were 2.114, 3.000, and 2.560, respectively. Among them, the estimated model of Cd was an excellent prediction. (3) Compared with PLSR, the model prediction results established by the MEA-BP neural network had higher estimation accuracy. In summary, the use of diffuse reflectance spectroscopy to predict heavy metal content provides a theoretical basis for further study of the large-scale monitoring of soil heavy-metal pollution and its remediation evaluation in the polluted area, which is of great significance.

1. Introduction

Soil provides basic habitats and natural resources for the survival of humans, plants, and animals [1,2]. However, as a result of metal mining and the associated dumping of tailings, soil contamination by heavy metals has become a serious problem [3,4]. Heavy metals become enriched in soil, leading to deleterious effects on food safety via the food chain [5], as well as indirect concertation and other harmful impacts on human health [6]. At present, many reclamation methods have been applied in areas affected by mining, of which phytoremediation is a common method for the remediation of heavy-metal-contaminated soil. Some research has shown that local plants could achieve ecological rehabilitation of soils contaminated by heavy metals [7,8]. Therefore, the rapid determination of heavy metals in soil is crucial for the monitoring and assessment of remediation.

In terms of the detection of heavy metal concentration, the traditional method relies on laboratory analysis and field sampling to determine the concentration and distribution of heavy metals. This is a time-consuming and expensive process [9,10,11,12]. In order to quickly determine heavy metal contamination in soil, the spectroscopic method has been widely used in soil monitoring [13,14]. Compared with traditional methods, the spectroscopic method has more advantages, such as less damage to the ecological environment, lower cost, and accurate reflection of the biological, chemical, and physical characteristics of objects [15]. This technology is convenient for obtaining the content of heavy metals and its dynamic evolution in soil [15,16,17,18]. Many researchers have found that soil organic matter (SOM), iron oxide, water, and other active ingredients of soil have a good correlation with visible-near-infrared spectroscopy [19] and have used reflectance spectroscopy to effectively estimate soil properties. Based on the relationship between heavy metals and soil active components, soil heavy metal concentrations can be indirectly predicted by spectroscopy [20,21]. Based on the significant relationship between soil composition and spectral characteristics, Kooistra successfully predicted the Cd and Zn contents in the soil of the Rhine Basin [22].

Recently, the study of estimating heavy metal concentration by visible and near-infrared spectroscopy (VNIRS) has mainly focused on soils with medium and high heavy metal concentrations [23,24,25,26]. For the detection of changes in pollutant content, more studies focused on measuring vegetation reflectance as a proxy to soil contamination [27,28,29]. However, in the process of soil reclamation, vegetation is sometimes absent. Therefore, in this study, we still directly used the soil reflectance to evaluate and monitor soil remediation effects. Because the content of heavy metals in soil under remediation varies widely and considerably, the relationship between heavy metal content and spectral reflectance is ambiguous, and it requires a more suitable model for estimation. In terms of the methods developed to detect heavy metals based on hyperspectral spectroscopy, various models have been established in combination with statistical theory to predict the concentration of heavy metals in soil. Partial least squares regression (PLSR) [30], stepwise multiple linear regression (SMLR) [31], support vector machine (SVM) [32], the random forest method (RF) [33], and the back propagation neural network (BPNN) [34,35,36,37] are some commonly used technical methods. A promising approach is to use the mind evolutionary algorithm back propagation neural network (MEA-BPNN) to estimate heavy metal concentrations in soil. In 1998, Sun Chengyi and others proposed the multilateral environmental agreement [38] to explain the principles of the MEA. The MEA-BPNN model uses an evolutionary algorithm to optimize back propagation (BP) neural networks, combining BPNN’s local search capabilities with MEA’s global search capabilities. According to our literature review, this model has often been used for the prediction of wave height [39], weather, and wind, and achieved good predictions. The results show that the MEA-BP takes less time and has higher prediction accuracy. However, the capacity of this method in the prediction of soil heavy-metal (Cd, Cr and Pb) concentration still lacks assessment.

This study analyzed the reflectance spectra of heavy-metal-contaminated soils that were undergoing phytoremediation, and investigated the feasibility of employing the MEA-BPNN model to predict the concentrations of Cd, Cr, and Pb in these soils. We compared the accuracy of the BPNN and MEA-BPNN models with the PLSR model commonly used in soil heavy metal inversion in order to determine the best estimation model. By exploring the more convenient and effective detection method of heavy metal content in soil, this study provides a basis for monitoring and evaluating the restoration of contaminated soil rapidly.

2. Materials and Methods

2.1. Experimental Design and Samples

The soil heavy metal remediation experiment was based on the heavy metal pollution area of the dumping field in the vanadium-titanium iron ore district of Panzhihua City, Sichuan Province, China. The study area has the world’s third-largest vanadium-titanium magnetite resource and has been plagued by potential and toxic metal pollution for many years. Solid waste contains many harmful substances such as Zn, Cd, Pb, Hg, As, V, Co, Ti, Cu, Ni, Cr, Mn, and so on. In our experiment, Cd, Cr, and Pb—the moderate and severe pollutants in the soil of this area—were used as stress metals. At present, the Panzhihua mining area is undergoing ecological restoration of the soil environment that has been polluted by mining. We used the local reed Neyraudia reynaudiana in the area as our experimental restoration plant. Based on these conditions, this study simulated the heavy metal pollution, ecological remediation plants, and climate in the Panzhihua mining area in order to provide a basis for the remediation and monitoring of heavy metals in soil.

The medium in each culture pot was a 3:1 ratio of clean river sand to loess soil. The metal-adsorbing N. reynaudiana reeds were planted in these pots. CdCl₂·2.5H₂O, K₂Cr₂O₇, and Pb(NO₃)₂ solutions containing Cd²⁺, Cr⁶⁺, and Pb²⁺ were selected as heavy metal stress sources. As shown in Table 1, based on the determination of soil heavy-metal content in the Panzhihua vanadium-titanium ore district [40], along with the Risk Control Standard for Soil Contamination (GB 15618-2018), 6 concentration gradients of Cd, Cr, and Pb in the simulated polluted soil samples were examined. Three sets of parallel experiments were set up for each gradient, for a total of 54 potted plants (as listed in Table 2). The pH of the soil was maintained at 6.5 (± 0.1) and continuously ventilated for 24 h.

The experiment began in a greenhouse in May 2018. After 14 days of slow seedling development, the plants with consistent growth were selected and transplanted into potted soil with heavy metal stress. During the experiment, the greenhouse was kept ventilated, and the potted plants were irrigated regularly (every 2 days). Soil samples were collected at different stages of the growth of Neyraudia reynaudiana. The specific experimental process is shown in Figure 1.

In this study, six stages of data were collected at 30, 60, 90, 120, 150, and 180 days after heavy metal stress. From the pots containing soils with different concentrations of heavy metals, 36 samples of surface soil under Cd, Cr, and Pb stress were collected, totaling 108 samples. The soil samples were placed in sealed plastic bags and labeled. In order to eliminate the interference of soil moisture on the spectrum, all soil samples were naturally air-dried at room temperature after removing small stones and plant roots. The soil was ground in an agate bowl and passed through a 100-mesh nylon sieve (sieve size: 0.150 mm) for laboratory analysis. The soil samples were divided in half and sent to the laboratory for hyperspectral measurements and soil property analysis.

2.2. Chemical Analysis

Data on the heavy metal content in the soil samples were obtained from laboratory tests. After sample preparation, the contents of Cd, Cr, and Pb in the soil samples were measured with an inductively coupled plasma mass spectrometry (ICP-MS) instrument according to the standard HJ766-2015 method [26].

2.3. Spectral Measurement and Preprocessing

Soil spectra were measured with an ASD Field Spec^®3 field spectrometer (Analytical Spectral Devices, Inc., Colorado, USA), which covers a spectral range of 350–2500 nm. Spectral measurements were collected in a dark room. A 50 W halogen lamp was installed on a tripod as the light source and placed 30 cm from each soil sample at a 15° angle to the vertical. An optical fiber probe was located directly above each sample and 10 cm from it to obtain the spectral information of the soil. Before spectral measurement, the instrument was preheated for 30 min and calibrated with a white spectral reference panel. The soil samples were shaken evenly after drying and grinding to make the soil more homogeneous. Then the soil sample was uniformly tiled in glassware with a diameter of 10 cm to make its surface smooth. Each sample was tested five times. ViewSpecPro, a type of spectral processing software, was used to view and eliminate the spectral curves with obvious errors, and the arithmetic average value was obtained as the final spectral reflectance data.

Two edge bands of 350–399 nm and 2400–2500 nm with obvious noise were removed from all sample spectral data. During the acquisition process, spectral data will be affected by random factors, resulting in noise, making it necessary to process the signal and smooth the reflectance spectrum. In this study, multiplicative scatter correction (MSC) was used to eliminate the scattering [41] and then spectra were subjected to the Savitzky–Golay smoothing method with a second-order polynomial fit and a window size of 7 data points. The first derivative (FD), the second derivative (SD), and the reciprocal logarithm (RL) were used to transform the original spectrum after pretreatment with SG smoothing and MSC. Derivative transformations are less sensitive to noise effects, eliminate interference from baselines and other backgrounds, and minimize the baseline drift of the original soil spectrum [42]. The reciprocal logarithm not only enhances the difference in the visible spectrum, but also reduces the influence of multiplication factors under illumination conditions [36]. The spectral preprocessing was performed using ViewSpecPro and MATLAB.

2.4. Model Calculation and Accuracy Evaluation

2.4.1. The Mind Evolutionary Algorithm Back Propagation Neural Network

The BP neural network (BPNN) is a multi-layer feedforward network with error back propagation, primarily consisting of an input layer, a hidden layer, and an output layer. The number of hidden layers is uncertain and can be 1 or more. The BPNN exhibits a good approximation ability for any nonlinear continuous function and is one of the most widely used neural network models [43]. Its structure is shown in Figure 2.

The mind evolutionary algorithm (MEA) is used to optimize the initial weights and thresholds of the BP neural network, and to solve the BPNN shortcoming of being only partially optimal due to its randomly generated weights and thresholds, which further affects the prediction results. The MEA-BPNN proposes similartaxis and dissimilation strategies, which are coordinated with one another and not opposed. The similartaxis strategy is aimed at one segment of the population, while the dissimilation strategy is aimed at the entire population, such that it can improve the selection efficiency, yield better prediction ability, and work more quickly [39].

The specific steps of the MEA are as follows [44]:

(1) Network initialization

The BP neural network in this study has a three-layer network structure. The nodes in the hidden layer are determined by the empirical model. The hidden layer defaults to the ‘tansig’ function, and the output layer defaults to the ‘purelin’ function. The learning rate was 0.1, the training frequency was 1000, and the expected error was 0.0001. The initial weights $w$ and thresholds $\theta$ of the BP neural network are optimized by the mind evolutionary algorithm, which can significantly reduce the number of network iterations and model errors. Based on the BP network structure, the solution space is mapped to the encoding space, and each solution corresponds to the encoding, where the encoding length s is:

(1)$s=s_1\times s_2+s_2\times s_3+s_2+s_{3,}$

(2)$s_3=\sqrt{s_1+s_2}+a$

where $s_1, s_2, s_3$ are the number of input, output, and hidden layer nodes, $\mathrm{a}\in \left[1,10\right],$ respectively.

(2) Generation of individuals and subpopulations

The size of the mind evolutionary algorithm is set to 8, the number of best subpopulations and temporary subpopulations are each set to 5, and the number of iterations is set to 20. The strategy of the MEA-BPNN is to generate individuals for a certain number (N), sort the scores, and search for the superior individual, that is, the individual with the highest score. It also searches for other temporary individuals, that is, those individuals that have the highest score. If M and T are respectively the number of superior and temporary subpopulations, then the size of each subpopulation m is as follows:

(3)$\mathrm{m}=\frac{\mathrm{N}}{\mathrm{M}+\mathrm{T}},$

The derivative of the mean square error of the training set is selected as the score function val of each individual and subpopulation:

(4)$\mathrm{val}=\frac{1}{\mathrm{SE}},$

(5)$\mathrm{SE}=\mathrm{mse}\left(\mathrm{T}-\mathrm{A}2\right),$

where val is the output value of the output layer after each iteration of the A2 bit; T is the expected output; and SE is the mean square error.

(3) Similartaxis operation

Constantly comparing the scores of each individual in a subpopulation, and searching for the individual with the highest score, the similartaxis operation then selects this optimal individual and uses it as the center for the generation of a new subpopulation. This process continues until a new optimal individual is no longer generated.

(4) Dissimilation operation

Comparing the score of the superior subpopulations with those of the temporary subpopulations, if there is a temporary subpopulation with a score that is higher than that of a superior subpopulation, the superior subpopulation is released, and is replaced by this temporary subpopulation. If the requirements are not met, the model returns to the network and repeats steps (2) and (3). After the evolution is complete, the optimal individual is decoded as the initial weight and threshold of the network structure.

The specific process of the MEA-BPNN model is shown in Figure 3.

2.4.2. Partial Least Squares Regression

PLSR is now the most frequently used method for the estimation of soil heavy-metal concentrations with VNIRS (Table 2). PLSR has the same general structure as principal component regression (PCR), but takes the dependent variables into account when calculating the principal components [45]. Therefore, PLSR is able to retain groups of collinear variables, and handle multivariate statistical analysis as well as situations in which the number of variables considerably exceeds the number of available samples. The optimum number of latent variables (LVs) was determined using a leave-one-out cross-validation procedure.

2.4.3. Model Calibration and Validation

First, the 36 soil samples contaminated by each heavy metal were sorted in ascending order of concentration. Then, they were divided into 6 groups according to the concentration. The samples in each group were randomly divided into a calibration set and a validation set at a ratio of 2:1. Thus, the calibration dataset included 24 samples, and the validation set comprised the other 12 samples.

Regression models need to be evaluated prior to their application. The R² value, root mean square error (RMSE) and residual prediction deviation (RPD) are the important statistical parameters used to evaluate the estimation accuracies of calibration models. The evaluation criteria are: RPD < 1.5 and R² < 0.50 indicate unsuccessful models, which are not recommended; 1.5 ≤ RPD < 2.0 and 0.50 ≤ R² < 0.65 indicate fair models; 2.0 ≤ RPD < 2.5 and 0.65 ≤ R² < 0.81 indicate good models; 2.5 ≤ RPD < 3.0 and 0.81 ≤ R² < 0.90 indicate very good quantitative models; and RPD ≥ 3.0 and R² ≥ 0.90 indicate excellent models [15]. The smaller the RMSE value is, the higher the prediction accuracy of the model is. The calculation formulas for the three indicators are as follows:

(6)$r^2=1-\frac{{\displaystyle \sum }{\left(y_m-y_p\right)}^2}{{\displaystyle \sum }{\left(y_m-\overline{y}\right)}^2},$

(7)$\mathrm{RMSE}=\sqrt{\frac{{\displaystyle \sum }{\left(y_m-y_p\right)}^2}{n}},$

(8)$\mathrm{RPD}=\frac{\mathrm{SD}}{\mathrm{RMSE}},$

where y_m is the measured value of heavy metal content; $\overline{y}$ is the average value of the measured value of heavy metal content; y_p is the estimated value of heavy metal content; n is the number of samples; and SD is the standard deviation of the measured value of the test set.

3. Results

3.1. Analysis of the Effect of Rehabilitation Experiments on Heavy Metal Pollution

The rehabilitation experiments were designed to investigate and analyze the remediation effect of the Neyraudia reynaudiana reeds and the changes of heavy metals in the soil. To accomplish this, the heavy metal contents of soils during the six growth stages of Neyraudia reynaudiana were collected (Figure 4). Because of the adsorption of heavy metal ions by Neyraudia reynaudiana, we found that the heavy metal content in the soil decreased with time, and the concentrations of Cr and Pb changed significantly. After 60 days of treatment with 1.0 mg/kg Cd and 1000 mg/kg Pb, the metal content in the soil reached the lowest values, 0.101 mg/kg and 17.430 mg/kg, respectively, while after 90 days of treatment with 1000 mg/kg Cr, the Cr content reached a minimum of 8.450 mg/kg. At that time, the maximum absorption amounts of Cd, Cr, and Pb in the soil were 29.199 mg/kg, 804.550 mg/kg, and 395.457 mg/kg, respectively. After 180 days, the content of heavy metals in the soil was lower than the screening value of soil pollution risk for agricultural land. At this time, the soil could be used for agricultural production.

Mathematical statistical analysis was carried out on the measured values of soil heavy-metal content, including the extreme values, standard deviation, skewness, kurtosis, and coefficient of variation. It can be seen from Table 3 that the variation ranges of Cd, Cr, and Pb in the soil samples were 0.13–29.30 mg/kg, 6.90–813.00 mg/kg, and 0.84–397.20 mg/kg, respectively. The coefficients of variation of the three heavy metals were relatively high, all >1.00. This is because the concentration gradients of the experimental design were quite different, resulting in correspondingly large variation coefficients. The skewness coefficients and kurtosis coefficients of the three heavy metals were all >0, indicating that their contents were all right-skewed distributions. Among them, the peak coefficient and skewness coefficient of Cd were the largest, indicating that the deviation of Cd content from a normal distribution was the greatest.

3.2. Analysis of Spectral Characteristics of Heavy metal Soil

Differences in soil spectra are mainly indicated by differences in reflectance. The characteristic absorption bands exhibited by the different soil types differed in location and strength. In the experiment, the spectral intensity of soil under heavy metal Cr stress was lower than that of the other two heavy metal stress soils. This may be due to the dark soil color, with the reflectance absorbing additional light energy [46]. Figure 5 shows the reflectance for the soil with different heavy metal contents in the experiment.

Generally speaking, the original spectral curves of the three heavy metal soils were relatively smooth. As shown in Figure 5, the spectral reflectance increased gently with wavelength from 400–760 nm (the visible wavelength band). The upward trend of reflectance was obvious between 850–900 nm. The reflectance reached the first reflection peak between 1350 nm and 1400 nm. Three more significant absorption peaks appeared at 1440 nm, 1900 nm, and 2200 nm, which were due to soil moisture. At 1400, 1800, 2000, and 2250 nm, prominent reflection peaks formed, while the rest of the bands were steady. The spectral reflectance range of Cd and Pb increased with the change of their heavy metal content, while the spectral reflectance of Cr exhibited a smaller response to the change of its content. In addition, we can also find the spectral reflectance of Cd- and Pb-contaminated soil generally decreases with time, and the spectral reflectance of Cr increased with time. Combining with Figure 4, heavy metal content in soil decreased over time, so, we could roughly judge that the spectral reflectance of soil was positively correlated with the content of Cd and Pb, and inversely correlated with the content of Cr.

After the differential transformation, the spectral characteristics of the three heavy metals were similar (Figure 6). Several absorption features were highlighted across all heavy metal classes, such as the significant absorption peaks and valleys at approximately 950 nm, 1300 nm and 1900 nm. The spectral reflectance after the second derivative transform was closer to 0, indicating that the baseline drift and mixed overlapping peaks had been gradually removed, and the spectral intensity had also been reduced. The reciprocal logarithmic transformation can stretch and enhance the reflectance of the original spectrum, but the information enhancement effect was less significant than that of the derivative transformation.

In this study, the Spearman correlation analysis was carried out with the soil heavy-metal contents and the spectral reflectance after the first-order differential transformation, second-order differential transformation and logarithmic transformation of the reciprocal. Compared with the original reflectance spectra, the correlation values of the derivative transformed reflectance spectra were significantly improved. As shown in Table 4, the first derivative effect on the reflection spectrum is best, followed by the second derivative. The correlation between the original spectral reflectance with Cd and Pb content was weak, and the correlation with Cr was more significant. The reflectance spectra of soil could be distinguished from the absorption characteristics after derivative transformation, which enhanced the correlation between heavy metal elements and reflectance. This is because the derivative transform determines the bending point of the spectrum and the wavelength of the maximum and minimum values by mathematical simulation of the reflection spectrum, and also sharpens peaks and slightly resolves overlapping bands. However, there was little difference between the spectral reflectance after logarithmic transformation of the reciprocal and the original spectral reflectance in correlation with heavy metal content.

The 2150–2170 nm band was most sensitive to Cd, 420–433 nm was most sensitive to Cr, and the spectra at 534 and 1333 nm were most sensitive to changes in Pb content. It can be seen from Figure 7 that the absolute value of the correlation coefficients between the spectral reflectance (the FD spectral reflectance is in the wavelength range of 420–495 nm and 1310–1330 nm, the SD spectral reflectance is in the range of 420–981 nm and 2110–2170 nm) and the soil Cd content were greater than 0.430 (significance level of p = 0.01). The original spectral reflectance in the wavelength range of 2000–2010 nm and 1899–1900 nm, the FD spectral reflectance in the range of 433–526 nm and 2250–2340 nm, the SD spectral reflectance in the range of 408–979 nm and 2290–2340 nm, and the RL spectral reflectance in the range of 2000–2020 nm and the soil Cr content were correlated significantly. The position ranges of larger correlation coefficient values with the soil Pb content were 1330–1340 nm and 1106–1130 nm of the FD spectra, 442–719 nm and 1645–1743 nm of the SD spectra. From each of the transformed reflection spectra, we selected 8 bands with the highest absolute value of the correlation coefficient and no collinearity to predict the soil metal contents.

3.3. Model Estimations and Comparisons

The original reflectance spectra were used as the control group. The corresponding PLSR, BPNN, and MEA-BPNN heavy metal estimation models were established using the reflectance spectra after first derivative, second derivative, and reciprocal logarithmic transformations. Finally, R², RMSE, and RPD were used for evaluation.

3.3.1. PLSR Model Prediction Results

The estimation model of three heavy metals was established by the partial least square method. The scatter plot of the measured value and the estimated value is shown in the Figure 8 below. The latent variables (LVs) of the model were determined by leaving-one cross-validation. The LVs of the optimal models of the three metals Cd, Cr, and Pb are 2, 3 and 3. Pb estimation was the best for all PLSR optimal models with FD reflectance, as indicated by the following statistical results: R² = 0.853, and RPD = 2.193. Furthermore, the estimated value of Pb is closer to 1:1.

3.3.2. BPNN Model Prediction Results

The selected characteristic bands were used to predict the heavy metal content using a three-layer BP neural network. The numbers of input and output nodes of the neural network were 8 and 1, respectively. Here we used trial and error method to determine the number of hidden layer nodes, which was 8. The scatter plot of the estimated and measured values of heavy metal content is shown in Figure 8. It can be seen that among all of the optimal BPNN models, Cr and Cd exhibit the best modeling effect, as reflected by: R² = 0.728, 0.737 and RPD = 1.950, 1.944. Then is Pb, with prediction accuracies (R²) of 0.675.

3.3.3. MEA-BPNN Model Prediction Results

The threshold and weight optimization of BPNN were performed by using MEA. The parameters of BPNN were consistent with those in Section 3.3.2 to enable comparison. The prediction accuracy of metal Cr was the best, and its scatter distribution was closest to the 1:1 line. The prediction accuracies R² of Cr was 0.884 and the RPD value was 3.000. The relatively poor modelling results were detected in Cd and Pb, where R² values were 0.873, 0.857, and the RPD values were 2.114 and 2.560, respectively.

3.3.4. Comparison of Models

In terms of the estimated value scatter plots (Figure 8, Figure 9 and Figure 10), compared with the PLSR and the BPNN, the MEA-BPNN yielded better modeling results, and its estimated values were more evenly distributed on the 1:1 line.

Table 5 lists the best estimates for the three heavy metal elements. It can be seen from the table that all R² was greater than 0.5, and the RPD values are greater than 1.5, which indicated that it was effective in estimating the heavy metals in soil by using the spectroscopic method. Among all the optimal models, the MEA-BPNN model has the best performance in modeling stability and prediction accuracy. It could be considered as a good prediction model. Compared with the BPNN model, R² was improved by 0.145, 0.147, and 0.182, respectively. It can be seen from Table 5 that the prediction ability of PLSR was relatively weaker than that of MEA-BPNN model, especially in terms of model stability. Among all of the optimal MEA-BPNN models, Cr had the best modeling effect, as indicated by: R² = 0.884, RMSE = 33.415, and RPD = 3.000.

4. Discussion

Heavy metal pollution in soil has become a major problem in the ecological restoration of contaminated areas. The use of plants to adsorb heavy metals in soil is one of the green restoration techniques that is inexpensive and causes little damage to soil. Previous studies have revealed that the plant Neyraudia reynaudiana has a strong sensitivity to Pb, and its adsorption capacity is relatively high, while its absorption rate of Cd is also significant, and it exhibits certain remediation potential. The change of Cd found in this study was not large, however, which may be related to the concentration range of Cd in the soil. The overall Cd concentration was lower than that of the other two metals, so the change was not obvious. In this study, it was also discovered that the concentration of the heavy metal Cr in the soil also decreased significantly during the growth of Neyraudia reynaudiana, indicating that this plant has the potential to absorb Cr. After 120 days of restoration (in November 2018), however, the heavy metal content in the soil increased slightly. Due to the decrease of temperature at that time of year in Xuzhou, Jiangsu Province, Neyraudia reynaudiana could not grow normally. Therefore, the activities of antioxidant enzymes and the ability of cells to absorb complex heavy metal ions in Neyraudia reynaudiana were weakened, and then some heavy metal ions were transferred to the soil. However, at the end of the experiment, the heavy metal contents in the soil were still lower than the screening values of soil pollution risk for agricultural land.

Due to the direct spectral response of C–H, N–H, O–H, and C–O covalent bonds in the VNIR spectrum, VNIR spectroscopy has proven to be a feasible method for detecting total soil properties such as SOM, TOC and clay [15]. However, VNIRS cannot directly estimate heavy metal concentrations. The only way VNIRS detects heavy metals is through their association with proxy (organic, iron, and clay minerals). Referring to previous studies, due to the strong adsorption capacity of SOM and iron oxides, heavy metal content is highly correlated with SOM and iron oxide content [17,20]. The overlap of important wavelengths of soil properties also indicated that Cd and Pb were closely related to SOM, while Cr was closely related to Fe. The important bands of SOM are 400–500, 1350–1400, and 1900–2400 nm [10,32,37,47,48,49], while the effective bands of Fe are 400–750, 1520–1680, and 1800–2400 nm [50]. In this study, the effective bands of Cd and SOM primarily overlapped at 450–460, 1310–1330 nm and 2140–2150 nm, while the characteristic wavelengths of SOM and Pb mostly overlapped at 560–680, 1330–1340 nm, 1900–1980 and 2320–2390 nm. The important wavelengths of Cr and SOM and Fe mostly overlapped in the ranges of 400–500, 510–670, and 2000–2020 nm. In summary, these preliminary analyses of spectral reflectance provided the basis for a reference for the selection of characteristic bands.

In this study, we employed the PLSR, BPNN, and MEA-BPNN models to estimate heavy metal concentrations in the soils from VNIR spectra. Many previous studies have analyzed models for the estimation of heavy metal concentrations from VNIR spectra and reported R² values of 0.102–0.95 for Cd, 0.010–0.990 for Cr, and 0.200–0.820 for Pb [10,17,20,37,45,46,47,51,52,53]. We have also summarized some studies on heavy metal estimation using limited samples (less than 50), which reported R² values of 0.102–0.880 for Cd, 0.460–0.860 for Cr, and 0.940 for Pb [17,26,42]. In this study, the optimal R² of the PLSR model for Cd, Cr and Pb contents were 0.759, 0.664 and 0.853, the BPNN model yielded corresponding optimal R² values of 0.728, 0.737, and 0.675, and the MEA-BPNN model yielded corresponding optimal R² values of 0.873, 0.884, and 0.857. These results were consistent with those of previous studies. In this study, the MEA-BPNN model usually produced better results. Compared with the BPNN model, the MEA can optimize the initial weights and thresholds of the BP neural network and it is easier to get the optimal result. For remediated soil, since there may be a nonlinear relationship between heavy metal content and spectral reflectance, the BPNN model would be superior to PLSR. In general, the generalization ability of MEA-BPNN model is better, and the MEA-BPNN model can be better applied to remediated soil with heavy metal content changes. Thus, VNIR spectroscopy coupled with the MEA-BPNN model is a promising method for estimating the concentration of some heavy metals (e.g., Cd, Cr, and Pb) in soils.

The RMSE values in this study were relatively large compared to those of previous investigations [37,46,52], mainly because the accuracy of the estimated model is related to the average of the measured values of soil heavy metal content (Table 3). Due to the concentration gradient setting in the heavy-metal pollution remediation soil experiments of this study, the heavy metal contents of the collected soils were too large during the early stages of plant growth. The arithmetic average of various heavy metal contents was thus higher than the measured values in other investigations, and therefore the predicted root mean square error also increased. However, one of the purposes of this study was to test the feasibility of the VNIR model for estimating heavy metals in soil samples that are undergoing pollution remediation. Therefore, samples with extreme values of soil heavy metal content were not excluded from the original dataset.

In addition, these models also had certain limitations in the study. The small number of soil samples under each heavy metal stress in the study will cause model instability and uncertainty in calibration, which will affect the accuracy and precision of model estimation. Furthermore, the composition and structure of the soil in the field is more complex, and correlations between different heavy metal contents will interfere with the accuracy of the model. Therefore, in future experimental research, optimization of the calibration dataset can be examined and a sufficient number of soil samples will be collected in the field to verify the feasibility of the model.

5. Conclusions

This study demonstrated that VNIR spectroscopy can be used to estimate Cd, Cr, and Pb content in soils under phytoremediation. First-order differential, second-order differential, and reciprocal logarithmic transformations were performed on the pre-processed original spectra and utilized to reduce the influence of external factors. Among them, differential transformation helped to improve the correlation between the heavy metal content and spectral reflectance of the soil, and the use of higher correlation bands significantly improved the stability and predictability of the model. The optimal spectral estimation models for Cd, Cr, and Pb heavy metal elements were established by the BPNN, MEA-BPNN, and PLSR models. The best spectral model established by the MEA-BPNN model was better and had higher precision than those of the other models, and could realize the rapid prediction of three heavy metal elements. Cr exhibited the highest modeling accuracy and a perfect prediction effect, reflected by statistical analysis values of R² = 0.884 and RPD = 3.000. The modeling effects of Cd and Pb were good, with R² values of 0.873 and 0.857, respectively, and RPD values of 2.114 and 2.560, respectively. Therefore, in the future, the use of spectroscopy to estimate heavy metal content will be the main method for the rapid monitoring of soil heavy metal content and soil remediation assessment.

Author Contributions

X.W. made substantial contributions toward the collection, processing, analysis, and interpretation of data and assisted with manuscript preparation; S.A. and Y.X. collected the data and assisted in the revision of the manuscript; Y.Y. provided his experience and was responsible for the proofreading of the manuscript; S.Z. and H.H. were mainly responsible for the content verification; F.C. designed and performed the experiments; and R.L. was responsible for the estimation model. All authors have read and agreed to the published version of the manuscript.

Funding

This research was supported by the National Natural Science Foundation of China (51874307).

Conflicts of Interest

The authors declare no conflict of interest.

Figures and Tables

Generate image description

Figure 1. Experimental process for the remediation of heavy metal contaminated soil. (a) Neyraudia reynaudiana under Cd stress; (b) Neyraudia reynaudiana under Cr stress; (c) Neyraudia reynaudiana under Pb stress.

Generate image description

Figure 2. Structure of the back propagation (BP) neural network.

Generate image description

Figure 3. Flowchart of the mind evolutionary algorithm back propagation neural network (MEA-BPNN) model.

Generate image description

Figure 4. Changes in soil heavy metal concentration with time. (a) The change in Cd concentration; (b) the change in Cr concentration; (c) the change in Pb concentration.

Generate image description

Figure 5. Original spectral reflectance of soils contaminated by (a) Cd, (b) Cr and (c) Pb.

Generate image description

Figure 6. Soil spectra with different spectra pretreated. (a) Spectra after first derivative transformation; (b) spectra after second derivative transformation; (c) spectra after logarithm transformation of reciprocal.

Generate image description

Figure 7. Correlation coefficient between soil spectral reflectance and (a) Cd, (b) Cr and (c) Pb content.

Generate image description

Figure 8. Individual points of heavy metal concentrations in soils predicted by PLSR model compared with measured values. (a) Individual points of Cd concentrations; (b) Individual points of Cr concentrations; (c) Individual points of Pb concentrations.

Generate image description

Figure 9. Individual points of heavy metal concentrations in soils predicted by BPNN model compared with measured values. (a) Individual points of Cd concentrations; (b) Individual points of Cr concentrations; (c) Individual points of Pb concentrations.

Generate image description

Figure 10. Individual points of heavy metal concentrations in soils predicted by MEA-BPNN model compared with measured values. (a) Individual points of Cd concentrations; (b) Individual points of Cr concentrations; (c) Individual points of Pb concentrations.

Table 1

Reference criteria for the experimental design (soil pH = 6.5).

Element	Content in the Panzhihua Mining Area	Risk Screening Value	Risk Intervention Value
Cd (mg/kg)	0.10–0.40	0.3	2.0
Cr (mg/kg)	19.00–495.00	150	850
Pb (mg/kg)	1.60–15.60	90	500

Table 2

Experimental soil concentration gradient settings.

Element	1	2	3	4	5	6
Cd (mg/kg)	0	0.2	0.4	0.6	0.8	1.0
Cr (mg/kg)	0	200	400	600	800	1000
Pb (mg/kg)	0	200	400	600	800	1000

Table 3

Statistical results of heavy metal concentrations in soil samples.

Element	Number of Samples	Maximum (mg/kg)	Minimum (mg/kg)	Standard Deviation (mg/kg)	Skewness	Kurtosis	Coefficient of Variation
Cd	36	29.30	0.13	5.20	4.32	20.49	2.16
Cr	36	813.00	6.90	202.02	2.45	5.83	1.65
Pb	36	397. 20	0.84	112.75	1.72	2.19	1.39

Table 4

Most sensitive spectral bands after different spectra pretreated of each metal pollution stress of soil.

Elements	Transformation	Wavelength (nm)	Correlation Coefficient
Cd	OR	2380	−0.139
	FD	1710	−0.540 **
	SD	2151	−0.602 **
	RL	462	0.150
Cr	OR	2009	−0.559 **
	FD	433	−0.742 **
	SD	480	0.731 **
	RL	2009	0.577 **
Pb	OR	2051	0.384 *
	FD	1333	0.691 **
	SD	534	−0.610 **
	RL	2051	−0.384 *

* Represents significance level p < 0.05, ** represents p < 0.01.

Table 5

Results of heavy metal estimations for three models.

Elements	Models	Transformation	Training Set		Test Set
			R2CV	RMSECV	R2	RMSE	RPD
Cd	PLSR	OR	0.685	3.779	0.759	2.484	1.654
	BPNN	RL	0.560	6.245	0.728	0.994	1.950
	MEA-BPNN	SD	0.767	1.462	0.873	4.316	2.114
Cr	PLSR	FD	0.539	116.997	0.664	141.516	1.730
	BPNN	SD	0.727	295.678	0.737	120.985	1.944
	MEA-BPNN	OR	0.679	103.845	0.884	33.415	3.000
Pb	PLSR	FD	0.502	102.271	0.852	60.585	2.193
	BPNN	FD	0.518	102.081	0.675	69.222	1.792
	MEA-BPNN	SD	0.711	67.169	0.857	50.926	2.560