Karst farmland heavy metal pollution migration mechanism identification and ecological risk assessment method and system

Abstract

The application discloses a kind of karst farmland heavy metal pollution migration mechanism identification and ecological risk assessment method and system, comprising: constructing standardized multidimensional space-time dataset;Establish a multi-medium coupling migration mechanism model describing soil-water-plant system, and generate event trigger connected network based on karst hydrological characteristics;The parameter field of mechanism model is constrained and inverted jointly using the network, and the initial prediction value is obtained;Based on the connected network, the prediction residual is decomposed into connected propagation component and local chemical component, closed-loop parameter correction and physical constraint projection are performed;Controlled disturbance is applied to the corrected model to extract counterfactual mechanism fingerprint, identify the dominant migration type and generate control strategy.The application solves the problem of parameter unidentifiable and complex migration path in karst area through connected constraint inversion and residual closed-loop return, realizes the explainable identification of pollution mechanism and the dynamic control of risk.

Description

Technical Field

[0001] This invention belongs to the interdisciplinary field of environmental science and information technology, especially the method and system for identifying the migration mechanism of heavy metal pollution in karst farmland and assessing ecological risks. Background Technology

[0002] In karst regions, the soil layer is shallow and karst fissures are well-developed, making cadmium and arsenic combined pollution in paddy fields highly insidious and sudden, threatening regional food security. Analyzing the nonlinear migration mechanism of heavy metals in dual media and quantifying the spatiotemporal dynamic ecological risk accordingly is a prerequisite for achieving safe utilization and zoned remediation of contaminated farmland.

[0003] Currently, the assessment and simulation of heavy metal pollution in farmland mainly rely on two types of methods: one is numerical simulation based on physicochemical processes, which usually uses convection-dispersion equations or geochemical models to describe the transport and transformation of solutes in the soil; the other is statistical or machine learning methods based on environmental monitoring data, such as using geostatistical interpolation, regression analysis or neural network models to analyze the spatial distribution characteristics of pollutants.

[0004] However, existing technologies face core challenges in dealing with the highly heterogeneous karst environment, namely, physical topological distortion and parameter-mechanism separation. Specifically, traditional modeling methods based on Euclidean distance or static geographical adjacency cannot characterize the instantaneous connectivity of fracture channels triggered by extreme rainfall events unique to karst regions, resulting in the absence of rapid preferential flow paths for pollutants across land parcels in the model. Furthermore, existing methods often use mechanistic models and data-driven algorithms in isolation. Simple residual correction only addresses numerical compensation at the output end, failing to feed back error signals from monitoring data into the parameter space of the multi-media coupled model. This prevents physical parameters from adaptively correcting under complex hydrological conditions, making it difficult to accurately identify dominant migration mechanisms such as preferential fracture flow and redox release in the absence of prior knowledge. Summary of the Invention

[0005] The purpose of this invention is to provide a method and system for identifying the migration mechanism of heavy metal pollution in karst farmland and assessing ecological risks, in order to solve one of the problems mentioned above in the existing technology.

[0006] Technical solution: A method for identifying the migration mechanism of heavy metal pollution in karst farmland and assessing ecological risks, including:

[0007] Acquire multi-source environmental monitoring data of the target area and construct a standardized multidimensional spatiotemporal dataset;

[0008] A multi-media coupled migration mechanism model describing the soil-water-plant system was established, and an event-triggered connectivity network was generated based on karst hydrological characteristics.

[0009] By using event-triggered connected networks, joint constraint inversion is performed on the parameter field of the multi-media coupled migration mechanism model to obtain the initial state prediction value;

[0010] Calculate the prediction residual between the initial state prediction value and the measured data, decompose the prediction residual based on the event-triggered connected network and perform closed-loop parameter correction to obtain the final prediction result that satisfies the physical constraints;

[0011] Controlled perturbations are applied to the corrected multi-media coupled migration mechanism model to extract counterfactual perturbation response fingerprints, and the dominant migration mechanism type of each evaluation unit is identified accordingly.

[0012] Based on the final prediction results and the dominant migration mechanism type, an ecological risk classification and management strategy is generated.

[0013] In one exemplary embodiment, a standardized multidimensional spatiotemporal dataset is constructed, including:

[0014] Based on the hydrogeological characteristics in multi-source environmental monitoring data, the minimum cumulative impedance path between each monitoring point and the point to be interpolated is calculated to obtain the connectivity impedance distance.

[0015] Using connectivity impedance distance as a measure of spatial correlation, co-kriging interpolation is performed on discrete sampling points in multi-source environmental monitoring data to generate a standardized multidimensional spatiotemporal dataset covering the target area.

[0016] In an exemplary embodiment, the multi-media coupled migration mechanism model includes mutually coupled matrix domain modules, fracture domain modules, and inter-domain mass exchange modules:

[0017] The matrix domain module is used to characterize Darcy flow and solute dispersion processes in soil porous media;

[0018] The fracture domain module is used to characterize non-Darcy preferential flow and solute convection processes in karst conduits and large fractures;

[0019] The interdomain mass exchange module is used to characterize the dynamic mass transfer process of solute between the matrix domain module and the fracture domain module based on concentration gradient.

[0020] In an exemplary embodiment, the matrix domain module and the fracture domain module respectively follow the convection-diffusion-reaction equation and are coupled through a mass exchange flux defined by the interdomain mass exchange module:

[0021] The governing equations of the matrix domain module include second-order partial derivatives describing the dispersion effect of porous media and source-sink terms describing matrix adsorption reactions; the governing equations of the fissure domain module include first-order partial derivatives describing the convection effect of preferred flow paths.

[0022] The interdomain mass exchange module calculates the mass exchange flux based on a first-order kinetic model. The mass exchange flux is proportional to the solute concentration difference between the matrix domain and the fracture domain and the interdomain mass exchange coefficient.

[0023] In an exemplary embodiment, the multi-media coupled migration mechanism model further includes a chemical reaction module for calculating reaction source and sink terms within each domain, specifically including:

[0024] An extended Langmuir competitive adsorption model was used to calculate the dynamic equilibrium between solid-phase adsorption capacity and liquid-phase concentration based on the competitive mechanism between cadmium ions and arsenate ions on the surface of iron and manganese oxides in soil.

[0025] The Fe-As coupled release model was used to establish the functional relationship between redox potential and iron oxide reduction rate, and the flux of adsorbed arsenic released during the reduction and dissolution of iron oxide was calculated accordingly.

[0026] A redox-driven morphological transformation model was adopted, and the reversible conversion rate between trivalent and pentavalent arsenic was calculated based on real-time redox potential.

[0027] In one exemplary embodiment, generating an event-triggered connected network based on karst hydrological features includes:

[0028] Construct the basic topology between evaluation units and define the edge impedance of the connecting edges, which is negatively correlated with the density of geological fractures;

[0029] An event-driven mechanism is introduced to dynamically update the edge impedance based on the cumulative rainfall and groundwater level fluctuation characteristics in the standardized multidimensional spatiotemporal dataset. When the rainfall or water level exceeds a preset threshold, the impedance value of the relevant connected edges is reduced.

[0030] The minimum cumulative impedance distance between evaluation units is calculated based on the updated edge impedance, and a time-varying connectivity weight matrix is ​​generated accordingly to construct an event-triggered connectivity network.

[0031] In an exemplary embodiment, joint constraint inversion is performed on the parameter field of a multi-media coupled migration mechanism model using an event-triggered connected network, including:

[0032] Construct a joint objective function that includes a data fitting term, a connectivity regularization term, and a prior constraint term;

[0033] The data fitting term is used to characterize the deviation between the model's simulated values ​​and the measured data;

[0034] The connectivity regularization term uses the time-varying connectivity weight matrix to impose a penalty constraint on the parameter differences of spatially adjacent evaluation units, forcing the parameters of highly connected regions to tend to be smooth, while allowing the parameters of weakly connected regions to remain abrupt.

[0035] Prior constraints are used to limit the inversion parameters to a physically reasonable range of values;

[0036] By minimizing the joint objective function, the model parameters of all evaluation units are solved simultaneously to obtain the initial state prediction values.

[0037] In an exemplary embodiment, decomposing the prediction residuals based on an event-triggered connected network and performing closed-loop parameter correction includes:

[0038] The prediction residuals are decomposed into connectivity propagation components caused by spatial network linkages and local chemical components caused by local processes.

[0039] Construct a parameter sensitivity matrix to inversely map the connectivity propagation component and the local chemical component into parameter correction increments in the model parameter space;

[0040] The parameter field of the multi-media coupled migration mechanism model is updated using parameter correction increments, and the updated model is used to recalculate to obtain the final prediction results.

[0041] In one exemplary embodiment, performing closed-loop parameter correction further includes performing physical constraint projection on the model output, specifically including:

[0042] By applying non-negative concentration constraints, negative concentration values ​​that violate physical meaning in the prediction results are forcibly truncated to zero or extremely small positive numbers;

[0043] Applying the regional mass conservation constraint, the total mass change of heavy metals within the evaluation area is calculated and compared with the boundary input and output fluxes. The residual deviation that violates the conservation law is evenly distributed to each evaluation unit within the area.

[0044] By applying the upper limit constraint of adsorption capacity, the upper limit of the predicted value of solid phase adsorption is truncated based on the maximum adsorption capacity parameter in the multi-media coupling migration mechanism model.

[0045] In an exemplary embodiment, a controlled perturbation is applied to the corrected multi-media coupled migration mechanism model to extract a counterfactual perturbation response fingerprint, including:

[0046] Key driving factors, including cumulative rainfall, redox potential, and preferential flow ratio, were selected.

[0047] While keeping other variables constant, small numerical perturbations are applied to the selected key driving factors, and the corrected multi-media coupled migration mechanism model is run to calculate the slope of the model output response to each perturbation.

[0048] The response slope combination is used to construct a counterfactual perturbation response fingerprint vector, which is then matched with a pre-defined mechanism template library to determine whether the evaluation unit belongs to the fracture preferential flow-dominated, rhizosphere enrichment, or redox interface release type.

[0049] Beneficial effects: This invention solves the problems of unidentifiable parameters and complex migration paths in karst areas by using connectivity constraint inversion and residual closed-loop backpropagation, thereby achieving interpretable identification of pollution mechanisms and dynamic and precise control of risks. Attached Figure Description

[0050] Figure 1 This is a diagram of a closed-loop assessment system that integrates multi-source sensing, physical-connectivity bidirectional coupling modeling, closed-loop parameter correction, counterfactual mechanism identification, and dynamic risk management in the embodiments of this application.

