A method and system for identifying ecological risk of heavy metal pollution in soil

Abstract

The application discloses a soil heavy metal pollution ecological risk identification method and system, relates to the technical field of soil pollution identification, and couples a module to construct a dynamic model of a physical field and fuse to obtain a coupling model, takes output of the coupling model as dynamic coupling data; an optimization module constructs a double-layer optimization framework of upper-layer space layout parameter screening and lower-layer strength judgment parameter optimization based on the coupling model output, adopts a probability agent model in iteration of the double-layer optimization framework, constructs a Pareto front solution set to screen candidate schemes based on multiple groups of risk identification schemes generated by iteration, and outputs a soil heavy metal pollution ecological risk dynamic identification atlas and hierarchical management and control suggestions in combination with a dynamic checking mechanism. The identification system can output the soil heavy metal pollution ecological risk dynamic identification atlas and the hierarchical management and control suggestions, not only realizes accurate positioning and hierarchical management of risks in space, but also achieves balance between ecological safety guarantee and economic cost control at a decision-making level.

Description

Technical Field

[0001] This invention relates to the field of soil pollution identification technology, specifically to a method and system for identifying ecological risks of heavy metal pollution in soil. Background Technology

[0002] With the acceleration of industrialization and agricultural intensification, soil heavy metal pollution has become an important environmental problem restricting the sustainable use of land resources and the quality and safety of agricultural products. Its risk depends not only on the total amount of heavy metals, but also on the comprehensive influence of soil physicochemical properties, environmental driving factors and ecological receptor responses, showing complex characteristics of multi-process coupling and multi-scale interaction.

[0003] Traditional risk identification methods lack simultaneous characterization of dynamic temperature-humidity stress, heavy metal speciation and spatial migration, and ecological function decline processes, making it difficult to reveal the causal chain and spatiotemporal evolution patterns between pollution, environment, and ecology.

[0004] Therefore, there is an urgent need for a soil heavy metal pollution ecological risk identification system that can integrate multi-source data, construct cross-scale parameter systems, couple multi-physics dynamic processes within a spatial explicit framework, and improve identification accuracy and decision applicability through intelligent optimization and dynamic verification, so as to meet the actual needs of precise pollution control and differentiated management. Summary of the Invention

[0005] The purpose of this invention is to provide a method and system for identifying ecological risks of heavy metal pollution in soil, so as to solve the problems in the prior art.

[0006] To achieve the above objectives, the present invention provides the following technical solution: a soil heavy metal pollution ecological risk identification system, comprising a construction module, a coupling module, and an optimization module:

[0007] Construction module: Collect multi-source data of soil region to form a cross-scale parameter system, discretize the soil region space into several homogeneous micro-elements, and construct a spatial correlation matrix based on the spatial orientation of micro-elements and the connectivity of ecological processes;

[0008] Coupling module: Based on the parameter system and spatial correlation matrix, construct a dynamic model of the physical field and fuse it to obtain a coupling model, and use the output of the coupling model as dynamic coupling data;

[0009] Optimization Module: Based on the output of the coupled model, a two-layer optimization framework is constructed, which involves screening upper-layer spatial layout parameters and optimizing lower-layer intensity determination parameters. In the iteration of the two-layer optimization framework, a probabilistic surrogate model and TPE strategy are adopted. Based on the multiple risk identification schemes generated by the iteration, a Pareto front solution set is constructed to screen candidate schemes. Combined with a dynamic verification mechanism, a dynamic identification map of soil heavy metal pollution ecological risk and hierarchical management suggestions are output.

[0010] Preferably, the optimization module constructs a two-layer optimization framework based on the output of the coupled model, which involves filtering upper-layer spatial layout parameters and optimizing lower-layer intensity determination parameters.

[0011] The upper-level model uses spatial layout parameters as decision variables to select key factor combinations from the micro-element.

[0012] Under the constraints of the upper-level screening results, the lower-level model optimizes the intensity judgment parameters to quantify the risk level.

[0013] Preferably, the optimization module employs a probabilistic surrogate model combined with a TPE strategy in the two-layer optimization framework iteration:

[0014] Randomly generate initial parameter schemes and use a surrogate model to fit the parameter-risk identification result mapping;

[0015] The TPE strategy is used to iteratively select solutions and establish feedback between upper and lower layers until the objective function converges.

[0016] Preferably, the optimization module employs a probabilistic surrogate model combined with a TPE strategy in the iterative two-layer optimization framework. Based on multiple risk identification schemes generated through iteration, it constructs a Pareto front solution set to screen candidate schemes. Combined with a dynamic verification mechanism, it outputs a dynamic identification map of soil heavy metal pollution ecological risks and hierarchical management recommendations.

[0017] The proxy model takes the parameter scheme as input and the spatial consistency index between the spatial distribution of risk level and the measured damage as output. It uses nonlinear fitting method to evaluate the performance of the parameter scheme of the coupled model that has not been actually run, and applies the TPE strategy for iterative optimization.

[0018] The TPE strategy divides existing samples into two groups: good performance and poor performance. It fits probability density functions to the two groups respectively and guides parameter sampling based on Bayesian inference principles.

[0019] In each iteration, the intensity parameters obtained from the lower-level optimization are fed back to the upper level to re-evaluate the explanatory power and directionality of the spatial layout rules;

[0020] All solutions are evaluated based on the dual objectives of control effectiveness and management cost according to the risk space clustering, and the non-dominated solution set is identified, which is the Pareto front solution set;

[0021] Several candidate schemes were selected from the Pareto frontier solution set based on decision preferences for analysis. The risk identification results of each candidate scheme in the current period were spatially overlaid and compared with the measured ecological damage data of the same period.

[0022] The deviation distribution between predicted and observed damage is calculated, and the weights of factors with significant deviations are readjusted. The output includes a dynamic identification map of ecological risks of soil heavy metal pollution and corresponding hierarchical control recommendations.

[0023] Preferably, the optimization module constructs a two-layer optimization framework based on the output of the coupled model, which involves filtering upper-layer spatial layout parameters and optimizing lower-layer intensity determination parameters.

[0024] The upper-level model constructs risk indicators based on the ecological function damage degree, bioavailability index and environmental stress field data output by the coupled model, including spatial gradient variation coefficient, risk contribution weight and proximity measure to known pollution sources;

[0025] The risk indicators are logically combined with the key environmental and social factors in the original parameter system and thresholded, and Boolean intersection and neighborhood superposition operations are used to generate candidate factor combination rules.

[0026] Calculate the proportion of explained variance of each candidate factor combination rule in the risk space distribution, and test its spatial clustering.

[0027] The lower-level model extracts the corresponding simulated values ​​of environmental stress, pollution load and ecological response from the set of high-risk micro-elements output by the upper-level model, and constructs a multi-dimensional feature space.

[0028] Candidate schemes for initializing risk thresholds and weights in a multidimensional feature space are proposed, and the risk level distribution is recalculated through a coupled model to compare its consistency with historical damage records.

[0029] Preferably, the coupling module constructs a dynamic model of the physical field based on the parameter system and spatial correlation matrix, and fuses them to obtain a coupling model, using the output of the coupling model as dynamic coupling data.

[0030] The outputs of the environment-driven field model, the pollution migration field model, and the ecological response field model are synchronized and aligned within the coupling module through a shared spatiotemporal grid index and data interface to form a dynamically coupled dataset. The dynamically coupled dataset contains independent state variables of each field and derived indicators generated by cross-field interactions.

[0031] Preferably, the coupling module constructs a dynamic model of the physical field, including an environment-driven field model, a pollution migration field model, and an ecological response field model.

[0032] Preferably, the environmental driving field model uses annual mean temperature and precipitation parameters as inputs to establish a temperature-humidity seasonal variation model, defines the nonlinear relationship between temperature and acid soluble conversion rate, humidity and leaching migration, and outputs dynamic environmental stress data;

