Agricultural land soil heavy metal priority control evaluation method

Abstract

This invention discloses a method for prioritizing the control of heavy metals in agricultural land soil, comprising: collecting heavy metal concentration data and deep background values ​​from surface soil of agricultural land; fitting the marginal distribution of each heavy metal, and constructing a multidimensional joint probability distribution model based on the AIC / BIC criterion to select a Copula function; performing Monte Carlo sampling to generate concentration simulation samples that maintain the original dependency structure; inputting the samples in parallel into a health risk assessment module and an ecological risk assessment module to calculate the hazard index, total carcinogenic risk, and potential ecological risk index, respectively; statistically analyzing the risk probability distribution and level probability; identifying and ranking priority control factors that simultaneously meet the health and ecological contribution conditions through Spearman rank correlation and standardized regression coefficients; identifying the primary control source by combining APCS-MLR source apportionment; and outputting control recommendations. This invention quantifies the correlation between heavy metals, achieves dual-risk collaborative assessment, identifies priority control factors and sources, and provides a scientific basis for the classification management and source control of agricultural land.

Description

Technical Field

[0001] This invention relates to the field of soil environmental risk assessment and pollution control technology, and in particular to a method for prioritizing the control of heavy metals in agricultural land soil, especially a method for prioritizing the control of heavy metals in agricultural land soil based on Copula-Monte Carlo dual risk coupling. Background Technology

[0002] Heavy metal pollution in agricultural land soil is characterized by its insidious nature, persistence, and ecotoxicological effects. It not only impacts soil ecological functions but can also enter the human body through the food chain, posing a threat to public health. With the rapid advancement of industrialization and urbanization in my country, heavy metal pollution in agricultural land soil has become an increasingly prominent problem, constraining the quality and safety of agricultural products and the sustainable development of agriculture. Therefore, scientifically assessing the ecological and human health risks of heavy metals in agricultural land soil and identifying priority control factors are crucial for implementing classified management of agricultural land and controlling pollution sources.

[0003] Currently, commonly used methods for assessing soil heavy metal risk include the geoaccumulation index (I_geo), the potential ecological risk index (RI), and the human health risk assessment model recommended by the US Environmental Protection Agency (USEPA). The geoaccumulation index determines the degree of pollution by comparing the heavy metal concentration with background levels; it is simple to calculate but only reflects the degree of pollution and cannot characterize ecological effects. The potential ecological risk index, proposed by Swedish scholar Hakanson, comprehensively considers factors such as the synergistic effects of multiple elements, toxicity levels, pollution concentrations, and environmental sensitivity to heavy metal pollution, and has become a widely used method for assessing the quality of heavy metal pollution in sediments and soils. The USEPA-recommended human health risk assessment model calculates hazard quotients and carcinogenic risks through exposure routes such as oral ingestion, skin contact, and inhalation, and is widely used in the field of soil pollution health risk assessment. However, existing methods have the following shortcomings when applied to heavy metal risk assessment in agricultural land soils: First, it ignores the correlations between variables. Traditional Monte Carlo simulations, when assessing risk uncertainty, typically assume that the concentrations of each heavy metal are independent and sampled independently from their respective marginal distributions. In reality, soil heavy metals often exhibit homology or co-migration; for example, Cd and Zn, and Pb and Hg often show significant positive correlations or tail dependencies. Ignoring this dependency structure leads to distorted risk assessment results, especially underestimating the combined risks under extreme pollution scenarios. Existing research has shown that the Copula-MCS integrated framework can significantly enhance the reliability of health risk assessments under multi-metal pollution scenarios.

[0004] Second, risk assessments often lack a single dimension. Existing research tends to focus on a single risk type, namely ecological risk or health risk, lacking a systematic approach that integrates both. Agricultural land, as a land type that simultaneously serves ecological functions and provides agricultural products, requires consideration of its dual impacts on soil organisms and human health. A single-dimensional risk assessment cannot fully reflect the actual severity of heavy metal pollution in agricultural land soil.

[0005] Third, the output results lack decision-making guidance. Existing risk assessment methods typically only output risk levels or risk probabilities, failing to further identify key heavy metal factors that lead to the risk and their main sources, making it difficult to directly guide priority control actions for agricultural land. When faced with large amounts of risk assessment data, decision-makers often struggle to quickly determine the heavy metal elements and pollution sources that should be prioritized for control.

[0006] Fourth, there is a lack of specific methods for agricultural land scenarios. Agricultural land has unique characteristics such as multi-media transport between soil, crops, and the human body, spatial variability in concentration caused by tillage disturbance, and frequent outdoor activities by children. Existing general risk assessment methods do not fully consider these characteristics. For example, exposure parameters often use default values ​​for built-up land or general scenarios, and fail to specifically set parameters such as soil intake rates for children (under 6 years old) on agricultural land.

[0007] In recent years, some studies have attempted to combine Copula functions with Monte Carlo simulations (Copula-MCS) to quantify uncertainties in soil heavy metal health risk assessment. However, this method remains limited to single health risks and is not coupled with ecological risks, nor does it output priority control factors. Furthermore, existing methods also have shortcomings in source identification, failing to effectively link risk contribution with pollution source analysis, making it difficult to achieve full-chain management from "risk identification" to "source control."

[0008] Therefore, developing a method for agricultural land scenarios that can simultaneously assess both ecological and health risks, quantify the correlation of variables, and output priority control factors has significant practical implications and industrial value. Summary of the Invention

[0009] To address the shortcomings of existing technologies, the purpose of this invention is to provide a method for prioritizing the control of heavy metals in agricultural land soil. This method solves the problems of neglecting the correlation between heavy metals, having a single dimension of risk assessment, weak decision-oriented output results, and lack of adaptability to agricultural land scenarios. This invention captures the dependency structure between heavy metal concentrations using a Copula function and performs Monte Carlo sampling based on joint distributions. It incorporates health and ecological risks into the same probabilistic assessment framework for parallel calculation. Through sensitivity analysis and contribution rate decomposition, it quantitatively identifies priority control factors and combines them with a source apportionment model to identify primary control sources. This achieves accurate, efficient, and comprehensive priority control assessment of heavy metals in agricultural land soil.

