Quantitative analysis method for heavy metals in soil

Abstract

The invention discloses a quantitative analysis method for heavy metals in soil, which comprises the following steps: dividing a high-difference region and a low-difference region according to the difference significance of the heavy metal concentration between in-situ spectral data and a pollution source distance, and carrying out scattering correction on the characteristic peak of the in-situ spectral data in the high-difference region by using the scattering characteristic of a heavy metal spectrum; reversely simulating heavy metal scattering artifacts according to the scattering correction spectrum and the in-situ spectrum data of the low-difference region, performing artifact repair through Mie scattering constraint based on the simulated heavy metal scattering artifacts, and performing bilinear interpolation based on spectral characteristics of an artifact repair spectrum in a detection region. And inputting the interpolated heavy metal spectrum distribution grid and heavy metal concentration into a multilayer perceptron quantitative model to obtain a heavy metal quantitative result and a heavy metal distribution result. According to the method, the accuracy of heavy metal detection and the interpretability of heavy metal distribution are improved by carrying out scattering correction and artifact repair on spectral data of a detection area.

Description

Technical Field

[0001] The present invention relates to the technical field of quantitative analysis, and in particular to a method for quantitative analysis of heavy metals in soil. Background Art

[0002] Heavy metal contamination of soil is a global environmental problem, posing a serious threat to ecosystems, human health, and sustainable agricultural development. With the rapid development of industrialization and urbanization, large amounts of heavy metals such as cadmium (Cd), lead (Pb), mercury (Hg), and arsenic (As) have entered the soil environment through various pathways, including industrial waste discharge, exhaust gas deposition from coal-fired power plants, agricultural activities such as fertilizer and pesticide application, and traffic exhaust. These heavy metals have poor mobility and are difficult to degrade in soil, easily accumulating and persisting in the soil for long periods of time. They then enter the human body through bioaccumulation in the food chain, causing serious harm to human health, including neurological diseases, kidney damage, and cancer.

[0003] However, traditional soil heavy metal detection technologies, such as atomic absorption spectroscopy (AAS), inductively coupled plasma optical emission spectroscopy (ICP-OES), and X-ray fluorescence (XRF), have the disadvantages of being complex, time-consuming, costly, and unable to achieve real-time on-site detection. Laser-induced breakdown spectroscopy (LIBS), as an emerging in-situ spectral analysis technology, can quickly and in-situ detect heavy metals in soil without complex pretreatment. LIBS uses high-energy laser pulses to excite soil samples to generate plasma, and uses a spectrometer to analyze the spectrum emitted by the plasma to determine the elemental composition and content. However, the accuracy of LIBS quantitative analysis is easily affected by factors such as soil matrix effects, particle size, and humidity. The quantitative analysis results only represent the content concentration at the sampling point and cannot effectively determine the diffusion distribution of heavy metals. Therefore, developing a method that can effectively correct for these interfering factors and improve the accuracy of LIBS quantitative analysis has become an urgent problem to be solved. Summary of the Invention

[0004] The purpose of the present invention is to provide a method for quantitative analysis of heavy metals in soil.

[0005] To achieve the above object, the present invention is implemented according to the following technical solutions:

[0006] A first aspect of the present invention provides a method for quantitative analysis of heavy metals in soil, comprising:

[0007] Collecting soil sample data for preprocessing, the soil sample data including in-situ spectral data, auxiliary data, spatial parameters, heavy metal concentrations, and pollution source coordinates;

[0008] The sampling coordinates of the in-situ spectral data are mapped to the detection area grid through spatial parameters, the significance of the difference in heavy metal concentration is obtained according to the distance between the detection area grid and the pollution source coordinates, and high difference areas and low difference areas are obtained according to the significance of the difference;

[0009] performing particle size analysis on the particle size distribution and moisture content of the auxiliary data in the high-difference region to obtain heavy metal spectral scattering characteristics of the soil, and performing scattering correction on characteristic peaks of the in-situ spectral data using the heavy metal spectral scattering characteristics;

[0010] Reversely simulate heavy metal scattering artifacts based on the scattering correction spectrum and the in-situ spectral data of the low-difference area, perform artifact repair on the scattering correction spectrum and the in-situ spectral data based on the simulated heavy metal scattering artifacts through Mie scattering constraints, extract heavy metal spectral features of the artifact repaired spectrum based on the detection area, perform bilinear interpolation on the heavy metal spectral features, and obtain a heavy metal spectral distribution grid;

[0011] The heavy metal spectrum distribution grid and heavy metal concentration are input into a multi-layer perceptron quantitative model to obtain heavy metal quantitative results and heavy metal distribution results in the detection area.

[0012] As a further method, the pretreatment method comprises:

[0013] LIBS spectral data of soil samples around heavy metal pollution sources are collected as in-situ spectral data. The heavy metal concentration corresponding to the in-situ spectral data is collected. The particle size distribution and humidity of the soil samples are collected as auxiliary data. The coordinates of the pollution source are collected, and the grid units are obtained according to the latitude and longitude of the sampling point coordinates. The sampling point coordinates and grid units are used as spatial parameters.

[0014] As a further method, a method for mapping the sampling coordinates of the in-situ spectral data to a detection area grid through spatial parameters includes:

[0015] The Euclidean distance between the center coordinates of the grid cells and the coordinates of the pollution source is calculated based on the spatial parameters. The in-situ spectral data of 3 to 5 soil samples are extracted as grid samples using the plum blossom method. The grid cells are interpolated using co-kriging based on the spectral characteristic peak intensity and Euclidean distance of the grid samples to obtain the continuous spectral distribution within the grid cells.

[0016] If the mean absolute error between the continuous spectral distribution within a grid cell and the measured value of the grid sample is greater than 10%, 1 to 2 grid samples are added to the grid and re-interpolated. If the mean absolute error is less than or equal to 10%, the continuous spectral distribution within the grid cell is retained to obtain the detection area grid.

[0017] As a further method, a method for obtaining the high-difference region and the low-difference region includes:

[0018] Based on the detection area grid, the Euclidean distance between the sampling point coordinates and the pollution source coordinates and the spectral characteristic peak intensity of the sampling point are obtained. The data pairs consisting of the Euclidean distance and the spectral characteristic peak intensity are clustered to obtain the near-source cluster and the far-source cluster. The middle value of the near-source cluster and the far-source cluster is used as the distance threshold.

[0019] The detection area grid is divided into several intervals according to the Euclidean distance. The spectral characteristic peak intensity within the interval is subjected to a one-way variance analysis to obtain the difference significance level values ​​of different intervals. If the difference significance level value is greater than or equal to 7%, no regional division is performed.

[0020] If the significance level of the difference is less than 7%, the reciprocal of the normalized Euclidean distance from the sampling point to the pollution source is used as the spatial weight. The multi-scale feature fusion value is calculated based on the spatial weight and the spectral characteristic peak intensity and half-height width of the spectral characteristic peak extracted from the in-situ spectral data of the corresponding sampling point. The dynamic threshold is obtained by using the multi-scale feature fusion value and the distance threshold through the dynamic threshold formula. The dynamic threshold formula is:

[0021]

[0022] Where τ(D) is the multi-scale feature fusion value F i The dynamic threshold of μ is the multi-scale feature fusion value F obtained by calculating the in-situ spectral data and spatial weight of the sub-sample in the grid unit. i The statistical mean of σ is the multi-scale feature fusion value F i The statistical standard deviation of the data, Φ is the cumulative distribution function of the standard normal distribution, D is the Euclidean distance from the sampling point to the pollution source, μ D is the mean of all D, σ D is the standard deviation of all D, is the average intensity of the spectral characteristic peak of the in-situ spectral data of the near-source region sub-sample, peak is the characteristic emission peak of heavy metal elements, is the mean intensity of the spectral characteristic peak of the in-situ spectral data of the far-source region sub-sample, D1 is the distance threshold between the near-source region cluster and the far-source region cluster, is the mean half-maximum width of the spectral characteristic peak of the in-situ spectral data of the near-source area sample, L is the side length of the grid unit, D2 is the sum of D1 and the side length of the grid unit, and is the transition area;

[0023] The grid cells whose multi-scale feature fusion values ​​of sampling points in the interval are greater than the dynamic threshold are divided into high-difference areas, and the rest are divided into low-difference areas.

[0024] As a further method, the scatter correction method includes:

[0025] The particle size distribution and moisture content were extracted based on auxiliary data from high-difference areas. The particle size distribution was discretized into an equivalent spherical diameter sequence. The equivalent spherical diameter sequence was classified according to the International Soil Texture Triangle Diagram standard to obtain the median particle size and soil type coefficient.

[0026] Based on the median particle size and soil type coefficient, the humidity is converted into a dielectric constant using the GWC model. The dielectric constant is used to calculate the real part of the complex refractive index. The scale parameter is obtained using the relative scale of the median particle size and the laser wavelength of the in-situ spectral data of the corresponding sampling point. The Mie scattering amplitude coefficient is calculated through Anning recursion based on the median particle size, dielectric constant, real part approximation of the complex refractive index, and scale parameter, and the spectral scattering characteristics of heavy metals including the median particle size, scale parameter, and amplitude coefficient are obtained.

[0027] The spectral data of high-difference areas are corrected using the scattering correction formula using the spectral scattering characteristics of heavy metals. The scattering correction formula is:

[0028]

[0029] Among them I corr (λ) is the spectrum after scattering correction, I raw (λ) is the spectral data of the first region, I0(λ) is the light intensity of the laser directly irradiating the Spectralon standard white plate after dark field subtraction, I ref (λ) is the spectral response of the standard white plate at the soil measurement position, The equivalent projected area of ​​soil particles obtained based on the median particle size, α is the scale parameter of soil particle size and laser wavelength λ, l is the multipole order of Mie scattering, N is the truncation term, N = α + 4α 1 / 3 +2, a l is the amplitude coefficient of the Mie scattering electric dipole, b l is the amplitude coefficient of the magnetic dipole of Mie scattering.

[0030] As a further method, a method for obtaining the artifact repair spectrum includes:

[0031] Scattering-corrected spectral data are extracted based on the high-difference area. The scattering-corrected spectral data and the in-situ spectral data of the low-difference area are used as the reference spectrum. The scattering correction formula is used to reversely simulate the heavy metal scattering artifacts through the intensity scaling method based on the reference spectrum, median particle size, and humidity. The in-situ spectral data of the high-difference area and the low-difference area are mixed with the simulated heavy metal scattering artifacts and divided into training set, validation set, and test set;

[0032] The training set is used as input to a Mie scattering constrained autoencoder constructed based on a convolutional encoder, a Mie scattering physical constraint layer, and a deconvolution decoder. The scattering correction formula is inserted into the Mie scattering physical constraint layer according to the differential constraints of automatic differentiation. Reconstruction loss and KL divergence are used as a joint loss for training. A validation set is added during training to control accuracy convergence. The trained autoencoder model is obtained, and the test set is used to optimize the model performance by the mean squared error.

[0033] The scattering-corrected spectra in the high-difference area and the in-situ spectral data in the low-difference area are input into the autoencoder model. The noise artifacts in the high-difference area and the low-difference area are repaired by forward propagation based on the autoencoder model to obtain the artifact-repaired spectrum.

[0034] As a further method, a method for obtaining the heavy metal spectral characteristics includes:

[0035] Based on the detection area, the artifact repair spectrum is extracted and peak detection is performed through continuous wavelet transform to obtain the characteristic peak, the total number of characteristic peaks and the central wavelength of the characteristic peak. The spectrum mean of the area far from the characteristic peak is used as the baseline initial value. According to the baseline, the spectral intensity at the center of the Lorentz peak is used as the initial value of the amplitude of the Lorentz peak. The typical half-maximum width of the LIBS plasma peak is used as the initial value of the half-maximum width of the Lorentz peak. The spectral intensity at the center of the Gaussian peak is used as the initial value of the amplitude of the Gaussian peak. The standard deviation of the baseline noise near the center of the Gaussian peak is used as the initial value of the broadening of the Gaussian peak.

[0036] The initial value of the baseline, the initial value of the amplitude of the Lorentz peak, the initial value of the half-height width of the Lorentz peak, the initial value of the amplitude of the Gaussian peak, and the initial value of the broadening of the Gaussian peak are used to perform Lorentz-Gaussian mixture fitting through the loss function formula. The loss function formula is:

[0037]

[0038] in is the loss function, which measures the deviation between the fitted model and the true repair spectrum, I recon (λ) is the artifact repair spectrum, M is the total number of characteristic peaks, A i is the amplitude of the i-th Lorentz peak, Γ i is the half-height width of the i-th Lorentz peak, λ 0i is the central wavelength of the i-th characteristic peak, B i is the amplitude of the i-th Gaussian peak, σ i is the standard deviation of the i-th Gaussian peak, C is the baseline offset, which is obtained from the humidity through the GWC model;

[0039] The loss function is optimized by least squares iteration until the residual change rate is less than 10 -6The iteration is terminated and the baseline, Lorentz amplitude, Lorentz half-maximum width, Gaussian amplitude and Gaussian broadening are obtained as the heavy metal spectral characteristics.

[0040] As a further method, a method for obtaining the heavy metal spectrum distribution grid includes:

[0041] Based on the detection area grid corresponding to the heavy metal spectral characteristics, the grid coordinates and sampling timestamps are extracted to form a feature data set. The feature data set is spatially bilinearly interpolated according to the side length of the grid unit. The minimum sampling time interval of the feature data set is used as the time axis, and the feature data set is temporally bilinearly interpolated according to the time axis. The mean absolute error and root mean square error between the interpolation result and the original sampling point are calculated. If the mean absolute error is greater than 4%, the bilinear interpolation coefficient is optimized using the covariance function of Kriging interpolation and then recalculated. If the mean absolute error is less than or equal to 4%, the interpolation result is retained to obtain the heavy metal spectral distribution grid.

[0042] As a further method, a method for obtaining the heavy metal spectrum distribution grid includes:

[0043] Based on the detection area grid corresponding to the heavy metal spectral characteristics, the grid coordinates and sampling timestamps are extracted to form a feature data set. The feature data set is spatially bilinearly interpolated according to the side length of the grid unit. The minimum sampling time interval of the feature data set is used as the time axis, and the feature data set is temporally bilinearly interpolated according to the time axis. The mean absolute error and root mean square error between the interpolation result and the original sampling point are calculated. If the mean absolute error is greater than 4%, the bilinear interpolation coefficient is optimized using the covariance function of Kriging interpolation and then recalculated. If the mean absolute error is less than or equal to 4%, the interpolation result is retained to obtain the heavy metal spectral distribution grid.

[0044] A second aspect of the present invention provides a system for quantitative analysis of heavy metals in soil, comprising:

[0045] Data acquisition module: collects soil sample data for preprocessing, including in-situ spectral data, auxiliary data, spatial parameters, heavy metal concentrations, and pollution source coordinates;

[0046] Region division module: used to map the sampling coordinates of the in-situ spectral data to the detection area grid through spatial parameters, obtain the significance of the difference in heavy metal concentration based on the distance between the detection area grid and the pollution source coordinates, and divide the high-difference area and low-difference area according to the significance of the difference;

[0047] Scattering correction module: used to perform particle size analysis on the particle size distribution and humidity of the auxiliary data in the high-difference area, obtain the heavy metal spectrum scattering characteristics of the soil, and use the heavy metal spectrum scattering characteristics to perform scattering correction on the characteristic peaks of the in-situ spectrum data;

[0048] A heavy metal spectrum distribution grid acquisition module is used to reversely simulate heavy metal scattering artifacts based on the scattering correction spectrum and the in-situ spectrum data of the low-difference area, perform artifact repair on the scattering correction spectrum and the in-situ spectrum data based on the simulated heavy metal scattering artifacts through Mie scattering constraints, extract heavy metal spectrum features of the artifact repair spectrum based on the detection area, perform bilinear interpolation on the heavy metal spectrum features, and obtain a heavy metal spectrum distribution grid;

[0049] Quantitative analysis module: used to input the heavy metal spectrum distribution grid and heavy metal concentration into the multi-layer perceptron quantitative model to obtain the heavy metal quantitative results and heavy metal distribution results of the detection area.

[0050] Compared with the prior art, the embodiments of the present invention have at least the following advantages or beneficial effects:

[0051] The present invention achieves differentiated processing of pollution sources through regional division and spatial analysis, and can obtain accurate data on heavy metal content in soil after a pollution incident occurs. Scatter correction and artifact repair are performed on the spectral data of soil samples to effectively eliminate the influence of factors such as soil particle size and humidity on LIBS spectral signals, thereby improving the quality of spectral data and enhancing the accuracy and reliability of quantitative analysis of heavy metals in soil. Feature extraction through Lorentz-Gaussian mixture fitting and combined with multi-layer perceptron quantitative model analysis improves the model's adaptability to different soil conditions, and the distribution, diffusion and trend of heavy metals in soil are obtained, making it widely applicable in different soil types and pollution scenarios. BRIEF DESCRIPTION OF THE DRAWINGS

[0052] Figure 1 The figure is a flowchart of the steps of a method for quantitative analysis of heavy metals in soil according to an embodiment of the present invention. DETAILED DESCRIPTION

[0053] The technical solutions in the embodiments of the present invention will be clearly and completely described below in conjunction with the drawings in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, rather than all the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without making creative work are within the scope of protection of the present invention.

[0054] Reference Figure 1 As shown, the present invention provides a method for quantitative analysis of heavy metals in soil, comprising:

[0055] Collecting soil sample data for preprocessing, the soil sample data including in-situ spectral data, auxiliary data, spatial parameters, heavy metal concentrations, and pollution source coordinates;

[0056] In the actual assessment, one hectare of farmland around a retired electroplating plant was selected. The electroplating plant had a history of heavy metal emissions of Cd and Pb. The soil type was clay loam. Sampling points were arranged in a 5-meter square grid, and a total of 200 sampling points were collected. The soil samples were air-dried and ground. A laser diffractometer was used to obtain the particle size distribution. The median particle size of the soil was 45 μm. A TDR soil moisture meter was used to obtain a moisture content of 14.6% as humidity. In situ spectral data were obtained using a LIBS system with a LIBS wavelength of 1064 nm, a pulse energy of 100 mJ, and a repetition rate of 10 Hz. The spectrometer wavelength range was 200–1000 nm, and the resolution was 0.1 nm. The heavy metal concentrations of the soil samples were collected, with Cd ranging from 0.2 to 4.8 mg / kg and Pb ranging from 12 to 198 mg / kg. The longitude and latitude of the sampling points in WGS84 format and the grid unit were used as spatial parameters to collect the coordinates of the electroplating plant, the pollution source.

[0057] The sampling coordinates of the in-situ spectral data are mapped to the detection area grid through spatial parameters, the significance of the difference in heavy metal concentration is obtained according to the distance between the detection area grid and the pollution source coordinates, and high difference areas and low difference areas are obtained according to the significance of the difference;

[0058] In the actual assessment, the Euclidean distance between the center coordinates of the grid cells and the coordinates of the pollution source was calculated based on spatial parameters. Five points (four vertices plus one center) were selected for each 5-meter square grid using the plum blossom point method. In-situ spectral data of four soil samples were extracted, totaling 800 grid samples. Based on the spectral characteristic peak intensities of Cd (228.8 nm) and Pb (405.7 nm) in the grid samples and the Euclidean distance, the grid cells were interpolated using co-kriging to obtain a continuous spectral distribution within the grid cells.

[0059] If the average absolute error of the continuous spectral distribution within a grid cell and the measured value of the grid sample is greater than 10%, one grid sample is added to the grid and re-interpolated. The average absolute error after re-interpolation is 7.8%, and a 50 by 50 grid covering 1 km is obtained. 2 , as the detection area grid, based on the detection area grid, the Euclidean distance between the sampling point coordinates and the pollution source coordinates and the spectral characteristic peak intensity of the sampling point are obtained, and the data pairs composed of the Euclidean distance and the spectral characteristic peak intensity are clustered by K-means++ to obtain the near-source cluster and the far-source cluster. The near-source cluster is an area less than 200m, and its spectral characteristic peak is high. The Euclidean distance between the sampling point coordinates and the pollution source coordinates in the far-source zone is greater than 500m, and the spectral characteristic peak is weak. The middle value of 350m, which is the near-source cluster and the far-source cluster, is used as the distance threshold;

[0060] The detection area grid was divided into several intervals according to the Euclidean distance, and the spectral characteristic peak intensity within the interval was subjected to a one-way analysis of variance to obtain the difference significance level values ​​of different intervals; the difference significance level value of 3.2% between the near-source area and the middle-source area was less than 7%. The reciprocal of the normalized Euclidean distance from the sampling point to the pollution source was used as the spatial weight. The multi-scale feature fusion value was calculated based on the spatial weight and the spectral characteristic peak intensity and half-height width of the spectral characteristic peak extracted from the in-situ spectral data of the corresponding sampling point. The dynamic threshold value of 0.873 was obtained by using the multi-scale feature fusion value and the distance threshold through the dynamic threshold formula. The grid cells whose multi-scale feature fusion values ​​of the sampling points in the interval were greater than the dynamic threshold value of 0.873 were divided into high-difference areas, totaling 120 grids, and the rest were divided into low-difference areas, totaling 80 grids.