[0033] The pollution migration field model combines porosity and connectivity parameters to adjust the pollution diffusion coefficient and water migration path weights, simulates the spatiotemporal distribution of different forms of heavy metals through the transformation potential of heavy metal forms, and outputs dynamic distribution data of pollution in the soil-crop system.

[0034] The ecological response field model uses ecological elastic strength and damage softening coefficient as state variables, coupled with environmental driving field and pollution migration field data, calculates bioavailability index and ecological function damage degree through ecotoxicological control equations, and introduces nonlinear regression relationship between heavy metal-environment interaction performance parameters and ecological damage characteristic parameters, and generates corrected elastoplastic damage constitutive parameters by combining structural parameters.

[0035] Preferably, the building blocks form a cross-scale parameter system:

[0036] A cross-scale parameter system is formed by combining heavy metal-environment interaction performance parameters, ecological receptor response characteristic parameters, and environmental-anthropogenic driving factor parameters.

[0037] The module collects multi-source data on soil regions, including total and available concentrations of heavy metals in the soil, soil physicochemical properties, ecological receptor indicators, environmental driving factors, and pollution control intensity.

[0038] This application also provides a method for identifying the ecological risk of heavy metal pollution in soil, the method comprising the following steps:

[0039] Multi-source data were collected from soil regions to form a cross-scale parameter system;

[0040] The soil region is spatially discretized into several homogeneous micro-elements, and a spatial correlation matrix is ​​constructed based on the spatial orientation of the micro-elements and the connectivity of ecological processes.

[0041] Based on the parameter system and spatial correlation matrix, dynamic models of three physical fields are constructed and fused to obtain a coupled model, and the output of the coupled model is used as dynamic coupling data.

[0042] Based on the output of the coupled model, a two-layer optimization framework is constructed, which involves filtering the upper-layer spatial layout parameters and optimizing the lower-layer intensity determination parameters. In the iteration of the two-layer optimization framework, a probabilistic surrogate model plus the TPE strategy is adopted.

[0043] Based on the multiple risk identification schemes generated through iteration, a Pareto front solution set is constructed to screen candidate schemes. Combined with a dynamic verification mechanism, a dynamic identification map of soil heavy metal pollution ecological risks and hierarchical management recommendations are output.

[0044] The technical effects and advantages provided by the present invention in the above technical solution are as follows:

[0045] 1. This invention constructs a dynamic model of three physical fields—environmental driving field, pollution migration field, and ecological response field—based on a parameter system and spatial correlation matrix. Through mechanism fusion, a unified coupled model is formed, realizing the synchronous dynamic simulation of temperature-humidity stress, heavy metal speciation and spatial migration, and ecological function decline processes. This overcomes the limitations of traditional single-field models in reflecting the interaction of multiple processes, making risk evolution prediction more mechanistic and spatiotemporally continuous.

[0046] 2. This invention introduces a two-layer framework of upper-layer spatial deployment parameter screening and lower-layer intensity judgment parameter optimization. It employs a probabilistic surrogate model combined with a TPE strategy for efficient iterative search, which can quickly approximate the global optimal solution in a high-dimensional parameter space. At the same time, it continuously corrects the spatial targeting direction and intensity judgment criteria through upper and lower layer feedback mechanisms, effectively improving the targeting and regional adaptability of risk identification. Furthermore, by constructing a Pareto front solution set to screen candidate schemes that take into account both spatial risk cluster control and management cost-effectiveness, and combining a dynamic verification mechanism to correct model parameters or weights with historical measured damage data, the reliability and timeliness of the identification results are significantly enhanced.

[0047] 3. This invention systematically collects multi-dimensional data covering soil heavy metal concentration, physicochemical properties, ecological receptor responses, environmental driving factors, and pollution control intensity, and forms a cross-scale parameter system to ensure that subsequent modeling has a comprehensive and consistent physical, chemical, and biological basis. Furthermore, the study area is spatially discretized into homogeneous micro-elements, and a spatial correlation matrix is ​​constructed by combining the spatial orientation of the micro-elements with the connectivity of ecological processes. This enables the model to characterize pollutant migration paths and ecological effect propagation mechanisms within a spatially explicit framework, significantly improving spatial resolution and the realism of process coupling. Attached Figure Description

[0048] To more clearly illustrate the technical solutions in the embodiments of this application or the prior art, the drawings used in the embodiments will be briefly introduced below. Obviously, the drawings described below are only some embodiments recorded in this invention. For those skilled in the art, other drawings can be obtained based on these drawings.

[0049] Figure 1 This is a framework diagram of the identification system of the present invention.

[0050] Figure 2 This is a flowchart of the identification method of the present invention. Detailed Implementation

[0051] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0052] Example: This example provides a soil heavy metal pollution ecological risk identification system. Please refer to [link / reference]. Figure 1 As shown, it includes:

[0053] The construction module collects multi-source data on soil regions, including total and available concentrations of heavy metals in the soil, soil physicochemical properties (pH, organic matter, clay content, porosity, aggregate stability), ecological receptor indicators (heavy metal content in edible parts of crops, abundance of microbial functional genes, survival rate of indicator organisms), environmental driving factors (annual mean temperature, precipitation, cropping system), and pollution control intensity (application of passivating agents, proportion of planting structure adjustment). It acquires heavy metal-environment interaction performance parameters (such as acid-soluble conversion rate, bioavailability coefficient, environmental buffer capacity), ecological receptor response characteristic parameters (tolerance threshold, toxicity effect kinetic parameters, functional group sensitivity weights), and environmental-anthropogenic driving factor parameters, forming a cross-scale parameter system. The soil region is spatially discretized into several homogeneous micro-elements, allowing each micro-element to independently characterize heavy metal concentration, soil structure, and ecological state. Based on the spatial orientation and ecological process connectivity of the micro-elements (such as pollutant migration paths and biological migration corridors), a spatial correlation matrix is ​​constructed to quantify the transmission coefficient of risk signals between grids. The parameter system and spatial correlation matrix are then sent to the coupling module.

[0054] Coupling Module: Based on the parameter system and spatial correlation matrix, dynamic models of three physical fields are constructed and fused to obtain a coupled model, including an environmental driving field model, a pollution migration field model, and an ecological response field model.

[0055] The environmental driving field model uses parameters such as annual average temperature and precipitation as inputs to establish a temperature-humidity seasonal variation model, defines the nonlinear relationship between temperature and acid soluble conversion rate, and humidity and leaching migration, and outputs dynamic environmental stress data.

[0056] The pollution migration field model combines porosity and connectivity parameters to adjust the pollution diffusion coefficient and water migration path weights. It simulates the spatiotemporal distribution of different forms (acid-soluble, organically bound, etc.) of heavy metals through the potential for heavy metal form transformation, and outputs dynamic distribution data of pollution in the soil-crop system.

[0057] The ecological response field model uses ecological elasticity and damage softening coefficient as state variables, coupling environmental driving field and pollution migration field data. It calculates the bioavailability index and ecological function damage degree through ecotoxicological control equations, and introduces nonlinear regression relationships between heavy metal-environment interaction performance parameters and ecological damage characteristic parameters. Combined with microstructural parameters (with weighted coefficients added when porosity is below a threshold), it generates corrected elastoplastic damage constitutive parameters, thus accurately characterizing the ecological function decline process under pollution accumulation. The output of the coupled model is used as dynamic coupling data.

[0058] Optimization Module: Based on the output of the coupled model, a two-layer optimization framework is constructed, consisting of upper-layer spatial layout parameter selection and lower-layer intensity determination parameter optimization.

