Method for analyzing and evaluating heavy metal pollution of sediment

Abstract

This invention discloses a method for analyzing and evaluating heavy metal pollution in sediments, belonging to the field of environmental statistical modeling technology. The method includes: acquiring sediment physicochemical property indicators; calculating the effective adsorption capacity of each component using a preset saturation response relationship; inputting the physicochemical indicators into a Bayesian network constructed with hierarchical mechanism constraints to infer the latent variable of heavy metal adsorption affinity; constructing a competitive mapping function containing a vacancy constant term based on a competitive Langmuir mechanism to calculate the capacity occupancy rate under multi-metal competition conditions; and finally calculating the predicted heavy metal content from the effective adsorption capacity and occupancy rate. This invention introduces a physical capacity upper limit and a multi-metal competition occupancy mechanism into the Bayesian network, solving the problems of prediction distortion and lack of competition mechanism in traditional models under high loads, and achieving accurate evaluation consistent with physicochemical laws.

Description

Technical Field

[0001] This invention belongs to the field of environmental statistical modeling technology, and in particular to the method for analyzing and evaluating heavy metal pollution in sediments. Background Technology

[0002] Sediments are important sinks and sources of heavy metals in aquatic ecosystems. Active sites on their surfaces, such as clay minerals, organic matter, and oxides, regulate the migration and transformation of heavy metals through mechanisms like adsorption and complexation. Accurately assessing the adsorption potential and pollution risk of heavy metals in sediments is of significant scientific importance and practical value for revealing the fate mechanisms of pollutants at the water-soil interface and guiding watershed environmental governance.

[0003] Currently, the analysis and evaluation of heavy metals in sediments mainly rely on environmental statistical indicators such as geoaccumulation index and potential ecological risk index, or use data-driven methods such as multiple regression, principal component analysis, and traditional Bayesian networks to establish quantitative relationships between environmental factors and metal content. Existing studies mostly focus on using statistical methods to explore correlations between data, attempting to simulate the enrichment characteristics of heavy metals in sediments through linear fitting or unconstrained probabilistic inference.

[0004] However, existing evaluation models generally lack physical mechanism constraints, mainly suffering from two core problems: overestimation of capacity due to neglecting the limited number of adsorption sites and the lack of multi-metal competitive occupancy mechanisms. Specifically, traditional linear or Gaussian models do not consider the saturation effect of active sites on sediment surfaces, often predicting adsorption amounts exceeding physical limits under high pollution loads. Simultaneously, existing methods typically evaluate single metals in isolation or rely solely on statistical correlation modeling, failing to characterize the competitive crowding-out behavior of different heavy metals on the same type of adsorption sites (such as cation exchange sites or coordination sites), leading to predictions deviating from the true chemical equilibrium state in multi-metal complex pollution scenarios. Furthermore, purely data-driven structure learning is prone to generating spurious correlation structures that violate geochemical causal timelines. Summary of the Invention

[0005] The purpose of this invention is to provide a method for analyzing and evaluating heavy metal pollution in sediments, in order to solve the aforementioned problems in the prior art.

[0006] Technical solution: A method for analyzing and evaluating heavy metal pollution in sediments, comprising:

[0007] Obtain the physicochemical properties of sediment samples from the target water body;

[0008] Based on physicochemical properties, the effective adsorption capacity of various adsorption components for heavy metals is calculated using a preset saturation response relationship.

[0009] The physicochemical properties are input into a pre-constructed Bayesian network model to infer the latent adsorption affinity of each heavy metal element, where the latent adsorption affinity follows a conditional linear Gaussian distribution.

[0010] Based on the latent variable of adsorption affinity, the capacity occupancy rate of each heavy metal element for adsorption sites is calculated using a preset competitive mapping function. The competitive mapping function introduces a vacancy competition mechanism through the denominator term to constrain the capacity occupancy rate within a preset range.

[0011] The predicted heavy metal content is calculated based on the effective adsorption capacity and capacity occupancy rate, and the heavy metal pollution assessment results are output accordingly.

[0012] Beneficial effects: This invention solves the problems of prediction distortion and lack of competition mechanism in traditional models under high load, and achieves accurate evaluation in accordance with the laws of physicochemical processes. Attached Figure Description

[0013] Figure 1 This is a complete technical flowchart of the embodiments of this application, from data acquisition and physical constraint calculation to pollution result output.

[0014] Figure 2 This is a causal network diagram showing the relationship between sediment properties and heavy metal content in the embodiments of this application.

[0015] Figure 3 This is a diagram of the Bayesian network structure in an embodiment of this application.

[0016] Figure 4 This is a scatter plot of the model predictions and measured values ​​in the embodiments of this application.

[0017] Figure 5 This is a diagram showing the Bootstrap analysis results of Zn in the embodiments of this application. Detailed Implementation

[0018] Example 1 describes the complete technical process from data acquisition and physical constraint calculation to pollution result output, such as... Figure 1 As shown in the figure. This embodiment focuses on demonstrating how to combine physicochemical mechanisms with Bayesian networks to solve the problem that traditional statistical models may violate physical common sense (such as predicted values ​​exceeding the upper limit of adsorption capacity) when predicting heavy metal pollution.

[0019] Step 101: Obtain the physicochemical properties of sediment samples from the target water body.

[0020] In this embodiment, the target water body is selected as a typical estuarine area, such as the waters adjacent to the North Channel of the Yangtze River Estuary. Sediment samples were collected from the top 0 to 5 cm of sediment using a grab sampler. The obtained physicochemical properties include, but are not limited to, the particle size distribution (e.g., percentage content of clay, silt, and sand), total organic carbon (TOC), total nitrogen (TN), total phosphorus (TP), specific surface area (SSA), and the content of major minerals (e.g., kaolinite, illite, quartz, feldspar, calcite, and hematite). In addition, the measured contents of various heavy metal elements (e.g., arsenic (As), chromium (Cr), copper (Cu), zinc (Zn), lead (Pb), and cadmium (Cd)) in the samples were simultaneously determined for model training or validation. Specifically, particle size distribution can be determined using a Bettersize 2600 laser particle size analyzer, heavy metal content can be determined using an ICP-optima 8000 inductively coupled plasma atomic emission spectrometer, and mineral composition can be determined using X-ray diffraction (XRD).

[0021] Based on this, the original measurement data of physicochemical properties and heavy metal content are subjected to a logarithmic transformation of log(1+x) to satisfy the normal distribution assumption of the Bayesian network model, i.e., x'=log(1+x), where x' is the transformed data value, log is the natural logarithm function (or the logarithm to base 10, which needs to be uniformly defined), x is the original measured physicochemical property or heavy metal content value, and the constant 1 is used to avoid the logarithm being meaningless when x=0. This preprocessing method can effectively reduce the skewness of the data distribution and reduce the impact of dimensional differences. Furthermore, for some highly correlated physicochemical properties, such as total organic carbon (TOC) and total phosphorus (TP), principal component analysis (PCA) can be performed. PCA comprehensive factor formula: F PC1 =w1*Z TOC +w2*Z TP Among them, F PC1 The first principal component extracted (i.e., the comprehensive factor TOC) TP_PC1 The values ​​of ), w1 and w2 are the loading coefficients (weights) calculated by principal component analysis, and Z. TOC Z represents the standardized total organic carbon content. TP This represents the standardized total phosphorus content.

[0022] Specifically, the first principal components of TOC and TP were extracted using PCA and named the comprehensive factor TOC. TP_PC1 This involves performing principal component analysis on highly correlated physicochemical property indicators, synthesizing multiple highly correlated indicators into a comprehensive factor, and using the comprehensive factor as the input node of a Bayesian network model to eliminate the interference of multicollinearity on the estimation of model parameters.

[0023] Step 102: Based on the physicochemical properties, calculate the effective adsorption capacity of various adsorption components for heavy metals using the preset saturation response relationship.

[0024] In this embodiment, the effective adsorption capacity S tot This refers to the maximum theoretical limit that sediments can adsorb a certain heavy metal under specific environmental conditions. This capacity is not a fixed constant, but a physical quantity that dynamically changes with the composition of the sediment. The calculation process identifies key components with adsorption capacity for heavy metals, such as clay minerals, organic matter, and iron and manganese oxides. For each component, its individual contribution capacity is calculated using a pre-defined saturation response relationship. The contributions of each component are summed to obtain the total effective adsorption capacity.

[0025] Step 103: Input the physicochemical property indicators into the pre-constructed Bayesian network model to infer the latent variables of adsorption affinity of each heavy metal element, wherein the latent variables of adsorption affinity follow a conditional linear Gaussian distribution.