[0010] The above-mentioned objective of this invention is achieved through the following technical solutions: A method for prioritizing the control of heavy metals in agricultural land soil includes the following steps: Step S1: Data Acquisition and Preprocessing. Obtain surface soil samples from the farmland plot to be evaluated and determine the concentration data of at least eight heavy metals, including As, Cd, Cu, Pb, Hg, Ni, Cr, and Zn, C=c1,c2,…,c n n is the number of samples; simultaneously, the average concentrations of various heavy metals in the deep soil of this plot are obtained as background values ​​B=b1,b2,…,b m m is the number of heavy metal species; Step S2: Construct a joint probability distribution model for heavy metal concentrations. Fit marginal distributions to the concentration data of each heavy metal. Use the Kolmogorov-Smirnov test to select the optimal marginal distribution type Fⱼ(cⱼ) from normal, log-normal, gamma, and Weibull distributions, and estimate the distribution parameters. Based on the AIC and BIC criteria, select the optimal Copula function C_θ from multiple Copula function families. Combine the marginal distributions of each heavy metal with the optimal Copula function to construct a multidimensional joint probability distribution model. Step S3: Monte Carlo sampling based on joint distribution, set the number of Monte Carlo simulations N, N≥10000; according to the joint probability distribution model, independently generate N groups of concentration simulation samples C⁽ᵏ⁾=(c1⁽ᵏ⁾,c2⁽ᵏ⁾,…,c8⁽ᵏ⁾),k=1,2,…,N, and the heavy metal concentrations in each group of samples maintain the same dependency structure as the original data; Step S4: Parallel calculation of dual risks. The N sets of concentration simulation samples generated in step S3 are input into the health risk assessment module and the ecological risk assessment module in parallel to calculate the health risk indicators and ecological risk indicators corresponding to each set of concentration samples. Step S5: Risk probability statistics and characterization. Statistically calculate the mean, median, standard deviation, and quantile of each risk indicator in N simulations, generate the probability density distribution map and cumulative distribution curve of each risk indicator, and calculate the probability of each risk indicator being at different risk levels according to the preset risk level classification standard. Step S6: Priority control factor identification. Using the Spearman rank correlation coefficient method and standardized regression coefficient method, the correlation coefficient and contribution ratio between the concentration of each heavy metal and each risk indicator are calculated, and a risk contribution matrix is ​​constructed. Heavy metals that simultaneously meet the health risk contribution conditions and ecological risk contribution conditions are identified as priority control factors and sorted according to the comprehensive contribution index. Step S7: Priority source identification. Based on the contribution rate of each pollution source output by the source apportionment model, the main source that leads to the priority control factor is identified as the primary control source. Step S8: Output control recommendations, output a list of priority control factors and their order, output a comparison chart of the health risk contribution rate and ecological risk contribution rate of each priority control factor, and output risk control recommendations for the agricultural land plot based on the agricultural land soil pollution risk screening value.

[0011] As a further technical solution of the present invention: In step S2, the various Copula function families include Gaussian Copula, t-Copula, Clayton Copula, Gumbel Copula, and Frank Copula; the parameter θ of each Copula function is estimated using the maximum likelihood estimation method, the AIC value and BIC value of each Copula are calculated, and the Copula with the smallest AIC value and BIC value is selected as the optimal Copula function; the multidimensional joint probability distribution model is: F(c1,c2,…,c8)=Cθ(F1(c1),F2(c2),…,F8(c8)).

[0012] As a further technical solution of the present invention: In step S4, the health risk assessment module adopts the USEPA-recommended human health risk assessment model, considering at least three exposure routes: oral ingestion, skin contact, and inhalation; for each concentration sample C⁽ᵏ⁾, the non-carcinogenic hazard index HI⁽ᵏ⁾ and total carcinogenic risk TCR⁽ᵏ⁾ for children and adults are calculated: Wherein, IRsoil is soil ingestion rate, EF is exposure frequency, ED is exposure duration, BW is body weight, AT is average exposure time, SA is skin surface area, AF is skin adhesion coefficient, ABS is skin absorption factor, InhR is respiration rate, PEF is particulate matter emission factor, RfD is reference dose, and SF is carcinogenic slope factor; the above parameters adopt exposure parameters specific to agricultural land scenarios. As a further technical solution of the present invention: in step S4, the ecological risk assessment module adopts the Hakanson potential ecological risk index model; for each group of concentration samples C⁽ᵏ⁾, the potential ecological risk coefficient of each element is calculated. In the formula Let be the toxicity response factor for the i-th heavy metal. The background value obtained in step S1 is used to calculate the comprehensive potential ecological risk index. The values ​​of the toxicity response factors are: Hg=40, Cd=30, As=10, Pb=5, Cu=5, Ni=5, Cr=2, Zn=1.

[0013] As a further technical solution of the present invention: in step S5, the risk level classification standard is as follows: For the hazard index HI: HI<1 is acceptable risk, and HI≥1 is unacceptable risk; For total carcinogenic risk As a negligible risk, As an acceptable risk, This is an unacceptable risk; For the comprehensive potential ecological risk index RI: RI<150 is low risk, 150≤RI<300 is medium risk, 300≤RI<600 is high risk, and RI≥600 is extremely high risk.

[0014] As a further technical solution of the present invention: in step S6, the health risk contribution condition is a standardized regression coefficient. Or Spearman rank correlation coefficient The ecological risk contribution condition is a standardized regression coefficient. Or Spearman rank correlation coefficient The comprehensive contribution index Where α is a preset health risk weighting coefficient, 0≤α≤1. As a further technical solution of the present invention: the health risk weighting coefficient α=0.5.

[0015] As a further technical solution of the present invention: In step S7, the source apportionment model is the APCS-MLR model; specifically, the APCS-MLR model is used to perform source apportionment on the original concentration data to identify x pollution sources and the contribution rate matrix of each source to each heavy metal; for the priority control factor determined in step S6, the pollution source with the highest contribution rate is found and determined as the primary control source.

[0016] As a further technical solution of the present invention: in step S1, the sampling depth of the deep soil is 2.0m, and the sampling depth is determined according to the actual stratum distribution of the plot.

[0017] As a further technical solution of the present invention: in step S8, the risk management recommendations include: recommended priority management factors for encrypted monitoring, recommended types of remediation measures, and recommended categories of pollution sources for priority management.

[0018] In summary, compared with the prior art, the present invention has at least one of the following beneficial technical effects: Quantifying the correlation of variables improves the accuracy of risk assessment. By capturing the dependency structure (especially tail dependency) among heavy metal concentrations using a Copula function and performing Monte Carlo sampling based on the joint distribution, the risk underestimation caused by the traditional independence assumption is resolved. Comparative experiments show that the method identified by this invention has a mean HI value that is approximately 22% higher and a P95 value that is approximately 25% higher than that of the traditional independent sampling method, particularly improving the accuracy of risk assessment under extreme pollution scenarios.

[0019] This invention employs a dual-risk synergistic assessment to comprehensively reflect the impact on farmland. For the first time, it integrates health and ecological risks into the same probabilistic assessment framework, utilizing the same set of simulation samples to calculate HI, TCR, and RI in parallel. This achieves coupled analysis and synergistic comparison of the two risks, filling the gap in existing technologies that only focus on a single risk. The invention can simultaneously output the probability distributions of health and ecological risks, comprehensively reflecting the actual degree of harm caused by heavy metal pollution in farmland soil.

[0020] This system outputs priority control factors to enhance decision-making guidance. Through Spearman rank correlation coefficient analysis and standardized regression coefficient decomposition, it quantitatively identifies the heavy metal elements that contribute most to dual risks and outputs a ranking list. This directly guides priority control actions for farmland, overcoming the limitation of existing technologies that only output risk levels. Decision-makers can directly formulate targeted control measures based on the priority control factor list.