[0059] The upper-level model uses spatial layout parameters as decision variables to select key factor combinations with strong explanatory power for ecological risks and clear spatial orientation from micro-elements (such as grids with pH≤5.5 and industrial source distance <1km) to ensure the spatial targeting of risk identification.

[0060] Under the constraints of the upper-level screening results, the lower-level model optimizes the intensity judgment parameters (such as the risk threshold and weight of different regions) so that the risk level quantification is both in line with the actual situation of the region and reflects the interaction effect of factors.

[0061] The two-layer optimization framework iteratively employs a probabilistic surrogate model + TPE strategy:

[0062] First, an initial parameter scheme is randomly generated. Then, a proxy model is used to quickly fit the parameter-risk identification result mapping. Next, a TPE strategy is used to iteratively screen high-potential schemes and establish feedback between upper and lower layers (the optimization results of the lower layer correct the spatial screening direction of the upper layer) until the objective function (spatial equilibrium prevention and control and cost-effectiveness adaptation) converges.

[0063] Based on the multiple risk identification schemes generated through iteration, a Pareto front solution set is constructed to screen out candidate schemes that can effectively control the spatial aggregation of risks and reduce the cost of redundant management and control. Combined with a dynamic verification mechanism (comparing historical measured damage data with predicted values ​​and correcting parameters or weights), a dynamic identification map of soil heavy metal pollution ecological risks and hierarchical management and control recommendations are output (such as implementing engineering remediation in primary control areas and strengthening agronomic regulation in secondary control areas).

[0064] This embodiment also provides a method for identifying the ecological risk of heavy metal pollution in soil. Please refer to [link / reference]. Figure 2 As shown, the identification method includes the following steps:

[0065] S1: Collect multi-source data on soil regions, including total and available concentrations of heavy metals in the soil, soil physicochemical properties (pH, organic matter, clay content, porosity, aggregate stability), ecological receptor indicators (heavy metal content in edible parts of crops, abundance of microbial functional genes, survival rate of indicator organisms), environmental driving factors (annual mean temperature, precipitation, cropping system), and pollution control intensity (application of passivating agents, proportion of planting structure adjustments). Obtain parameters of heavy metal-environment interaction performance (such as acid-soluble conversion rate, bioavailability coefficient, environmental buffer capacity), ecological receptor response characteristic parameters (tolerance threshold, toxicity effect kinetic parameters, functional group sensitivity weights), and environmental-anthropogenic driving factor parameters to form a cross-scale parameter system.

[0066] The soil region is spatially discretized into several homogeneous micro-elements, so that each micro-element can independently characterize heavy metal concentration, soil structure and ecological status; and a spatial correlation matrix is ​​constructed based on the spatial orientation of the micro-elements and the connectivity of ecological processes (such as pollutant migration paths and biological migration corridors) to quantify the transmission coefficient of risk signals between grids.

[0067] S2: Based on the parameter system and spatial correlation matrix, dynamic models of three physical fields are constructed and fused to obtain coupled models, including an environmental driving field model, a pollution migration field model, and an ecological response field model:

[0068] The environmental driving field model uses parameters such as annual average temperature and precipitation as inputs to establish a temperature-humidity seasonal variation model, defines the nonlinear relationship between temperature and acid soluble conversion rate, and humidity and leaching migration, and outputs dynamic environmental stress data.

[0069] The pollution migration field model combines porosity and connectivity parameters to adjust the pollution diffusion coefficient and water migration path weights. It simulates the spatiotemporal distribution of different forms (acid-soluble, organically bound, etc.) of heavy metals through the potential for heavy metal form transformation, and outputs dynamic distribution data of pollution in the soil-crop system.

[0070] The ecological response field model uses ecological elasticity and damage softening coefficient as state variables, coupling environmental driving field and pollution migration field data. It calculates the bioavailability index and ecological function damage degree through ecotoxicological control equations, and introduces nonlinear regression relationships between heavy metal-environment interaction performance parameters and ecological damage characteristic parameters. Combined with microstructural parameters (with weighted coefficients added when porosity is below a threshold), it generates corrected elastoplastic damage constitutive parameters, thus accurately characterizing the ecological function decline process under pollution accumulation. The output of the coupled model is used as dynamic coupling data.

[0071] S3: Based on the output of the coupled model, a two-layer optimization framework is constructed, consisting of upper-layer spatial layout parameter selection and lower-layer intensity determination parameter optimization.

[0072] The upper-level model uses spatial layout parameters as decision variables to select key factor combinations with strong explanatory power for ecological risks and clear spatial orientation from micro-elements (such as grids with pH≤5.5 and industrial source distance <1km) to ensure the spatial targeting of risk identification.

[0073] Under the constraints of the upper-level screening results, the lower-level model optimizes the intensity judgment parameters (such as the risk threshold and weight of different regions) so that the risk level quantification is both in line with the actual situation of the region and reflects the interaction effect of factors.

[0074] The two-layer optimization framework iteratively employs a probabilistic surrogate model + TPE strategy:

[0075] First, an initial parameter scheme is randomly generated. Then, a proxy model is used to quickly fit the parameter-risk identification result mapping. Next, a TPE strategy is used to iteratively screen high-potential schemes and establish feedback between upper and lower layers (the optimization results of the lower layer correct the spatial screening direction of the upper layer) until the objective function (spatial equilibrium prevention and control and cost-effectiveness adaptation) converges.

[0076] Based on the multiple risk identification schemes generated through iteration, a Pareto front solution set is constructed to screen out candidate schemes that can effectively control the spatial aggregation of risks and reduce the cost of redundant management and control. Combined with a dynamic verification mechanism (comparing historical measured damage data with predicted values ​​and correcting parameters or weights), a dynamic identification map of soil heavy metal pollution ecological risks and hierarchical management and control recommendations are output (such as implementing engineering remediation in primary control areas and strengthening agronomic regulation in secondary control areas).

[0077] In this embodiment, the functions of each module of this application are also described in detail, as follows:

[0078] The construction module collects multidimensional heterogeneous data sources covering the target soil area and forms a dataset with spatiotemporal continuity and ecological significance.

[0079] In one embodiment disclosed in this application, the data acquisition framework covers five core categories:

[0080] The total concentration of heavy metals in soil and the content of their available components were obtained through standardized sampling design and chemical stepwise extraction methods to reflect the potential migration and biological effects of heavy metals in different occurrence forms in soil. Soil physicochemical properties, including but not limited to pH, organic matter content, clay content, soil porosity, and water-stable aggregate stability indicators, provide a foundation for characterizing the physicochemical regulation capacity of the soil matrix for heavy metal fixation and release. Ecological receptor response indicators included the accumulation level of heavy metals in edible parts of crops and key microbial functional genes (such as the arsC ars reductase gene and the mer mercury resistance gene). The abundance distribution of heavy metals (such as A, B, C, and D) and the survival rates of typical indicator organisms (such as earthworms and nematodes) under exposure conditions are used to quantitatively describe the actual response status of ecosystems to pollution stress. Environmental driving factors, mainly including multi-year average temperature, annual precipitation, and agricultural farming systems (such as the frequency and pattern of continuous cropping, crop rotation, and fallow) at the regional scale, serve as background variables affecting the transformation of heavy metal forms and ecological effects. Pollution control intervention intensity indicators cover the types, amounts, and frequency of application of passivating agents (such as lime, biochar, and phosphate rock powder), as well as the proportion of crop variety substitution or structural adjustment based on policy guidance, to characterize the potential contribution of human regulation to pollution mitigation. Three types of key performance and response parameters are extracted and prepared to form a cross-scale parameter system.

[0081] The first category consists of heavy metal-environment interaction performance parameters, including acid-soluble conversion rate (used to characterize the speed at which heavy metals are converted from a stable state to a bioavailable state under acidification conditions), bioavailability coefficient (characterizing the tendency of a specific heavy metal form to cross biofilms or enter the food chain), and environmental buffer capacity (reflecting the ability of the soil system to resist external disturbances and maintain its original chemical state).