[0026] In this embodiment, the Bayesian network model is a pre-trained directed acyclic graph model, whose nodes include the aforementioned physicochemical properties (as parent nodes) and the latent adsorption affinity variables of heavy metals (as child nodes). The latent adsorption affinity variable θ is an intermediate variable at the statistical inference level, reflecting the strength of the sediment's adsorption tendency for a certain heavy metal under the current combination of environmental factors. This variable numerically follows a conditionally linear Gaussian distribution, meaning its mean is determined by a linear combination of parent nodes, and its variance is constant. By inputting measured physicochemical properties (such as clay content and pH value) into the model, the θ value corresponding to each heavy metal in each sample can be calculated using the Bayesian inference algorithm.

[0027] Step 104: Based on the latent variable of adsorption affinity, calculate the capacity occupancy rate of each heavy metal element for adsorption sites using a preset competitive mapping function. The competitive mapping function introduces a vacancy competition mechanism through the denominator term to constrain the capacity occupancy rate within a preset range.

[0028] In this embodiment, the capacity occupancy rate r is a dimensionless ratio between 0 and 1, representing the proportion of the effective adsorption capacity actually occupied by heavy metals. The competition mapping function constructs a nonlinear mapping bridge from the unconstrained latent variable θ to the constrained occupancy rate r. Specifically, the denominator of this function includes a constant term representing vacancies, allowing the model to accommodate unoccupied adsorption sites, consistent with real-world physicochemical adsorption scenarios.

[0029] Step 105: Calculate the predicted value of heavy metal content based on the effective adsorption capacity and capacity occupancy rate, and output the heavy metal pollution assessment result accordingly.

[0030] Specifically, the predicted heavy metal content is calculated based on the effective adsorption capacity and capacity occupancy rate. This involves calculating the product of the effective adsorption capacity and the capacity occupancy rate, ensuring that the predicted heavy metal content automatically meets the physical constraint of not exceeding the effective adsorption capacity. The formula is expressed as q=S tot ×r. Where q is the predicted value of heavy metal content, S tot S represents the total effective adsorption capacity of the sediment for this heavy metal, and r represents the occupancy rate of the effective adsorption capacity for this heavy metal (dimensionless, ranging from 0 to 1). Since r is strictly constrained to be between 0 and 1, the calculated q must be less than or equal to S. tot This mechanism eliminates the unreasonable phenomenon that traditional linear regression models might predict unlimited adsorption under extremely high concentration conditions. Ultimately, the system outputs this predicted value q, which can be compared with environmental quality standards to provide a heavy metal pollution assessment result, such as determining the pollution level or identifying major pollution sources.

[0031] Example 2 describes how to physically decompose adsorption parameters based on surface complexation theory, and how to use the Michaelis-Menten saturation model to calculate the nonlinear effective adsorption capacity, so as to more accurately characterize the adsorption behavior under high content of sediment components.

[0032] Step 201: Based on physicochemical properties, calculate the effective adsorption capacity of various adsorption components for heavy metals using a preset saturation response relationship. This includes: identifying multiple sediment components that provide the same type of adsorption sites from the physicochemical properties; calculating the component contribution capacity of each sediment component based on the saturation response relationship; and summing the component contribution capacities of multiple sediment components to obtain the effective adsorption capacity.

[0033] Total effective adsorption capacity summation formula: S tot,i =∑ g=1 G S g,i Among them, S tot,i S represents the total effective adsorption capacity of the i-th sample, G is the total number of sediment components involved in adsorption (e.g., clay, organic matter, iron and manganese oxides), and S is the total effective adsorption capacity of the i-th sample. g,i The capacity contributed to each component.

[0034] In this embodiment, sediment components are classified according to the chemical properties of the adsorption sites. For example, for cation exchange type adsorption sites, the components providing these sites typically include clay minerals and organic matter; for coordination adsorption type sites, the components providing these sites typically include iron oxides and manganese oxides. After the system identifies the specific components, it calculates their respective capacity contributions and then sums them up. This summation mechanism reflects the additive effect of the different components' adsorption capacities for heavy metals.

[0035] Step 202: Calculate the component contribution capacity of each sediment component, depending on the unit capacity coefficient of the sediment component; the unit capacity coefficient is decomposed into the product of the surface adsorption site density and the metal-specific binding factor; wherein, the surface adsorption site density is calculated based on the specific surface area of ​​the sediment component and the number of active sites per unit area; the metal-specific binding factor is determined based on the ratio of the binding constant of the heavy metal to the sediment component to the binding constant of the reference metal.

[0036] In a preferred embodiment, the unit capacity coefficient β obtained from physical decomposition g =ρ g ×φ m,g It can be used in conjunction with the Michaelis-Menten saturation function. In this case, the maximum contributing capacity S... max,g Determined by site density and binding factor: S max,g =ρ g ×φ m,g ×B g,sat ;

[0037] Among them, B g,sat This represents the upper limit of the effective content of the g-th sediment component under aggregate saturation. When the component content is much lower than the half-saturation constant, the saturation function degenerates into a linear model, and the linear proportionality coefficient is approximately equal to S. max,g / K half,g =ρ g ×φ m,g ×B g,sat / K half,g Therefore, the physical decomposition scheme can be regarded as a degenerate form of the saturation response relationship in the low content range, and the two remain consistent in the theoretical framework.

[0038] Specifically, the unit capacity factor β g β is a proportionality coefficient connecting the content of sediment components to the adsorption capacity, expressed in mol / kg or mg / kg, depending on the conversion coefficient. In this embodiment, β g Instead of being an empirically fitted value, it is decomposed into two parameters with clear physical meaning: the density of surface adsorption sites ρ. g and metal-specific binding factor φ m,g That is, β. g =ρ g *φ m,g .

[0039] Among them, the surface adsorption site density ρ g This represents the number of moles of active sites on the surface of a unit mass of a g-th sediment component (such as clay minerals), expressed in mol / kg. The formula for calculation is:

[0040] ρ g =(SSAg ×σ g ) / N A ;

[0041] In the formula, SSA g The specific surface area of ​​the g-th component is expressed in m². 2 / kg; σ g The number of active sites per unit area, expressed in sites / m². 2 N A is Avogadro's constant, approximately 6.022 × 10⁻⁶. 23 sites / mol.

[0042] For example, referring to a typical sediment composition parameter table, the specific surface area (SSA) of clay minerals can be taken as 50 × 10⁻⁶. 3 Up to 800×10 3 m 2 Between / kg, its surface site density σ is approximately 2.3×10 18 sites / m 2 The SSA value for iron oxides can be in the range of 100 × 10⁻⁶. 3 Up to 300×10 3 m 2 Between / kg, σ is approximately 2.0 × 10 18 sites / m 2 By substituting these physical parameters, ρ can be calculated accurately. g .

[0043] Metal-specific binding factor φ m,g This reflects the difference in affinity of different metals for the same site. The calculation formula is:

[0044] φ m,g =K m,g / K ref,g ;

[0045] In the formula, φ m,g K is the dimensionless binding factor (K) of metal m to the g-th class component. m,g K is the apparent binding constant between metal m and the g-th site. ref,g Reference metal (e.g., Cu) 2+ The binding constant of ). This type of constant can be obtained from standard coordination chemistry databases (such as the NIST database).

[0046] Step 203: The saturation response relationship is the Michaelis-Menten saturation function; calculate the component contribution capacity for each sediment component, specifically using the pre-configured maximum contribution capacity and half-saturation constant of that sediment component, according to the following formula:

[0047] S g,i =S max,g *B g,i / K half,g +B g,i ;

[0048] Among them, S g,i B is the component contribution capacity of the g-th sediment component in the i-th sample. g,i S represents the measured content of this sediment component. max,g K represents the maximum contributing capacity of this sediment component. half,g This is the half-saturation constant of the sediment component, i.e., when the capacity reaches S... max,g The component content corresponding to half.

[0049] In this specification, the measured content of sediment components B g,i It is uniformly expressed as mass fraction, with units of g / kg or percentage content %. When entering as percentage content, the conversion relationship is: B g,i (g / kg)=B g,i (%)×10. In the Michaelis-Menten saturation function, B g,i With half-saturation constant K half,g Use the same units.

[0050] In real-world environments, when the content of sediment components (such as organic matter) is too high, the growth rate of effective adsorption sites gradually slows down and tends towards saturation due to aggregation effects or site shading. To characterize this nonlinear phenomenon, this embodiment employs the Michaelis-Menten saturation function. Where S... max,g K represents the maximum capacity limit that can be provided when the content of this component is infinitely large; half,g It is the half-saturation constant, that is, when the content of component B is... g,i equals K half,g At that time, the capacity contributed by this component reaches exactly half of its maximum value.