[0021] Optimized for agricultural land scenarios, this approach is highly applicable. It employs exposure parameters specific to agricultural land (such as the distribution of children's soil ingestion rates), background value acquisition methods, and risk classification standards, fully considering multi-media transmission from soil to crops to the human body and the characteristics of children's exposure. The method is more closely aligned with the actual management needs of agricultural land.

[0022] By combining source identification with risk contribution analysis, the APCS-MLR source apportionment model is coupled with risk contribution analysis to identify the main pollution sources leading to priority control factors, providing technical support for the entire "source-risk" chain of control. Complete tracing from priority control factors to priority control sources enables control measures to shift from end-of-pipe treatment to source reduction.

[0023] The results are visualized, facilitating decision-making. Risk probability distribution maps, contribution rate comparison charts, and priority control factor ranking tables are generated. The results are intuitive, enabling ecological and environmental departments and land management departments to quickly understand and apply them, thus lowering the barrier to transforming technological achievements into management decisions.

[0024] In summary, this invention provides a scientific, comprehensive, and practical assessment method for priority control of heavy metals in agricultural land soil, which has significant technological advancements and broad prospects for widespread application. Attached Figure Description

[0025] Figure 1This is a flowchart of the overall process of the method of the present invention. Detailed Implementation

[0026] The technical solutions in the embodiments of this application will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of this application, and not all embodiments. All other embodiments obtained by those skilled in the art based on the embodiments of this application without creative effort are within the scope of protection of this application. Example

[0027] Reference Figure 1 This invention provides a method for prioritizing the management and assessment of heavy metals in agricultural land soil. The specific implementation methods of each step are described in detail below.

[0028] Step S1: Data Acquisition and Preprocessing In step S1, the first step is to determine the scope of the agricultural land parcel to be evaluated. The scope of the parcel is determined based on administrative boundaries, watershed boundaries, or soil type units. There is no limit to the parcel area; it can be a local parcel of less than 1 mu (approximately 0.16 acres) or a contiguous area of ​​agricultural land of more than 1,000 mu (approximately 667 acres).

[0029] For topsoil sampling, either a grid sampling method or a systematic random sampling method is used. The grid sampling method is suitable for plots with regular shapes and large areas; the grid size is determined based on the plot area and the required number of samples, generally set to 50m×50m to 200m×200m. The systematic random sampling method is suitable for plots with irregular shapes or significant differences in the distribution of pollution sources. Sampling depth is 0-20cm for the topsoil layer, with at least 1kg of soil sample collected at each sampling point. After collection, the samples are air-dried, ground, and sieved (using a 100-mesh or 200-mesh nylon sieve) before use.

[0030] Methods for determining heavy metal concentrations include, but are not limited to, inductively coupled plasma mass spectrometry (ICP-MS), inductively coupled plasma optical emission spectrometry (ICP-OES), atomic absorption spectrometry (AAS), or atomic fluorescence spectrometry (AFS). Specific determination methods shall be performed in accordance with the "Technical Specification for Soil Environmental Monitoring" (HJ / T 166) and the corresponding standard analytical methods for heavy metals. The heavy metals to be measured must include at least eight: As, Cd, Cu, Pb, Hg, Ni, Cr, and Zn. Other heavy metals (such as Co, Mn, Sb, etc.) may be added depending on the pollution characteristics of the site.

[0031] There are two ways to obtain the background value of deep soil: (1) Collect deep soil samples in the area around the plot that has not been disturbed by human activities. The sampling depth is generally 2.0m, and the specific depth is determined according to the actual stratum distribution of the plot. The sample should pass through the topsoil layer and the cultivated layer and reach the undisturbed original soil layer. (2) Refer to the regional soil background value survey data and adopt the background value that is consistent with the parent material and soil type of the plot to be evaluated.

[0032] Data preprocessing includes outlier identification and handling, and missing value imputation of the raw concentration data. Outlier identification uses the Grubbs test or Dixon test. For data identified as outliers, the original records must be checked before retention, removal, or correction. Missing value imputation uses multiple imputation or mean substitution to ensure the integrity and reliability of the dataset.

[0033] Step S2: Construct a joint probability distribution model of heavy metal concentrations Step S2 is one of the core steps of this invention, and its purpose is to construct a joint probability distribution model that can characterize the dependence structure between heavy metal concentrations through the Copula function.

[0034] S2.1 Marginal Distribution Fitting Marginal distribution fitting was performed on the concentration data of each heavy metal. Specifically, for the j-th heavy metal (j=1,2,…,m, where m is the number of heavy metal species, m=8 in this embodiment), its n concentration observations cj,1,cj,2,…,cj,n were used as samples, and parameter estimation and goodness-of-fit tests were performed using four theoretical distributions: normal distribution, log-normal distribution, gamma distribution, and Weibull distribution.

[0035] The probability density function of the normal distribution is: ; in For location parameters (mean). This is the scale parameter (standard deviation).

[0036] The probability density function of the log-normal distribution is: ; Where μ is the logarithmic mean and σ>0 is the logarithmic standard deviation.

[0037] The probability density function of the gamma distribution is: ; Where α>0 is the shape parameter, β>0 is the scale parameter, and Γ(α) is the gamma function.

[0038] The probability density function of the Weibull distribution is: ; Where k>0 is the shape parameter and λ>0 is the scale parameter.

[0039] The maximum likelihood estimation method is used to estimate the distribution parameters. Maximum likelihood estimation is achieved by maximizing the likelihood function. To estimate the parameter θ, where Let be the probability density function.

[0040] The goodness-of-fit test was performed using the Kolmogorov-Smirnov test (KS test). The KS test statistic Dn is defined as follows: in Let be the empirical distribution function. Let be the theoretical distribution function. For a significance level α = 0.05, when When the value is less than the critical value, the sample is considered to come from the theoretical distribution. The distribution type with the largest p-value and passing the KS test (p>0.05) is selected as the optimal marginal distribution for the j-th heavy metal. .

[0041] S2.2 Copula Function Optimization The Copula function is a function that connects the joint distribution of a multidimensional random variable with its respective marginal distributions. According to Sklar's theorem, for an m-dimensional random variable... If its marginal distribution function is Then there exists a Copula function. Make: ; In this embodiment, the optimal Copula is selected from the following five commonly used Copula function families: (1) Gaussian Copula: ; Where Φ is the standard normal distribution function, and Φθ is the m-dimensional normal distribution function with a mean of 0 and a covariance matrix of θ. Gaussian Copula is suitable for characterizing symmetric dependencies between variables, but it does not have tail dependencies.

[0042] (2) t-Copula: ; in For degrees of freedom of Distribution function, It is an m-dimensional t-distribution function with a mean of 0, a covariance matrix of θ, and ν degrees of freedom. t-Copula can capture the symmetric tail dependence between variables and is suitable for situations where variables may exhibit cooperative changes at both ends. (3) Clayton Copula: ; Clayton Copula has strong lower tail dependency and is suitable for situations where variables show stronger correlation in the low-value region.

[0043] (4) Gumbel Copula: ; Gumbel Copula has strong upper tail dependency and is suitable for situations where variables show stronger correlation in high-value regions.

[0044] (5) Frank Copula: ; Frank Copula is suitable for characterizing symmetrical relationships with weak tail dependence.

[0045] The parameters of the Copula function are estimated using the maximum likelihood estimation method (a two-stage method). The first stage estimates the parameters of each marginal distribution, and the second stage estimates the parameters θ of the Copula function based on the given marginal distribution parameters. Specifically, the estimation of the Copula parameters is achieved by maximizing the following log-likelihood function: ; The selection of Copula functions is based on the AIC and BIC criteria. AIC (Akaike Information Criterion) is defined as follows: ; Where ℓ(θhat) is the maximum log-likelihood value, and k is the number of parameters in the Copula function. The Bayesian Information Criterion (BIC) is defined as: BIC = −2·ℓ(θhat) + ln(n)·k. For the five Copula functions, calculate their AIC and BIC values ​​respectively, and select the Copula function with the smallest AIC and BIC values ​​as the optimal Copula function. When the optimal Copula selected by the two criteria are inconsistent, prioritize the Copula function with the smallest BIC value, because BIC imposes a stricter penalty on the number of parameters, which helps prevent overfitting.

[0046] S2.3 Constructing a joint probability distribution model The optimal edge distribution of each heavy metal determined in S2.1 The optimal Copula function selected in S2.2 By combining these methods, a multidimensional joint probability distribution model can be constructed: ; This joint probability distribution model fully characterizes the marginal distribution features of each heavy metal concentration variable and the dependency structure between variables, including linear and nonlinear correlations, symmetric and asymmetric dependencies, and upper and lower tail dependencies.

[0047] Step S3: Monte Carlo sampling based on joint distribution The purpose of step S3 is to generate a large number of concentration simulation samples that maintain the original data dependency structure based on the joint probability distribution model constructed in step S2.

[0048] S3.1 Set the number of simulations Set the number of Monte Carlo simulations, N, where N ≥ 10000. The value of the number of simulations, N, determines the stability and accuracy of the risk assessment results. A larger N results in smaller statistical fluctuations in the simulation results, but also increases computation time. Depending on the application scenario, N can be set to 10000, 20000, 50000, or 100000. For regional-scale assessments (involving multiple plots), N=10000 is recommended to balance computational efficiency; for detailed assessments of a single plot, N=50000 is recommended to obtain a more accurate risk probability estimate.

[0049] S3.2 Generate simulated samples N sets of simulated concentration samples are generated independently based on the joint probability distribution model. The specific sampling method used is either conditional sampling or Gibbs sampling. For a Copula-based joint distribution, the sampling steps are as follows: (1) From the Copula function Generate m-dimensional uniformly distributed random vectors For t-Copula, the sampling steps are as follows: ① From a multivariate t-distribution with ν degrees of freedom... ① Generate sample vector y; ② Calculate This yields uniformly distributed samples.

[0050] (2) The simulated values ​​of heavy metal concentrations are obtained by inverse function transformation of each marginal distribution: .

[0051] Repeat the above steps N times to obtain N sets of concentration simulation samples. Because the Copula function captures the complete dependency structure between the original variables, the generated simulated samples are consistent with the original data in terms of the correlations and tail dependencies between variables.

[0052] Step S4: Dual-risk parallel computing Step S4 involves inputting the N sets of concentration simulation samples generated in step S3 into the health risk assessment module and the ecological risk assessment module in parallel.

[0053] S4.1 Health Risk Assessment Module Health risk assessment employed the USEPA-recommended human health risk assessment model, considering three exposure routes: oral ingestion, skin contact, and inhalation. For the recipient population, risks were calculated separately for children (under 6 years old) and adults (over 18 years old), with children being the most vulnerable exposure group in agricultural settings.

[0054] For each group of concentration samples Calculate the following exposure doses: SF is the carcinogenic slope factor (kg·d / mg), and the SF values ​​for different heavy metals are referenced from the USEPA IRIS database. Carcinogenic risk is calculated for heavy metals with carcinogenic effects such as As and Cd, while only non-carcinogenic hazard indices are calculated for non-carcinogenic substances such as Cu and Zn. The exposure parameters are specifically designed for agricultural land scenarios. Specifically, the distribution values ​​of the soil intake rate parameter for children in agricultural land fully consider the characteristics of children's outdoor activities on agricultural land (such as playing in fields and direct contact with soil), and are higher than the default values ​​for built-up land, to reflect the actual exposure level in agricultural land scenarios.

[0055] S4.2 Ecological Risk Assessment Module The ecological risk assessment adopted the Hakanson potential ecological risk index model. This model was proposed by Swedish scholar Hakanson in 1980 and comprehensively considers the toxicity level of heavy metals, the degree of pollution, and environmental sensitivity.

[0056] For each group of concentration samples Calculate the potential ecological risk coefficient of each element. : ; in: For the first The toxicity response factors of heavy metals reflect their toxicity levels, environmental mobility, and bioaccumulation. For the first Simulated concentrations (mg / kg) of several heavy metals; For the first Background values ​​(mg / kg) for each heavy metal.

[0057] The values ​​of the toxicity response factors were determined according to Hakanson's principle: Hg=40, Cd=30, As=10, Pb=5, Cu=5, Ni=5, Cr=2, Zn=1.

[0058] Calculate the comprehensive potential ecological risk index RI (dimensionless): ; For each heavy metal, when hour, This indicates that no significant ecological risk has been generated; when When it is large, A corresponding increase indicates a rise in ecological risk.

[0059] Step S5: Risk Probability Statistics and Characterization The purpose of step S5 is to perform statistical analysis on the risk indicators obtained from N simulations and output the probability distribution characteristics of the risk.

[0060] S5.1 Calculation of Statistical Characteristic Values In the Nth simulation The following statistical characteristics: Mean: ; Median: The middle value among N simulated values ​​after sorting them. When N is even, the median is the average of the two middle values.

[0061] Standard deviation (SD): ; Quantile values: After sorting the N simulated values ​​from smallest to largest, P5 is the 5th quantile, P50 is the 50th quantile (i.e., the median), and P95 is the 95th quantile.

[0062] S5.2 Probability Distribution Generation: This section generates probability density distribution plots and cumulative distribution curves for each risk indicator. The probability density distribution is generated using kernel density estimation, with a Gaussian kernel as the kernel function and the bandwidth determined by the Silverman rule. The cumulative distribution curve is plotted directly from the sorted simulated values.

[0063] S5.3 Risk Level Probability Statistics Based on the preset risk level classification criteria, calculate the probability that each risk indicator falls into a different risk level: For the hazard index HI: HI < 1 indicates acceptable risk, meaning the non-carcinogenic effects have not exceeded the threshold; HI ≥ 1 indicates unacceptable risk, meaning there may be non-carcinogenic health effects. The percentage of occurrences of HI ≥ 1 is used as the probability of unacceptable risk.

[0064] For total carcinogenic risk TCR: TCR < 1 × 10-6 Negligible risk means a carcinogenic risk of less than one in a million, which is generally considered acceptable; 1×10 -6 ≤TCR<1×10 -4 "Acceptable risk" indicates that the carcinogenic risk is within an acceptable range but requires attention; TCR ≥ 1 × 10 -4 An unacceptable risk indicates a carcinogenic risk exceeding one in ten thousand, requiring risk management measures. The percentage of occurrence for each risk level is statistically analyzed. For the Comprehensive Potential Ecological Risk Index (RI): RI < 150 is low risk; 150 ≤ RI < 300 is medium risk; 300 ≤ RI < 600 is high risk; and RI ≥ 600 is extremely high risk. The percentage of occurrence for each risk level is statistically analyzed.

[0065] Step S6: Identification of Priority Control Factors Step S6 is another core step of the present invention, the purpose of which is to quantitatively identify the heavy metal factors that contribute the most to the dual risk.

[0066] S6.1 Spearman's Rank Correlation Coefficient Method The Spearman rank correlation coefficient method was used to calculate the correlation between the concentration of each heavy metal and the risk index. The Spearman rank correlation coefficient is a nonparametric statistic used to measure the strength of the monotonic correlation between two variables, and does not require the variables to follow a normal distribution.

[0067] For the concentration of the i-th heavy metal and the HI output, the Spearman rank correlation coefficient is... The calculation formula is: ; in For the first The rank difference between the heavy metal concentration and HI in the group samples after sorting.

[0068] Similarly, the Spearman rank correlation coefficients between each heavy metal and RI were calculated. .

[0069] S6.2 Standardized Regression Coefficient Method The contribution ratio of each heavy metal to the total variation of risk indicators was estimated using the standardized regression coefficient method. The specific steps are as follows: (1) Establish a multiple linear regression model with the concentration of each heavy metal as the independent variable and the risk indicator as the dependent variable: ; Where X is a risk indicator (HI or RI). ε is the regression coefficient, and ε is the residual. (2) Standardize the independent and dependent variables respectively (subtract the mean and divide by the standard deviation) to obtain the standardized regression coefficients. .

[0070] (3) The contribution ratio of each heavy metal to the total variation of risk indicators is calculated as follows: ; S6.3 Priority Control Factor Determination Construct health risk contribution matrices and ecological risk contribution matrices. For each heavy metal i, examine its contribution to both health and ecological risks: A heavy metal is designated as a priority control factor only if it simultaneously meets both of the above conditions. The thresholds are chosen based on the following considerations: a contribution ratio exceeding 10% indicates that the heavy metal has a significant impact on risk; a correlation coefficient exceeding 0.5 indicates a moderate to strong correlation between the heavy metal concentration and the risk indicator. The dual conditions (contribution ratio or correlation coefficient) ensure the robustness of priority control factor identification, considering both linear contributions and non-linear correlations.

[0071] S6.4 Ranking by Comprehensive Contribution Index The selected priority control factors are ranked according to their comprehensive contribution index. (Comprehensive Contribution Index) Defined as the weighted sum of health risk contribution and ecological risk contribution: ; Here, α is a preset health risk weighting coefficient, with a value ranging from 0 ≤ α ≤ 1. The value of α reflects the relative importance that decision-makers attach to health risks and ecological risks: α > 0.5 indicates a greater emphasis on health risks, and α < 0.5 indicates a greater emphasis on ecological risks. In the context of agricultural land, since agricultural land simultaneously undertakes agricultural product supply and ecological functions, α = 0.5 is recommended, meaning that health risks and ecological risks are equally important.

[0072] Step S7: Prioritize source identification Step S7 connects risk contribution with pollution source analysis to achieve complete traceability from "priority control factors" to "priority control sources".

[0073] S7.1 APCS-MLR Source Resolution The APCS-MLR (Absolute Principal Component Score-Multiple Linear Regression) model was used for source apportionment of the original concentration data. The APCS-MLR model is a source apportionment method developed based on principal component analysis (PCA). It does not require prior assumptions about the number and spectrum of pollution sources and can objectively identify pollution sources and their contributions.

[0074] The specific steps of the APCS-MLR model are as follows: (1) Perform principal component analysis (PCA) on the heavy metal concentration data, extract principal components with eigenvalues ​​greater than 1, and obtain the factor loading matrix and factor score of each principal component.

[0075] (2) Calculate the absolute principal component score (APCS) by converting the factor scores into absolute scores, eliminating the negative values ​​introduced by PCA rotation and centering, so that each principal component score represents the actual pollution source contribution.

[0076] (3) A multiple linear regression model was established with the concentration of each heavy metal as the dependent variable and each APCS as the independent variable. The regression coefficients are the contribution rates of each pollution source to the heavy metal. The intercept term represents the contribution of unidentified sources (natural background).

[0077] The APCS-MLR model was used to identify x pollution sources and the contribution rate matrix of each source to each heavy metal.

[0078] S7.2 Identification of Primary Control Sources For the priority control factors determined in step S6, the pollution source with the highest contribution rate is identified. Specifically, for each heavy metal in the priority control factor list, the contribution rate of each source for that heavy metal is found in the contribution rate matrix obtained in S7.1, and the pollution source with the highest contribution rate is identified as the main source of that factor. When the main sources of multiple priority control factors point to the same pollution source, that pollution source is identified as the primary control source.

[0079] Step S8: Output Control Recommendations S8.1 Output priority control factor list The list of priority control factors, ranked from highest to lowest based on their comprehensive contribution index, is output, along with their health risk contribution rate. .

[0080] S8.2 Output Contribution Rate Comparison Chart Output a comparison chart of the health risk contribution rate and ecological risk contribution rate of each priority control factor, presented in the form of a two-sided bar chart or a grouped bar chart.

[0081] S8.3 Output risk management recommendations Based on the agricultural land soil pollution risk screening values ​​(GB 15618-2018), risk management recommendations are provided for this agricultural land plot. Specifically, these include: (1) Recommended priority control factors for intensified monitoring: For heavy metals in the priority control factor list, it is recommended to increase the sampling density within the plot for detailed investigation in order to accurately characterize their spatial distribution characteristics.

[0082] (2) Recommended remediation measures: Based on the soil chemical behavior of the priority control factors and the suitability of remediation technologies, appropriate remediation measures are recommended respectively. For example, passivation remediation or phytoremediation is recommended for Cd pollution, stabilization remediation or thermal desorption remediation is recommended for Hg pollution, and phosphate stabilization remediation is recommended for Pb pollution.

[0083] (3) Recommended pollution source categories for priority control: The primary control sources identified in step S7 shall be designated as priority pollution source categories for priority control, and targeted source control recommendations shall be proposed, including pollution source investigation, emission supervision, process improvement and other measures. Example

[0084] This embodiment uses a typical agricultural plot in southern my country as an example to apply the method of this invention for priority control assessment of heavy metals in soil. The plot covers an area of ​​50 mu (approximately 3.3 hectares), mainly used for rice cultivation, and is surrounded by small industrial and mining enterprises.

[0085] Data acquisition and preprocessing (step S1) Following a systematic random sampling method, 20 sampling points were established within the site to collect surface soil samples (0-20 cm). After air-drying, grinding, and sieving (100 mesh), the concentrations of eight heavy metals (As, Cd, Cu, Pb, Hg, Ni, Cr, and Zn) were determined using inductively coupled plasma mass spectrometry (ICP-MS). Simultaneously, deep soil samples (2.0 m depth) were collected from undisturbed areas surrounding the site to determine the background values ​​of each heavy metal. The results are shown in Table 1.

[0086] heavy metal Concentration range average concentration Standard deviation Background value As 8.5-25.6 15.2 4.3 10.0 Cd 0.15-1.80 0.62 0.42 0.12 Cu 18.5-85.3 42.6 16.8 22.0 Pb 25.6-120.5 58.3 22.5 26.0 Hg 0.05-0.85 0.28 0.19 0.06 Ni 12.5-45.6 26.8 8.2 25.0 Cr 45.6-128.3 72.5 20.1 65.0 Zn 45.8-280.5 120.6 58.4 70.0 As can be seen from Table 1, the average concentrations of heavy metals such as Cd, Pb, Hg, and Zn are significantly higher than the background values, indicating that there is a certain degree of heavy metal enrichment in this plot.

[0087] Construct a joint probability distribution model for heavy metal concentrations (step S2) S2.1 Marginal Distribution Fitting For the concentration data of each heavy metal, normal, log-normal, gamma, and Weibull distributions were used for fitting, and the optimal distribution was selected by Kolmogorov-Smirnov test. The goodness-of-fit test results of the marginal distributions of each heavy metal are shown in Table 2.

[0088] heavy metal Optimal marginal distribution type KS test p value Goodness-of-fit evaluation As Log-normal distribution 0.325 Passed (p>0.05) Cd Log-normal distribution 0.287 Passed (p>0.05) Cu Gamma distribution 0.412 Passed (p>0.05) Pb Weibull distribution 0.198 Passed (p>0.05) Hg Log-normal distribution 0.356 Passed (p>0.05) Ni normal distribution 0.452 Passed (p>0.05) Cr normal distribution 0.398 Passed (p>0.05) Zn Gamma distribution 0.265 Passed (p>0.05) As shown in Table 2, the optimal marginal distributions of all heavy metals passed the KS test (p values ​​were all greater than 0.05), indicating that the selected distribution types can fit the empirical distributions of each heavy metal concentration well.

[0089] S2.2 Copula Function Optimization The parameters of Gaussian Copula, t-Copula, Clayton Copula, Gumbel Copula, and Frank Copula were estimated using the maximum likelihood estimation method, and the AIC and BIC values ​​of each Copula were calculated. The calculation results are shown in Table 3.

[0090] Copula function AIC value BIC value Optimal Results Gaussian Copula -156.3 -142.8 - t-Copula -168.7 -152.6 Optimal Clayton Copula -98.5 -88.3 - Gumbel Copula -112.4 -102.2 - Frank Copula -105.6 -95.4 - As shown in Table 3, the t-Copula has the smallest AIC and BIC values, therefore, the t-Copula is selected as the optimal Copula function. The t-Copula can simultaneously capture the symmetric tail dependence between variables, making it suitable for situations where soil heavy metal concentrations may exhibit synergistic changes at both high and low ends.

[0091] S2.3 Constructing a joint probability distribution model By combining the marginal distributions of each heavy metal with t-Copula, an 8-dimensional joint probability distribution model is constructed: ; Monte Carlo sampling based on joint distribution (step S3) The number of Monte Carlo simulations was set to N=50000. Based on the joint probability distribution model described above, 50000 sets of concentration simulation samples were generated. To verify the ability of the joint distribution model to capture the dependency structure between variables, the Spearman rank correlation coefficients among the key heavy metals in the simulation samples were calculated, and the results are shown in Table 4.

[0092] Variable pairs correlation coefficient of raw data Simulated sample correlation coefficient deviation Cd-Zn 0.68 0.66 -0.02 Pb-Hg 0.72 0.70 -0.02 As-Cu 0.45 0.43 -0.02 Cd-Pb 0.52 0.50 -0.02 Hg-Zn 0.38 0.36 -0.02 As can be seen from Table 4, the correlation coefficients between the simulated samples and the original data on each key variable pair have very small deviations (all less than 0.03), indicating that the constructed joint distribution model captures the dependencies between variables well, especially preserving the significant positive correlation structure between Cd and Zn, and Pb and Hg.

[0093] Dual-risk parallel computation (step S4) S4.1 Health Risk Assessment The USEPA-recommended human health risk assessment model was adopted, considering three exposure routes: oral ingestion, skin contact, and inhalation. Exposure parameters specific to agricultural settings were used: soil ingestion rate for children was 200 mg / day (range 100-400 mg / day), exposure frequency was 350 days / year, exposure duration was 6 years, body weight was 15 kg, and mean exposure time was ED × 365 days for non-carcinogens and 70 years × 365 days for carcinogens. Other parameters were based on USEPA's recommended values ​​in the exposure parameter manual.

[0094] For each concentration sample, the non-carcinogenic hazard index (HI) and total carcinogenic risk (TCR) for children were calculated. The probability distributions of HI and TCR were obtained through 50,000 simulations.

[0095] S4.2 Ecological Risk Assessment The Hakanson potential ecorisk index model was used. The toxicity response factors were set as follows: Hg=40, Cd=30, As=10, Pb=5, Cu=5, Ni=5, Cr=2, Zn=1. For each concentration sample, the potential ecorisk coefficient for each element was calculated. And the comprehensive potential ecological risk index RI.

[0096] Risk probability statistics and characterization (step S5) The statistical characteristic values ​​of HI, TCR and RI in 50,000 simulations were statistically analyzed, and the results are shown in Table 5.

[0097] Risk indicators mean median Standard deviation P5 P50 P95 Skewness coefficient Distribution pattern HI (Children) 1.85 1.62 0.92 0.68 1.62 3.56 1.25 Significant right deviation TCR (Children) 2.35 x 10 -4 ]]> 1.98 x 10 -4 ]] <![CDATA[1.56×10 -4 ]]> <![CDATA[5.62×10 -5 ]]> <![CDATA[1.98×10 -4 ]]> <![CDATA[5.28×10 -4 ]]> 1.38 Significant right deviation RI 285.6 256.3 142.5 98.5 256.3 568.2 1.42 Significant right deviation As shown in Table 5, the skewness coefficients of the three risk indicators are all greater than 1, exhibiting a significantly right-skewed distribution, and the mean values ​​are all greater than the median. This indicates that there are a few high-pollution scenarios that lead to a higher risk, highlighting the necessity of using probabilistic risk assessment rather than deterministic assessment.