[0051] Figure 2 This is a flowchart illustrating the steps involved in constructing a standardized multidimensional spatiotemporal dataset in an embodiment of this application.

[0052] Figure 3 This is a diagram of the internal physical architecture of the multi-media coupling migration mechanism model in the embodiments of this application.

[0053] Figure 4 This is a flowchart illustrating the steps of generating an event-triggered connected network based on karst hydrological features in an embodiment of this application.

[0054] Figure 5 This is a flowchart illustrating the steps of performing joint constraint inversion on the parameter field of a multi-media coupled migration mechanism model using an event-triggered connected network in this embodiment of the application. Detailed Implementation

[0055] Example 1, such as Figure 1 As shown, this embodiment constructs a closed-loop assessment system that integrates multi-source sensing, physical-connectivity bidirectional coupling modeling, closed-loop parameter correction, counterfactual mechanism identification, and dynamic risk management, solving problems such as complex pollution migration paths, difficulty in identifying model parameters, and delayed risk assessment caused by the development of underground karst conduits in karst areas.

[0056] Step 101: Obtain multi-source environmental monitoring data of the target area and construct a standardized multidimensional spatiotemporal dataset.

[0057] In this embodiment, the target area specifically refers to karst farmland with a carbonate rock geological background, which typically exhibits dual-medium hydrological characteristics with frequent transformation between surface water and groundwater. Acquiring multi-source environmental monitoring data is the foundation of the entire assessment process, specifically including multi-dimensional information such as soil physicochemical properties, groundwater level and water quality, and crop growth indicators collected through ground sensor networks, as well as land cover and vegetation indices obtained through remote sensing. To eliminate differences in temporal frequency and spatial resolution between different data sources, the system performs spatiotemporal alignment operations, mapping discrete monitoring point data to a unified spatiotemporal grid, and performing outlier removal and missing value imputation, ultimately generating a standardized multi-dimensional spatiotemporal dataset.

[0058] Step 102: Establish a multi-media coupled migration mechanism model describing the soil-water-plant system, and generate an event-triggered connectivity network based on karst hydrological characteristics.

[0059] In this embodiment, the multi-media coupled migration mechanism model considers slow seepage and solute dispersion in the soil matrix, incorporates rapid preferential flow and convection processes in fissures or karst conduits, and the physicochemical behavior of pollutants at interfaces such as the soil solid-liquid interface and rhizosphere microzones. Meanwhile, considering the unique hydrological connectivity of karst regions—connected during the rainy season and disconnected during the dry season—this invention constructs an event-triggered connectivity network. This network is not a static geographic adjacency graph, but rather a dynamically evolving topology that changes with rainfall or groundwater levels. When a hydrological event occurs, disconnected distant plots may be reconnected through underground conduits, and the edge weights in the network are dynamically updated accordingly.

[0060] Step 103: Use the event-triggered connected network to perform joint constraint inversion on the parameter field of the multi-media coupled migration mechanism model to obtain the initial state prediction value.

[0061] In this embodiment, a joint constraint inversion mechanism is introduced. Instead of calibrating parameters for each evaluation unit in isolation, this mechanism uses the event-triggered connected network as a spatial constraint operator to construct a joint objective function containing a connectivity regularization term. Specifically, the system forces highly connected regions in the network to share similar parameter distributions, while allowing parameter independence in disconnected regions. By minimizing this joint objective function, the optimal parameter field across the entire domain can be solved simultaneously, and the mechanistic model can be run accordingly to output the initial state predictions of heavy metal concentration and flux for each evaluation unit. This process solves the problems of parameter overfitting and spatial non-transferability.

[0062] Step 104: Calculate the prediction residual between the initial state prediction value and the measured data. Decompose the prediction residual based on the event-triggered connected network and perform closed-loop parameter correction to obtain the final prediction result that satisfies the physical constraints.

[0063] This embodiment constructs a prediction-feedback closed-loop control loop. The system calculates the deviation between the model's predicted values ​​and the actual observed values, i.e., the prediction residual. Utilizing the topology of the event-triggered connected network, the total residual is decomposed into two parts: one part is the external connectivity error originating from upstream plots through fractured conduits, and the other part is the local error caused by deviations in local chemical reaction processes. Based on this decomposition, the system calculates targeted parameter correction increments, updates the physical parameters of the mechanistic model in reverse, and applies physical constraints such as nonnegativity and mass conservation to the updated output, obtaining a final prediction result that conforms to both statistical data patterns and physical laws.

[0064] Step 105: Apply controlled perturbation to the corrected multi-media coupled migration mechanism model to extract counterfactual perturbation response fingerprints, and identify the dominant migration mechanism type of each evaluation unit accordingly.

[0065] In this embodiment, a counterfactual inference method is employed to enhance model interpretability. The system conducts controlled variable experiments on the calibrated model; that is, while keeping other conditions constant, small perturbations are applied to key driving factors such as rainfall and redox potential, and the changes in the model output response are observed. The slope of this response constitutes a fingerprint vector reflecting the dynamic characteristics of the land parcel. By matching this fingerprint with preset mechanism templates (such as preferential flow-dominated and redox release-dominated), the system can automatically determine the dominant migration mechanism type of each assessment unit at the current moment.

[0066] Step 106: Based on the final prediction results and the dominant migration mechanism type, generate an ecological risk classification and control strategy.

[0067] Optionally, the system integrates quantitative prediction results (such as crop exceedance probability and groundwater leaching flux) with qualitative mechanism labels to conduct a multi-dimensional ecological risk assessment. Specifically, the system dynamically adjusts assessment weights based on climate conditions, for example, increasing the weight of groundwater risk under extreme rainfall events. Based on the risk level (such as high, medium, and low risk) and its dominant causes derived from the assessment, the system automatically generates tiered management strategies, such as recommending drainage interception measures in high-risk areas with priority flow patterns, or optimizing water management in redox risk areas, achieving site-specific prevention and control measures.

[0068] Example 2: This example describes the acquisition of multi-source data and a spatial interpolation method based on connectivity impedance. According to one aspect of this application, a standardized multidimensional spatiotemporal dataset is constructed, such as... Figure 2 As shown, it specifically includes:

[0069] Step 201: Based on the hydrogeological characteristics in the multi-source environmental monitoring data, calculate the minimum cumulative impedance path between each monitoring point and the point to be interpolated, and obtain the connectivity impedance distance.

[0070] In this embodiment, the acquisition of multi-source environmental monitoring data relies on a three-dimensional sensing network. Specifically, the ground monitoring system includes soil sensors deployed at different depths (e.g., layered burial at 0-20cm, 20-40cm, and 40-60cm), which transmit soil moisture content, conductivity, and redox potential data in real time using LoRa or NB-IoT wireless protocols; and water level and water quality sensors deployed in groundwater observation wells. Spatial data acquisition includes retrieving the Normalized Difference Vegetation Index (NDVI) from Sentinel-2 satellite multispectral imagery, and acquiring high-resolution surface temperature fields using a UAV equipped with a thermal infrared camera. When performing spatial interpolation, traditional straight-line distances cannot characterize the complex paths in karst regions that are blocked by rocks or connected by fissures. Therefore, this invention defines a connectivity impedance distance. The system constructs an impedance field based on fissure distribution maps and topographic slopes obtained from geological surveys, assigning low impedance values ​​to fissure-developed areas and high impedance values ​​to intact bedrock areas. For any two spatial points i and j, the system searches for all possible paths connecting the two points and calculates the minimum value of the integral of the impedance of each infinitesimal element along the path, which is the path with the minimum cumulative impedance. This distance d... ij Physically, it represents the ease or difficulty of a solute migrating from one point to another, rather than the geometric length of a straight line.

[0071] The formula for defining the foundation impedance is: R base_ij =L ij / (1+k frac *ρ frac_ij );

[0072] Among them, R base_ij To evaluate the base-side impedance between units i and j; L ij k is the Euclidean linear distance between the center points of the two units. frac ρ is the gap connectivity gain coefficient, used to adjust the degree to which the gap reduces impedance; frac_ij This represents the average crack line density along the path connecting the two units.

[0073] The formula for defining the impedance distance is: d imp_ij =min path (Σ(R e (t)));

[0074] Where, d imp_ij To evaluate the minimum cumulative connection impedance distance between units i and j; min path Σ represents the minimum value among all possible connection paths; Σ represents the summation of impedances over all segments along a given path; R e (t) represents the dynamic impedance of a line element e on the path at time t.

[0075] Step 202: Using the connectivity impedance distance as a measure of spatial correlation, perform co-kriging interpolation on discrete sampling points in the multi-source environmental monitoring data to generate a standardized multidimensional spatiotemporal dataset covering the target area.

[0076] In this embodiment, the system employs the Co-Kriging method to expand sparse sampling point data into grid data covering the entire domain. Optionally, when constructing the variation function, its independent variable is replaced from the traditional Euclidean distance h with the connectivity impedance distance d calculated above. ij The co-kriging estimation formula based on connectivity impedance is: Z0 * =Σ i=1 N λ i *Z i ;

[0077] Among them, Z0 * Z represents the attribute estimate of the point to be interpolated (0); N is the number of observation points in the neighborhood; i λ represents the measured attribute value of observation point i; i The weighting coefficients assigned to observation point i are obtained by solving the Kriging equations. The independent variable h in the semivariogram γ(h) of the equations is replaced by the aforementioned connectivity impedance distance d. imp_i0 In interpolation calculations, points that are spatially distant but connected by low-impedance fissures may have a higher correlation weight than points that are spatially close but blocked by high-impedance bedrock. For example, when using soil heavy metal point data to infer the concentration of unsampled points, the system will refer more to sample points with hydraulic connectivity. Furthermore, the system uses high-density auxiliary variables (such as vegetation indices retrieved from remote sensing) and low-density target variables (such as measured heavy metal content) to establish a covariance relationship, further improving interpolation accuracy. The final standardized multidimensional spatiotemporal dataset D contains a series of aligned spatiotemporal slices. Each slice D(t) is a matrix or tensor containing the attribute values ​​of each grid cell at time t. Specific fields include, but are not limited to, plot ID, sampling time, soil pH, redox potential Eh, iron and manganese oxide content, and 7-day cumulative rainfall Rain7, providing a data foundation for subsequent modeling.

[0078] Example 3: This example describes the internal physical architecture of the multi-media coupled migration mechanism model, such as... Figure 3 The diagram shows the key biochemical processes.