[0051] As an example, suppose that for a certain sample, the clay mineral content is B. clay It is 20% (i.e., 200g / kg). If the maximum contribution capacity S of this type of clay is preset... max,clay The half-saturation constant K is 100 mg / kg. half,clay The value is 100 g / kg. Substituting this into the formula, the volume contributed by clay in this sample is calculated as follows:

[0052] S clay =(100×200) / (100+200)=66.67mg / kg;

[0053] In contrast, if a simple linear model is used (assuming β=S)max / K half If the concentration is 1, the predicted value is 1×200=200mg / kg. The linear model will significantly overestimate the capacity at high concentrations, while the saturation model in this embodiment gives an estimate that is more in line with physical laws.

[0054] In some alternative implementations, when the content of sediment components is relatively low, i.e., B g,i <<K half,g When sufficient saturation parameter data is lacking, the above saturation model can be simplified to a linear adsorption capacity model. In this case, the formula approximates:

[0055] S g,i ≈(S max,g / K half,g )*B g,i ;

[0056] This provides a computationally simple alternative suitable for sediment scenarios with low pollution or low organic matter content.

[0057] Example 3 describes how to construct a competitive occupancy mechanism for a limited number of adsorption sites for multi-metal pairs at the mathematical model level by introducing latent variables of adsorption affinity, competitive grouping strategies, and improved competitive Langmuir mapping functions, thereby solving the technical problem that traditional statistical models cannot characterize the physical competition process.

[0058] Step 301: The latent adsorption affinity variable follows a conditional linear Gaussian distribution; the latent adsorption affinity variable of each heavy metal element is inferred, specifically based on the parent node variable pointing to the heavy metal element in the Bayesian network structure, and calculated according to the following linear Gaussian model:

[0059] θ m,i =a m,0 +∑ k=1 K (a m , k *X k,i )+ε m,i ;

[0060] Where, θ m,i Let a be the latent variable of the adsorption affinity of metal m in the i-th sample. m,0 X is the intercept term. k,i Let a be the observed value of the k-th parent node (physicochemical property index). m,k For the corresponding regression coefficient, ε m,i The noise term follows a Gaussian distribution.

[0061] In this embodiment, the latent adsorption affinity θ is defined as a dimensionless intermediate variable used to characterize the driving force for heavy metals to occupy adsorption sites. Specifically, θm =ln(K m *C m ); where θ m Let K be the latent variable for adsorption affinity, ln be the natural logarithm, and K be the latent variable for adsorption affinity. m C is the adsorption equilibrium constant. m The concentration of free metal ions in the liquid phase is θ. Since Bayesian networks employ the conditional Gaussian distribution assumption, the value of θ is determined by a linear combination of its parent nodes (i.e., environmental factors influencing the enrichment of the heavy metal). For example, for the heavy metal zinc (Zn), its parent nodes might include clay content, total nitrogen (TN) content, and zeta potential. In this case, the latent variable θ representing the adsorption affinity of zinc is... Zn It equals the constant term plus the sum of the products of the observations of each parent node and their corresponding regression coefficients. This linear structure ensures the efficiency of Bayesian network inference.

[0062] Correspondingly, the latent variable θ of adsorption affinity can be physically interpreted as the logarithm of the product of the adsorption equilibrium constant and the liquid phase concentration. This interpretation provides a theoretical basis for the functional form of the competition mapping function. In actual calculations, θ serves as an intermediate latent variable for statistical inference. Its value is approximately inferred by a Bayesian network based on observable physicochemical properties, rather than by directly measuring the liquid phase concentration and equilibrium constant. Physicochemical properties (such as clay content, organic carbon content, and pH value) indirectly determine the levels of the adsorption equilibrium constant and the liquid phase metal concentration by regulating the number of adsorption sites and adsorption environmental conditions, and therefore can serve as effective predictors of θ.

[0063] Step 302: Calculate the capacity occupancy rate of each heavy metal element for adsorption sites using a preset competition mapping function. This includes: dividing multiple heavy metal elements into several competition groups based on the chemical properties of adsorption sites on the sediment surface. Heavy metal elements belonging to the same competition group share the same type of adsorption sites and have a competitive relationship, while heavy metal elements belonging to different competition groups correspond to different types of adsorption sites. For each competition group, the capacity occupancy rate of each heavy metal element in the group is jointly calculated based on the latent adsorption affinity of all heavy metal elements in the group using the competition mapping function.

[0064] Step 303: Divide the various heavy metal elements into several competing groups, specifically including: dividing copper, zinc, lead, and cadmium into a cation competing group, and setting the cation competing group to share cation exchange sites on the surface of clay minerals and organic matter; dividing arsenic and chromium into an oxygen-containing anion competing group, and setting the oxygen-containing anion competing group to share coordination adsorption sites on the surface of iron-manganese oxides.

[0065] In this embodiment, to conform to the environmental chemistry mechanism, heavy metals with different properties must be grouped. For example, copper (Cu), zinc (Zn), lead (Pb), and cadmium (Cd) usually exist as divalent cations. They mainly compete for negatively charged sites on the surface of clay minerals or carboxyl and phenolic hydroxyl sites on the surface of organic matter. Therefore, they are classified into the cation competition group, denoted as G1. Arsenic (As) and chromium (Cr) usually exist as oxyanions (such as arsenate and chromate). They mainly undergo coordination exchange reactions with hydroxyl groups on the surface of iron and manganese oxides. Therefore, they are classified into the oxyanion competition group, denoted as G2. In the calculation, the metals in group G1 compete with each other, and the metals in group G2 compete with each other, but there is no interference between groups G1 and G2. Each group occupies its corresponding effective adsorption capacity.

[0066] Step 304: The competitive mapping function is constructed based on the competitive Langmuir adsorption mechanism; the capacity occupancy rate of each heavy metal element in the group is jointly calculated, specifically according to the following formula: r m,i =exp(θ m,i ) / 1+∑ m'∈Gk exp(θ m',i );

[0067] Where, r m,i For the i-th sample, belonging to the k-th competitive group G k The capacity utilization rate of heavy metal m, θ m,i Let ∑ be the latent variable of the adsorption affinity of the heavy metal m. m'∈Gk exp(θ m',i ) represents the exponential sum of the adsorption affinity potential variables of all heavy metal elements in this competition group, and the constant term 1 in the denominator is used to characterize the vacancies present in this type of adsorption site.

[0068] Unlike the traditional Softmax normalization function, this formula explicitly includes a constant term 1 in the denominator. This constant term has a clear physical meaning, representing the vacant sites in the adsorption sites that are not occupied by any metal. This makes the sum of the capacity occupancy rates of all metals in the group ∑r m,i It is always strictly less than 1, rather than being forced to equal 1. This aligns with actual physical processes, meaning that some sites on the sediment surface will always be vacant and cannot always be completely occupied. When the adsorption affinity θ of all metals is very low (i.e., a large negative value), exp(θ) approaches 0, and the occupancy rate r approaches 0, indicating that most sites are vacant. When the θ of a certain metal is large, its occupancy rate r will increase, while the summation term in the denominator suppresses the occupancy rate of other metals in the same group, reflecting a competitive effect of one increasing while the other decreases.