[0098] Calculate the probability of each risk level based on the risk level classification criteria: The probability of HI≥1 (unacceptable risk) is 72.5%.

[0099] TCR≥1×10 -4 The probability of (unacceptable risk) is 58.3%; 1×10 -6 ≤TCR<1×10 -4 The probability of (acceptable risk) is 35.2%; TCR < 1×10 -6 The probability of (negligible risk) is 6.5%.

[0100] The probability of RI≥600 (extremely high risk) is 8.2%; the probability of 300≤RI<600 (high risk) is 32.6%; the probability of 150≤RI<300 (medium risk) is 42.8%; and the probability of RI<150 (low risk) is 16.4%.

[0101] The above results indicate that the site poses high health and ecological risks and requires priority control measures.

[0102] Priority control factor identification (step S6) The Spearman rank correlation coefficient method and the standardized regression coefficient method were used to calculate the correlation coefficient and contribution ratio of each heavy metal with the risk indicator, and the results are shown in Table 6.

[0103] Table 6. Contribution ratio of Spearman correlation coefficient and standardized regression coefficient of each heavy metal to risk indicators Based on the criteria for priority control factors (β>10% or ρ>0.5), Cd, Pb, and Hg are the heavy metals that simultaneously meet both the health risk contribution and ecological risk contribution criteria. Therefore, Cd, Pb, and Hg are determined to be the priority control factors for this site.