[0079] Step 301: The multi-media coupled migration mechanism model includes mutually coupled matrix domain modules, fracture domain modules, and inter-domain mass exchange modules.

[0080] In this embodiment, considering the dual-medium characteristics of karst soil, the model is designed as a dual-domain structure. The matrix domain module represents the microporous system in the soil, occupying most of the soil volume. Water flows slowly here, and chemical reactions such as adsorption and desorption mainly occur. The fissure domain module represents the macropores, fissures, or karst conduits in the soil. Although it occupies a small volume, it has strong water conductivity and is the main channel for the rapid migration of pollutants. The inter-domain mass exchange module is responsible for calculating the mass exchange flux between these two overlapping continuums, coupling the two independent flow fields.

[0081] Step 302: The matrix domain module is used to characterize Darcy flow and solute dispersion processes in the soil porous media; the fracture domain module is used to characterize non-Darcy preferential flow and solute convection processes in karst conduits and large fractures; the interdomain mass exchange module is used to characterize the dynamic mass transfer process of solute driven by concentration gradient between the matrix domain module and the fracture domain module.

[0082] As an optional implementation, the matrix domain module and the fracture domain module each follow the convection-dispersion-reaction equation and are coupled through the mass exchange flux defined by the inter-domain mass exchange module:

[0083] The governing equations of the matrix domain module include second-order partial derivatives describing the dispersion effect of porous media and source-sink terms describing matrix adsorption reactions; the governing equations of the fissure domain module include first-order partial derivatives describing the convection effect of preferred flow paths.

[0084] The interdomain mass exchange module calculates the mass exchange flux based on a first-order kinetic model. The mass exchange flux is proportional to the solute concentration difference between the matrix domain and the fracture domain and the interdomain mass exchange coefficient.

[0085] Specifically, the matrix domain module follows the convection-dispersion-reaction equation, and its calculation formula is as follows: ;

[0086] Where, θ m C represents the volumetric water content of the matrix domain. m The concentration of the solute in the liquid phase within the matrix domain; t is time; z is vertical depth; D m v is the hydrodynamic dispersion coefficient of the matrix domain; m Γ represents the average pore flow velocity in the matrix domain. m f is the inter-domain mass exchange flux per unit volume; R m These are the sources and sinks of biochemical reactions within the matrix domain.

[0087] The governing equations for the fracture domain module focus on rapid convection, and their calculation formulas are as follows: ;

[0088] Where, θ f C represents the volumetric water content of the fractured domain. fD represents the concentration of the solute in the liquid phase within the fracture domain. f v is the dispersion coefficient of the fractured domain; f The preferred flow velocity in the fractured domain is typically much greater than v. m ;R f This refers to the reaction source and sink terms within the fracture domain.

[0089] The exchange flux calculated by the inter-domain quality exchange module adopts a first-order dynamic model, as shown in the following formula: Γ m f=α w ×(C m -C f );

[0090] Among them, Γ m f is the mass exchange flux; α w This is the first-order mass exchange coefficient, expressed as the reciprocal of the daily rate. This coefficient is positively correlated with fracture development density, as shown in the following formula: α w =α0×(ρ f / ρ ref ) n ;

[0091] Where α0 is the reference commutation coefficient, ρ f ρ represents the fracture density of the current evaluation element. ref For reference fracture density, n is an empirical exponent, typically ranging from 1.5 to 2.5; C m C represents the matrix domain concentration. f This represents the concentration in the fractured domain.

[0092] Optionally, the multi-media coupled migration mechanism model also includes a chemical reaction module, which is used to calculate the source and sink terms of the reactions in each domain. Specifically, it includes: using an extended Langmuir competitive adsorption model, based on the competition mechanism between cadmium ions and arsenate ions on the surface of iron and manganese oxides in soil, to calculate the dynamic equilibrium between the amount of adsorption in the solid phase and the concentration in the liquid phase; using an Fe-As coupled release model, to establish a functional relationship between the redox potential and the reduction rate of iron oxides, and to calculate the flux of adsorbed arsenic released along with the reduction and dissolution of iron oxides; and using a redox-driven speciation model, based on the real-time redox potential, to calculate the reversible conversion rate between trivalent and pentavalent arsenic.

[0093] In this embodiment, the model incorporates a complex chemical mechanism to reflect the interaction between cadmium and arsenic. For the competitive adsorption process, the formula for calculating the solid-phase adsorption amount is: q Cd =(Q max ×K Cd ×C Cd ) / (1+K Cd ×C Cd +K As ×CAs );

[0094] Where, q Cd Q represents the amount of cadmium adsorbed by soil solids; max K represents the maximum adsorption capacity of the soil. Cd C is the adsorption affinity constant for cadmium ions; Cd This refers to the concentration of cadmium in the liquid phase; K As C is the adsorption affinity constant for arsenate ions; As The concentration of arsenic in the liquid phase is given. The denominator of this formula reflects the competitive effect of arsenic on cadmium adsorption sites.

[0095] For the Fe-As coupled release process, the model defines the arsenic release rate caused by the reduction of iron oxides, as follows: R rel =-k Fe ×f(Eh)×q AsFe ;

[0096] Among them, R rel k is the release rate of adsorbed arsenic. Fe q is the maximum reduction rate constant of iron oxides; AsFe The content of arsenic bound to iron oxides; f(Eh) is a switching function controlled by redox potential, specifically: f(Eh) = 1 / (1 + exp((Eh - Eh)) crit ) / S Eh ));

[0097] Where Eh is the real-time redox potential, in millivolts (mV); crit The critical potential threshold for iron reduction, for example, +100 millivolts; S Eh The smoothing parameter, which controls the steepness of the transition band, is in millivolts, the same as Eh. The smaller the SEh value, the steeper the transition of the function near the critical potential Ehcrit, and the closer the reduction reaction triggering is to a step characteristic. The specific value of SEh can be determined by nonlinear least squares fitting of f(Eh) based on the redox experimental data of the iron reduction process, and the preferred value range is 10~50mV.

[0098] Optionally, the multi-media coupled migration mechanism model also includes a root absorption module, which is used to calculate the flux of heavy metals from the rhizosphere soil into the plant. Specifically, this includes: using a two-substrate Michaelis kinetic model to describe the competitive absorption process of cadmium and arsenic by roots; introducing an interaction inhibition constant to characterize the competitive effect of cadmium ions and arsenic on root transport proteins, so that an increase in the rhizosphere concentration of one heavy metal will inhibit the absorption rate of the other heavy metal; and introducing a rhizosphere enrichment factor to correct the rhizosphere solution concentration based on the oxidative microenvironment characteristics of the rhizosphere microregion, which serves as the input to the two-substrate Michaelis kinetic model.

[0099] To predict the cumulative risk to crops, the model incorporates a two-substrate Michaelis equation at the biological interface. Taking cadmium uptake flux as an example, the calculation formula is: J Cd =(V max ×C rz_Cd ) / (K m_Cd +C rz_Cd ×(1+C rz_As / K i ));

[0100] Among them, J Cd V represents the cadmium uptake flux per unit surface area of ​​the root system; max C is the maximum absorption rate. rz_Cd Cadmium concentration in rhizosphere soil solution; K m_Cd C is the Michaelis constant, representing the concentration at which half of the maximum absorption rate is reached; rz_As The concentration of arsenic in the rhizosphere soil solution; K i is the interaction inhibition constant of arsenic to cadmium absorption.

[0101] This formula indicates that when the rhizosphere arsenic concentration C rz_As As the concentration increases, the denominator increases, leading to an increase in the cadmium absorption flux J. Cd The concentration is reduced. Meanwhile, considering the enrichment or blocking effect of the iron film formed by oxygen secretion from rice roots on heavy metals, a correction factor is introduced into the model when calculating the rhizosphere concentration. The formula is: C rz =C bulk ×ξ;

[0102] Among them, C bulk ξ represents the concentration of the non-rhizosphere soil solution; ξ is the rhizosphere enrichment factor, which is the rhizosphere redox potential Eh. rz The function completes the final segment of the migration path from soil to plant.

[0103] Example 4: This example describes the construction of a spatiotemporal topology that adapts to the connectivity of karst during the rainy season and the disconnection during the dry season, as well as the sensitivity screening of complex physical model parameters beforehand.

[0104] Specifically, before performing joint constraint inversion, sensitivity screening of the initial parameter set of the multi-media coupled migration mechanism model is also included.

[0105] In this embodiment, given that the multi-media coupling model contains dozens of physicochemical parameters, inverting the entire set would lead to overfitting and computational divergence; therefore, it is necessary to first lock the core parameters. The system calculates the normalized error sequence e(t) based on historical simulation error data. Robust statistical methods are introduced to calculate the median(e) and median absolute deviation (MAD)(e) of the normalized error sequence. Based on these two statistics, the system defines a first screening threshold ε1 and a second screening threshold ε2, with the specific calculation formulas being: ε1 = a1 * median(e), ε2 = a2 * MAD(e);

[0106] Where a1 and a2 are preset multiplier coefficients, for example, 0.5 and 2.0 respectively. For each undetermined parameter θ in the model... k The system calculates its sensitivity index S. k Only when S k >ε1 and S k Only when the fluctuation amplitude is less than ε2 is this parameter retained in the subset Ω of the core parameters to be inverted. opt This dual filtering mechanism removes redundant parameters that have little impact on the output or cause extreme instability in the model, thus reducing the dimensionality of the inversion problem.

[0107] Step 401, as follows Figure 4 As shown, an event-triggered connected network is generated based on karst hydrological characteristics. Specifically, this includes: constructing the basic topology between evaluation units and defining the edge impedance of the connecting edges, which is negatively correlated with the density of geological fissures.

[0108] In this embodiment, the system constructs a static basic topology graph G. base =(V,E base ). Where node V represents a farmland assessment unit, and edge E base The existence of a connection is determined by both geographical adjacency and the orientation of potential underground fissures. For each connecting edge (i,j), the system defines its basic edge impedance R. ij The impedance is not a simple Euclidean distance, but a physical measure based on the degree of fracture development. Specifically, the impedance R... ij ,base and the crack density ρ along this path f They are negatively correlated, and the calculation formula can be expressed as: R ij ,base=L ij / (1+k*ρ f );

[0109] Among them, L ij Let be the straight-line distance between the centers of the two plots, and k be the fracture connectivity coefficient. In areas with dense fracture development, even if the physical distance is far, the connectivity impedance may be very low, presupposing potential fast migration channels.