[0069] Step 305: The competition mapping function also introduces a competition intensity coefficient to characterize the difference in affinity of different heavy metal elements for the same type of adsorption site; the capacity occupancy rate of each heavy metal element in the group is jointly calculated, specifically according to the following formula: r m,i =(α m *exp(θ m,i )) / 1+∑ m'∈Gk α m' *exp(θ m',i );

[0070] Where, α m α is the competition intensity coefficient for heavy metal m, used to correct for differences in intrinsic affinity between different metal pairs at specific sites. m' is the competition intensity coefficient of heavy metal m' in the same group.

[0071] To more accurately characterize the differences in chemical properties among different metals, this embodiment introduces a competition intensity coefficient α into the basic Langmuir formula. This coefficient is used to correct for adsorption affinity and reflects the inherent strength of the metal's adsorption capacity. In some preferred embodiments, for the cation competition group, based on the Irving-Williams sequence or relevant adsorption experimental data, lead is typically set to have the highest competition intensity, followed by copper, then zinc, and cadmium the lowest. Specifically, α can be set... Pb >α Cu >α Zn >α Cd For example, let α be... Pb =5.0, α Cu =2.0, α Zn =1.0, α Cd =0.5.

[0072] As a specific numerical example, suppose that in a certain sample, the cation group only considers two metals, copper and zinc. After inferring the latent variable θ of copper's adsorption affinity using a Bayesian network... Cu =0.5 (i.e., exp(0.5)≈1.648), the latent variable θ of zinc adsorption affinity. Zn =0.2 (i.e., exp(0.2)≈1.221).

[0073] If the competition intensity coefficient is not introduced (i.e., α=1), then:

[0074] Denominator term = 1 + 1.648 + 1.221 = 3.869;

[0075] Copper utilization rate r Cu =1.648 / 3.869≈0.426;

[0076] Zinc utilization rate r Zn=1.221 / 3.869≈0.316;

[0077] The combined occupancy rate of the two is 0.742, with 0.258 remaining empty slots.

[0078] If we introduce a competition intensity coefficient, let α Cu =2.0, α Zn =1.0.

[0079] The weighted term for copper is 2.0 × 1.648 = 3.296, and the weighted term for zinc is 1.0 × 1.221 = 1.221.

[0080] Denominator term = 1 + 3.296 + 1.221 = 5.517;

[0081] Copper utilization rate r Cu =3.296 / 5.517≈0.597;

[0082] Zinc utilization rate r Zn =1.221 / 5.517≈0.221;

[0083] It can be seen that copper, with its higher competition intensity coefficient α, occupies a larger share in the competition (from 0.426 to 0.597), while zinc is more strongly suppressed (from 0.316 to 0.221). This calculation result is more consistent with the chemical fact that strongly adsorbed metals preferentially occupy sites.

[0084] Example 4 describes a simplified alternative for evaluating a single heavy metal. When the practical application focuses only on the pollution status of a specific heavy metal, or when the content of other competing metals in the sediment is so low as to be negligible, the complex competition model can degenerate into an independent Logit model.

[0085] Step 401: When analyzing and evaluating a single heavy metal element, the preset competition mapping function degenerates into the inverse Logit transform; the capacity utilization rate is calculated according to the following formula: r m,i =1 / 1+exp(-θ m,i ); constrain the capacity utilization rate between 0 and 1.

[0086] When analyzing and evaluating a single heavy metal element, the competition group contains only one metal. If a competition intensity coefficient is not introduced or is set to 1, the competition mapping function degenerates into the standard Sigmoid function, i.e., r m,i =1 / (1+exp(-θ m,i If a competition intensity coefficient α is introduced... m The degenerate form is r m,i =α m ×exp(θm,i ) / (1+α m ×exp(θ m,i In the case of a single metal, α m The effect is equivalent to the effect on the intercept term a m,0 Apply offset ln(α) m This can be absorbed by the parameter learning of the Bayesian network, therefore α in single-metal evaluation m It can be set to 1 by default.

[0087] In this embodiment, the inverse transform of the Logit transform is the Sigmoid function. This formula is the mathematical degenerate form of the competing Langmuir formula in the case of a single component. The specific derivation is as follows: When there is only one metal m in the competing group, the summation term in the denominator is reduced to only one term: exp(θ). m,i The formula then becomes:

[0088] r m,i = exp(θ m,i ) / (1+exp(θ m,i ));

[0089] Divide both the numerator and denominator by exp(θ) m,i That is, we get 1 / 1+exp(-θ) m,i ).

[0090] The physical implication of this alternative is that it assumes heavy metals compete only with vacancies, not with other metals. In this case, the larger the latent variable of adsorption affinity θ, the closer the occupancy rate r approaches 1; conversely, it approaches 0. This approach is computationally simpler and still retains the physical constraint that the result is between 0 and 1, making it suitable for single-factor pollution assessment or preliminary screening scenarios.

[0091] Example 5 describes the offline construction process of the Bayesian network model. To obtain a Bayesian network with a reasonable structure and accurate parameters, this example proposes a hierarchical structure learning method that integrates prior geochemical knowledge, effectively avoiding the spurious correlations or reverse causality problems that may arise from purely data-driven algorithms.

[0092] Step 501: The pre-constructed Bayesian network model is trained through the following steps: historical observation data containing sedimentary physicochemical properties and heavy metal content are obtained as a training dataset; based on the training dataset, a heuristic search algorithm is used to learn the network structure and determine the directed acyclic graph structure of the Bayesian network; based on the determined directed acyclic graph structure, the maximum likelihood estimation method is used to learn the conditional probability parameters in the network to obtain the linear Gaussian model parameters used to infer the latent variables of adsorption affinity.

[0093] Step 502 involves learning the network structure using a heuristic search algorithm, including a hierarchical constraint process based on geochemical mechanisms: dividing the node variables in the Bayesian network into three levels, namely the material source variable layer, the environmental state variable layer, and the heavy metal content variable layer; applying directed edge constraints during the network structure learning process, including: allowing points from the material source variable layer to the environmental state variable layer or the heavy metal content variable layer, allowing points from the environmental state variable layer to the heavy metal content variable layer, and prohibiting directed edges originating from the heavy metal content variable layer.

[0094] In this embodiment, the hierarchical division of variables is based on the source-sink mechanism and causal time sequence.

[0095] Specifically, the material source variable layer (Layer 1) includes the basic materials that make up the sedimentary skeleton, such as clay, sand, quartz, feldspar, iron (Fe), and aluminum (Al). These variables are usually determined by the geological background and are exogenous variables.

[0096] The environmental state variable layer (Layer 2) includes environmentally modulated physicochemical indicators such as pH, redox potential (Eh), zeta potential, and the organic matter-phosphate complex (TOC). TP_PC1 This variable plays a role in regulating the adsorption environment.

[0097] The heavy metal content variable layer (Layer 3) includes target pollutants such as As, Cr, Cu, Zn, Pb, and Cd, which are the result variables of the adsorption process.

[0098] Step 503, applying directed edge constraints during network structure learning, also includes setting a whitelist and a blacklist: the whitelist contains directed edges that must exist based on geochemical mechanisms, including directed edges pointing from pH value, redox potential, organic matter content or iron and manganese oxide content to specific heavy metal element variables; the blacklist contains directed edges that are prohibited from existing based on causal time sequence, including directed edges pointing from heavy metal content variable layer to material source variable layer or environmental state variable layer.

[0099] Specifically, this embodiment employs the constrained hill-climbing algorithm for structure learning. The algorithm's input includes a training dataset D, a whitelist W, and a blacklist B. Initially, the network structure G0 contains all edges from the whitelist W. During iteration, the algorithm attempts to add, delete, or reverse edges and calculates the change in the Bayesian Information Criterion (BIC) score. Crucially, any operation belonging to the blacklist B (e.g., attempting to add an edge from Cd to pH) is forcibly prohibited, as is any operation involving deleting edges from the whitelist W. For example, the whitelist can be configured to require the presence of Fe->As edges because iron oxides are the primary carriers of arsenic; the blacklist can be configured to prohibit Cd->Clay edges because heavy metal pollution does not alter the particle size composition of sediments.

[0100] The Bayesian network structure obtained through the above-mentioned constraint training is as follows: Figure 3 As shown in the figure, the causal path from material sources (such as muscovite and hematite) to heavy metals (such as As and Zn) is clearly illustrated, with no reverse arrows, verifying the effectiveness of the hierarchical constraints. After determining the structure, maximum likelihood estimation is used to determine the conditional probability distribution parameters (i.e., linear regression coefficients a) of the MLE for each node. m,k The model parameters are fitted using the variance and variance to obtain the complete model parameters.

[0101] Example 6 describes how Bootstrap robustness analysis was used to verify the reliability of the model and identify key driving factors affecting heavy metal enrichment. This process not only verifies the model's predictive ability but is also an important tool for analyzing the causes of environmental pollution.

[0102] Step 601, the method further includes performing a robustness analysis of the regression coefficients of the Bayesian network model: using the bootstrap method to perform multiple resamplings with replacement on the observed data of physicochemical properties and heavy metal content, generating multiple resampling datasets; while keeping the Bayesian network structure unchanged, re-estimating the regression coefficients of the latent variable of adsorption affinity using each resampling dataset; statistically analyzing the distribution characteristics and confidence intervals of the regression coefficients in the multiple resampling processes, and assessing the significance of the influence of physicochemical properties on heavy metal enrichment based on this.

[0103] In this embodiment, the Bootstrap method is used to evaluate the stability of the model parameters. Specifically, the following steps are taken: Random sampling with replacement is performed on the original observation dataset to construct 1000 resampled datasets with the same sample size as the original dataset. For these 1000 datasets, while maintaining the previously determined Bayesian network structure, the parameters are relearned to obtain 1000 sets of regression coefficients.

[0104] By statistically analyzing the distribution of these 1000 sets of regression coefficients, the mean and 95% confidence interval of the influence coefficient of each environmental factor on heavy metals can be calculated. For example... Figure 5 As shown, the Bootstrap analysis results for Zn are presented. For the heavy metal zinc (Zn), the confidence intervals of the regression coefficients for its parent nodes, clay and total nitrogen (TN), are above the zero mark. This indicates that these two factors have a significant positive promoting effect on zinc enrichment, and this conclusion is statistically robust. Conversely, the confidence intervals of the regression coefficients for some factors (such as carboxyl COO) cross zero, indicating that their influence direction is uncertain and they are not key control factors. Furthermore, this embodiment also uses leave-one-out cross-validation (LOOCV) to evaluate the predictive performance of the model. Figure 4 As shown, the scatter plot of the model predictions and measured values ​​closely follows the y=x line, and the coefficient of determination R0 2 The values ​​are generally between 0.6 and 0.9, which proves that the method of the present invention still has excellent prediction accuracy and generalization ability under limited sample conditions.

[0105] Describe the specific implementation process of one of the above embodiments using a specific scenario:

[0106] Step 1: Data Acquisition: Select representative sediment samples, determine the content of each heavy metal in the sediments, and determine the set of material indicators used to characterize the supply of exchangeable sites, including mineralogical composition, surface functional group characteristics, and specific surface area.

[0107] Step 2: Preprocessing the obtained data: Remove irrelevant fields such as sample numbers, and use logarithmic transformation (e.g., taking log(1+x) for non-negative real numbers to reduce the influence of dimensional differences and skewed distributions). For highly correlated indicators, dimensionality reduction methods can be used to synthesize them into comprehensive factors (e.g., in this invention, principal component analysis is used to combine TOC and TP into a comprehensive factor TOC). TP_PC1 (This represents the combined effects of organic matter and phosphate). The preprocessed factor and heavy metal content data constitute the input data table for the Bayesian network model. This invention is based on surface complexation theory, using the unit capacity coefficient β... g It is decomposed into the product of site density and binding factor, giving it a clear physicochemical meaning.

[0108] S2.1, Site density calculation:

[0109] For the g-th type of sediment component, its surface adsorption site density ρ g (Unit: mol / kg) Calculated by the following formula:

[0110] ρ g =(SSA g ·σ g ) / N A ;

[0111] Among them, SSA g σ is the specific surface area (m² / kg) of the g-th component. g N represents the number of active sites per unit area (sites / m²). A is Avogadro's constant.

[0112] Table 1: Parameter values ​​of typical sediment components

[0113]

[0114] In this specification, the specific surface area SSA g The unit is uniformly m² / kg, and the number of active sites σ g All units are uniformly expressed as sites / m². For ease of understanding, the parameter table also separately indicates the commonly used units m² / g and sites / nm², with the conversion relationships as follows: 1 m² / g = 1000 m² / kg; 1 site / nm² = 10 18 sites / m². When substituting into the formula for calculation, the values ​​should be uniformly converted to m² / kg and sites / m².

[0115] S2.2, Metal-specific binding factor:

[0116] Considering the differences in affinity of different metals for the same type of site, a metal-specific binding factor φ is introduced. m,g , representing the effective occupancy ability of metal m at the g-th type site. This factor is determined by coordination chemistry theory:

[0117] φ m,g =K m,g / K ref,g ;

[0118] Among them, K m,g K is the apparent binding constant between metal m and the g-th site. ref,g Reference metal (e.g., Cu) 2+ The binding constant of ). This ratio can be obtained from the NIST coordination chemistry database or thermodynamic models such as WHAM / Model VII.

[0119] S2.3 Calculation of effective adsorption capacity:

[0120] Considering both site density and binding factor, the effective adsorption capacity of the g-th component for metal m is:

[0121] S m,g,i =ρ g ·φ m,g ·B g,i ;

[0122] Then the total effective capacity of the competing group k to which metal m belongs is:

[0123] S tot,i (k,m) =∑ g∈Ωk S m,g,i =∑ g∈Ωk ρ g ·φ m,g ·B g,i ;

[0124] At this time, β g Decomposed into ρ g ·φ m,g , where ρ g Derived from physical parameters, φ m,g All data are derived from thermodynamic databases and have clear physical origins.

[0125] Step 3: Gaussian Latent Variables and Network Structure Learning: To maintain the applicability of conditional Gaussian modeling in Bayesian networks, a logit transformation is performed on the capacity occupancy variable to construct Gaussian latent variable nodes:

[0126] y m,i =logr m,i / 1-r m,i ;

[0127] And assume y m,i Given its parent node, it satisfies a linear Gaussian relationship. Bayesian network structure learning is performed at the Gaussian latent variable level, employing a heuristic search algorithm combined with prior constraints to determine the dependencies between environmental factors, adsorbent material indicators, and each metal occupancy latent variable.

[0128] Step 4: Parameter Learning and Deterministic Nonlinear Mapping: After completing the network structure learning, parameter learning is performed on the Gaussian latent variable nodes to obtain the following linear Gaussian model:

[0129] y m,i =a m,0 +∑ k a m,k X k,i +ε m,i ;

[0130] Subsequently, capacity utilization and heavy metal content were calculated using deterministic nonlinear mapping:

[0131] r m,i =1 / 1+e -y(m,i) ,q m,i =S tot,i ·r m,i ;

[0132] It automatically ensures that the predicted heavy metal results do not exceed the upper limit of adsorption capacity.

[0133] Step 5: Multi-metal competition constraint: In the preferred embodiment, a competition constraint mechanism is introduced for multiple heavy metals sharing the same type of adsorption sites. This can be achieved by setting remaining capacity nodes or group capacity nodes and using a vacancy-containing mapping function based on the competitive Langmuir mechanism to allocate the group content, ensuring that the sum of the contents of each metal in the group does not exceed the corresponding adsorption capacity, thus realizing the structured expression of multi-metal competitive site occupancy.

[0134] Specifically, the competition mapping function in this scheme is similar in form to the Softmax function, but it includes an additional constant term 1 in the denominator to represent the vacancy state. Therefore, the sum of the capacity occupancy rates of each metal in the group is strictly less than 1, rather than equal to 1. This difference corresponds to different physical assumptions: Softmax assumes that adsorption sites are always occupied, while this scheme allows some sites to be in an vacant state, which is more in line with the actual sediment adsorption scenario.

[0135] Step 6: Model Prediction Performance Evaluation: To verify the prediction accuracy and generalization ability of the established model, this invention employs the Leave-One-Out Cross-Validation (LOOCV) method to evaluate the model. Specifically, data from each sample is sequentially reserved as a test sample, while the remaining samples are used to train the Bayesian network model. In each round of validation, the Gaussian latent variables are predicted, and the corresponding predicted heavy metal content is calculated using deterministic nonlinear mapping. This predicted value is then compared with the measured value of the reserved sample, and the prediction error is recorded. After iterating through all samples, the prediction error indices for each heavy metal are calculated, including the root mean square error (RMSE) and the coefficient of determination (R²). RMSE reflects the degree to which the predicted value deviates from the measured value, and R² characterizes the model's ability to explain the variance of actual observations. The closer R² is to 1, the higher the model's prediction accuracy. LOOCV evaluation effectively verifies the robustness of the model under small sample conditions, avoiding bias in the model performance evaluation caused by individual outliers. Actual results show that the model of this invention has high prediction accuracy for most heavy metals, and the predicted values ​​are generally close to a 1:1 relationship with the measured values.

[0136] Step 7: Robustness Analysis of Regression Coefficients (Bootstrap): To evaluate the stability and reliability of the quantitative relationships related to heavy metals in the Bayesian network model under sample perturbation conditions, this invention introduces the bootstrap method to conduct robustness analysis of the model results. With the network structure unchanged, the original sample data is repeatedly sampled with replacement, and the model parameters are re-estimated on each sample dataset to obtain the repeated sample distribution of the estimated parameters of key influencing factors. Statistical analysis of the multiple sampling results allows for a comprehensive evaluation of the variation characteristics of model parameters under different sample combinations, determining the sensitivity of the model results to sample fluctuations and its overall stability. Through this robustness analysis process, this invention can effectively verify the reliability of the established model under limited sample conditions, providing support for subsequent result interpretation and pollution control factor identification.

[0137] To verify the effectiveness of the method of the present invention, the technical solution of the present invention is specifically described below with reference to heavy metal pollution data of sediments from a typical water body. This embodiment is only used to illustrate the implementation process of the method of the present invention and is not limited to a specific watershed or region. The study area is the waters adjacent to the North Channel of the Yangtze River Estuary. Several surface sediment sampling points were set up in this area, and the collected sediment samples were analyzed in the laboratory to obtain the heavy metal content and multiple physicochemical properties. The six heavy metal indicators selected are arsenic (As), chromium (Cr), copper (Cu), zinc (Zn), lead (Pb), and cadmium (Cd).

[0138] Table 2: Heavy metal content in sediments at sampling points

[0139]

[0140] The sedimentary physicochemical indicators include grain size composition (percentage content of clay, silt, and sand, and average grain size Md), total organic carbon (TOC), total nitrogen (TN), total phosphorus (TP), content of major clay minerals (such as kaolinite, illite, and chlorite) and non-clay minerals (such as quartz and calcite), surface functional groups, and specific surface area. In data preprocessing, all heavy metals and physicochemical indicators underwent a unified logarithmic transformation to reduce the influence of biases between data of different dimensions. TOC and TP were highly correlated and were synthesized into a comprehensive index, TOC, using principal component analysis. TP_PC1 Its first principal component explains most of the common variation between the two, representing the organic matter-phosphate complex factor.

[0141] The preprocessed data was loaded into the RStudio environment, and Bayesian network structure learning was performed using the bnlearn package. A whitelist was created based on the general mechanism of heavy metal adsorption in sediments, allowing directed edges from clay, silt, sand, average grain size (Md), specific surface area, nutrients, major minerals (such as quartz, biotite, chlorite, hematite, calcite), surface charge, and functional groups to nodes of various heavy metals (As, Cr, Cu, Zn, Pb, Cd). This is equivalent to assuming that the above factors can directly affect the enrichment degree of a certain heavy metal. In the blacklist, any edge emanating from a heavy metal node pointing to other nodes (including pointing to other heavy metals or any factors) is prohibited, reflecting the causal assumption that cause precedes effect, i.e., heavy metal content will not conversely affect the physicochemical properties of sediments or the content of other metals. Furthermore, direct connections between heavy metals are also prohibited, eliminating the possibility of heavy metals acting as causally related to each other. After setting the above constraints, the Hill-Climbing algorithm is used to search for the network structure. The scoring function adopts the BIC criterion based on the data continuity (for Gaussian distributed variables), and the final Bayesian network is obtained.

[0142] like Figure 2 As shown, the nodes are named after heavy metal elements and environmental factors, with arrows pointing from cause to effect, clearly presenting a causal network between sediment properties and heavy metal content. For example, the diagram shows that Quartz points to Cu, indicating that quartz content affects Cu enrichment; Fe points to Zn, indicating that iron oxides promote Zn enrichment; and Hematite points to As, indicating that iron oxides such as hematite are important adsorption carriers for As. The overall network structure is consistent with existing environmental geochemical understanding, validating the rationality of the structure learning results.

[0143] Model Validation and Result Analysis: Based on the learned network structure, parameter learning is performed using all sample data to obtain a linear Gaussian model for each heavy metal under the condition of its parent nodes. For example, if the parent nodes of the Cu node are Md, sand, and Ni (nickel), the corresponding linear regression model is as follows:

[0144] Cu=β0+β1·Md+β2·sand+β3·Ni+ε;

[0145] Here, ε represents the Gaussian noise term, and the coefficients β1, β2, β3, etc., can be estimated through parameter learning. To evaluate the model's predictive performance, LOOCV was used to predict the heavy metal content of each sample, and the calculation error was compared with the observed values. The results show that the predicted R0 for each heavy metal is... 2 A value between 0.6 and 0.9 (varying depending on the metal) indicates that the model can explain most of the observed variance; for example, the predicted R-values ​​for Zn, Pb, and Cu are [missing value]. 2 They are all around 0.8.

[0146] Figure 4 A comparison of the predicted and measured logarithmic concentrations of Zn is presented, along with the residual distribution. It can be seen that the predicted and measured values ​​of Zn closely approximate the sloping line y=x, indicating that the model's prediction of Zn is quite accurate. The residual scatter plots are randomly distributed around zero without a clear trend, indicating that the model has no systematic errors. Similar prediction-measurement plots are shown for other heavy metals, with most data points distributed around y=x, demonstrating that the model has good overall predictive consistency.

[0147] Based on the acceptance of the model's predictive performance, the robustness of the regression relationships among various heavy metals is further examined using the Bootstrap method. Taking Zn as an example, its parent nodes include clay, TN, and TOC. TP_PC1 COO functional group. Perform 1000 bootstrapping resampling operations on the original sample, each time fitting a linear model Zn ~ clay + TN + TOC. TP_PC1 +COO was used to collect the regression coefficient distributions of various factors. The results showed that the coefficients for clay content and total nitrogen had a positive mean in 1000 experiments, with a 95% confidence interval greater than 0, meaning they maintained a positive correlation with Zn under all sampling conditions. This indicates that the promoting effect of clay minerals and nitrogen source abundance on Zn enrichment is highly significant and robust. In contrast, although the regression coefficient of the COO functional group had a positive mean, its confidence interval crossed zero, implying that the direction of this factor's effect might be uncertain under different sampling conditions, and its influence might be insignificant. This suggests that the carboxyl functional group is not the main controlling factor for Zn adsorption.

[0148] Similarly, bootstrap analyses were performed for each heavy metal, ultimately identifying a set of key sedimentary factors that significantly influence each heavy metal. For example, significant positive influencing factors for Cu include Md and Zeta potential, while Pb is mainly affected by clay particles and TOC. TP_PC1 The positive effects were observed, with As and Cr showing a significant correlation with iron oxides (Fe-O). These findings are consistent with the mechanisms revealed by traditional experimental studies. Therefore, the method of this invention not only provides quantitative predictions but also statistically identifies the most influential pollution control factors, offering new insights and tools for a deeper understanding of the causal mechanisms of heavy metal pollution.

[0149] The method of this invention has been described in the above embodiments using sediments from the Yangtze River estuary as an example, but it is not limited to applications in a specific watershed or region. Since the method relies on only conventional sediment indices and heavy metal content measurements, the general workflow of this invention can be used for heavy metal pollution analysis of sediments in any water body (including rivers, lakes, reservoirs, and bays). For different scenarios, the types of prior factors in the network whitelist can be adjusted as needed (e.g., environmental factors such as pH and Eh can be added to mine wastewater sediments, and salinity can be added to marine sediments), flexibly expanding the model to adapt to new environmental conditions. The open-source tools provided by the RStudio platform ensure the convenience and reproducibility of the method implementation: users only need to change the input data and run the R script according to the same steps to automatically complete the entire process of model training, validation, and visualization, obtaining an analysis report for that dataset.

[0150] According to one aspect of this application, the remaining capacity node or group capacity node is set, and the group capacity is allocated using a softmax method, as follows:

[0151] In multimetallic sediment systems, different metal ions compete for a limited number of adsorption sites. This invention constructs a grouped competitive site-occupancy model based on the mathematical form of the competing Langmuir equation, embedding the chemisorption mechanism into the nonlinear mapping layer of a Bayesian network.

[0152] Adsorption site grouping:

[0153] Based on the chemical properties of the adsorption sites on the sediment surface, the M heavy metals were divided into K competing groups G1, G2, ..., G... K Metals in the same group share the same type of adsorption sites and compete with each other, while metals in different groups occupy different types of sites and are independent of each other.

[0154] Typical grouping schemes are as follows:

[0155] Cation group: G1=Cu 2+ ,Zn 2+ ,Pb 2+ Cd 2+ They compete for cation exchange sites on the surfaces of clay minerals and organic matter.

[0156] Oxygen-containing anion group: G2=AsO4 3- CrO4 2- Competition for coordination sites on the surface of iron-manganese oxides

[0157] Define the latent adsorption affinity variable: For each heavy metal m, define the latent adsorption affinity variable θ. m,i , representing the adsorption driving force of metal m under sample i conditions:

[0158] θm,i =a m,0 +∑ k a m,k X k,i ;

[0159] Among them, X k,i Let a be the k-th environmental factor of sample i. m,k The regression coefficient is to be estimated. This variable satisfies a linear Gaussian relationship under the condition of the parent node, maintaining the applicability of the standard Bayesian network learning algorithm.

[0160] Constructing a grouped competition placeholder model:

[0161] For metals m∈G belonging to the k-th group k Its capacity utilization rate is calculated using the following competitive occupancy formula: r m,i =exp(θ m,i ) / 1+∑ m'∈G_k exp(θ m',i );

[0162] This formula originates from a mathematical transformation of the competing Langmuir equation. Let the competing Langmuir equation be: q m =q max ·K m C m / 1+∑ m' K m' C m' ;

[0163] Where K m C is the adsorption equilibrium constant. m Let θ be the liquid phase concentration. m =ln(K m C m ), then: q m / q max =exp(θ m ) / 1+∑ m' exp(θ m' ); This is consistent with the above occupancy rate formula.

[0164] When a physical decomposition scheme with a unit capacity coefficient is adopted, different heavy metals have their own effective adsorption capacities due to differences in binding factors. In this case, the physical constraints reflected by the competition mapping function include two levels: first, the predicted content of any heavy metal does not exceed its own upper limit of effective adsorption capacity; second, when the adsorption affinity of any heavy metal in the competing group increases, the capacity occupancy rate of other heavy metals in the same group is competitively suppressed. These two levels together ensure that the prediction results conform to the chemical mechanism of multi-metal competitive adsorption.

[0165] This formula has the following physical properties:

[0166] Single-metal constraint: 0 <r m,i <1, the occupancy rate of any metal does not exceed 100%.

[0167] Intra-group total constraint: ∑ m∈Gk r m,i <1, the total metal occupancy rate within the group is strictly less than 1, and empty slots are allowed.

[0168] Competition effect: When the θ of a certain metal m,i When it increases, its r m,i Increases, while the r of other metals in the same group increases. m',i Decrease.

[0169] Group independence: The occupancy rate calculations for different groups are independent of each other, and each group satisfies its own total constraint.

[0170] Grouping capacity allocation: For each type of adsorption site k, define its effective adsorption capacity S. tot,i (k) :S tot,i (k) =∑ g∈Ωk β g B g,i ;

[0171] Among them, Ω k B is the set of sediment indicators associated with the k-th site. g,i β represents the measured value of the corresponding indicator. g This is the unit capacity factor.

[0172] Calculate heavy metal content: For metal m belonging to group k, its predicted content is: q m,i =S tot,i (k) ·r m,i ;

[0173] The result automatically satisfies the following physical constraints:

[0174] Single metal not exceeding group capacity: q m,i tot,i (k) ;

[0175] The total amount of metals in a group does not exceed the group capacity: ∑ m∈Gk q m,i tot,i (k) 。

[0176] Furthermore, a competition intensity coefficient is introduced.

[0177] ​​The basic competitive occupancy formula assumes that all metals have the same competitive ability for the same site (each metal has a weight of 1 in the denominator). In reality, different metals have different affinities for the same site. Therefore, a competition intensity coefficient α is introduced. m , representing the ratio of metal m's competitiveness to that of the reference metal: r m,i =α m ·exp(θ m,i ) / (1+∑ m'∈Gk α m' ·exp(θ m',i );