[0104] Calculate the overall contribution index of each priority control factor (α=0.5): The ranking of the comprehensive contribution index of each priority control factor is shown in Table 7.

[0105] Ranking Priority control factors Health risk contribution rate (β_HI) Ecological risk contribution rate (β_RI) Overall Contribution Index (W) 1 Cd 35.2% 28.6% 31.9% 2 Hg 22.5% 38.2% 30.4% 3 Pb 12.6% 10.8% 11.7% The ranking results are: Cd (31.9%) > Hg (30.4%) > Pb (11.7%). Among them, the combined contribution index of Cd and Hg is significantly higher than that of Pb, and they should be the heavy metals with the highest priority for control.

[0106] Prioritize source identification for control (step S7) The APCS-MLR model was used for source apportionment of the raw concentration data. First, principal component analysis was performed on the concentration data of the eight heavy metals, extracting principal components with eigenvalues ​​greater than 1, resulting in three principal components with a cumulative variance contribution rate of 82.5%. Then, the absolute principal component scores were calculated, and multiple linear regression was performed with each heavy metal concentration as the dependent variable and the absolute principal component scores as the independent variables to obtain the contribution rate of each pollution source to each heavy metal. The source apportionment results are shown in Table 8.

[0107] Types of pollution sources Main contributing elements Contribution rate Cumulative contribution rate Source Feature Description Source 1 (Industrial Source) Cd (52.3%), Pb (45.6%), Hg (48.2%), Zn (38.5%) 46.2% 46.2% Non-ferrous metal smelting / electroplating wastewater and waste gas Source 2 (Agricultural Source) As (42.8%), Cu (35.2%), Zn (28.6%) 28.5% 74.7% Fertilizer and pesticide application / Livestock and poultry manure Source 3 (Natural Source) Cr (65.2%), Ni (58.3%), As (32.5%) 25.3% 100% Parent material / rock weathering For the priority control factors Cd, Pb, and Hg determined in step S6, identify the pollution sources with the highest contribution rates: The main source of Cd is Source 1 (industrial source), contributing 52.3%; The main source of Pb is Source 1 (industrial source), contributing 45.6%; The main source of Hg is source 1 (industrial source), contributing 48.2%.