[0082] The second category consists of ecological receptor response characteristic parameters, which specifically include the tolerance threshold (i.e., critical concentration value, exceeding which will trigger significant negative effects) of different species or functional groups to specific heavy metals, toxicity effect kinetic parameters (describing the time and intensity coupling characteristics in the dose-response relationship, such as the change pattern of median inhibitory concentration with exposure duration), and functional group sensitivity weights (used to express the differences in the sensitive roles played by different ecological functional groups in the overall system stability).

[0083] The third category comprises environmental and anthropogenic driving factors, covering the seasonal fluctuations of climate factors, the temporal patterns of cropping system changes, and the coverage and intensity of control measures, thereby establishing a full-spectrum driving logic from the natural background to human intervention. All parameters must be standardized by unit, matched by dimension, and verified for credibility before being archived and coded according to spatial entities or statistical units to form a structured parameter database.

[0084] The continuous soil region is divided into a limited number of homogeneous micro-cells (grid-cells or pixels) based on the consistency of dominant attributes, so as to ensure that each micro-cell has statistical homogeneity in terms of heavy metal concentration distribution, soil physical structure and ecological function status.

[0085] The partitioning strategy can be determined comprehensively based on spatial resolution requirements (such as 10m×10m or 1km×1km) and data availability. Within each micro-element, the aforementioned multi-source data and parameter system are aggregated and stored, enabling it to independently characterize the local pollution state and ecological response, providing basic units for subsequent simulations.

[0086] Based on spatial discretization, a spatial correlation matrix that reflects the interaction relationships between micro-elements needs to be constructed to support cross-unit coupled operations of pollutant migration and ecological effect diffusion processes. The construction of this matrix comprehensively considers the spatial orientation of micro-elements (such as topological relationships such as adjacency, oblique connection, and long-distance isolation) and connectivity elements of ecological processes. The latter includes the spatial migration trajectory of pollutants along hydrological pathways (such as preferential flow zones and groundwater flow directions) or soil media (such as topsoil mixing processes), as well as the activity range of organisms along habitat corridors or seasonal migration routes.

[0087] Based on a Geographic Information System (GIS) platform, spatial indexing and encoding of micro-elements are performed, and adjacency or distance weight tables are established. For each type of ecological process, a connectivity model is defined, and minimum-cost path analysis is used to extract the lowest-resistance pathways for biological migration. For each pair of interconnected micro-elements, a risk signal transmission coefficient is calculated, with the coefficient's value determined by a comprehensive evaluation of process type, medium properties, distance attenuation patterns, and barrier effects. The final output is a weighted directed or undirected spatial correlation matrix, whose element values ​​quantify the intensity of risk or effect transmission between different micro-elements driven by a specific process.

[0088] If the study area is divided into 1km×1km micro-elements, where micro-element A and micro-element B are adjacent on the hydrological path and are connected by a priority flow zone, the two are confirmed to be hydrologically connected based on the GIS spatial index and adjacency list.

[0089] The hydrological path is determined to have strong hydraulic conductivity, with sandy loam as the medium and a distance of approximately 1 km. Referring to existing distance attenuation laws (e.g., the transmission coefficient decreases by 20% for every additional 1 km) and the barrier effect correction table (sandy loam has no significant barrier, and the reduction factor is 1), the basic transmission coefficient can be set to 1.0. After distance attenuation calculation, it is 1.0 × (1 − 0.2) = 0.8. Multiplying this by the medium and barrier reduction factor of 1, the final risk signal transmission coefficient is 0.8. If the path is changed to red soil with high clay content and low permeability interlayers, the barrier reduction factor is taken as 0.6, and the transmission coefficient is 0.8 × 0.6 = 0.48, indicating that the risk transmission intensity of the red soil path is significantly reduced at the same distance.

[0090] For biomigration connectivity, minimum cost path analysis yields a minimum resistance channel length of 2.5 km from element C to D. The base transmission coefficient of 1.0, calculated with a 20% reduction per 1 km distance, is: 1.0 × (1-0.2)^2.5 ≈ 1.0 × 0.8^2.5 ≈ 1.0 × 0.455 ≈ 0.455. If this channel crosses a composite barrier of farmland and roads, with a resistance weight reduction factor of 0.7, the final transmission coefficient is 0.455 × 0.7 ≈ 0.319, reflecting the inhibitory effect of landscape patterns and human activities on the intensity of risk diffusion in ecological processes. This can be expressed as: Transmission Coefficient = Base Coefficient × Distance Attenuation Factor × Medium and Barrier Reduction Coefficient. The distance attenuation factor is determined by the power factor or percentage reduction based on the process type and measured patterns. The medium and barrier reduction coefficient is assigned based on soil texture, interlayer distribution, and ecological resistance, thereby quantifying the transmission intensity of risk or effect signals between elements under different scenarios.

[0091] Based on the parameter system and spatial correlation matrix established in the early stage, the coupling module develops environmental driving field model, pollution migration field model and ecological response field model respectively, and realizes dynamic simulation of the three fields in synergy through mechanism coupling and data linkage.

[0092] In one embodiment of this application, the construction of the environmental driving field model uses climate factors such as regional-scale annual average temperature and precipitation as input variables, while incorporating the temporal arrangement of farming systems to establish a seasonal dynamic variation model of temperature and humidity.

[0093] Climate sequences are generated based on historical meteorological observation data and trend fitting algorithms. Periodic perturbation functions are introduced according to cropping system nodes (such as sowing, irrigation, and harvesting periods) to characterize the moderating effect of agricultural management measures on local microclimates. A nonlinear response rule is defined between temperature and the rate of acid-soluble state conversion.

[0094] A segmented incremental response logic is adopted, that is, when the temperature is in the range of low biochemical activity, the conversion rate increases slowly with the increase of temperature, while after exceeding a certain activation energy threshold, the rate increases rapidly. At the same time, the diurnal temperature variation is introduced as a modulation factor for the reaction rate.

[0095] To investigate the relationship between humidity and leaching migration, a hierarchical response logic based on a saturation threshold is constructed. When soil volumetric water content is below field capacity, leaching is weak and migration increases approximately linearly. However, when the water content exceeds field capacity and approaches saturation, migration increases rapidly due to enhanced gravity drainage. The model also considers the amplification effect of rainfall duration and intensity on short-term peak migration. The model ultimately outputs dynamic environmental stress data for each time step and each spatial micro-element, including the equivalent temperature stress index and the comprehensive humidity leaching potential field, providing driving boundary conditions for subsequent pollution form transformation and ecological response calculations.

[0096] Under the dynamic stress provided by the environmental driving field, the pollution migration field model further integrates the migration path weights in soil porosity, connectivity, and spatial correlation matrix to construct a simulation framework for the spatiotemporal distribution of multi-form heavy metals in the soil-crop system.

[0097] Spatial heterogeneity correction is performed on the classical diffusion coefficient based on the porosity distribution and the connectivity level in the spatial correlation matrix:

[0098] In micro-elements with well-developed pore networks and high connectivity, higher diffusion coefficients and water transport weights are assigned; conversely, reduction coefficients are applied to regions with low porosity or fracture blockages. For different heavy metal occurrence forms (e.g., acid-soluble, exchangeable, organically bound, and residual forms), time-step morphological transformation simulations are performed based on previously obtained morphological transformation potential parameters and environmental stress data within each micro-element. The specific processing is as follows:

[0099] Within each time step, based on the current temperature stress index and humidity leaching potential of the microelement, the corresponding acid-soluble conversion rate and organic-bound stability coefficient are used to calculate the net conversion flux between different forms and update the form concentration distribution. Simultaneously, combined with water migration path weights, dissolved or suspended heavy metals are spatially pushed along preferred flow paths, considering path hindrance effects (such as clay adsorption and organic matter chelation) and dilution effects (such as infiltration recharge). The model's final output includes a spatiotemporal raster sequence of heavy metal concentrations in each form and the cumulative flux in the crop root layer and edible parts, providing pollution load input for the ecological response field.

[0100] The ecological response field model aims to couple information from both environmental drivers and pollution migration fields to quantitatively characterize the functional decline trajectory of ecosystems under continuous pollution stress. The model uses ecological resilience and damage softening coefficient as core state variables, with initial values ​​set based on regional ecological baseline surveys and literature benchmarks, and dynamically updated during the simulation process.

[0101] The equivalent temperature stress index and humidity leaching potential from the environmental driving field, as well as the spatiotemporal distribution of bioavailable heavy metal concentrations output from the pollution migration field, are received. The bioavailability index within each micro-element is then calculated using ecotoxicological control equations.

[0102] The calculation process is based on a dose-response relationship library of heavy metal types and ecological receptor types. It uses table lookup and interpolation logic to determine the effectiveness correction factor at different concentrations and introduces a hysteresis decay kernel function for long-term cumulative effects.

[0103] Combining the ecological function damage assessment logic, the bioavailability index is compared with the receptor tolerance threshold. If the threshold is exceeded, the damage increment is calculated based on the toxicity effect kinetic parameters, and each damage component is weighted and summed using functional group sensitivity weights to obtain the comprehensive ecological function damage degree.

[0104] To further accurately characterize the nonlinear degradation process under pollution accumulation, the model introduces a multivariate nonlinear regression relationship between heavy metal-environment interaction performance parameters and ecological damage characteristic parameters. This regression relationship is obtained through training on historical case data, and is adjusted for prediction based on variables such as the current proportion of acid-soluble substances in micro-elemental organisms, organic matter content, and abundance of microbial functional genes. Furthermore, the model incorporates a feedback mechanism for microstructural parameters.

[0105] When the porosity of a certain micro-element is detected to be lower than a preset threshold (such as when clay accumulation leads to deterioration of aeration), the weighting coefficient of ecological damage at that location is automatically increased to reflect the inhibitory effect of physical structure degradation on ecological resilience.

[0106] Based on the above calculations, the corrected elastoplastic damage constitutive parameters are generated to describe the irreversible decline characteristics of ecological functions under repeated stress, and output a dynamic dataset containing the spatiotemporal evolution of ecological elasticity, the distribution of functional damage degree, and the damage constitutive parameters.

[0107] Suppose that the annual average temperature of a certain micro-element at a certain time step is 28℃ and the daily average diurnal temperature range is 8℃. According to the piecewise increasing response rule of the environmental driving field model, the low biochemical activity range is assumed to be 0~20℃ (the conversion rate increases slowly with temperature, increasing by 0.01 for every 1℃ increase). After exceeding the activation energy threshold of 20℃, it enters the acceleration range (the rate increases by 0.03 for every 1℃ increase). Then, the basic conversion rate of the micro-element = 0.01×20 + 0.03×(28-20) = 0.2 + 0.24 = 0.44. Introducing the diurnal temperature range fluctuation modulation factor (assuming that the modulation factor increases by 0.005 for every 1℃ increase in temperature range, with a base value of 1.0), the modulation factor = 1.0 + 0.005×8 = 1.04, and the acid-soluble conversion rate = 0.44×1.04≈0.458.

[0108] If the soil volumetric water content of this micro-element is 32% and the field capacity is 30%, according to the saturation threshold graded response logic of the pollution migration field model, the water content enters a rapid increase phase when it surpasses the field capacity. Let the baseline leaching migration rate be 5 kg·ha. -1 ·d -1 For every 1% increase in water content exceeding the threshold, the migration rate increases by 0.8 kg·ha. -1 ·d -1 Therefore, the migration amount = 5 + 0.8 × (32 - 30) = 6.6 kg·ha -1 ·d -1 Considering the heavy rainfall lasted for 2 hours and was of high intensity, a short-term amplification factor of 1.5 was introduced. The final leaching migration amount was 6.6 × 1.5 = 9.9 kg·ha. -1 ·d -1 .

[0109] In the pollution migration field, the micro-element has a porosity of 0.45 (higher than the threshold of 0.35), high connectivity, a diffusion coefficient correction value of 1.2, and a water transport weight of 0.9. Combining the acid-soluble conversion rate and the organic-bound stability coefficient (assuming a stability coefficient of 0.85), the net conversion flux increases the acid-soluble concentration by 0.458 × Δt. Simultaneously, it is pushed along the preferred flow path. The adsorption retardation coefficient of clay particles is 0.7, and the chelation retardation coefficient of organic matter is 0.8. The overall retardation is 0.7 × 0.8 = 0.56. The dilution coefficient of the pushing loss is 0.9. Therefore, the effective pushing concentration is equal to the original concentration × 0.56 × 0.9.

[0110] In the ecological response field, assuming the concentration of bioavailable heavy metals is 120 mg·kg⁻¹ -1 Receptor tolerance threshold 100 mg / kg -1 From the table, the efficacy correction factor is 1.1, the hysteresis decay kernel function (with a cumulative weight of 0.9 for 3 time steps) yields a bioefficacy index of 120 × 1.1 × 0.9 = 118.8; the excess dose is 18.8 mg / kg.-1 Based on the toxicity effect kinetic parameters (damage increment of 0.05 per unit increment), the damage increment is calculated as 18.8 × 0.05 = 0.94. The sensitivity weights for functional groups are 0.5 for crops, 0.3 for microorganisms, and 0.2 for indicator organisms. Therefore, the comprehensive ecological function damage degree is 0.94 × (0.5 + 0.3 + 0.2) = 0.94. Since the porosity of 0.45 is greater than the threshold of 0.35, the weight increase is not triggered.

[0111] By correcting through multivariate nonlinear regression (assuming a correction coefficient of 1.02 in the regression output), the corrected damage degree is obtained as 0.94 × 1.02 ≈ 0.959. The ecological elastic strength reduction value and elastoplastic damage constitutive parameters of the micro-element at this time are generated to describe the irreversible functional degradation.

[0112] The outputs of the environment-driven field model, the pollution migration field model, and the ecological response field model are synchronized and aligned within the coupling module through a shared spatiotemporal grid index and data interface, forming a unified, dynamically coupled data set. This set not only includes the independent state variables of each field but also covers derived indicators generated by cross-field interactions (such as environmental stress-weighted pollution bioavailability and ecological vulnerability index under structural constraints). It can directly drive higher-level decision analysis modules, such as risk zoning, control priority ranking, and adaptive management scenario simulation.

[0113] If the shared spatiotemporal grid index of a certain micro-element in the study area is G(35,78) at a certain time step, the environmental driving field outputs an equivalent temperature stress index of 1.25 (baseline value 1.0, indicating that temperature conditions promote the transformation of acid-soluble states), and the comprehensive humidity leaching potential field value is 0.92 (baseline value 1.0, slightly lower than the strong leaching level); the pollution migration field outputs a bioavailable heavy metal concentration of 110 mg·kg⁻¹ for this micro-element. -1 The proportion of acid-soluble state in the morphological distribution is 0.40; the ecological response field output at this location is 0.82 (1.0 is the undamaged state), and the porosity is 0.38 (below the structural threshold of 0.40, triggering structural constraints).