[0178] Where, α m A value greater than 0 can be set as a priori based on adsorption experimental data or the adsorption affinity sequence reported in the literature. For example, for the adsorption of cationic metals with clay minerals, the common affinity sequence is Pb>Cu>Zn>Cd, and α can be set accordingly. Pb >α Cu >α Zn >α Cd .

[0179] This extension enables the model to characterize the physical mechanisms of preferential site occupancy by strongly adsorbed metals, and to more accurately predict the competitive allocation of each metal under high metal loading conditions.

[0180] According to one aspect of this application, the capacity factor β g The method for determining it is as follows:

[0181] Based on surface complexation theory, the unit capacity coefficient β g It is decomposed into the product of site density and binding factor, giving it a clear physicochemical meaning.

[0182] For the g-th type of sediment component, its surface adsorption site density ρ g =(SSA g ·σ g ) / N A .

[0183] Among them, SSA g σ is the specific surface area (m² / kg) of the g-th component. g N represents the number of active sites per unit area. A is Avogadro's constant.

[0184] Table 3: Parameter values ​​of typical sediment components

[0185]

[0186] Considering the differences in affinity of different metals for the same type of site, a metal-specific binding factor φ is introduced. m,gφ represents the effective occupancy capacity of metal m at the g-th site: m,g =K m,g / K ref,g ;