[0061] performing particle size analysis on the particle size distribution and moisture content of the auxiliary data in the high-difference region to obtain heavy metal spectral scattering characteristics of the soil, and performing scattering correction on characteristic peaks of the in-situ spectral data using the heavy metal spectral scattering characteristics;

[0062] In the actual assessment, particle size distribution and moisture content were extracted based on auxiliary data from high-variability areas. The particle size distribution was discretized into an equivalent spherical diameter sequence ranging from 1 to 100 μm. The equivalent spherical diameter sequence was classified according to the USDA International Soil Texture Triangle Diagram standard, resulting in a median particle size of 45 μm and a soil type coefficient of 0.85 for clay loam. Since the median particle size can represent most morphologies of soil samples, when the particles are approximately spherical, the calculation of the T matrix can be degenerated into an analytical solution for Mie scattering. In this case, the Anning recursive method is used to calculate the Mie scattering amplitude coefficient.

[0063] In the actual evaluation, humidity was converted into a dielectric constant of 35.6 through the GWC model based on the median particle size and soil type coefficient. The real part of the complex refractive index was approximately 1.49, which was calculated by taking the square root of the dielectric constant divided by 2. The scale parameter 263 was obtained by using the relative scale of the median particle size and the laser wavelength of the in-situ spectral data of the corresponding sampling point. The Mie scattering amplitude coefficient was calculated through Anning recursion based on the median particle size, dielectric constant, real part approximation of the complex refractive index, and scale parameter, and its truncation term was 275. The spectral scattering characteristics of heavy metals such as the median particle size, scale parameter, and amplitude coefficient were used to correct the spectral data in the high-difference area through the scattering correction formula, so that the half-maximum width of the Cd peak was reduced from 2.1nm to 0.78nm, and the half-maximum width of the Pb peak was reduced from 1.8nm to 0.65nm.

[0064] Reversely simulate heavy metal scattering artifacts based on the scattering correction spectrum and the in-situ spectral data of the low-difference area, and perform artifact repair on the scattering correction spectrum and the in-situ spectral data based on the simulated heavy metal scattering artifacts through Mie scattering constraints

[0065] In the actual evaluation, scattering correction spectral data were extracted based on the high-difference area, and the scattering correction spectral data and the in-situ spectral data of the low-difference area were used as the reference spectra. The scattering artifacts were reversely simulated by the intensity scaling method using the scattering correction formula according to the reference spectrum, median particle size, and humidity to obtain 2000 spectra of simulated heavy metal scattering artifacts. The simulated heavy metal scattering artifacts were mixed with the in-situ spectral data of the high-difference area and the low-difference area and divided into training set, validation set, and test set according to the ratio of 8 to 1 to 1. The training set was used as input to the Mie scattering constrained autoencoder constructed based on the convolution encoder, Mie scattering physical constraint layer, and deconvolution decoder. The convolution encoder has 3 layers and the number of channels is 8→16→32. The deconvolution decoder is also 3. The number of channels is increased from 32 to 16 to 8. The scattering correction formula is inserted into the Mie scattering physical constraint layer according to the differential constraint of automatic differentiation. Reconstruction loss and KL divergence are used as joint losses for training. A validation set is added during training to control the convergence of accuracy. After 120 rounds of iteration, the MAE is reduced from 0.12 to 0.03, and the trained autoencoder model is obtained. The model performance is optimized by mean square error using the test set. The corrected spectral data of the high-difference area and the in-situ spectral data of the low-difference area are input into the autoencoder model. The noise artifacts in the high-difference area and the low-difference area are repaired by forward propagation based on the autoencoder model, and the artifact repair spectrum is output, which reduces the noise standard deviation in the detection area by 67% and makes the peak shape sharper.

[0066] The heavy metal spectrum features of the artifact repair spectrum are extracted based on the detection area, and bilinear interpolation is performed on the heavy metal spectrum features to obtain the heavy metal spectrum distribution grid;

[0067] In the actual evaluation, peak detection was performed on the artifact repair spectrum using continuous wavelet transform, setting the Morlet basis and the scale to 1–30. Two characteristic peaks, Cd at 228.8 nm and Pb at 405.7 nm, were detected. The initial values ​​of the Lorentz-Gaussian mixture were set according to the characteristic peaks: the baseline was the background mean of 0.08 in the 250–300 nm region away from the peak, the Lorentz peak amplitude was the peak center intensity A1 = 0.85 and A2 = 0.72, the Lorentz peak half-maximum width was the typical LIBS value, where Γ1 = 0.8 nm and Γ2 = 0.7 nm, the Gaussian peak amplitude was the peak center noise, where B1 = 0.12 and B2 = 0.11, and the Gaussian peak broadening was the standard deviation of the baseline noise, where σ1 = 0.3 nm and σ2 = 0.25 nm. The loss function was fitted and optimized using the Levenberg-Marquardt algorithm until the residual change rate was less than 10 -6The iteration was terminated, and the final fitting residual MAE was 0.012. Seven types of features were extracted: baseline, Cd Lorentz amplitude, Cd Lorentz half-width, Cd Gaussian amplitude, Cd Gaussian broadening, Pb Lorentz amplitude, and Pb Lorentz half-width. Although M is 2, the Gaussian parameters of Pb can be ignored due to weak pollution, and Pb Gaussian amplitude and Pb Gaussian broadening are close to 0. Therefore, the effective feature is 6-dimensional as the heavy metal spectral feature. The grid coordinates and sampling timestamps are extracted based on the detection area grid corresponding to the heavy metal spectral feature to form a feature data set. There are 200 sampling points in the feature data set multiplied by 24 hours, with a time interval of 1 hour, for a total of 4,800 records. The feature data set is spatially bilinearly interpolated according to the grid unit side length of 5m, and temporally bilinearly interpolated according to the time axis of 1 hour interval. The mean absolute error and root mean square error of the interpolation result and the original sampling point are calculated. The initial interpolation mean absolute error MAE is 4.2%, which is greater than 4%. By introducing the covariance function of Kriging interpolation d is the spatial distance, 3×5m=15m. After optimization, the mean absolute error MAE is 3.7%, and a 50×50×24 three-dimensional space-time heavy metal spectrum distribution grid is obtained, each grid contains 6-dimensional spectral features.

[0068] The heavy metal spectrum distribution grid and heavy metal concentration are input into a multi-layer perceptron quantitative model to obtain heavy metal quantitative results and heavy metal distribution results in the detection area.