[0114] The coupling module first performs a weighted synthesis of environmental stress and pollution bioavailability, deriving an environmental stress-weighted pollution bioavailability index: Weighted bioavailability = bioavailable concentration × acid-soluble proportion × equivalent temperature stress index × humidity leaching potential field value. Substituting the values, we get 110 × 0.40 × 1.25 × 0.92 = 110 × 0.40 = 44, 44 × 1.25 = 55, 55 × 0.92 = 50.6 mg·kg. -1This indicates that the bioavailability of this micro-element under environmental stress is significantly higher than that of simple concentration measurements. Next, combining ecological resilience and structural constraints, an ecological vulnerability index under structural constraints was generated, calculated as: Vulnerability Index = (1 - Ecological Resilience) × Structural Constraint Weight. The structural constraint weight is 1.2 when the porosity is below a threshold, and 1.0 otherwise. Here, the porosity is 0.38 < 0.40, so the weight is 1.2. The calculated vulnerability index is (1 - 0.82) × 1.2 = 0.18 × 1.2 = 0.216, reflecting that this location is more prone to functional decline under conditions of impaired ecological resilience and unfavorable physical structure.

[0115] The two derived indicators mentioned above, along with the original state variables of each field, are synchronously aligned under a unified grid index to form a dynamically coupled dataset. For example, risk zoning can be based on weighted bioavailability >45 mg·kg. -1 Furthermore, a vulnerability index > 0.20 is classified as a Level 1 risk zone; the control priority ranking can prioritize micro-elements that simultaneously meet both conditions; adaptive management scenario simulation can test the response changes of the vulnerability index under scenarios of reduced temperature stress index or increased porosity, thereby assessing the effectiveness of engineering or agronomic measures.

[0116] Based on the dynamic risk information output by the coupled model, the optimization module establishes a two-layer optimization framework for precise spatial identification and scientific intensity determination, thereby achieving targeted identification and intelligent generation of hierarchical management strategies for the ecological risks of soil heavy metal pollution. This framework adopts a hierarchical structure of upper-layer spatial layout parameter screening and lower-layer intensity determination parameter optimization. Furthermore, it introduces a probabilistic surrogate model combined with a tree-structured-Parzen-Estimator (TPE) strategy in the iterative solution process to ensure efficient search for feasible solutions in a high-dimensional and complex parameter space.

[0117] In one embodiment disclosed in this application, the upper-level model uses spatial layout parameters as decision variables to focus on screening key factor combinations with strong explanatory power for ecological risks and clear spatial orientation from the global micro-element, thereby delineating risk target areas with priority attention value.

[0118] Based on the ecological function damage degree, bioavailability index and environmental stress field data output by the coupled model, risk indicators such as spatial gradient variation coefficient, risk contribution weight and proximity measure to known pollution sources are constructed.

[0119] These indicators are logically combined and thresholded with key environmental and social factors in the original parameter system (such as pH value, distance from industrial sources, land use type, and intensity of farming system), and Boolean intersection and neighborhood superposition operations are used to generate candidate factor combination rules.

[0120] For example, a rule can be set that the pH value is less than or equal to 5.5 and the straight-line distance from the industrial source is less than 1 kilometer, and micro-elemental cells that meet this rule are marked as high-risk potential grids. To evaluate the explanatory power of different factor combinations, a screening criterion based on information entropy and spatial autocorrelation is introduced:

[0121] The proportion of the variance explained by each candidate rule in the risk spatial distribution is calculated, and its spatial clustering significance (such as Moran's I index) is tested. Rules that exhibit both high explanatory power and significant spatial clustering are retained. The output of the upper-level model is a set of optimized spatial layout parameters and their corresponding set of high-risk micro-elements, providing a spatial constraint domain for the lower-level optimization.

[0122] Within the spatial range defined by the upper-level screening results, the lower-level model optimizes the intensity determination parameters to ensure that the quantitative results of the risk level not only match the regional ecological background and actual stress level, but also fully reflect the marginal effect of multi-factor interaction on risk amplification. The intensity determination parameters mainly include the thresholds for different risk level classifications (such as the specific numerical boundaries for classifying ecological function damage into three levels: mild, moderate, and severe) and the weight coefficients of each factor in the comprehensive risk index.

[0123] Based on the set of high-risk micro-elements output from the upper layer, the corresponding simulated values ​​of environmental stress, pollution load, and ecological response are extracted to construct a multi-dimensional feature space. Within this space, a set of candidate schemes with risk thresholds and weights are initialized, and the risk level distribution is recalculated through a coupled model, comparing its consistency with historical damage records. To improve search efficiency, a probabilistic surrogate model is used to approximate the mapping relationship between intensity determination parameters and risk identification results.

[0124] The surrogate model takes a parameter scheme as input and outputs spatial consistency indicators between the spatial distribution of risk levels and measured damage (such as the Kappa coefficient, the harmonic mean of hit rate and false alarm rate). It uses nonlinear fitting methods such as Gaussian processes or random forests to quickly evaluate the performance of parameter schemes that have not been actually run with the coupled model. Iterative optimization is performed using a TPE strategy.

[0125] TPE divides existing samples into two groups, one with better performance and one with worse performance. It then fits probability density functions to each group and uses Bayesian inference principles to guide the next parameter sampling step, prioritizing the exploration of potentially high-yield regions. In each iteration, the optimized intensity parameters obtained from the lower layer are fed back to the upper layer to reassess the explanatory power and directionality of the spatial deployment rules. If necessary, factor combinations are adjusted or supplemented, thus forming a two-way correction mechanism between the upper and lower layers. This drives the overall scheme towards the objective function of spatially balanced prevention and control and cost-effectiveness.

[0126] If the entire area is divided into 1km² micro-elements, the upper-level model first calculates the risk significance assessment index, such as the ecological function damage degree of a certain micro-element being 0.75 and the bioavailability index being 52 mg·kg. -1 The equivalent temperature stress index is 1.25. The spatial gradient coefficient of variation is the standard deviation of the damage degree of the micro-element and the average damage degree of the neighboring 8 grids divided by the mean. If the neighboring average is 0.60, then the standard deviation is 0.10, and the coefficient of variation = 0.10 ÷ 0.60 ≈ 0.167. The risk contribution weight is determined according to the proportion of the damage degree of the micro-element to the total damage of the whole area. Assuming the total damage of the whole area is 750, and the contribution of this grid is 0.75, then the weight = 0.75 ÷ 750 = 0.001. The proximity measurement with known industrial sources adopts inverse distance weighting. If the straight-line distance from the industrial source is 0.8km, the function is set as f(d) = max(0, 1-d / 1.0), and f(0.8) = 0.2 is obtained.

[0127] Combining the above indicators with the original parameter system's pH=5.2, distance from industrial sources 0.8km, and land use type of industrial and mining perimeter farmland, and applying the Boolean intersection rule (pH≤5.5 and distance from industrial sources<1km), the grid is marked as a high-risk potential grid. To assess explanatory power, the proportion of variance explained by the risk spatial distribution of the covered micro-element is calculated:

[0128] Assuming the rule's coverage area has an average damage score of 0.72 and the global average is 0.50, with variances of 0.015 and 0.020 respectively, the explained variance percentage is (0.020 − 0.015) ÷ 0.020 = 0.25, meaning the rule can explain 25% of the spatial variation in damage scores. The Moran's I index test yields 0.62 (p < 0.01), indicating significant spatial clustering, and therefore the rule is retained. The upper layer outputs this rule and the corresponding set of high-risk micro-elements for use by the lower layer.

[0129] The lower-level model optimizes the intensity judgment parameters within this constraint domain, setting initial ecological function damage thresholds: mild <0.4, moderate 0.4–0.7, and severe ≥0.7. The initial weighting coefficients are damage degree 0.5, bioavailability index 0.3, and environmental stress index 0.2. A cell within the high-risk set is selected with damage degree 0.75, bioavailability index 52, and environmental stress index 1.25. The comprehensive risk index is calculated as: 0.5 × normalized damage degree + 0.3 × normalized bioavailability + 0.2 × normalized stress index (assuming normalized values ​​are 0.75, 0.52, and 0.625 respectively). The comprehensive risk is then calculated as: 0.5 × 0.75 + 0.3 × 0.52 + 0.2 × 0.625 = 0.375 + 0.156 + 0.125 = 0.656. Based on the initial threshold, this cell is classified as medium risk, but historical measurements indicate severe damage, resulting in low consistency.

[0130] A probabilistic surrogate model is introduced, taking parameter schemes (thresholds and weights) as input and outputting consistency indicators such as the Kappa coefficient. Suppose that one scheme is evaluated by the surrogate model with a Kappa of 0.58, while another scheme has a Kappa of 0.72; the latter is clearly superior. The TPE strategy divides the samples into high-Kappa and low-Kappa groups, fitting probability densities to each group. Based on Bayesian inference, in the next iteration, more parameters in the high-Kappa range are sampled. For example, the damage threshold is adjusted to mild <0.35, moderate 0.35–0.65, and severe ≥0.65, with the weights changed to damage 0.6, bioavailability 0.25, and stress 0.15. The overall risk of this grid is recalculated as 0.6 × 0.75 + 0.25 × 0.52 + 0.15 × 0.625 = 0.45 + 0.13 + 0.094 = 0.674, falling into the severe range, thus improving the consistency.

[0131] This superior strength parameter is fed back to the upper layer to recalculate the rule's explained variance and Moran's I. If the explained variance increases to 35% and the I value increases to 0.68, the rule's explanatory power and directionality are enhanced. This bidirectional correction between the upper and lower layers drives the convergence of the objective function (spatial equilibrium prevention and control and cost-effectiveness adaptation), ultimately obtaining a feasible solution that balances risk identification accuracy and control costs.

[0132] In one embodiment disclosed in this application, during the multi-round iteration of the two-level optimization framework, several risk identification schemes with different performance in the objective function space are generated. To select the optimal solution that balances risk control effectiveness and economic rationality, a Pareto front solution set needs to be constructed:

[0133] All schemes are evaluated under two main optimization objectives: the control effectiveness of risk spatial agglomeration (such as the reduction rate of high-risk area area and the extent of connectivity reduction) and the economy of control costs (such as the spatial coverage and unit cost of required repair or control measures). The non-dominated solution set that does not have one side better without weakening the other is the Pareto front.

[0134] Based on decision-making preferences (e.g., favoring schemes with high clustering control efficiency if ecological security is prioritized, and favoring low-cost schemes if economic efficiency is prioritized), several candidate schemes are selected from the frontier solution set for detailed analysis. To ensure the reliability and timeliness of the schemes, a dynamic verification mechanism is introduced:

[0135] The risk identification results of each candidate scheme in the current period are spatially overlaid and compared with the measured ecological damage data of the same period. The deviation distribution between predicted damage and observed damage is calculated, and the weights of factors with significant deviations are readjusted, triggering local re-optimization when necessary. The final output includes a dynamic identification map of soil heavy metal pollution ecological risk and corresponding hierarchical management recommendations. The risk level is marked by different color levels in the map. In the first-level control area (usually a high-risk cluster with extremely low ecological resilience), engineering remediation (such as deep plowing and replacement, chemical passivation combined with phytoremediation) is recommended. In the second-level control area (medium risk with some self-recovery capacity), enhanced agronomic regulation (such as optimizing fertilization structure, rotating nitrogen-fixing plants, and adjusting irrigation system to inhibit leaching) is recommended. In the third-level area, monitoring priority and preventive management are the main focus, thereby achieving differentiated and refined soil pollution risk management based on scientific simulation and multi-objective trade-offs.

[0136] If the iteration generates 5 risk identification schemes, each scheme is evaluated according to two optimization objectives. The first objective is the control effectiveness of risk spatial aggregation, which is calculated by weighted average of the area reduction rate and the decrease in connectivity of high-risk areas. The formula is: Control effectiveness = 0.6 × Area reduction rate + 0.4 × Decrease in connectivity (both area reduction rate and decrease in connectivity are expressed as decimals, such as 0.30 for a 30% reduction). The second objective is the economy of control costs, which is expressed as the reciprocal of the product of the spatial coverage of the required repair or control measures and the unit cost. The formula is: Economy score = 1 / (Coverage × Unit cost). Coverage and unit cost are also decimals or normalized values.

[0137] Option A has a high-risk area reduction rate of 0.30 and a connectivity reduction of 0.20, resulting in a control effectiveness of 0.6 × 0.30 + 0.4 × 0.20 = 0.18 + 0.08 = 0.26. With a coverage of 0.40 and a unit cost of 0.80, the economic score is 1 / (0.40 × 0.80) = 1 / 0.32 ≈ 3.13. Option B has a control effectiveness of 0.35 and an economic score of 2.50; Option C has a control effectiveness of 0.28 and an economic score of 3.80; Option D has a control effectiveness of 0.36 and an economic score of 2.40; and Option E has a control effectiveness of 0.25 and an economic score of 4.00.

[0138] Perform bi-objective non-dominated screening:

[0139] Option D has the highest control effectiveness (0.36), but its economic efficiency is lower than that of Option B; Option E has the best economic efficiency (4.00), but its control effectiveness is the worst. There is no option that is superior to the other in both objectives. Therefore, among A, B, C, D, and E, B, C, and D constitute the Pareto front (assuming A is dominated by B and E by C). If the decision preference leans towards ecological security, then D (0.36) with the highest control effectiveness and B (0.35) with the second highest are selected as candidates; if the preference leans towards economic efficiency, then C (3.80) and E (4.00) with the best economic efficiency are selected. Subsequent dynamic verification follows.

[0140] The risk identification results of the candidate schemes are spatially overlaid with the measured ecological damage data of the same period, and the mean square error between the predicted damage and the observed damage is calculated. For example, Scheme D predicts a damage degree of 0.70 and the measured damage degree of 0.55 in a certain micro-element, with an error of 0.15. Scheme B has an error of 0.08, and Scheme C has an error of 0.12. If the significant deviation threshold is set to 0.10, then Scheme D has a significant deviation and its factor weights or thresholds need to be readjusted. For example, the weight of ecological function damage degree can be reduced from 0.6 to 0.55, and the weight of bioavailability can be increased from 0.25 to 0.30, and the local optimization can be run again. Output a dynamic identification map of soil heavy metal pollution ecological risk, using color levels to distinguish between Level 1 (high risk clustering and ecological resilience <0.4, recommended engineering remediation such as deep tillage replacement + chemical passivation + phytoremediation), Level 2 (medium risk and resilience 0.4–0.7, recommended agronomic regulation such as optimized fertilization + crop rotation of nitrogen-fixing plants + controlled leaching irrigation), and Level 3 (low risk and resilience >0.7, monitoring priority and preventive management) control zones, to achieve differentiated and refined risk management based on scientific simulation and multi-objective trade-offs.

[0141] In the description of this specification, references to terms such as "an embodiment," "example," "specific example," etc., indicate that a specific feature, structure, material, or characteristic described in connection with that embodiment or example is included in at least one embodiment or example of the invention. In this specification, illustrative expressions of the above terms do not necessarily refer to the same embodiment or example. Furthermore, the specific features, structures, materials, or characteristics described may be combined in any suitable manner in one or more embodiments or examples.

[0142] The preferred embodiments of the present invention disclosed above are merely illustrative of the invention. These preferred embodiments do not exhaustively describe all details, nor do they limit the invention to any specific implementation. Clearly, many modifications and variations can be made based on the content of this specification. This specification selects and specifically describes these embodiments to better explain the principles and practical applications of the invention, thereby enabling those skilled in the art to better understand and utilize the invention. The invention is limited only by the claims and their full scope and equivalents.

Claims

1. A soil heavy metal pollution ecological risk identification system, characterized in that: This includes building modules, coupling modules, and optimization modules: Construction module: Collect multi-source data of soil region to form a cross-scale parameter system, discretize the soil region space into several homogeneous micro-elements, and construct a spatial correlation matrix based on the spatial orientation of micro-elements and the connectivity of ecological processes; Coupling module: Based on the parameter system and spatial correlation matrix, construct a dynamic model of the physical field and fuse it to obtain a coupling model, and use the output of the coupling model as dynamic coupling data; The coupling module constructs a dynamic model of the physical field, including an environmental driving field model, a pollution migration field model, and an ecological response field model. The environmental driving field model takes the annual average temperature and precipitation parameters as inputs, establishes a temperature-humidity seasonal variation model, defines the nonlinear relationship between temperature and acid soluble conversion rate, humidity and leaching migration, and outputs dynamic environmental stress data. The pollution migration field model combines porosity and connectivity parameters to adjust the pollution diffusion coefficient and water migration path weights, simulates the spatiotemporal distribution of different forms of heavy metals through heavy metal form transformation potential, and outputs dynamic distribution data of pollution in the soil-crop system. The ecological response field model uses ecological elastic strength and damage softening coefficient as state variables, couples environmental driving field and pollution migration field data, calculates bioavailability index and ecological function damage degree through ecotoxicological control equation, and introduces nonlinear regression relationship between heavy metal-environment interaction performance parameters and ecological damage characteristic parameters, and generates corrected elastoplastic damage constitutive parameters by combining structural parameters. Optimization Module: Based on the output of the coupled model, a two-layer optimization framework is constructed, which involves screening upper-layer spatial layout parameters and optimizing lower-layer intensity determination parameters. In the iteration of the two-layer optimization framework, a probabilistic surrogate model and TPE strategy are adopted. Based on the multiple risk identification schemes generated by the iteration, a Pareto front solution set is constructed to screen candidate schemes. Combined with a dynamic verification mechanism, a dynamic identification map of soil heavy metal pollution ecological risk and hierarchical management suggestions are output.

2. The soil heavy metal pollution ecological risk identification system according to claim 1, characterized in that: The optimization module, based on the output of the coupled model, constructs a two-layer optimization framework: upper-layer spatial layout parameter selection and lower-layer intensity determination parameter optimization. The upper-level model uses spatial layout parameters as decision variables to select key factor combinations from the micro-element. Under the constraints of the upper-level screening results, the lower-level model optimizes the intensity judgment parameters to quantify the risk level.

3. The soil heavy metal pollution ecological risk identification system according to claim 2, characterized in that: The optimization module employs a probabilistic surrogate model plus a TPE strategy in the two-layer optimization framework iteration: Randomly generate initial parameter schemes and use a surrogate model to fit the parameter-risk identification result mapping; The TPE strategy is used to iteratively select solutions and establish feedback between upper and lower layers until the objective function converges.

4. The soil heavy metal pollution ecological risk identification system according to claim 3, characterized in that: The optimization module employs a probabilistic surrogate model and a TPE strategy in the iterative two-layer optimization framework. Based on multiple risk identification schemes generated through iteration, it constructs a Pareto front solution set to screen candidate schemes. Combined with a dynamic verification mechanism, it outputs a dynamic identification map of soil heavy metal pollution ecological risks and hierarchical management recommendations. The proxy model takes the parameter scheme as input and the spatial consistency index between the spatial distribution of risk level and the measured damage as output. It uses nonlinear fitting method to evaluate the performance of the parameter scheme of the coupled model that has not been actually run, and applies the TPE strategy for iterative optimization. The TPE strategy divides existing samples into two groups: good performance and poor performance. It fits probability density functions to the two groups respectively and guides parameter sampling based on Bayesian inference principles. In each iteration, the intensity parameters obtained from the lower-level optimization are fed back to the upper level to re-evaluate the explanatory power and directionality of the spatial layout rules; All solutions are evaluated based on the dual objectives of control effectiveness and management cost according to the risk space clustering, and the non-dominated solution set is identified, which is the Pareto front solution set; Several candidate schemes were selected from the Pareto frontier solution set based on decision preferences for analysis. The risk identification results of each candidate scheme in the current period were spatially overlaid and compared with the measured ecological damage data of the same period. The deviation distribution between predicted and observed damage is calculated, and the weights of factors with significant deviations are readjusted. The output includes a dynamic identification map of ecological risks of soil heavy metal pollution and corresponding hierarchical control recommendations.

5. The soil heavy metal pollution ecological risk identification system according to claim 3, characterized in that: The optimization module, based on the output of the coupled model, constructs a two-layer optimization framework: upper-layer spatial layout parameter selection and lower-layer intensity determination parameter optimization. The upper-level model constructs risk indicators based on the ecological function damage degree, bioavailability index and environmental stress field data output by the coupled model, including spatial gradient variation coefficient, risk contribution weight and proximity measure to known pollution sources; The risk indicators are logically combined with the key environmental and social factors in the original parameter system and thresholded. Boolean intersection and neighborhood superposition operations are used to generate candidate factor combination rules. Calculate the proportion of explained variance of each candidate factor combination rule in the risk space distribution, and test its spatial clustering. The lower-level model extracts the corresponding simulated values ​​of environmental stress, pollution load and ecological response from the set of high-risk micro-elements output by the upper-level model, and constructs a multi-dimensional feature space. Candidate schemes for initializing risk thresholds and weights in a multidimensional feature space are proposed, and the risk level distribution is recalculated through a coupled model to compare its consistency with historical damage records.

6. The soil heavy metal pollution ecological risk identification system according to claim 2, characterized in that: The coupling module constructs a dynamic model of the physical field based on the parameter system and spatial correlation matrix, and fuses them to obtain a coupling model. The output of the coupling model is used as dynamic coupling data. The outputs of the environment-driven field model, the pollution migration field model, and the ecological response field model are synchronized and aligned within the coupling module through a shared spatiotemporal grid index and data interface to form a dynamically coupled dataset. The dynamically coupled dataset contains independent state variables of each field and derived indicators generated by cross-field interactions.

7. The soil heavy metal pollution ecological risk identification system according to claim 6, characterized in that: The building blocks form a cross-scale parameter system: A cross-scale parameter system is formed by combining heavy metal-environment interaction performance parameters, ecological receptor response characteristic parameters, and environmental-anthropogenic driving factor parameters. The construction module collects multi-source data on soil regions, including total and available concentrations of heavy metals in the soil, soil physicochemical properties, ecological receptor indicators, environmental driving factors, and pollution control intensity.

8. A method for identifying ecological risks of heavy metal pollution in soil, implemented using the identification system described in any one of claims 1-7, characterized in that: The identification method includes the following steps: Multi-source data were collected from soil regions to form a cross-scale parameter system; The soil region is spatially discretized into several homogeneous micro-elements, and a spatial correlation matrix is ​​constructed based on the spatial orientation of the micro-elements and the connectivity of ecological processes. Based on the parameter system and spatial correlation matrix, dynamic models of three physical fields are constructed and fused to obtain a coupled model, and the output of the coupled model is used as dynamic coupling data. Based on the output of the coupled model, a two-layer optimization framework is constructed, which involves filtering the upper-layer spatial layout parameters and optimizing the lower-layer intensity determination parameters. In the iteration of the two-layer optimization framework, a probabilistic surrogate model plus the TPE strategy is adopted. Based on the multiple risk identification schemes generated through iteration, a Pareto front solution set is constructed to screen candidate schemes. Combined with a dynamic verification mechanism, a dynamic identification map of soil heavy metal pollution ecological risks and hierarchical management recommendations are output.

Images