[0187] Among them, K m,g K is the apparent binding constant between metal m and the g-th site. ref,g Reference metal (e.g., Cu) 2+ The binding constant of ). This ratio can be obtained from the NIST coordination chemistry database or thermodynamic models such as WHAM / Model VII.

[0188] Considering both site density and binding factor, the effective adsorption capacity of the g-th component for metal m is:

[0189] S m,g,i =ρ g ·φ m,g ·B g,i ;

[0190] Then the total effective capacity of the competing group k to which metal m belongs is:

[0191] S tot,i (k,m) =∑ g∈Ωk S m,g,i =∑ g∈Ωk ρ g ·φ m,g ·B g,i ;

[0192] At this time, β g Decomposed into ρ g ·φ m,g , where ρ g Derived from physical parameters, φ m,g All data are derived from thermodynamic databases and have clear physical origins.

[0193] Furthermore, a nonlinear model of capacity saturation effect can be introduced.

[0194] When the sediment composition is high, the growth of effective adsorption sites exhibits a saturation effect (e.g., organic matter aggregation leads to site shading). To characterize this phenomenon, a Michaelis-Menten type saturation function is used instead of a linear model:

[0195] S tot,i (k) =∑ g∈Ωk S max,g ·B g,i / K half,g +B g,i ;