[0110] Step 402 introduces an event-driven mechanism to dynamically update the edge impedance based on the cumulative rainfall and groundwater level fluctuation characteristics in the standardized multidimensional spatiotemporal dataset. When the rainfall or water level exceeds a preset threshold, the impedance value of the relevant connected edges is reduced.

[0111] Furthermore, this embodiment implements event-triggered characteristics. The system monitors driving factors in the standardized dataset in real time and employs different activation function mechanisms for matrix flow connection edges and fissure flow connection edges. For matrix flow connection edges, the system uses a power function activation mechanism based on soil moisture content. As the soil moisture content θ approaches the saturation moisture content θ... s Edge weight w ij m The reciprocal of impedance increases exponentially, and its calculation formula is: w ij m (t)=w max *((θ(t)-θ r ) / (θ s -θ r )) β ;

[0112] Where, θ r θ represents the residual soil moisture content. s Let β represent the soil saturation water content, and β be the pore connectivity index, typically taken as 3-4. For fracture flow connection edges, the system employs a threshold-based step activation mechanism. Activation only occurs when the 7-day cumulative rainfall (Rain7) or the groundwater level (Level) exceeds a preset trigger threshold (P). th Only when the high conduction weight of the fracture flow connection edge is activated will it remain in a low conduction state or an open state; otherwise, it will remain in a low conduction state. The calculation formula for this mechanism is: w ij f (t)=w high *H(Rain7(t)-P th )+w low ;

[0113] Where w ij f (t) represents the dynamic weight of the fracture flow connection edge; w high This represents the high-weight value when the fracture is conducting; w low The low-weight value is the value at which the fracture closes; H(·) is the Herveyd step function, which takes the value 1 when the independent variable is greater than 0, and 0 otherwise; Rain7(t) is the cumulative rainfall over 7 days; P th The rainfall threshold used to trigger fissure flow is defined. This type of design captures the nonlinear hydrological response of karst regions where heavy rains result in high flow and light rains result in low flow.

[0114] Step 403: Calculate the minimum cumulative impedance distance between evaluation units based on the updated edge impedance, and generate a time-varying connectivity weight matrix accordingly to construct an event-triggered connectivity network.

[0115] Accordingly, based on the dynamically updated impedance field described above, the system calculates the minimum cumulative impedance distance d between any two points at each time step t using either the Dijkstra algorithm or the Floyd-Warshall algorithm. ij (t). This distance is transformed into a time-varying connectivity weight matrix W(t) using a negative exponential decay function, as shown in the formula: W ij (t)=exp(-d ij (t) / σ);

[0116] Where σ is the connectivity decay feature scale. The generated W(t) matrix is ​​not only sparse but also highly dynamic, with high-weight elements W ij (t) indicates the pairs of land parcels with strong hydraulic connections under the current hydrological conditions, which constitute the spatial constraint operator for subsequent joint inversion.

[0117] Example 5: This example describes how to use dynamic networks to resolve the contradiction between parameter space heterogeneity and unidentifiability.

[0118] Step 501, as follows Figure 5 As shown, the parameter field of the multi-media coupled migration mechanism model is jointly constrained and inverted using an event-triggered connected network. Specifically, this includes: constructing a joint objective function that includes a data fitting term, a connectivity regularization term, and a priori constraint terms; the data fitting term is used to characterize the deviation between the model's simulated values ​​and the measured data.

[0119] In this embodiment, in order to simultaneously solve the model parameter field Θ={θ1,θ2,...,θ} for all evaluation units in the global domain, N The system constructs a global joint objective function J(Θ). This function consists of three parts, and its calculation formula is: J(Θ) = J data +ρ*J reg +η*J prior In calculating the connectivity regularization term J reg and prior constraint term J prior At that time, each parameter subvector θ i The different physical quantity parameters included must undergo dimensionless preprocessing before participating in norm calculations. Specifically, for the k-th parameter component θ... k The following normalization transformation is adopted: When θ prior,k ≠0 o'clock or When θ prior,k When θ = 0, scale,k The dimensional characteristic scale of this parameter can be taken as the maximum value within its physically permissible range; the normalized dimensionless parameter vector Used for norm calculations to ensure dimensional consistency.

[0120] Among them, J data This is the data fitting term, used to characterize the simulated value y of the model. pred Compared with measured data y obs The deviation between them is usually expressed in the form of mean square error (MSE): J data =∑ i ||y obs,i -y pred,i (θ i )|| 2 ρ and η are the connectivity regularization coefficient and prior constraint coefficient, respectively, used to balance the accuracy of data fitting with the rationality of the parameter space structure.

[0121] Step 502: The connectivity regularization term applies a penalty constraint to the parameter differences between spatially adjacent evaluation units using a time-varying connectivity weight matrix, forcing the parameters of highly connected regions to tend to be smooth, while allowing the parameters of weakly connected regions to remain abrupt.

[0122] Specifically, the connectivity regularization term J reg The construction directly utilizes the time-varying connectivity weight matrix W(t). Its calculation formula is: J reg =∑t∑{i,j}W ij (t)*||θ i -θ j || 2 The physical meaning of this formula is that when the connectivity weight W between two plots i and j... ij When (t) is large, indicating the existence of active fracture channels, minimizing the objective function will force the physical parameters θ of both to change. i and θ j The parameters tend to be consistent, achieving a smooth transition within the connected region; conversely, when W... ij When (t) approaches zero, that is, when the hydraulic connection between the plots is broken, ||θ i -θ j || 2 The penalty weights of the terms disappear, allowing for differences in the parameters between the two. This mechanism adapts to the homogeneous geological features and abrupt boundary changes within karst formations, avoiding excessive parameter homogenization caused by traditional Laplace smoothing.

[0123] Step 503: Prior constraints are used to limit the inversion parameters to a physically reasonable range; by minimizing the joint objective function, the model parameters of all evaluation units are solved simultaneously to obtain the initial state prediction values.

[0124] As an optional implementation, the prior constraint J prior The Tikhonov regularization form is usually used, i.e. To prevent the parameters from deviating from the prior mean θ given by the geological survey. prior Too far. The system employs gradient-based optimization algorithms (such as L-BFGS-B) or intelligent evolutionary algorithms (such as CMA-ES) to minimize the overall objective function J(Θ). During the iterative solution process, the selected subset of core parameters is continuously updated until convergence. The final parameter field Θ is obtained. opt Substituting the values ​​into the multi-media coupling mechanism model, the initial state prediction values ​​y of each evaluation unit in the uncorrected state are obtained. init .

[0125] The objective function formula for the joint inversion is: J(Θ)=Σ i ||y obs_i -y pred_i (θ i )|| 2 +ρ*Σ i,j W ij (t)*||θ i -θ j || 2 +η*||Θ-Θ prior || 2 ;