[0069] In the actual evaluation, the heavy metal concentration of the corresponding coordinate is used as a label based on the heavy metal spectral distribution grid, and the heavy metal spectral distribution grid is flattened from a three-dimensional tensor to a two-dimensional feature matrix to obtain a spatiotemporal feature concentration dataset. The Z-score is normalized according to the mean and standard deviation of the spatiotemporal feature concentration dataset. The standardized spatiotemporal feature concentration dataset is divided into training set, validation set, and test set at a ratio of 7:2:1, with 3360 training sets, 960 validation sets, and 480 test sets. The training set is input into a 4-layer multi-layer perceptron, and the spatiotemporal feature dimension of the training set data is used as the number of input layer neurons, with channels of 6→128→64→1, 128 neurons in hidden layer 1, and ReLU activation is used, with Dropout=0.2, 64 neurons in hidden layer 2, and ReLU activation plus BatchNormalization. The output layer 1 neurons output the regression concentration, and the mean square error is used as the loss function for training. The Adam optimizer is used, and the loss precision lr is set to 10 -3During training, the mean square error of the validation set did not decrease for 10 consecutive rounds and the training was completed at the 85th round. The multi-layer perceptron quantitative model was obtained. According to the predicted concentration, the mean of the true concentration of the training set and the total number of continuous spatiotemporal grids of the test set, the determination coefficient of Cd was 0.92 and the determination coefficient of Pb was 0.91. The generalization ability of the multi-layer perceptron quantitative model was optimized based on the determination coefficient. The multi-layer perceptron quantitative model was forward propagated on the complete continuous spatiotemporal grid, and a three-dimensional heavy metal concentration prediction grid was output to obtain the heavy metal quantitative analysis results. The Cd in one grid was 3.25 mg / kg and the Pb was 140.25 mg / kg. According to the quantitative results, a heavy metal distribution heat map was generated according to the predicted grid to show the changes in heavy metal distribution around the pollution source.

[0070] In this embodiment, the preprocessing method includes:

[0071] LIBS spectral data of soil samples around heavy metal pollution sources are collected as in-situ spectral data. The heavy metal concentration corresponding to the in-situ spectral data is collected. The particle size distribution and humidity of the soil samples are collected as auxiliary data. The coordinates of the pollution source are collected, and the grid units are obtained according to the latitude and longitude of the sampling point coordinates. The sampling point coordinates and grid units are used as spatial parameters.

[0072] In this embodiment, the method of mapping the in-situ spectral data to the detection area grid according to the spatial parameters includes:

[0073] The Euclidean distance between the center coordinates of the grid cells and the coordinates of the pollution source is calculated based on the spatial parameters. The in-situ spectral data of 3 to 5 soil samples are extracted as grid samples using the plum blossom method. The grid cells are interpolated using co-kriging based on the spectral characteristic peak intensity and Euclidean distance of the grid samples to obtain the continuous spectral distribution within the grid cells.

[0074] If the mean absolute error between the continuous spectral distribution within a grid cell and the measured value of the grid sample is greater than 10%, 1 to 2 grid samples are added to the grid and re-interpolated. If the mean absolute error is less than or equal to 10%, the continuous spectral distribution within the grid cell is retained to obtain the detection area grid.

[0075] In this embodiment, the method for obtaining the high-difference region and the low-difference region includes:

[0076] Based on the detection area grid, the Euclidean distance between the sampling point coordinates and the pollution source coordinates and the spectral characteristic peak intensity of the sampling point are obtained. The data pairs consisting of the Euclidean distance and the spectral characteristic peak intensity are clustered to obtain the near-source cluster and the far-source cluster. The middle value of the near-source cluster and the far-source cluster is used as the distance threshold.

[0077] The detection area grid is divided into several intervals according to the Euclidean distance. The spectral characteristic peak intensity within the interval is subjected to a one-way variance analysis to obtain the difference significance level values ​​of different intervals. If the difference significance level value is greater than or equal to 7%, no regional division is performed.

[0078] If the significance level of the difference is less than 7%, the reciprocal of the normalized Euclidean distance from the sampling point to the pollution source is used as the spatial weight. The multi-scale feature fusion value is calculated based on the spatial weight and the spectral characteristic peak intensity and half-height width of the spectral characteristic peak extracted from the in-situ spectral data of the corresponding sampling point. The dynamic threshold is obtained by using the multi-scale feature fusion value and the distance threshold through the dynamic threshold formula. The dynamic threshold formula is:

[0079]

[0080] Where τ(D) is the multi-scale feature fusion value F i The dynamic threshold of μ is the multi-scale feature fusion value F obtained by calculating the in-situ spectral data and spatial weight of the sub-sample in the grid unit. i The statistical mean of σ is the multi-scale feature fusion value F i The statistical standard deviation of the data, Φ is the cumulative distribution function of the standard normal distribution, D is the Euclidean distance from the sampling point to the pollution source, μ D is the mean of all D, σ D is the standard deviation of all D, is the average intensity of the spectral characteristic peak of the in-situ spectral data of the near-source region sub-sample, peak is the characteristic emission peak of heavy metal elements, is the mean intensity of the spectral characteristic peak of the in-situ spectral data of the far-source region sub-sample, D1 is the distance threshold between the near-source region cluster and the far-source region cluster, is the mean half-maximum width of the spectral characteristic peak of the in-situ spectral data of the near-source area sample, L is the side length of the grid unit, D2 is the sum of D1 and the side length of the grid unit, and is the transition area;

[0081] The grid cells whose multi-scale feature fusion values ​​of sampling points in the interval are greater than the dynamic threshold are divided into high-difference areas, and the rest are divided into low-difference areas.

[0082] In this embodiment, the scatter correction method includes:

[0083] The particle size distribution and moisture content were extracted based on auxiliary data from high-difference areas. The particle size distribution was discretized into an equivalent spherical diameter sequence. The equivalent spherical diameter sequence was classified according to the International Soil Texture Triangle Diagram standard to obtain the median particle size and soil type coefficient.

[0084] Based on the median particle size and soil type coefficient, the humidity is converted into a dielectric constant using the GWC model. The dielectric constant is used to calculate the real part of the complex refractive index. The scale parameter is obtained using the relative scale of the median particle size and the laser wavelength of the in-situ spectral data of the corresponding sampling point. The Mie scattering amplitude coefficient is calculated through Anning recursion based on the median particle size, dielectric constant, real part approximation of the complex refractive index, and scale parameter, and the spectral scattering characteristics of heavy metals including the median particle size, scale parameter, and amplitude coefficient are obtained.

[0085] The spectral data of high-difference areas are corrected using the scattering correction formula using the spectral scattering characteristics of heavy metals. The scattering correction formula is:

[0086]

[0087] Among them I corr (λ) is the spectrum after scattering correction, I raw (λ) is the spectral data of the first region, I0(λ) is the light intensity of the laser directly irradiating the Spectralon standard white plate after dark field subtraction, I ref (λ) is the spectral response of the standard white plate at the soil measurement position, The equivalent projected area of ​​soil particles obtained based on the median particle size, α is the scale parameter of soil particle size and laser wavelength λ, l is the multipole order of Mie scattering, N is the truncation term, N = α + 4α 1 / 3 +2, a l is the amplitude coefficient of the Mie scattering electric dipole, b l is the amplitude coefficient of the magnetic dipole of Mie scattering.

[0088] In this embodiment, the method for obtaining the artifact repair spectrum includes:

[0089] Scattering-corrected spectral data are extracted based on the high-difference area. The scattering-corrected spectral data and the in-situ spectral data of the low-difference area are used as the reference spectrum. The scattering correction formula is used to reversely simulate the heavy metal scattering artifacts through the intensity scaling method based on the reference spectrum, median particle size, and humidity. The in-situ spectral data of the high-difference area and the low-difference area are mixed with the simulated heavy metal scattering artifacts and divided into training set, validation set, and test set;

[0090] The training set is used as input to a Mie scattering constrained autoencoder constructed based on a convolutional encoder, a Mie scattering physical constraint layer, and a deconvolution decoder. The scattering correction formula is inserted into the Mie scattering physical constraint layer according to the differential constraints of automatic differentiation. Reconstruction loss and KL divergence are used as a joint loss for training. A validation set is added during training to control accuracy convergence. The trained autoencoder model is obtained, and the test set is used to optimize the model performance by the mean squared error.

[0091] The scattering-corrected spectra in the high-difference area and the in-situ spectral data in the low-difference area are input into the autoencoder model. The noise artifacts in the high-difference area and the low-difference area are repaired by forward propagation based on the autoencoder model to obtain the artifact-repaired spectrum.

[0092] In this embodiment, the method for obtaining the heavy metal spectral characteristics includes:

[0093] Based on the detection area, the artifact repair spectrum is extracted and peak detection is performed through continuous wavelet transform to obtain the characteristic peak, the total number of characteristic peaks and the central wavelength of the characteristic peak. The spectrum mean of the area far from the characteristic peak is used as the baseline initial value. According to the baseline, the spectral intensity at the center of the Lorentz peak is used as the initial value of the amplitude of the Lorentz peak. The typical half-maximum width of the LIBS plasma peak is used as the initial value of the half-maximum width of the Lorentz peak. The spectral intensity at the center of the Gaussian peak is used as the initial value of the amplitude of the Gaussian peak. The standard deviation of the baseline noise near the center of the Gaussian peak is used as the initial value of the broadening of the Gaussian peak.

[0094] The initial value of the baseline, the initial value of the amplitude of the Lorentz peak, the initial value of the half-height width of the Lorentz peak, the initial value of the amplitude of the Gaussian peak, and the initial value of the broadening of the Gaussian peak are used to perform Lorentz-Gaussian mixture fitting through the loss function formula. The loss function formula is:

[0095]

[0096] in is the loss function, which measures the deviation between the fitted model and the true repair spectrum, I recon (λ) is the artifact repair spectrum, M is the total number of characteristic peaks, A i is the amplitude of the i-th Lorentz peak, Γ i is the half-height width of the i-th Lorentz peak, λ 0i is the central wavelength of the i-th characteristic peak, B i is the amplitude of the i-th Gaussian peak, σ i is the standard deviation of the i-th Gaussian peak, C is the baseline offset, which is obtained from the humidity through the GWC model;

[0097] The loss function is optimized by least squares iteration until the residual change rate is less than 10 -6 The iteration is terminated and the baseline, Lorentz amplitude, Lorentz half-maximum width, Gaussian amplitude and Gaussian broadening are obtained as the heavy metal spectral characteristics.

[0098] In this embodiment, the method for obtaining the heavy metal spectrum distribution grid includes:

[0099] Based on the detection area grid corresponding to the heavy metal spectral characteristics, the grid coordinates and sampling timestamps are extracted to form a feature data set. The feature data set is spatially bilinearly interpolated according to the side length of the grid unit. The minimum sampling time interval of the feature data set is used as the time axis, and the feature data set is temporally bilinearly interpolated according to the time axis. The mean absolute error and root mean square error between the interpolation result and the original sampling point are calculated. If the mean absolute error is greater than 4%, the bilinear interpolation coefficient is optimized using the covariance function of Kriging interpolation and then recalculated. If the mean absolute error is less than or equal to 4%, the interpolation result is retained to obtain the heavy metal spectral distribution grid.

[0100] In this embodiment, the method for obtaining the heavy metal quantitative results and heavy metal distribution results includes:

[0101] Based on the heavy metal spectral distribution grid, the heavy metal concentration of the corresponding coordinate is used as a label. The heavy metal spectral distribution grid is flattened from a three-dimensional tensor to a two-dimensional feature matrix to obtain a spatiotemporal characteristic concentration dataset. The spatiotemporal characteristic concentration dataset is Z-score normalized according to its mean and standard deviation. The normalized spatiotemporal characteristic concentration dataset is divided into a training set, a validation set, and a test set with a ratio of 7:2:1.

[0102] The training set was input into a 4-layer multilayer perceptron. The spatiotemporal feature dimensions of the training set data were used as the number of neurons in the input layer. The mean square error was used as the loss function for training. If the mean square error of the validation set did not decrease for 10 consecutive rounds during training, the training was completed in advance. The multilayer perceptron quantitative model was obtained. The determination coefficient was obtained based on the predicted concentration, the mean of the true concentration of the training set, and the total number of continuous spatiotemporal grids of the test set. The generalization ability of the multilayer perceptron quantitative model was optimized based on the determination coefficient. The complete continuous spatiotemporal grid was forward propagated according to the multilayer perceptron quantitative model, and a three-dimensional heavy metal concentration prediction grid was output to obtain the heavy metal quantitative results and heavy metal distribution results of the detection area.

[0103] The second aspect of the present invention further provides a system for quantitative analysis of heavy metals in soil, comprising:

[0104] Data acquisition module: collects soil sample data for preprocessing, including in-situ spectral data, auxiliary data, spatial parameters, heavy metal concentrations, and pollution source coordinates;

[0105] Region division module: used to map the sampling coordinates of the in-situ spectral data to the detection area grid through spatial parameters, obtain the significance of the difference in heavy metal concentration based on the distance between the detection area grid and the pollution source coordinates, and divide the high-difference area and low-difference area according to the significance of the difference;

[0106] Scattering correction module: used to perform particle size analysis on the particle size distribution and humidity of the auxiliary data in the high-difference area, obtain the heavy metal spectrum scattering characteristics of the soil, and use the heavy metal spectrum scattering characteristics to perform scattering correction on the characteristic peaks of the in-situ spectrum data;

[0107] A heavy metal spectrum distribution grid acquisition module is used to reversely simulate heavy metal scattering artifacts based on the scattering correction spectrum and the in-situ spectrum data of the low-difference area, perform artifact repair on the scattering correction spectrum and the in-situ spectrum data based on the simulated heavy metal scattering artifacts through Mie scattering constraints, extract heavy metal spectrum features of the artifact repair spectrum based on the detection area, perform bilinear interpolation on the heavy metal spectrum features, and obtain a heavy metal spectrum distribution grid;

[0108] Quantitative analysis module: used to input the heavy metal spectrum distribution grid and heavy metal concentration into the multi-layer perceptron quantitative model to obtain the heavy metal quantitative results and heavy metal distribution results of the detection area.

[0109] The above content is merely an example and explanation of the structure of the present invention. Those skilled in the art may make various modifications or additions to the described specific embodiments or replace them in a similar manner. As long as they do not deviate from the structure of the invention or exceed the scope defined by the claims, they should all fall within the scope of protection of the present invention.