[0108] Therefore, Source 1 (industrial source) was identified as the primary source for control.

[0109] Output control recommendations (step S8) Based on the agricultural land soil pollution risk screening values ​​(GB 15618-2018), risk management recommendations for this plot of land are output and summarized in Table 9.

[0110] Control Dimensions Specific content Basis / Explanation Priority control factors ranking 1. Cd (31.9%); 2. Hg (30.4%); 3. Pb (11.7%) Ranking by Comprehensive Contribution Index W Recommended encryption monitoring factors Cd, Hg, Pb Based on the aforementioned 20 sampling points, sampling points were further densified in areas with high Cd and Hg content. Recommended repair measures type Cd: Apply passivating materials such as lime and biochar to reduce bioavailability, or plant black nightshade / sedum for phytoremediation; Hg: Add stabilizing materials such as sulfides for stabilization remediation; Pb: Take phosphate stabilization remediation measures. Soil chemical behavior and remediation technology adaptability based on different heavy metals It is recommended to prioritize source control. Primary source to be controlled: Industrial sources (Source 1) It has the highest contribution rate (>45%) to Cd, Pb, and Hg. Specific actions for source control Investigate the pollution discharge situation of surrounding industrial and mining enterprises, urge them to improve wastewater treatment facilities and exhaust gas dust removal facilities; carry out soil remediation in areas where pollution has already occurred. Cut off the input of heavy metals from industrial sources Long-term planning Establish a soil environmental quality monitoring and early warning mechanism and conduct regular risk assessments. Achieve dynamic control and risk early warning To verify the technical effectiveness of this invention, a comparative experiment was conducted. Risk assessments were performed on the same plot of land using the following three methods: Method A (the method of this invention): a priority control assessment method based on Copula-Monte Carlo dual risk coupling.