[0126] Where J(Θ) is the total objective function value; Θ is the global parameter set; the first term is the data fitting term, y obs_i y is the measured value. pred_i (θ i The first term is the model prediction; the second term is the connectivity regularization term, where ρ is the regularization coefficient and W is the model prediction. ij (t) is the time-varying connectivity weight matrix, θ i and θ j These are the parameter vectors of adjacent elements; the third term is the prior constraint term, where η is the prior coefficient and Θ is the parameter vector of the adjacent element. prior Let be the prior mean of the parameter.

[0127] Example 6: This example extends error compensation to mechanism correction by converting the residual signal into a correction value for physical parameters. According to one aspect of this application, the predicted residual is decomposed based on an event-triggered connected network, and closed-loop parameter correction is performed, specifically including:

[0128] Step 601: Calculate the prediction residual between the initial state prediction value and the measured data. Decompose the prediction residual based on the event-triggered connected network, and decompose the prediction residual into the connectivity propagation component caused by spatial network linkage and the local chemical component caused by local processes.

[0129] In this embodiment, the system calculates the total prediction residual r. total (t)=y obs (t)-y init(t). To trace the cause of the error, the system uses the time-varying connectivity weight matrix W(t) to spatially decouple the residual.

[0130] Specifically, the elements W in the connection weight matrix W(t) ij Defined as the hydraulic connectivity weight from unit i to unit j, i.e., with i as the source and j as the target. For any evaluation unit p, its upstream neighborhood set N(p) = {q|W qp The set of cells (t)>0 represents the set of cells from q to p with a valid connected path. The formula for calculating the connected propagation component of p is:

[0131] r conn,p (t)=Σ q∈N(p) W qp (t)*r total,q (t);

[0132] Where, r conn,p (t) represents the connectivity propagation residual component of the evaluation unit p; W qp (t) represents the connectivity weight from upstream cell q to current cell p in the time-varying connectivity weight matrix; r total,q (t) represents the total prediction residual of the upstream unit q. For example, for unit B, its upstream neighborhood N(B) = {A}, W AB =0.8, then r conn,B =W AB *r total,A =0.8×2.0=1.6.

[0133] Step 602: Construct a parameter sensitivity matrix, and inversely map the connectivity propagation component and the local chemical component into parameter correction increments in the model parameter space; use the parameter correction increments to update the parameter field of the multi-media coupling migration mechanism model, and recalculate using the updated model to obtain the final prediction result.

[0134] Optionally, the prediction residual is decomposed into connectivity propagation components and local chemical components, specifically including:

[0135] Based on the time-varying connectivity weight matrix in the event-triggered connected network, the prediction residuals of the neighboring evaluation units are weighted and summed to calculate the connectivity propagation component of the current evaluation unit, which represents the external input error transmitted through the gap or pipe.

[0136] A nonlinear mapping function is constructed based on the local state characteristics of the current assessment unit, or the local chemical component is calculated by removing the connectivity propagation component from the total prediction residual, which characterizes the internal reaction error caused by redox potential fluctuations or adsorption / desorption processes.

[0137] To achieve closed-loop control, the system needs to transform the residual *r* in the data space into a correction variable *Δθ* in the parameter space. The system constructs a Jacobian matrix (sensitivity matrix) *S* based on the Taylor expansion principle, where the elements... This represents the effect of the change in the nth parameter on the mth output variable. Using the iterative approach of the Gauss-Newton method, the formula for calculating the parameter correction increment Δθ is: Δθ = (S... T *S+λ*I) -1 *S T *r local ;

[0138] Where Δθ is the parameter correction increment vector; S is the Jacobian sensitivity matrix, i.e., the partial derivative matrix of the output with respect to the parameters; S T Let S be the transpose of S; λ is the damping coefficient, used to ensure the stability of matrix inversion; I is the identity matrix; r _local This represents the local chemical residual component. A damping factor λ*I is introduced to ensure the stability of matrix inversion. Through this mapping, the system calculates the value to eliminate the local error r. local The required parameter adjustment amount, and its addition to the original parameter: θ new =θ opt +Δθ. Using the updated parameter θ new By running the mechanistic model again, we can obtain the prediction results that have been corrected by the physical layer.

[0139] Optionally, performing closed-loop parameter correction also includes performing physical constraint projection on the model output results, specifically including: applying non-negative concentration constraints to forcibly truncate negative concentration values ​​in the prediction results that violate physical meaning to zero or a very small positive number; applying regional mass conservation constraints to calculate the total mass change of heavy metals in the evaluation region and verify it with the boundary input-output flux, and evenly distribute the residual deviations that violate the conservation law to each evaluation unit in the region; and applying an upper limit constraint on adsorption capacity to truncate the predicted value of solid phase adsorption based on the maximum adsorption capacity parameter in the multi-media coupled migration mechanism model.

[0140] As the final step in closed-loop calibration, to prevent over-tuning of parameters from causing the model output to violate physical laws (e.g., negative concentration), the system introduces a Physics Projection Layer at the output. For non-negative concentration constraints, the system performs a truncation operation: C final =max(0,C updated For the regional mass conservation constraint, the system calculates the change in total mass ΔM within the region. totalThe difference between the net flux and the boundary flux (NetFlux) is weighted and distributed back to each cell if a mass imbalance exists. For the upper limit constraint on adsorption capacity, the system checks whether the solid-phase adsorption capacity (q) exceeds the maximum adsorption capacity (Q) for that soil type. max , if q>Q max Then force q=Q max The excess heavy metals are returned to the liquid phase. Through this series of hard-constraint projections, the final output prediction results are ensured to be numerically accurate and physically reliable.

[0141] Example 7: This example describes the extraction of physically meaningful mechanism features from a corrected model using counterfactual inference. According to one aspect of this application, a controlled perturbation is applied to the corrected multi-medium coupled migration mechanism model to extract a counterfactual perturbation response fingerprint, specifically including:

[0142] Step 701: Select key driving factors including cumulative rainfall, redox potential, and preferential flow ratio coefficient.

[0143] In this embodiment, key driving factors refer to the set of variables in the input data vector or model parameters that influence the migration and transformation process of heavy metals. Specifically, the system pre-defines a set of driving factors U={u1,u2,...,u...} k} Where u1 can be selected as the 7-day cumulative rainfall Rain7, used to detect the hydrodynamically driven rapid migration mechanism; u2 can be selected as the soil redox potential Eh, used to detect the chemically driven dissolution and release mechanism; u3 can be selected as the preferential flow ratio coefficient w in the model. f This system is used to detect the contribution of fractured conduits. Furthermore, it allows users to add custom factors based on prior knowledge, such as soil pH or iron and manganese oxide content. These factors were selected because different migration mechanisms often exhibit extremely high sensitivity to specific environmental factors; for example, preferential flow mechanisms are highly sensitive to heavy rainfall, while reduction-release mechanisms are highly sensitive to a decrease in Eh.

[0144] Step 702: While keeping other variables constant, apply small numerical perturbations to the selected key driving factors, run the corrected multi-media coupled migration mechanism model, and calculate the slope of the model output response to each perturbation.

[0145] In this embodiment, the system performs a control variable experiment, i.e., a counterfactual perturbation operation. For each evaluation unit i and each selected driving factor u... k The system constructs a counterfactual input scenario: while keeping u as the remainder of the input... k With all other input variables and model parameters remaining unchanged, give u k A tiny increment δuk (For example, 1% of the original value or 10% of one standard deviation). Furthermore, using the already parameter-corrected multi-media coupled migration mechanism model, simulations are run under both the original input and counterfactual input conditions to obtain the original output y. i Output y after perturbation i The system uses the finite difference method to calculate the slope s of the model's response to this factor. i k is calculated using the formula: s i k=(y i '-y i ) / δu k ;

[0146] Among them, s i k is the slope of the response of the i-th unit to the k-th driving factor; y i ' is the output function of the corrected model; u k δ represents the original value of the k-th driving factor; δ is the applied small perturbation. The slope s i k quantitatively characterizes the instantaneous sensitivity of heavy metal migration flux or concentration to the driving factor under the current state. For example, if a plot of land has a very large response slope to Rain7, it indicates that rainfall is the main controlling switch triggering the spread of pollution in that plot.

[0147] Step 703: Construct the counterfactual perturbation response fingerprint vector by combining the response slopes, and match it with the preset mechanism template library to determine whether the evaluation unit belongs to the fracture preferential flow-dominated type, the rhizosphere enrichment type, or the redox interface release type.

[0148] In this embodiment, for each evaluation unit i, the system combines its response slopes to all key driving factors into a feature vector, namely the counterfactual perturbation response fingerprint vector m. i =[s i 1,s i 2,...,s i k]. To identify the dominant mechanism, the system pre-configures a mechanism template library T=t. type1 ,t type2 ...。 Among them, t type1 The template vector for (fracture-preferred flow-dominated) types has a high value in the rainfall factor dimension; t type2 The template vector for (redox interface release type) has a high negative value in the Eh factor dimension, indicating that a decrease in Eh leads to an increase in release; t type3 The template vector (rhizosphere enriched type) has high values ​​in the root biomass or transpiration rate factor dimensions. The system uses the cosine similarity algorithm to calculate the measured fingerprint m. i With each template vector t j Similarity between Sim ijThe calculation formula is: Sim ij =(Σ k (m ik *t jk )) / (sqrt(Σ k m ik 2 )*sqrt(Σ k t jk 2 ));

[0149] Among them, Sim ij For actual fingerprint m i With the j-th type mechanism template t j Similarity; m ik t represents the k-th component of the measured fingerprint vector; jk This is the k-th component of the mechanism template vector. The system will use the similarity Sim... ij The largest template type is determined as the dominant migration mechanism type of the evaluation unit. This process transforms data response features into mechanism labels with physical meaning, achieving interpretability of mechanism identification.

[0150] Example 8: This example describes dynamic risk assessment and intelligent governance recommendations based on prediction results and mechanism labels. According to one aspect of this application, an ecological risk classification and management strategy is generated based on the final prediction results and the dominant migration mechanism type, specifically including:

[0151] Optionally, based on the final prediction results, multi-dimensional risk indicators including the probability of heavy metal contamination in crops, groundwater leaching flux, and health hazard quotients in the population can be calculated.

[0152] In this embodiment, the system calculates risk indicators in three dimensions based on the final prediction results that satisfy physical constraints. The probability P of heavy metal contamination in crops is one of these indicators. crop This refers to the predicted heavy metal concentration C in crop grains. grain Exceeding the national food safety standard limit C limit The probability distribution of groundwater leaching flux F is typically determined using Monte Carlo simulations or confidence intervals of predicted values. gw This refers to the mass of heavy metals entering the groundwater system through the bottom boundary of a soil profile per unit time. The formula for its calculation is: F gw =C leach *q perc Among them, C leach q represents the concentration of the leached water. perc This refers to the volume of water seepage from deep layers.

[0153] The Health Hazard Quotient (HQ) is calculated using the model recommended by the USEPA (U.S. Environmental Protection Agency), with the formula: HQ = (C grain*IR*EF*ED) / (BW*AT*RfD); where C grain The values ​​represent heavy metal concentrations in crop grains, IR (intake rate), EF (exposure frequency), ED (exposure duration), BW (body weight), AT (average time), and RfD (reference dose). These indicators form the basis for a comprehensive assessment of ecohealth risks.

[0154] Optionally, the basic objective weights of each risk indicator can be calculated using the entropy weight method, and extreme climate events in the standardized multidimensional spatiotemporal dataset can be monitored in real time.

[0155] In this embodiment, the system uses the entropy weight method to determine the basic weight W of each risk indicator under normal conditions. base =[w crop ,w gw ,w hq The entropy weighting method objectively assigns values ​​to each indicator based on its spatial variability, with indicators exhibiting greater variability receiving higher weights. Simultaneously, the system monitors meteorological data streams in real time to identify extreme weather events. Specifically, the system sets an extreme rainfall threshold P. extreme This threshold can be defined as the 95th percentile of daily rainfall over the past 30 years. When the monitored 7-day cumulative rainfall Rain7 exceeds this threshold, the system determines that the system is currently under extreme hydrological stress.

[0156] Optionally, when the cumulative rainfall exceeds the historical high percentile threshold, a weight adjustment mechanism is automatically triggered to increase the weight ratio of the groundwater leaching flux index in the comprehensive assessment and generate dynamic risk weights. Based on the dynamic risk weights, a fuzzy comprehensive evaluation of multi-dimensional risk indicators is performed to determine the risk level of the assessment unit, and a corresponding ecological risk classification and control strategy is generated accordingly.

[0157] This embodiment achieves dynamic adaptation in risk assessment. In karst regions, extreme rainfall can sometimes cause fissure conduits to open instantaneously, drastically increasing the risk of groundwater contamination, while the risk of crop absorption is relatively delayed. Therefore, once an extreme event is triggered, the system automatically adjusts the weight vector, using a preset adjustment coefficient γ (γ>1) to increase the weight w of groundwater leaching flux. gw Simultaneously, the weight vector is normalized to generate dynamic risk weight W. dynamic (t). For example, the adjusted groundwater weight is calculated using the following formula: w gw_new =(w gw_base *γ) / (w crop +w gw_base *γ+w hq Correspondingly, the normalized values ​​for the other two weights are: w crop_new =w crop / (w crop +wgw_base *γ+w hq ), w hq_new =w hq / (w crop +w gw_base *γ+w hq ), where w gw_new w crop_new w hq_new These are the normalized weights corresponding to the adjusted groundwater leaching flux, crop exceedance probability, and human health hazard quotient; w gw_base w crop w hq The weights are the base weights determined by the entropy weight method before adjustment; γ is the enhancement coefficient triggered by extreme climate, γ>1. The sum of the three adjusted weights is always equal to 1, satisfying the normalization constraint.

[0158] This mechanism ensures that risk assessment results can capture acute risks brought about by sudden environmental changes, avoiding the lag of static assessments. Based on this dynamic weighting, the system uses the fuzzy comprehensive evaluation method to calculate the comprehensive risk index ERI and divides the assessment units into four levels: safe, low-risk, medium-risk, and high-risk.

[0159] Optionally, a risk-mechanism knowledge graph can be constructed to encapsulate the land parcel attributes, dominant migration mechanism type, risk level, and governance measures of the assessment unit into structured triples for storage.

[0160] In this embodiment, the system establishes a professional graph database. This graph stores historical case knowledge in the form of triples of evaluation unit-attribute-value. Entity nodes include specific land parcel objects, environmental characteristics (such as pH value, organic matter), dominant mechanisms (such as preferred flow patterns), risk levels (such as Level IV), and verified remediation measures (such as passivating agent application, water control). Relationship edges include logical connections such as belonging, causing, and applicable.

[0161] Optionally, a graph embedding algorithm is used to map the current evaluation unit into a low-dimensional vector, and the historical case node with the highest cosine similarity is retrieved in the risk-mechanism knowledge graph; the governance measures suggestions associated with the retrieved historical case nodes are output as a supplementary recommendation scheme for the ecological risk classification and control strategy.

[0162] To achieve intelligent recommendation, the system employs graph embedding algorithms (such as TransE or GraphSAGE) to map all nodes in the knowledge graph to a low-dimensional vector space. For the current assessment unit to be addressed, the system generates a query vector V based on its attribute features (such as mechanism labels, risk levels, and environmental parameters). query Furthermore, the system searches in the vector space for V. queryThe system identifies the top-N historical case nodes with the highest cosine similarity. It then extracts successful remediation measures associated with these similar historical cases, weights them according to their similarity scores, and generates a recommendation list. For example, if the current site is identified as acidic soil with redox release and high risk, the system might recommend a combined strategy of applying an alkaline passivating agent and intermittent flooding management.

[0163] Example 9 illustrates the specific computational process of residual decomposition, closed-loop correction, and connectivity construction in this example through a specific numerical calculation scenario, in order to further demonstrate the feasibility of the solution.

[0164] Suppose there are three spatially adjacent evaluation units in a microkarst farmland system, labeled as unit A, unit B, and unit C. Unit A is located at a higher elevation, unit B is downstream of unit A, and unit C is lateral to unit B.

[0165] In step S1, during the data preparation phase, the system obtains the measured heavy metal concentration vector y at a certain time t. obs =[10.0,12.0,5.0] (unit: mg / kg).

[0166] In step S2, the modeling and initial value prediction stage, the system uses the uncorrected mechanistic model to calculate the initial prediction vector y. init =[8.0,8.0,4.5]. At this point, the total predicted residual vector r total =y obs -y init =[2.0,4.0,0.5].

[0167] During the network construction phase, assuming it is currently the rainy season, the time-varying connectivity weight matrix W(t) calculated by the system based on rainfall and terrain fissure data is as follows: W AA =0,W AB =0.8,W AC =0.1,W BA =0,W BB =0,W BC =0.2,W CA =0,W CB =0,W CC The matrix =0 indicates that cell A has strong connectivity to cell B (weight 0.8, possibly indicating the existence of gap channels), while the connectivity between A and C, and between B and C, is weak.

[0168] During the residual decomposition stage, the connectivity propagation component r of system computation unit B conn B. According to the formula, r conn B≈W AB *r totalA = 0.8 * 2.0 = 1.6. This means that of the total residual of 4.0 in element B, 1.6 is due to the error (2.0) from upstream element A propagating through the fracture network. The system calculates the local chemical component r of element B. local B=r total ,Br conn B = 4.0 - 1.6 = 2.4.

[0169] During the parameter calibration phase, for unit B, the system identifies a mixed source of error: one part originates from external input, requiring inspection of upstream A, and the other part (2.4) originates from the local process. The system utilizes the sensitivity matrix S B The maximum adsorption capacity Q in the local parameters was found. max It has the highest sensitivity to output. Through inverse mapping, the system calculates the Q value of unit B that needs to be adjusted. max The parameters were lowered to account for the higher concentration in this part, with a predicted value of 8.0 < measured value of 12.0. Meanwhile, for unit A, the error is entirely localized, and as the source of this error, the system will significantly adjust the parameters of unit A.

[0170] This numerical example demonstrates that our method can separate external input errors from local model errors. Traditional methods might incorrectly attribute all errors in unit B (4.0) to its own parameter issues, leading to erroneous parameter correction, i.e., overfitting. Our method, through connected network decomposition, identifies a portion of the error originating upstream, achieving a correction more consistent with physical facts.

[0171] Example 10: This example describes the use of Physically Constrained Graph Residual Network (PC-GRN) to achieve nonlinear residual fitting and physical constraint projection, which is suitable for scenarios with large amounts of data and nonlinear characteristics.

[0172] In this alternative approach, performing closed-loop parameter correction specifically involves constructing a deep neural network containing graph convolutional layers and physical projection layers. The network utilizes the GraphSAGE architecture to aggregate features from the spatiotemporal residuals. For each node v, its (l+1)th layer hidden state h... v (l+1) The update formula is: h v (l+1) =σ(W l ·concat(h v l ,h agg l ));

[0173] Where σ is a nonlinear activation function, such as ReLU, W l Here, h represents the learnable weight matrix, concat denotes the vector concatenation operation, and h is the weight matrix. aggl This is the aggregation vector of neighborhood node features. The aggregation operation utilizes the time-varying connectivity weight matrix W(t), with the formula: h agg l =mean(h u l |u∈N(v)), or use an attention mechanism to weight the neighboring nodes.

[0174] Furthermore, to ensure that the neural network output conforms to physical laws, this embodiment designs a composite loss function L at the network output layer. total Used for training. The loss function is defined as: L total =L mSE +λ mass ×L mass +λ neg ×L neg ;

[0175] Among them, L total Total loss; L mSE The mean square error of the data fit is the standard data; L neg For non-negative penalty terms, the formula is: L neg =∑max(0,-C pred ) 2 Used to penalize negative concentration predictions; L mass The mass conservation penalty term is used to constrain the consistency between the total mass change rate and the boundary flux within the constraint region; λ mass and λ neg The weights are hyperparameters. By minimizing this composite loss function, the trained PC-GRN network can output residual compensation values ​​or parameter corrections that satisfy physical constraints, achieving a deep integration of deep learning and physical mechanisms.

[0176] Example 11: This example describes alternative implementations of data cleaning, standardization, and parameter sensitivity analysis commonly used in engineering practice.

[0177] As an optional implementation, during the construction of the standardized multidimensional spatiotemporal dataset, to eliminate outlier noise caused by sensor malfunctions or transmission errors, the system employs the Laida criterion (3σ criterion) to clean the raw monitoring data. Specifically, for each monitoring variable sequence x, the system calculates its arithmetic mean μ and standard deviation σ. The system iterates through each data point x in the sequence. i If its absolute deviation |x iIf -μ|>3σ, the data point is considered an outlier and is either removed or replaced with an interpolated value from a nearby time period. After outlier handling, to eliminate the significant differences in dimensions and orders of magnitude between different environmental variables (such as pH, rainfall, and heavy metal concentration), the system uses Z-score standardization to make the data dimensionless. The standardized variable z... i The calculation formula is: z i =(x i -μ) / σ. This operation ensures that the subsequent neural network or inversion algorithm is not dominated by a large number of variables during training, thus accelerating convergence.

[0178] As an alternative or auxiliary scheme to parameter screening based on error statistics, this system can also use the Morris global sensitivity analysis method to initially screen the parameters of the multi-media coupling migration mechanism model. The Morris method, based on the principle of the Elementary Effect (EE), qualitatively assesses the importance of a parameter by randomly sampling within the parameter space and calculating the rate of change in output caused by small perturbations in the parameter. Specifically, for the k-th parameter θ... k Its basic effect EE k The calculation formula is: EE k =(y(θ1,...,θ k +Δ,..,θ m )-y(θ)) / Δ; where Δ is the preset perturbation step size. The system calculates EE through multiple random samplings. k The mean μ and standard deviation σ k If μ is large, it indicates that the parameter has a linear effect on the output; if σ k A large value indicates a strong interaction or nonlinear effect between this parameter and other parameters. Based on the ranking of μ*, the system selects the top-ranked parameters as the core parameter set to be inverted.

[0179] Furthermore, in the ecological risk assessment phase, in addition to dynamic risk indicators, this system also integrates the geoaccumulation index Igeo as a supplementary indicator for static assessment. The geoaccumulation index reflects the enrichment degree of heavy metals in the soil relative to background levels, and its calculation formula is: Igeo = log2(C n / (1.5*B n ));

[0180] Among them, C n The measured concentration of heavy metal n in the soil, B nis the soil geochemical background value for this area, and 1.5 is the correction factor used to account for the background value fluctuations caused by diagenesis. Based on the calculated Igeo value, the system outputs a static pollution level according to a preset grading standard (such as Igeo < 0 indicating no pollution, 0 < Igeo < 1 indicating mild pollution, etc.), which corroborates with the dynamic risk index to provide a more comprehensive environmental quality profile.

[0181] Example 12. According to one aspect of the present application, a method and system for identifying the migration mechanism of heavy metal pollution and assessing ecological risks in karst farmland are provided, including the following steps:

[0182] Step S1: Construct a multi-source heterogeneous environmental data fusion system. Through in-situ sensing, laboratory analysis, satellite / UAV remote sensing, and geological surveys, multi-dimensional spatio-temporal data such as the Cd / As speciation in the soil vertical profile, the dynamic groundwater quality, the accumulation in rice organs, meteorological precipitation, and fracture development are obtained. After spatio-temporal alignment, quality control, and feature fusion, a standardized cadmium and arsenic multi-media migration database is generated.

[0183] Step S11: Layout an in-situ sensing and sampling network. In a typical karst paddy field area, soil multi-layer profile monitoring points, groundwater observation wells, automatic surface drainage outlet water quality stations, and rice growth monitoring plots are laid out according to grids or hydrographic units; soil samples are collected at depths of 0 - 20 cm, 20 - 40 cm, and 40 - 60 cm at the soil profile monitoring points to measure total Cd, total As, and speciation classification (such as exchangeable state, carbonate-bound state, iron-manganese oxide-bound state, organic-bound state, residual state); the groundwater observation wells are equipped with pH, Eh, conductivity, dissolved oxygen, and Cd / As on-line sensors, and the data is uploaded in real time through the LoRa / NB-IoT protocol; the rice plots are sampled by organ (such as roots, stems, leaves, grains) at the tillering stage, heading stage, and maturity stage, and after microwave digestion, the Cd / As content is measured by ICP-MS;

[0184] Step S12: Obtain remote sensing and geospatial data: Use Sentinel-2 or high-resolution satellite multi-spectral images to retrieve the normalized difference vegetation index NDVI, enhanced vegetation index EVI, and red-edge band characteristics to identify crop stress areas; based on UAV-mounted hyperspectral cameras or thermal infrared sensors, obtain the vegetation physiological state and soil moisture distribution at the sub-meter field scale; integrate regional geological maps, karst fracture development maps, digital elevation models DEM, and land use status maps with a scale of 1:50,000 or higher to construct a karst hydrogeological spatial base;

[0185] Step S13: Perform spatiotemporal alignment and data standardization: Unify all observation data to the UTC+8 time base and align events according to key nodes in the rice growth cycle (such as transplanting, tillering, booting, and maturity); use Kriging or co-kriging to interpolate soil / water quality data from discrete points to a regular grid (such as 30m×30m) and match it with the remote sensing pixel center; perform correlation screening between the vegetation stress index retrieved from remote sensing and the measured heavy metal content, and remove abnormal pixels dominated by non-pollution factors such as pests, diseases, and drought;

[0186] Step S14: Implement multi-level data fusion and quality control. In data-level fusion, high-precision laboratory analysis data is used as the true value to correct the systematic bias of portable XRF or online sensors, and a dynamic offset correction model is adopted. In feature-level fusion, a ternary feature vector of soil geochemistry-hydrological response-vegetation characterization is constructed and input into an autoencoder for noise reduction and redundancy removal. Finally, a standardized multidimensional data table with spatiotemporal index is generated, with fields including but not limited to: plot ID, sampling time, soil depth, pH, Eh, Fe / Mn content, Cd. total As total Cd 可交换态 As 铁锰结合态 The data includes groundwater Cd, rice grain As, NDVI, fracture density index, and cumulative rainfall; this data table serves as a unified driving dataset for subsequent migration modeling and risk assessment.

[0187] Step S2: Establish a coupled mechanism model of cadmium and arsenic migration-transformation-enrichment, integrating karst dual-domain water flow, interfacial adsorption-desorption kinetics and root absorption process, and introduce graph neural network (GNN) to compensate for the spatial dependence of model residuals; at the same time, Morris screening and Sobol global sensitivity analysis are used to identify key uncertain parameters.

[0188] Step S21: Based on karst hydrogeochemical processes and paddy field cultivation management practices, construct an initial parameter set for the cadmium-arsenic migration-transformation-enrichment coupling mechanism model, and set a physically reasonable range for each parameter; the parameters include, but are not limited to: the soil-fracture interface preferential flow ratio coefficient α (value range [0.1, 0.6]), the maximum adsorption capacity Qmax (mg / kg) of Cd / As on iron-manganese oxides, and the first-order desorption rate constant kd (d - ¹), root absorption half-saturation constant Km (mg / L), irrigation water infiltration distribution coefficient β ([0,1]), and organic matter complexation stability constant logK for As(III); the initial values ​​of the parameters were determined based on the regional soil survey report, literature review and small-scale tracer experiment.

[0189] Step S22: Perform preliminary sensitivity analysis on the initial parameter set using the improved Morris screening method: Generate multiple random trajectories in the parameter space, and perturb a single parameter sequentially on each trajectory and run the coupled model; calculate the mean absolute initial effect μ* and effect standard deviation σ of each parameter on key output variables (such as grain Cd concentration, groundwater As flux, and rhizosphere exchangeable As proportion); Parameters that satisfy μ*<ε1 and σ<ε2 are judged as low-sensitivity parameters and removed, where ε1 and ε2 are thresholds set based on historical simulation errors;

[0190] Step S23: For the subset of highly sensitive parameters retained after Morris screening, refine the quantification using the Sobol global sensitivity analysis method: Generate a high-coverage Quasi-Monte Carlo sample set based on the Sobol sequence to ensure uniform sampling in the parameter space; run the coupled model to obtain the corresponding output response, and calculate the first-order sensitivity index Si (main effect) and the total-order sensitivity index STi (including interaction) of each parameter through variance decomposition; identify key parameters (such as the Eh–Fe content–As speciation coupling term) that have significant nonlinear or strong interaction effects on the output of compound pollution, and mark them as core parameters to be corrected;

[0191] Step S24: Introduce a graph neural network (GNN) to perform data-driven compensation for the residuals of the mechanistic model: Using field plots as graph nodes, construct an adjacency matrix based on hydrological connectivity, lithological similarity, and irrigation canal adjacency relationships; use the difference between the simulated and measured values ​​of the mechanistic model as node labels to train the GNN to learn spatial dependency residual patterns; the correction term output by the GNN is superimposed with the predicted value of the mechanistic model to form the final mixed prediction result, improving the ability to capture local abnormal migration paths (such as pollution from hidden karst caves). The learnable weight parameters in the GNN are also included in the parameter set to be optimized and participate in subsequent adaptive correction.

[0192] Step S3: Construct a migration mechanism identification framework, define typical scenarios such as fracture-preferred flow, rhizosphere enrichment, and redox release, integrate Sobol sensitivity index, SHAP feature attribution and structural equation model (SEM), quantify multi-factor causal paths, automatically label and generate mechanism explanation reports.

[0193] Step S31: Define multiple potential migration path scenarios, including: fracture preferential flow-dominated type, i.e., pollutants rapidly infiltrate down along karst conduits; rhizosphere micro-enrichment type, i.e., Cd / As locally accumulate under the influence of iron film or exudates on rice roots; redox interface release type, i.e., flooding-drying alternation causes Eh fluctuations, triggering As(III) desorption or Cd²⁺ activation; groundwater recharge and upwelling type, i.e., the rise of groundwater level during the rainy season carries deep pollutants into the cultivated layer; each scenario corresponds to a set of characteristic index combinations (such as fracture density > 0.3 fractures / m, rhizosphere Fe content > 25 g / kg, Eh amplitude > 200 mV, etc.).

[0194] Step S32: Integrate global and local sensitivity analysis results to quantify the contribution of each parameter to the output of different migration paths: Using the Sobol total order sensitivity index STi obtained in step S2, identify parameters with high nonlinear response to grain As concentration or Cd flux at 60cm depth; combine the attribution weights of each input feature (such as NDVI anomaly, high As value of neighboring plots) in the SHAP (SHapley Additive exPlanations) value parsing graph neural network GNN compensation module to the local prediction bias.

[0195] Step S33: Construct a structural equation model (SEM) to reveal multi-factor causal pathways: Rainfall intensity, irrigation regime, soil pH, Fe / Mn oxide content, and fissure development degree are set as exogenous latent variables; the proportion of exchangeable Cd, the As(III) / As(V) ratio, and root absorption flux are set as endogenous latent variables; based on measured data, fit the path coefficients to identify significant causal chains (e.g., high rainfall - enhanced fissure flow - deep As upwelling - increased rhizosphere As(III) - excessive As in grains).

[0196] Step S34: Based on the above analysis results, each field is automatically tagged with a migration mechanism and a mechanism explanation report is generated. The tags are used for dynamic weight adjustment and intelligent recommendation of remediation measures in subsequent ecological risk assessments.

[0197] Step S4: Construct a dynamic ecological risk classification and assessment model, which integrates crop exceedance probability, groundwater leaching flux, and human health risk (HQ / CR). The comprehensive risk index is calculated using the entropy weight-fuzzy comprehensive evaluation method, and the model supports adaptive adjustment of weights under extreme climate events. The model outputs a Level I spatial risk map.

[0198] Step S41: Quantify multi-dimensional risk indicators: Crop safety risk: Based on the measured concentration of Cd / As in rice grains or the model prediction, calculate the probability of exceeding the limit Pexceed=P(C>Climit), where Climit is the limit of the National Food Safety Standard for Limits of Contaminants in Food (GB2762) (Cd: 0.2mg / kg, As: 0.3mg / kg); Monte Carlo simulation is used to consider parameter uncertainty and generate a probability distribution; Groundwater environmental risk: Calculate the annual average Cd / As downward leaching flux per unit area F=∫J(t)dt, where J(t) is the solute flux in the soil layer below 60cm output by the coupled model; Human health risk: Calculate the hazard quotient HQ=EDI / RfD and the carcinogenic risk CR=EDI×SF based on the daily intake EDI, where RfD (reference dose) and SF (carcinogenic slope factor) adopt the USEPA or WHO recommended values;

[0199] Step S42: Construct a dynamic weight allocation mechanism: Introduce the risk entropy weight method to determine the objective weights of each indicator: Construct a decision matrix for m risk indicators of n plots, calculate the information entropy Ej, and then obtain the weight ωj=(1-Ej) / Σ(1-Ek); At the same time, support manual intervention to adjust the weights (such as increasing the weight of groundwater in drinking water source protection areas and increasing the weight of crop safety in major grain-producing areas); When the region encounters extreme rainfall events (such as 7-day cumulative rainfall > historical 90th percentile), automatically increase the weight of leaching flux by 20%~50% to achieve adaptive risk response under climate disturbances;

[0200] Step S43: Calculate the comprehensive ecological risk index ERI using the fuzzy comprehensive evaluation method: Define the evaluation set V={low risk, medium risk, high risk, extremely high risk} and establish a membership function (such as Gaussian or trapezoidal function) for each risk indicator to map the quantitative value to a fuzzy vector; calculate ERI=Σ(ωj×μj), where μj is the membership degree of the j-th indicator; classify the risk level according to the preset threshold: ERI∈[0,0.3)-Level I (low), [0.3,0.5)-Level II (medium), [0.5,0.7)-Level III (high), [0.7,1.0]-Level IV (extremely high).

[0201] Step S5: Establish a risk-mechanism knowledge graph database, which stores land parcel attributes, mechanism tags, risk levels, and governance suggestions in a structured manner in the graph database. Based on graph embedding, it enables similar scenario retrieval and intelligent recommendation of parameters / measures, and supports expert collaboration and continuous learning.

[0202] Step S51: Structure and encapsulate each completed cadmium and arsenic pollution assessment case to generate a knowledge graph triple record, including: Entity nodes: plot ID, watershed, soil type (red soil / calcareous soil), rice variety, year; Attribute information: dominant migration mechanism tag (e.g., redox interface release type), comprehensive risk level (Level I~IV), key sensitive parameter values ​​(e.g., Eh amplitude, Fe content, fracture density), remediation measures recommendations (e.g., application of silicon-calcium fertilizer, intermittent irrigation, planting of low-accumulation varieties); Relationship edges: define semantic relationships such as belonging to, controlled by, recommended to be adopted, similar to, etc.

[0203] Step S52: Use a graph database (such as Neo4j or Nebula Graph) to store the knowledge graph, supporting efficient graph traversal and subgraph matching; establish a unique hash fingerprint for each case, generated from input data features (such as pH mean, As morphological distribution entropy, rainfall variation coefficient), for fast deduplication and version management; build an index to accelerate multi-dimensional retrieval based on mechanism type + risk level + geographical proximity.

[0204] Step S53: Implement intelligent recommendation of similar scenarios based on graph embedding: Map land parcel nodes to low-dimensional vector representations using GraphSAGE or TransR algorithms; when a new land parcel is input, calculate the cosine similarity between its vector and the historical case library, and return the Top-K most similar cases; Automatic recommendation: ① Initial model parameters (such as priority flow ratio α, adsorption capacity Qmax); ② High-probability effective governance measures; ③ Risk indicators that need to be monitored.

[0205] Step S54: Support continuous learning and evolution of the knowledge graph: Each newly corrected or successfully verified case is automatically added to the database, triggering incremental updates to the graph; a lightweight parameter recommendation model (such as random forest or XGBoost) is trained regularly, using rainfall features and soil properties as inputs to directly predict the optimal parameter set, reducing the number of times the full model is rerun; an expert review interface is set up to allow domain experts to correct the automatic recommendation results, and the correction records are incorporated into the reinforcement learning feedback loop.

[0206] Step S6: Visualize pollution distribution, migration path heatmap, risk zoning, and historical evolution through the WebGIS platform, and automatically push early warning and remediation suggestion lists to high-risk areas.

[0207] Accordingly, the present invention also provides a system for implementing the above method, including a multi-source environmental data fusion module, a coupled modeling and parameter identification module, a migration mechanism intelligent identification module, an ecological risk dynamic assessment module, a risk-mechanism knowledge graph library module, and a visualization and intelligent decision support module. Each module is deployed on a cloud platform in a microservice architecture to support efficient collaboration and expansion.

[0208] This invention constructs an event-triggered connectivity network based on hydrogeological characteristics to address the problem of physical topological distortion in karst regions. By defining edge impedances negatively correlated with fracture density and introducing rainfall and water level thresholds as dynamic activation functions, the nonlinear hydrological response of karst regions—characterized by connectivity during the rainy season and disconnection during the dry season—is depicted. This mechanism enables the model to capture the instantaneous connectivity characteristics of long-distance fractured conduits, solving the problem of the lack of preferred flow paths in traditional models.

[0209] Furthermore, to address the issues of parameter-mechanism disconnect and parameter unidentification, this invention constructs a closed-loop correction system based on a mechanistic model and data-driven approaches. Utilizing event-triggered connected networks as spatial constraint operators, a connectivity regularization term is used in the inversion phase to force smoothing of parameters in connected regions, solving the parameter calibration problem in heterogeneous environments. In the prediction phase, by decomposing the residuals into connectivity propagation components and local chemical components, and inversely mapping the error into correction increments for physical parameters, adaptive correction of model parameters is achieved. Combined with counterfactual perturbation fingerprinting, in the absence of prior knowledge, dominant migration mechanisms such as fracture-preferred flow and redox release are distinguished and identified, achieving a leap from prediction to interpretable mechanism identification.

[0210] The above description is merely a specific embodiment of this application, but the scope of protection of this application is not limited thereto. Any variations or substitutions that can be easily conceived by those skilled in the art within the scope of the technology disclosed in this application should be included within the scope of protection of this application. Therefore, the scope of protection of this application should be determined by the scope of the claims.

Claims

1. A method for identifying the migration mechanism of heavy metal pollution in karst farmland and assessing its ecological risk, characterized in that, include: Acquire multi-source environmental monitoring data of the target area and construct a standardized multidimensional spatiotemporal dataset; A multi-media coupled migration mechanism model describing the soil-water-plant system was established, and an event-triggered connectivity network was generated based on karst hydrological characteristics. By using event-triggered connected networks, joint constraint inversion is performed on the parameter field of the multi-media coupled migration mechanism model to obtain the initial state prediction value; Calculate the prediction residual between the initial state prediction value and the measured data, decompose the prediction residual based on the event-triggered connected network and perform closed-loop parameter correction to obtain the final prediction result that satisfies the physical constraints; Controlled perturbations are applied to the corrected multi-media coupled migration mechanism model to extract counterfactual perturbation response fingerprints, and the dominant migration mechanism type of each evaluation unit is identified accordingly. Based on the final prediction results and the dominant migration mechanism type, an ecological risk classification and management strategy is generated.

2. The method according to claim 1, characterized in that, Construct a standardized multidimensional spatiotemporal dataset, including: Based on the hydrogeological characteristics in multi-source environmental monitoring data, the minimum cumulative impedance path between each monitoring point and the point to be interpolated is calculated to obtain the connectivity impedance distance. Using connectivity impedance distance as a measure of spatial correlation, co-kriging interpolation is performed on discrete sampling points in multi-source environmental monitoring data to generate a standardized multidimensional spatiotemporal dataset covering the target area.

3. The method according to claim 1, characterized in that, The multi-media coupled migration mechanism model includes mutually coupled matrix domain modules, fracture domain modules, and inter-domain mass exchange modules: The matrix domain module is used to characterize Darcy flow and solute dispersion processes in soil porous media; The fracture domain module is used to characterize non-Darcy preferential flow and solute convection processes in karst conduits and large fractures; The interdomain mass exchange module is used to characterize the dynamic mass transfer process of solute between the matrix domain module and the fracture domain module based on concentration gradient.

4. The method according to claim 3, characterized in that, The matrix domain module and the fracture domain module follow the convection-dispersion-reaction equations respectively, and are coupled through the mass exchange flux defined by the inter-domain mass exchange module: The governing equations of the matrix domain module include second-order partial derivatives describing the dispersion effect of porous media and source-sink terms describing matrix adsorption reactions; the governing equations of the fissure domain module include first-order partial derivatives describing the convection effect of preferred flow paths. The interdomain mass exchange module calculates the mass exchange flux based on a first-order kinetic model. The mass exchange flux is proportional to the solute concentration difference between the matrix domain and the fracture domain and the interdomain mass exchange coefficient.

5. The method according to claim 3, characterized in that, The multi-media coupled migration mechanism model also includes a chemical reaction module, which is used to calculate the reaction source and sink terms in each domain, specifically including: An extended Langmuir competitive adsorption model was used to calculate the dynamic equilibrium between solid-phase adsorption capacity and liquid-phase concentration based on the competitive mechanism between cadmium ions and arsenate ions on the surface of iron and manganese oxides in soil. The Fe-As coupled release model was used to establish the functional relationship between redox potential and iron oxide reduction rate, and the flux of adsorbed arsenic released during the reduction and dissolution of iron oxide was calculated accordingly. A redox-driven morphological transformation model was adopted, and the reversible conversion rate between trivalent and pentavalent arsenic was calculated based on real-time redox potential.

6. The method according to claim 1, characterized in that, Event-triggered connected networks are generated based on karst hydrological features, including: Construct the basic topology between evaluation units and define the edge impedance of the connecting edges, which is negatively correlated with the density of geological fractures; An event-driven mechanism is introduced to dynamically update the edge impedance based on the cumulative rainfall and groundwater level fluctuation characteristics in the standardized multidimensional spatiotemporal dataset. When the rainfall or water level exceeds a preset threshold, the impedance value of the relevant connected edges is reduced. The minimum cumulative impedance distance between evaluation units is calculated based on the updated edge impedance, and a time-varying connectivity weight matrix is ​​generated accordingly to construct an event-triggered connectivity network.

7. The method according to claim 6, characterized in that, Using event-triggered connected networks, joint constraint inversion is performed on the parameter field of a multi-media coupled migration mechanism model, including: Construct a joint objective function that includes a data fitting term, a connectivity regularization term, and a prior constraint term; The data fitting term is used to characterize the deviation between the model's simulated values ​​and the measured data; The connectivity regularization term uses the time-varying connectivity weight matrix to impose a penalty constraint on the parameter differences of spatially adjacent evaluation units, forcing the parameters of highly connected regions to tend to be smooth, while allowing the parameters of weakly connected regions to remain abrupt. Prior constraints are used to limit the inversion parameters to a physically reasonable range of values; By minimizing the joint objective function, the model parameters of all evaluation units are solved simultaneously to obtain the initial state prediction values.

8. The method according to claim 1, characterized in that, The prediction residuals are decomposed based on event-triggered connected networks, and closed-loop parameter correction is performed, including: The prediction residuals are decomposed into connectivity propagation components caused by spatial network linkages and local chemical components caused by local processes. Construct a parameter sensitivity matrix to inversely map the connectivity propagation component and the local chemical component into parameter correction increments in the model parameter space; The parameter field of the multi-media coupled migration mechanism model is updated using parameter correction increments, and the updated model is used to recalculate to obtain the final prediction results.

9. The method according to claim 1, characterized in that, Performing closed-loop parameter correction also includes performing physical constraint projection on the model output, specifically including: By applying non-negative concentration constraints, negative concentration values ​​that violate physical meaning in the prediction results are forcibly truncated to zero or extremely small positive numbers; Applying the regional mass conservation constraint, the total mass change of heavy metals within the evaluation area is calculated and compared with the boundary input and output fluxes. The residual deviation that violates the conservation law is evenly distributed to each evaluation unit within the area. By applying the upper limit constraint of adsorption capacity, the upper limit of the predicted value of solid phase adsorption is truncated based on the maximum adsorption capacity parameter in the multi-media coupled migration mechanism model.

10. The method according to claim 1, characterized in that, Controlled perturbations are applied to the corrected multi-media coupled migration mechanism model to extract counterfactual perturbation response fingerprints, including: Key driving factors, including cumulative rainfall, redox potential, and preferential flow ratio, were selected. While keeping other variables constant, small numerical perturbations are applied to the selected key driving factors, and the corrected multi-media coupled migration mechanism model is run to calculate the slope of the model output response to each perturbation. The response slope combination is used to construct a counterfactual perturbation response fingerprint vector, which is then matched with a pre-defined mechanism template library to determine whether the evaluation unit belongs to the fracture preferential flow-dominated, rhizosphere enrichment, or redox interface release type.

Images