Claims

1. A method for quantitative analysis of heavy metals in soil, characterized in that: The following steps are involved: Collecting soil sample data for preprocessing, the soil sample data including in-situ spectral data, auxiliary data, spatial parameters, heavy metal concentrations, and pollution source coordinates; The sampling coordinates of the in-situ spectral data are mapped to the detection area grid through spatial parameters, the significance of the difference in heavy metal concentration is obtained according to the distance between the detection area grid and the pollution source coordinates, and high difference areas and low difference areas are obtained according to the significance of the difference; performing particle size analysis on the particle size distribution and moisture content of the auxiliary data in the high-difference region to obtain heavy metal spectral scattering characteristics of the soil, and performing scattering correction on characteristic peaks of the in-situ spectral data using the heavy metal spectral scattering characteristics; Reversely simulate heavy metal scattering artifacts based on the scattering correction spectrum and the in-situ spectral data of the low-difference area, perform artifact repair on the scattering correction spectrum and the in-situ spectral data based on the simulated heavy metal scattering artifacts through Mie scattering constraints, extract heavy metal spectral features of the artifact repaired spectrum based on the detection area, perform bilinear interpolation on the heavy metal spectral features, and obtain a heavy metal spectral distribution grid; The heavy metal spectrum distribution grid and heavy metal concentration are input into a multi-layer perceptron quantitative model to obtain heavy metal quantitative results and heavy metal distribution results in the detection area.

2. The method for quantitative analysis of heavy metals in soil according to claim 1, characterized in that: The pretreatment method comprises: LIBS spectral data of soil samples around heavy metal pollution sources are collected as in-situ spectral data. The heavy metal concentration corresponding to the in-situ spectral data is collected. The particle size distribution and humidity of the soil samples are collected as auxiliary data. The coordinates of the pollution source are collected, and the grid units are obtained according to the latitude and longitude of the sampling point coordinates. The sampling point coordinates and grid units are used as spatial parameters.

3. The method for quantitative analysis of heavy metals in soil according to claim 1, characterized in that: The method of mapping the sampling coordinates of the in-situ spectral data to the detection area grid through spatial parameters includes: The Euclidean distance between the center coordinates of the grid cells and the coordinates of the pollution source is calculated based on the spatial parameters. The in-situ spectral data of 3 to 5 soil samples are extracted as grid samples using the plum blossom method. The grid cells are interpolated using co-kriging based on the spectral characteristic peak intensity and Euclidean distance of the grid samples to obtain the continuous spectral distribution within the grid cells. If the mean absolute error between the continuous spectral distribution within a grid cell and the measured value of the grid sample is greater than 10%, 1 to 2 grid samples are added to the grid and re-interpolated. If the mean absolute error is less than or equal to 10%, the continuous spectral distribution within the grid cell is retained to obtain the detection area grid.

4. The method for quantitative analysis of heavy metals in soil according to claim 1, wherein: The method for obtaining the high-difference area and the low-difference area comprises: Based on the detection area grid, the Euclidean distance between the sampling point coordinates and the pollution source coordinates and the spectral characteristic peak intensity of the sampling point are obtained. The data pairs consisting of the Euclidean distance and the spectral characteristic peak intensity are clustered to obtain the near-source cluster and the far-source cluster. The middle value of the near-source cluster and the far-source cluster is used as the distance threshold. The detection area grid is divided into several intervals according to the Euclidean distance. The spectral characteristic peak intensity within the interval is subjected to a one-way variance analysis to obtain the difference significance level values ​​of different intervals. If the difference significance level value is greater than or equal to 7%, no regional division is performed. If the significance level of the difference is less than 7%, the reciprocal of the normalized Euclidean distance from the sampling point to the pollution source is used as the spatial weight. The multi-scale feature fusion value is calculated based on the spatial weight and the spectral characteristic peak intensity and half-height width of the spectral characteristic peak extracted from the in-situ spectral data of the corresponding sampling point. The dynamic threshold is obtained by using the multi-scale feature fusion value and the distance threshold through the dynamic threshold formula. The dynamic threshold formula is: Where τ(D) is the multi-scale feature fusion value F i The dynamic threshold of μ is the multi-scale feature fusion value F obtained by calculating the in-situ spectral data and spatial weight of the sub-sample in the grid unit. i The statistical mean of σ is the multi-scale feature fusion value F i The statistical standard deviation of the data, Φ is the cumulative distribution function of the standard normal distribution, D is the Euclidean distance from the sampling point to the pollution source, μ D is the mean of all D, σ D is the standard deviation of all D, is the average intensity of the spectral characteristic peak of the in-situ spectral data of the near-source region sub-sample, peak is the characteristic emission peak of heavy metal elements, is the mean intensity of the spectral characteristic peak of the in-situ spectral data of the far-source region sub-sample, D1 is the distance threshold between the near-source region cluster and the far-source region cluster, is the mean half-maximum width of the spectral characteristic peak of the in-situ spectral data of the near-source area sample, L is the side length of the grid unit, D2 is the sum of D1 and the side length of the grid unit, and is the transition area; The grid cells whose multi-scale feature fusion values ​​of sampling points in the interval are greater than the dynamic threshold are divided into high-difference areas, and the rest are divided into low-difference areas.

5. The method for quantitative analysis of heavy metals in soil according to claim 1, characterized in that: The scatter correction method comprises: The particle size distribution and moisture content were extracted based on auxiliary data from high-difference areas. The particle size distribution was discretized into an equivalent spherical diameter sequence. The equivalent spherical diameter sequence was classified according to the International Soil Texture Triangle Diagram standard to obtain the median particle size and soil type coefficient. Based on the median particle size and soil type coefficient, the humidity is converted into a dielectric constant using the GWC model. The dielectric constant is used to calculate the real part of the complex refractive index. The scale parameter is obtained using the relative scale of the median particle size and the laser wavelength of the in-situ spectral data of the corresponding sampling point. The Mie scattering amplitude coefficient is calculated through Anning recursion based on the median particle size, dielectric constant, real part approximation of the complex refractive index, and scale parameter, and the spectral scattering characteristics of heavy metals including the median particle size, scale parameter, and amplitude coefficient are obtained. The spectral data of high-difference areas are corrected using the scattering correction formula using the spectral scattering characteristics of heavy metals. The scattering correction formula is: Among them I corr (λ) is the spectrum after scattering correction, I raw (λ) is the spectral data of the first region, I0(λ) is the light intensity of the laser directly irradiating the Spectralon standard white plate after dark field subtraction, I ref (λ) is the spectral response of the standard white plate at the soil measurement position, The equivalent projected area of ​​soil particles obtained based on the median particle size, α is the scale parameter of soil particle size and laser wavelength λ, l is the multipole order of Mie scattering, N is the truncation term, N = α + 4α 1 / 3 +2, a l is the amplitude coefficient of the Mie scattering electric dipole, b l is the amplitude coefficient of the magnetic dipole of Mie scattering.