[0111] Method B (Traditional Monte Carlo Independent Sampling): Assuming that the concentrations of each heavy metal are independent of each other, Monte Carlo simulation is performed by independently sampling from their respective marginal distributions. Only health risks are calculated, without ecological risk assessment or identification of priority control factors.

[0112] Method C (Deterministic Risk Assessment): The average concentration of each heavy metal is substituted into the risk assessment model to calculate a single risk value without performing uncertainty analysis.

[0113] The comparison results are shown in Table 10.

[0114] Comparison indicators Method A (This invention) Method B (Traditional MCS) Method C (Deterministic Approach) Mean HI (children) 1.85 1.52 1.38 HI P95 value 3.56 2.85 — The probability of HI≥1 72.5% 58.6% — Should variable correlation be considered? Yes (Copula joint distribution) No (independent sampling) no Ecological risks assessed? Yes (RI probability distribution) no no Have priority control factors been identified? It is (Cd, Hg, Pb) no no Is the priority control source identified? It is (an industrial source). no no As can be seen from Table 10: First, the mean HI and P95 values ​​calculated by Method A are higher than those calculated by Method B. This is because Method A captures the positive correlations (especially tail dependence) between heavy metals through the Copula function, more realistically reflecting the probability of multiple heavy metals being at high concentrations simultaneously under extreme pollution scenarios, thus avoiding an underestimation of the combined risk. In contrast, Method B assumes that each heavy metal is independent, ignoring this synergistic relationship, leading to a lower risk assessment result.

[0115] Second, Method A assesses both health and ecological risks, while Methods B and C only assess health risks, failing to fully reflect the actual hazards of heavy metal pollution in agricultural land soil. In this embodiment, Method A identifies RI≥600 (extremely high risk) with a probability of 8.2%, information that is of significant reference value for the classification and management of agricultural land.

[0116] Third, Method A identified three priority control factors—Cd, Hg, and Pb—and their ranking, and further identified industrial sources as the primary control sources, providing decision-makers with clear control targets. Methods B and C, however, only output risk values ​​and cannot guide specific priority control actions.