[0196] Among them, S max,gK represents the maximum contributing capacity (mg / kg) of the g-th component. half,g It is the half-saturation constant, that is, when B g,i =K half,g At that time, this component contributes 50% of its maximum capacity.

[0197] This model has the following characteristics:

[0198] When B g,i <<K half,g When the model degenerates into a linear model: S≈S max,g / K half,g ·B g,i ;

[0199] When B g,i >>K half,g When the temperature approaches saturation: S≈S max,g ;

[0200] Avoid overestimating the capacity of samples with high organic matter or high iron content;

[0201] Parameter S max,g and K half,g It can be estimated jointly using Bayesian networks and other parameters, or pre-defined based on adsorption isotherm data reported in the literature. Specifically, before training the Bayesian network model, the parameters in the saturation response relationship need to be determined in advance, including the maximum contribution capacity S of each sediment component. max,g and the half-saturation constant K half,g This parameter is determined independently based on published adsorption isotherm experimental data or surface complexation model databases and does not participate in the parameter learning process of the Bayesian network. After the parameter is determined, the effective adsorption capacity S of each sample is calculated using the measured sediment composition content. tot,i The measured heavy metal content q m,i Divide by the corresponding effective adsorption capacity S tot,i Get the capacity utilization rate r m,i Then for r m,i A logit transformation is performed to obtain the target variable of the Bayesian network. This ensures that the calculation of the effective adsorption capacity is strictly separated from the structure learning and parameter learning of the Bayesian network in the process, avoiding circular dependencies.

[0202] Additionally, boundary cases need to be handled:

[0203] When calculating the capacity occupancy rate, if the measured heavy metal content is close to or exceeds the effective adsorption capacity, a truncation process is applied to the occupancy rate value, and r is reduced to a threshold value. m,i The value is limited to a preset open interval, such as (0.001, 0.999), to ensure the numerical stability of the logit transformation.

[0204] Furthermore, the hierarchical structure constraints are further defined as follows:

[0205] Based on the geochemical behavior of heavy metals in sediments, a hierarchical directed edge constraint system is constructed, including a whitelist (edges that must exist), a blacklist (edges that are prohibited from existing), and a graylist (edges to be learned).

[0206] Based on causal time sequence and chemical mechanism, Bayesian network variables are divided into three levels:

[0207] Table 4: Bayesian Network Variables

[0208]

[0209] Based on recognized geochemical mechanisms, the following whitelisted edges are set:

[0210] Table 5: Whitelisted Edges

[0211] Parent node child nodes Chemical basis pH Various cationic metals Hydrogen ion competition and hydroxide solubility control Eh As, Cr Redox Transformation iron oxide content As Specific adsorption of arsenic Organic carbon content Cu, Pb Organic complex clay content Cd, Zn Cation exchange

[0212] Whitelisted edges are forcibly retained in structure learning and are not included in the scoring deletion. Based on causal temporal order and chemical independence, the following blacklisted edges are defined:

[0213] Table 6: Blacklisted Sides

[0214] Forbidden border Reasons for constraint Metal → First-level variable Reverse causality (metal content cannot affect sediment composition) Metal → pH / Eh Except for extremely high levels of pollution, trace metals do not alter the overall acid-base redox state. Second floor → First floor Reverse cause and effect (pH cannot change clay content)

[0215] The above constraints are embedded into the Hill-Climbing algorithm, and the process is as follows:

[0216] Input dataset D, whitelist W, blacklist B;

[0217] Initialize the initial structure G0 containing all edges of W;

[0218] Perform the following procedure:

[0219] For each candidate edge :

[0220] Calculate the change in BIC score after adding / deleting / flipping e;

[0221] Retain the actions that maximize the improvement in score;

[0222] If an edge in W is selected for deletion, then skip that operation;

[0223] Continue until no action can improve the score; output the network structure G* that satisfies the constraints.