6. The method for quantitative analysis of heavy metals in soil according to claim 1, characterized in that: The method for obtaining the artifact repair spectrum comprises: Scattering-corrected spectral data are extracted based on the high-difference area. The scattering-corrected spectral data and the in-situ spectral data of the low-difference area are used as the reference spectrum. The scattering correction formula is used to reversely simulate the heavy metal scattering artifacts through the intensity scaling method based on the reference spectrum, median particle size, and humidity. The in-situ spectral data of the high-difference area and the low-difference area are mixed with the simulated heavy metal scattering artifacts and divided into training set, validation set, and test set; The training set is used as input to a Mie scattering constrained autoencoder constructed based on a convolutional encoder, a Mie scattering physical constraint layer, and a deconvolution decoder. The scattering correction formula is inserted into the Mie scattering physical constraint layer according to the differential constraints of automatic differentiation. Reconstruction loss and KL divergence are used as a joint loss for training. A validation set is added during training to control accuracy convergence. The trained autoencoder model is obtained, and the test set is used to optimize the model performance by the mean squared error. The scattering-corrected spectra in the high-difference area and the in-situ spectral data in the low-difference area are input into the autoencoder model. The noise artifacts in the high-difference area and the low-difference area are repaired by forward propagation based on the autoencoder model to obtain the artifact-repaired spectrum.

7. The method for quantitative analysis of heavy metals in soil according to claim 1, characterized in that: The method for obtaining the heavy metal spectral characteristics comprises: Based on the detection area, the artifact repair spectrum is extracted and peak detection is performed through continuous wavelet transform to obtain the characteristic peak, the total number of characteristic peaks and the central wavelength of the characteristic peak. The spectrum mean of the area far from the characteristic peak is used as the baseline initial value. According to the baseline, the spectral intensity at the center of the Lorentz peak is used as the initial value of the amplitude of the Lorentz peak. The typical half-maximum width of the LIBS plasma peak is used as the initial value of the half-maximum width of the Lorentz peak. The spectral intensity at the center of the Gaussian peak is used as the initial value of the amplitude of the Gaussian peak. The standard deviation of the baseline noise near the center of the Gaussian peak is used as the initial value of the broadening of the Gaussian peak. The initial value of the baseline, the initial value of the amplitude of the Lorentz peak, the initial value of the half-height width of the Lorentz peak, the initial value of the amplitude of the Gaussian peak, and the initial value of the broadening of the Gaussian peak are used to perform Lorentz-Gaussian mixture fitting through the loss function formula. The loss function formula is: in is the loss function, which measures the deviation between the fitted model and the true repair spectrum, I recon (λ) is the artifact repair spectrum, M is the total number of characteristic peaks, A i is the amplitude of the i-th Lorentz peak, Γ i is the half-height width of the i-th Lorentz peak, λ 0i is the central wavelength of the i-th characteristic peak, B i is the amplitude of the i-th Gaussian peak, σ i is the standard deviation of the i-th Gaussian peak, C is the baseline offset, which is obtained from the humidity through the GWC model; The loss function is optimized by least squares iteration until the residual change rate is less than 10 -6 The iteration is terminated and the baseline, Lorentz amplitude, Lorentz half-maximum width, Gaussian amplitude and Gaussian broadening are obtained as the heavy metal spectral characteristics.

8. The method for quantitative analysis of heavy metals in soil according to claim 1, characterized in that: The method for obtaining the heavy metal spectrum distribution grid comprises: Based on the detection area grid corresponding to the heavy metal spectral characteristics, the grid coordinates and sampling timestamps are extracted to form a feature data set. The feature data set is spatially bilinearly interpolated according to the side length of the grid unit. The minimum sampling time interval of the feature data set is used as the time axis, and the feature data set is temporally bilinearly interpolated according to the time axis. The mean absolute error and root mean square error between the interpolation result and the original sampling point are calculated. If the mean absolute error is greater than 4%, the bilinear interpolation coefficient is optimized using the covariance function of Kriging interpolation and then recalculated. If the mean absolute error is less than or equal to 4%, the interpolation result is retained to obtain the heavy metal spectral distribution grid.

9. The method for quantitative analysis of heavy metals in soil according to claim 1, characterized in that: The method for obtaining the heavy metal quantitative results and heavy metal distribution results comprises: Based on the heavy metal spectral distribution grid, the heavy metal concentration of the corresponding coordinate is used as a label. The heavy metal spectral distribution grid is flattened from a three-dimensional tensor to a two-dimensional feature matrix to obtain a spatiotemporal characteristic concentration dataset. The spatiotemporal characteristic concentration dataset is Z-score normalized according to its mean and standard deviation. The normalized spatiotemporal characteristic concentration dataset is divided into a training set, a validation set, and a test set with a ratio of 7:2:

1. The training set was input into a 4-layer multilayer perceptron. The spatiotemporal feature dimensions of the training set data were used as the number of neurons in the input layer. The mean square error was used as the loss function for training. If the mean square error of the validation set did not decrease for 10 consecutive rounds during training, the training was completed in advance. The multilayer perceptron quantitative model was obtained. The determination coefficient was obtained based on the predicted concentration, the mean of the true concentration of the training set, and the total number of continuous spatiotemporal grids of the test set. The generalization ability of the multilayer perceptron quantitative model was optimized based on the determination coefficient. The complete continuous spatiotemporal grid was forward propagated according to the multilayer perceptron quantitative model, and a three-dimensional heavy metal concentration prediction grid was output to obtain the heavy metal quantitative results and heavy metal distribution results of the detection area.

10. A quantitative analysis system for heavy metals in soil, used to execute the quantitative analysis method for heavy metals in soil according to any one of claims 1 to 9, characterized in that: The system comprises: Data acquisition module: collects soil sample data for preprocessing, including in-situ spectral data, auxiliary data, spatial parameters, heavy metal concentrations, and pollution source coordinates; Region division module: used to map the sampling coordinates of the in-situ spectral data to the detection area grid through spatial parameters, obtain the significance of the difference in heavy metal concentration based on the distance between the detection area grid and the pollution source coordinates, and divide the high-difference area and low-difference area according to the significance of the difference; Scattering correction module: used to perform particle size analysis on the particle size distribution and humidity of the auxiliary data in the high-difference area, obtain the heavy metal spectrum scattering characteristics of the soil, and use the heavy metal spectrum scattering characteristics to perform scattering correction on the characteristic peaks of the in-situ spectrum data; A heavy metal spectrum distribution grid acquisition module is used to reversely simulate heavy metal scattering artifacts based on the scattering correction spectrum and the in-situ spectrum data of the low-difference area, perform artifact repair on the scattering correction spectrum and the in-situ spectrum data based on the simulated heavy metal scattering artifacts through Mie scattering constraints, extract heavy metal spectrum features of the artifact repair spectrum based on the detection area, perform bilinear interpolation on the heavy metal spectrum features, and obtain a heavy metal spectrum distribution grid; Quantitative analysis module: used to input the heavy metal spectrum distribution grid and heavy metal concentration into the multi-layer perceptron quantitative model to obtain the heavy metal quantitative results and heavy metal distribution results of the detection area.

Images