[0117] Fourth, Method A provides information on the probability distribution of risk (e.g., the probability of HI≥1 is 72.5%), while Method C only gives a single definite value (HI=1.38), which cannot reflect the uncertainty of risk and may lead to misjudgment of the risk level.

[0118] In summary, the method of this invention is significantly superior to existing technologies in terms of the accuracy, comprehensiveness, and decision-making guidance of risk assessment.

[0119] The embodiments described herein are preferred embodiments of the present invention and are not intended to limit the scope of protection of the present invention. Therefore, all equivalent changes made in accordance with the structure, shape, and principle of the present invention should be covered within the scope of protection of the present invention.

Claims

1. A method for assessing priority control of heavy metals in agricultural land soil, characterized in that, Includes the following steps: Step S1: Data acquisition and preprocessing. Obtain surface soil samples from the farmland plot to be evaluated and determine the concentration data of at least eight heavy metals, including As, Cd, Cu, Pb, Hg, Ni, Cr, and Zn. n is the number of samples; simultaneously, the average concentrations of various heavy metals in the deep soil of this plot are obtained as background values. m is the number of heavy metal species; Step S2: Construct a joint probability distribution model for heavy metal concentrations. Fit marginal distributions to the concentration data of each heavy metal. Use the Kolmogorov-Smirnov test to select the optimal marginal distribution type Fⱼ(cⱼ) from normal, log-normal, gamma, and Weibull distributions, and estimate the distribution parameters. Based on the AIC and BIC criteria, select the optimal Copula function C_θ from multiple Copula function families. Combine the marginal distributions of each heavy metal with the optimal Copula function to construct a multidimensional joint probability distribution model. Step S3: Monte Carlo sampling based on joint distribution, setting the number of Monte Carlo simulations N, N≥10000; generating N sets of concentration simulation samples independently according to the joint probability distribution model. The concentrations of each heavy metal in each sample group maintained the same dependency structure as the original data; Step S4: Parallel calculation of dual risks. The N sets of concentration simulation samples generated in step S3 are input into the health risk assessment module and the ecological risk assessment module in parallel to calculate the health risk indicators and ecological risk indicators corresponding to each set of concentration samples. Step S5: Risk probability statistics and characterization. Statistically calculate the mean, median, standard deviation, and quantile of each risk indicator in N simulations, generate the probability density distribution map and cumulative distribution curve of each risk indicator, and calculate the probability of each risk indicator being at different risk levels according to the preset risk level classification standard. Step S6: Priority control factor identification. Using the Spearman rank correlation coefficient method and standardized regression coefficient method, the correlation coefficient and contribution ratio between the concentration of each heavy metal and each risk indicator are calculated, and a risk contribution matrix is ​​constructed. Heavy metals that simultaneously meet the health risk contribution conditions and ecological risk contribution conditions are identified as priority control factors and sorted according to the comprehensive contribution index. Step S7: Priority source identification. Based on the contribution rate of each pollution source output by the source apportionment model, the main source that leads to the priority control factor is identified as the primary control source. Step S8: Output control recommendations, output a list of priority control factors and their order, output a comparison chart of the health risk contribution rate and ecological risk contribution rate of each priority control factor, and output risk control recommendations for the agricultural land plot based on the agricultural land soil pollution risk screening value.

2. The method for priority control and assessment of heavy metals in agricultural land soil according to claim 1, characterized in that, In step S2, the various Copula function families include Gaussian Copula, t-Copula, Clayton Copula, Gumbel Copula, and Frank Copula; the parameter θ of each Copula function is estimated using the maximum likelihood estimation method, and the AIC and BIC values ​​of each Copula are calculated. The Copula function with the smallest AIC and BIC values ​​is selected as the optimal Copula function; the multidimensional joint probability distribution model is as follows: .

3. The method for priority control and assessment of heavy metals in agricultural land soil according to claim 1, characterized in that, In step S4, the health risk assessment module adopts the USEPA-recommended human health risk assessment model, considering at least three exposure routes: oral ingestion, skin contact, and inhalation; for each concentration sample... Calculate the non-carcinogenic hazard index for children and adults Total carcinogenic risk : in, For soil uptake rate, For exposure frequency, For the duration of exposure, For weight, The average exposure time, For skin surface area, The coefficient of skin adhesion. As a skin absorption factor, Respiratory rate, For particulate matter emission factors, For reference dosage, The above parameters are carcinogenic slope factors; they are exposure parameters specific to agricultural land scenarios.

4. The method for priority control assessment of heavy metals in agricultural land soil according to claim 1, characterized in that, In step S4, the ecological risk assessment module adopts the Hakanson potential ecological risk index model; for each group of concentration samples Calculate the potential ecological risk coefficient of each element. In the formula For the first Toxicity response factors of heavy metals, The background value obtained in step S1; Calculate the comprehensive potential ecological risk index The values ​​of the toxicity response factors are: Hg=40, Cd=30, As=10, Pb=5, Cu=5, Ni=5, Cr=2, Zn=1.

5. The method for priority control and assessment of heavy metals in agricultural land soil according to claim 1, characterized in that, In step S5, the risk level classification standard is as follows: For the hazard index HI: HI<1 is acceptable risk, and HI≥1 is unacceptable risk; For total carcinogenic risk TCR: TCR < 1 × 10 -6 For negligible risk, 1×10 -6 ≤TCR<1×10 -4 For acceptable risk, TCR ≥ 1 × 10 -4 This is an unacceptable risk; For the comprehensive potential ecological risk index RI: RI<150 is low risk, 150≤RI<300 is medium risk, 300≤RI<600 is high risk, and RI≥600 is extremely high risk.

6. The method for priority control and assessment of heavy metals in agricultural land soil according to claim 1, characterized in that, In step S6, the health risk contribution condition is a standardized regression coefficient. Or Spearman rank correlation coefficient The ecological risk contribution condition is a standardized regression coefficient. Or Spearman rank correlation coefficient The comprehensive contribution index ,in The preset health risk weighting coefficient, .

7. The method for priority control assessment of heavy metals in agricultural land soil according to claim 6, characterized in that, The health risk weighting coefficient .

8. The method for priority control assessment of heavy metals in agricultural land soil according to claim 1, characterized in that, In step S7, the source apportionment model is the APCS-MLR model. Specifically, the APCS-MLR model is used to perform source apportionment on the original concentration data to identify x pollution sources and the contribution rate matrix of each source to each heavy metal. For the priority control factors determined in step S6, the pollution source with the highest contribution rate is found and determined as the primary control source.

9. The method for priority control assessment of heavy metals in agricultural land soil according to claim 1, characterized in that, In step S1, the sampling depth of the deep soil is 2.0m, and the sampling depth is determined according to the actual stratum distribution of the plot.

10. The method for priority control assessment of heavy metals in agricultural land soil according to claim 1, characterized in that, In step S8, the risk management recommendations include: priority management factors for which intensive monitoring is recommended, types of remediation measures to be taken, and categories of pollution sources to be prioritized for management.

Images