[0224] To address the problem of overestimation of adsorption capacity, this solution employs the Michaelis-Menten saturated response model and physical parameter decomposition technology. By introducing the maximum contributing capacity and the half-saturation constant, the traditional linear growth model is modified into a nonlinear saturation model, forcibly limiting the capacity growth limit under high organic matter or high clay content. Simultaneously, by calculating the product of the effective adsorption capacity and the capacity occupancy rate as the final predicted value, the mathematical structure rigorously ensures that the predicted heavy metal content never exceeds the physical adsorption upper limit, effectively solving the prediction distortion problem under high pollution loads.

[0225] To address the lack of a multi-metal competition mechanism, this scheme constructs a grouped occupancy model based on a competitive Langmuir mechanism. By dividing different heavy metals into competing groups sharing specific sites (such as cation groups and anion groups), and introducing latent variables of adsorption affinity and constant terms representing vacancies into the capacity occupancy rate calculation formula, the microscopic chemical processes of strongly adsorbed metals crowding out weakly adsorbed metals and adsorption only possible when vacancies exist are explicitly simulated. This overcomes the shortcomings of traditional statistical models that isolate the evaluation of single metals or rely solely on data correlation, and realistically reproduces the chemical equilibrium under complex systems.

[0226] To address the issue of spurious correlations arising from structure learning, this approach implements a hierarchical mechanistic constraint strategy. By defining a three-layer variable structure of material source-environmental state-heavy metal, and combining prior geochemical knowledge to set up a whitelist (mandatory causal) and a blacklist (prohibition of reverse engineering), erroneous connections that violate causal time order are directly eliminated during the model training phase, ensuring the scientific validity and interpretability of the Bayesian network structure.

[0227] In summary, this invention, by constructing an analytical framework that combines physical constraints and statistical inference, not only achieves high-precision prediction of heavy metal pollution, but also reveals the microscopic mechanism of multi-metal competition, providing a scientific basis for sediment pollution control.

[0228] The preferred embodiments of the present invention have been described in detail above. However, the present invention is not limited to the specific details of the above embodiments. Within the scope of the technical concept of the present invention, various equivalent transformations can be made to the technical solutions of the present invention, and these equivalent transformations all fall within the protection scope of the present invention.

Claims

1. A method for analyzing and evaluating heavy metal pollution in sediments, characterized in that, include: Obtain the physicochemical properties of sediment samples from the target water body; Based on physicochemical properties, the effective adsorption capacity of various adsorption components for heavy metals is calculated using a preset saturation response relationship. The physicochemical properties are input into a pre-constructed Bayesian network model to infer the latent adsorption affinity of each heavy metal element, where the latent adsorption affinity follows a conditional linear Gaussian distribution. Based on the latent variable of adsorption affinity, the capacity occupancy rate of each heavy metal element for adsorption sites is calculated using a preset competitive mapping function. The competitive mapping function introduces a vacancy competition mechanism through the denominator term to constrain the capacity occupancy rate within a preset range. The predicted heavy metal content is calculated based on the effective adsorption capacity and capacity occupancy rate, and the heavy metal pollution assessment results are output accordingly. The latent variable of adsorption affinity follows a conditionally linear Gaussian distribution, specifically including: Based on the parent node variable pointing to the heavy metal element in the Bayesian network structure, the calculation is performed according to the following linear Gaussian model: i m,i =a m,0 +∑ k=1 K (a m , k *X k,i )+e m,i ; Where, θ m,i Let a be the latent variable of the adsorption affinity of metal m in the i-th sample. m,0 X is the intercept term. k,i Let a be the observation value of the k-th parent node. m,k For the corresponding regression coefficient, ε m,i The noise term follows a Gaussian distribution; The capacity occupancy rate of each heavy metal element for adsorption sites is calculated using a preset competition mapping function, including: Based on the chemical properties of adsorption sites on sediment surfaces, various heavy metal elements are divided into several competing groups. Heavy metal elements belonging to the same competing group share the same type of adsorption sites and compete with each other, while heavy metal elements belonging to different competing groups correspond to different types of adsorption sites. For each competition group, the capacity utilization rate of each heavy metal element in the group is jointly calculated based on the adsorption affinity latent variable of all heavy metal elements in the competition group using the competition mapping function. The capacity utilization rate of each heavy metal element within this group was calculated jointly, specifically including: The competitive mapping function is constructed based on the competitive Langmuir adsorption mechanism; Calculate the capacity utilization rate r of each heavy metal element in this group. m,i =exp(θ m,i ) / 1+∑ m'∈Gk exp(θ m',i ); Where, r m,i For the i-th sample, belonging to the k-th competitive group G k The capacity utilization rate of heavy metal m, θ m,i Let ∑ be the latent variable of the adsorption affinity of the heavy metal m. m'∈Gk exp(θ m',i ) represents the exponential sum of the adsorption affinity potential variables of all heavy metal elements in this competition group, and the constant term 1 in the denominator is used to characterize the vacancies present in this type of adsorption site.

2. The method according to claim 1, characterized in that, Based on physicochemical properties, the effective adsorption capacity of various adsorbent components for heavy metals is calculated using a pre-defined saturation response relationship, including: Multiple sediment components providing the same type of adsorption sites were identified from physicochemical properties. Based on the saturation response relationship, the component contribution capacity of each sediment component was calculated separately; The effective adsorption capacity is obtained by summing the component contribution capacities of various sediment components.

3. The method according to claim 2, characterized in that, The saturation response relationship is the Michaelis-Menten saturation function, which specifically includes: Calculate the component contribution capacity S for each sediment component separately. g,i =S max,g ·B g,i / K half,g +B g,i ; Among them, S g,i B is the component contribution capacity of the g-th sediment component in the i-th sample. g,i S represents the measured content of this sediment component. max,g K represents the maximum contributing capacity of this sediment component. half,g is the half-saturation constant of the sediment component.

4. The method according to claim 2, characterized in that, The component contribution capacity of each sedimentary fraction is calculated separately, depending on the unit capacity coefficient of that sedimentary fraction, specifically including: The unit capacity coefficient is decomposed into the product of the surface adsorption site density and the metal-specific binding factor; The density of surface adsorption sites is calculated based on the specific surface area of ​​the sediment component and the number of active sites per unit area. The metal-specific binding factor is determined based on the ratio of the binding constant of the heavy metal to the sediment component to the binding constant of a reference metal.

5. The method according to claim 1, characterized in that, The competition mapping function also introduces a competition intensity coefficient to characterize the differences in affinity of different heavy metal elements for the same type of adsorption site, specifically including: Calculate the capacity utilization rate r of each heavy metal element in this group. m,i =α m *exp(θ m,i ) / 1+∑ m'∈Gk α m' *exp(θ m',i ); Where, α m Let α be the competition intensity coefficient of heavy metal m. m' is the competition intensity coefficient of heavy metal m' in the same group.

6. The method according to claim 1, characterized in that, When analyzing and evaluating a single heavy metal element, the pre-defined competition mapping function degenerates into the inverse Logit transform; m,i =1 / 1+exp(-θ m,i ); Constrain capacity utilization between 0 and 1.

7. The method according to claim 1, characterized in that, The predicted heavy metal content is calculated based on the effective adsorption capacity and capacity occupancy rate. Specifically, the product of the effective adsorption capacity and capacity occupancy rate is calculated so that the predicted heavy metal content automatically meets the physical constraint that it does not exceed the effective adsorption capacity.

8. The method according to claim 1, characterized in that, The pre-built Bayesian network model is obtained through the following steps: Historical observational data containing sedimentary physicochemical properties and heavy metal content were used as the training dataset. Based on the training dataset, a heuristic search algorithm is used to learn the network structure and determine the directed acyclic graph structure of the Bayesian network. Based on a defined directed acyclic graph structure, the conditional probability parameters in the network are learned using the maximum likelihood estimation method, resulting in linear Gaussian model parameters for inferring latent variables of adsorption affinity.

9. The method according to claim 8, characterized in that, Network structure learning is performed using heuristic search algorithms, including hierarchical constraint processes based on geochemical mechanisms: The node variables in the Bayesian network are divided into three levels: the material source variable level, the environmental state variable level, and the heavy metal content variable level. During the network structure learning process, directed edge constraints are imposed, including: allowing points from the material source variable layer to the environmental state variable layer or the heavy metal content variable layer, allowing points from the environmental state variable layer to the heavy metal content variable layer, and prohibiting directed edges originating from the heavy metal content variable layer.

10. The method according to claim 9, characterized in that, Applying directed edge constraints during network structure learning also includes setting whitelists and blacklists: The whitelist includes directed edges that must exist based on geochemical mechanisms, including directed edges pointing from pH value, redox potential, organic matter content, or iron and manganese oxide content to specific heavy metal element variables. The blacklist contains directed edges that are prohibited based on causal time sequence, including directed edges from the heavy metal content variable layer to the material source variable layer or the environmental state variable layer.

Images