A graded remediation method using modified biochar for synergistic passivation in cadmium and arsenic contaminated paddy fields.

Abstract

This invention discloses a graded remediation method for cadmium and arsenic contaminated paddy fields using modified biochar for synergistic passivation. The method includes: collecting multi-source heterogeneous environmental data from the paddy field area; constructing coupled feature vectors and inputting them into a pre-constructed risk assessment model; calculating the composite risk probability of grid cells; and dividing the paddy field area into risk control zones of different levels based on this probability. The method involves: differentially applying functionally matched modified biochar materials, including composite-loaded biochar for high-risk areas; real-time monitoring of the rhizosphere microenvironment during the critical growth stages of rice; constructing an adaptive target control feasible domain boundary for each plot based on the coupled feature vectors; and implementing a multivariate synergistic control strategy to maintain the rhizosphere microenvironment state within a safe domain that inhibits cadmium and arsenic activation. This invention effectively solves the remediation failure problems caused by regional Eh-pH decoupling and fissure water disturbance, achieving simultaneous and long-term remediation of cadmium and arsenic.

Description

Technical Field

[0001] This invention belongs to the field of soil environmental remediation, and in particular to a method for graded remediation of cadmium and arsenic contaminated paddy fields by synergistic passivation with modified biochar. Background Technology

[0002] Currently, remediation technologies for paddy fields contaminated with cadmium and arsenic mainly include in-situ passivation and water management. In-situ passivation often uses materials such as lime, biochar, or iron-manganese oxides to reduce the activity of heavy metals through adsorption or precipitation. Water management widely employs wet-dry alternation (AWD) or intermittent irrigation patterns, attempting to balance cadmium fixation and arsenic adsorption by adjusting the soil redox potential (Eh). Existing research and engineering practices typically establish fixed water management systems or universal Eh control thresholds, combined with quantitative application of passivating agents, in order to achieve the goal of passivation remediation.

[0003] However, existing technologies suffer from deep-seated technical problems in specific geological habitats, such as regulatory failure due to the decoupling of the Eh-pH relationship and control lag caused by disturbances from karst fissure water. Specifically, traditional methods do not adequately consider the pH buffering effect of high concentrations of bicarbonate in the soil, leading to the breakdown of the classic coupling relationship between Eh and pH. Regulation based solely on the Eh threshold is insufficient in taking into account the stability of cadmium carbonate precipitation dominated by carbon dioxide partial pressure (pCO2), resulting in the possibility of cadmium reactivation even at the safe potential. At the same time, the rapid heterogeneous intrusion of karst fissure water can cause drastic fluctuations in the rhizosphere microenvironment. Static control based on fixed rules cannot adapt to such high-frequency hydrological variations, causing regulatory actions to always lag behind the speciation of pollutants, making it difficult to achieve simultaneous compliance with cadmium and arsenic standards. Summary of the Invention

[0004] The purpose of this invention is to provide a method for graded remediation of cadmium and arsenic contaminated paddy fields using modified biochar through synergistic passivation, in order to solve the aforementioned problems existing in the prior art.

[0005] Technical solution: A graded remediation method for cadmium and arsenic contaminated paddy fields using modified biochar with synergistic passivation, including:

[0006] Collect multi-source heterogeneous environmental data from paddy field areas and construct a coupled feature vector that combines pollution characteristics, hydrogeological parameters, and ecological indicators;

[0007] The coupled feature vector is input into a pre-built risk assessment model to calculate the composite risk probability of the grid cell, and the paddy field area is divided into risk control zones of different levels accordingly.

[0008] Based on the risk management zone level, differentiated modified biochar materials with matching functions are deployed, including composite supported biochar for high-risk areas.

[0009] During the critical growth period of rice, the rhizosphere microenvironment was monitored in real time, and the boundary of the target regulation feasible domain for plots was constructed based on coupled feature vectors.

[0010] Based on the feasible boundary of target regulation, a multivariate collaborative regulation strategy is implemented to maintain the rhizosphere microenvironment state within a safe range that inhibits cadmium and arsenic activation.

[0011] Beneficial effects: This invention effectively solves the problem of remediation failure caused by regional Eh-pH decoupling and fissure water disturbance, and achieves simultaneous and long-term treatment of cadmium and arsenic. Attached Figure Description

[0012] Figure 1 A flowchart illustrating the steps of a graded remediation method using modified biochar for synergistic passivation in cadmium and arsenic contaminated paddy fields, provided in this application embodiment.

[0013] Figure 2 A flowchart illustrating the steps of implementing a multivariate collaborative regulation strategy as provided in this application embodiment.

[0014] Figure 3 A flowchart illustrating the steps for entering an Eh-pCO2-based synergistic regulation mode, as provided in this application embodiment.

[0015] Figure 4 A flowchart illustrating the steps for entering a model predictive control-based regulation mode, as provided in this application embodiment. Detailed Implementation

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

[0017] It should be noted that the terms include and have, and any variations thereof, are intended to cover non-exclusive inclusion. For example, a process, method, system, product, or device that includes a series of steps or units is not necessarily limited to those steps or units that are explicitly listed, but may include other steps or units that are not explicitly listed or that are inherent to such process, method, product, or device.

[0018] First, the basic implementation method of the synergistic passivation and graded remediation method using modified biochar for cadmium and arsenic contaminated paddy fields proposed in this invention is described. This basic implementation method constructs the overall technical framework and general processing flow of this invention, covering the implementation details of data acquisition, basic model construction, and conventional remediation operations. Those skilled in the art can combine and upgrade the technical elements based on this, in conjunction with the subsequent preferred embodiments. The basic implementation method is specifically as follows:

[0019] Step S1: Conduct multi-source heterogeneous environmental data collection to obtain multi-dimensional information on available heavy metals in soil, karst hydrogeology and farmland management, construct pollution feature vectors, and establish a multi-dimensional pollution risk assessment index system based on pollution feature vectors, including available cadmium concentration, available arsenic concentration, pH, Eh, fissure density and groundwater fluctuation amplitude.

[0020] For example, after rice harvest and before transplanting the next season, sampling points were set up in a 50m×50m grid to collect soil samples from the 0-20cm topsoil layer. Simultaneously, miniature groundwater monitoring tubes were deployed within a 5m radius of each sampling point to continuously monitor groundwater level dynamics throughout a complete hydrological year. The difference in groundwater level between the dry and rainy seasons was calculated as the seasonal variation in groundwater level. Soil samples were deployed in situ or quasi-in situ using DGT (diffusion gradient film) technology to determine the concentrations of bioavailable cadmium (DGT-Cd) and bioavailable arsenic (DGT-As). Parallel samples were also collected to determine soil pH, redox potential (Eh), organic matter content, and calcium carbonate content. Based on calcium carbonate content, pH, and groundwater ion concentration, the calcium carbonate saturation index (SIc) was calculated. High-density resistivity dating was conducted in paddy fields to obtain a three-dimensional distribution map of karst fissures within a depth range of 0-8m. The number of fissures per unit area and the average aperture were extracted, and the fissure flux index (KFI = fissure density × average aperture) was calculated to characterize the vertical migration potential of pollutants. A multispectral camera mounted on a drone was used to acquire the field-scale normalized difference vegetation index (NDVI) to identify anomalous vegetation areas under heavy metal stress, and the NDVI value was incorporated as a spatial covariate into the assessment system. Integrating the data obtained above, an eight-dimensional pollution-hydrology-ecology coupled feature vector was constructed, including: DGT-Cd, DGT-As, seasonal variation of groundwater level, fracture flux index (KFI), soil pH, calcium carbonate saturation index (SIc), arsenic valence ratio (As(III) / As(V)), and NDVI vegetation stress index; among which, the As(III) / As(V) ratio was obtained by valence analysis of the diffusion gradient film (DGT) adsorbent phase eluent using synchrotron X-ray absorption near-edge structure analysis (XANES) or high performance liquid chromatography-inductively coupled plasma mass spectrometry (HPLC-ICP-MS).

[0021] Step S2: Based on the evaluation index system, the pollution risk level is divided using the weighted comprehensive evaluation method. The paddy fields are divided into high-risk areas, medium-risk areas and low-risk areas, and a dynamic risk zoning map with geographic coordinates is generated.

[0022] Specifically, the constructed pollution-hydrology-ecology coupled feature vector is input into a graph convolutional neural network (GCN-Attention) model based on an attention mechanism. This model uses sampling points as nodes and spatial adjacency and hydraulic connectivity as edge weights to learn the heterogeneous diffusion pattern of pollutants in the fracture-soil medium. The model outputs the cadmium-arsenic combined risk probability distribution for each grid cell and introduces the Monte Carlo Dropout method to quantify the prediction uncertainty, automatically triggering supplementary sampling instructions for areas with high uncertainty. Based on the risk probability and uncertainty interval, a fuzzy C-means clustering (FCM) algorithm is used to generate soft boundary risk zones, allowing transitional zones to simultaneously belong to two risk levels and assigning membership weights. Combining the root distribution depth of rice and the vertical flux direction of karst fractures, hydrogeochemical corrections are performed on the zone boundaries: if the dominant water flow in the fractures is downward, the vertical influence radius of the high-risk zone is expanded; if lateral runoff is significant, the high-risk buffer zone is extended along the groundwater flow direction. The final partitioning results are encoded in GeoJSON format and embedded into the remediation decision engine to drive the modified biochar variable application system to mix and apply different types of materials according to their membership ratio.

[0023] Step S3: Based on the risk zoning results, apply functionally matched modified biochar in a targeted manner: apply iron-manganese dioxide co-supported biochar in high-risk areas to simultaneously adsorb pentavalent arsenic (As(V)) and promote the adsorption of divalent cadmium ions (Cd). 2+ This forms cadmium hydroxide (Cd(OH)2) or cadmium ferrite (CdFe2O4) precipitates; in medium-risk areas, phosphate-impregnated modified biochar is applied, releasing phosphate ions (PO42-). 3- Inducing cadmium phosphate (Cd3(PO4)2) precipitation and inhibiting trivalent arsenic (As(III)) desorption; applying alkaline-thermal treated microporous biochar in low-risk areas enhances soil aggregate stability and controls the vertical leaching of cadmium / arsenic (Cd / As) along karst fissures.

[0024] For example, using waste corn stalks or fir branches from a certain region as raw materials, pyrolysis was performed at 500±20℃ for 2 hours under a nitrogen atmosphere to obtain primary biochar with a pH of 9.5-10.5 and a specific surface area of ​​150-250 m². 2 / g. For high-risk areas, raw biochar was impregnated in a mixed solution of ferric nitrate (Fe(NO3)3) and manganese sulfate (MnSO4) (iron:manganese (Fe:Mn) molar ratio 3:1), stirred in a 60℃ water bath for 6 h, and then sodium hydroxide (NaOH) solution was added dropwise until pH=10. After aging for 12 h, washing and drying, it was calcined at 300℃ for 1 h to obtain iron-manganese bioxide nanoclusters supported biochar. The Fe / Mn oxide particle size on its surface was 5–20 nm, and the maximum adsorption capacity for As(V) was ≥45 mg / g. Simultaneously, Cd was co-precipitated... 2+ A stable Cd(OH)₂ or CdFe₂O₄ phase is formed. For medium-risk areas, the raw biochar is impregnated in a potassium dihydrogen phosphate (KH₂PO₄) solution (1.0 mol / L) at a solid-liquid ratio of 1:5, ultrasonically treated for 30 min, dried at 80 °C, and then calcined at 350 °C under an inert atmosphere for 2 h to obtain phosphate-intercalated mesoporous biochar with a mesopore content ≥65%, capable of simultaneously releasing PO₄²⁻. 3- Cd3(PO4)2 precipitation was induced, and As(III) desorption was inhibited through surface protonation. For low-risk areas, the raw biochar was subjected to alkaline heat treatment at 90℃ for 3 hours in a 1.5 mol / L NaOH solution, washed with water until neutral, and dried at 60℃ to obtain alkaline-heat modified biochar rich in micropores (<2 nm) and with increased oxygen-containing functional groups on the surface. This modified biochar was used to enhance soil aggregate stability and control the vertical leaching of Cd / As with karst fissure water. Based on the risk zoning map, the corresponding type of modified biochar was precisely applied to each area using a variable-rate spreader at the following rates: 8–12 t / ha in high-risk areas, 5–8 t / ha in medium-risk areas, and 2–4 t / ha in low-risk areas. This was applied as a single basal application before tillage to ensure uniform mixing with the 0-20 cm topsoil layer.

[0025] Step S4: Before rice transplanting, apply the corresponding type of modified biochar according to the zoning variables. During the two key growth stages of tillering and heading, based on real-time soil Eh and pH monitoring data, regulate the field surface water depth and drainage frequency through the intelligent irrigation system to stably maintain the rhizosphere soil redox potential (Eh) within the synergistic passivation window of +100 mV to -150 mV.

[0026] Specifically, the first and second phases of Eh regulation windows were initiated 15-20 days after rice transplanting (early tillering stage) and 55-60 days after transplanting (early heading stage), respectively, with each window lasting 7-10 days. Within each regulation window, the soil redox status was collected in real time using a multi-node Eh / pH wireless sensor array deployed in the rhizosphere, and weather forecasts and groundwater level predictions were acquired simultaneously. The measured Eh value was compared with the target window (+100mV to -150mV). If Eh > -50mV, intermittent shallow irrigation was initiated: the water depth on the field surface was maintained at 2-3cm for 24 hours, followed by drainage and drying for 48 hours to promote the formation of Fe / Mn oxides to adsorb arsenic (As); if Eh < -150mV, drainage was immediately initiated and aeration was carried out for 36 hours to prevent the excessive release of As(III) and excessive dissolution of cadmium sulfide (CdS). During Eh regulation, the ratio of sulfate to nitrate in the irrigation water is dynamically adjusted: when the goal is to enhance Cd fixation, sulfate-containing (SO4) ions are introduced. 2- Water sources (such as gypsum dissolution water) promote CdS precipitation; when the goal is to inhibit As activation, nitrate ions (NO3) are added. - As an alternative electron acceptor, it inhibits the activity of iron-reducing bacteria and maintains the stable adsorption state of As(V). Based on the Eh response hysteresis characteristics, a soil Eh prediction model based on LSTM (Long Short-Term Memory Network) is constructed to predict the Eh change trend 12-24 hours in advance and automatically optimize the duration and intensity of the next irrigation-drainage cycle, realizing feedforward-feedback composite control.

[0027] Step S5: During the rice maturity period, in-situ DGT technology and portable X-ray fluorescence spectrometer were used to simultaneously measure the concentrations of bioavailable cadmium and arsenic in rhizosphere and non-rhizosphere soils, and grain samples were collected to detect the accumulation of heavy metals and assess the remediation status.

[0028] Specifically, during the milk-ripe to waxy-ripe stages of rice, DGT (diffusion gradient film) probe arrays were deployed in each risk zone at two depths: 0-10 cm (rhizosphere zone) and 20-30 cm (non-rhizosphere / fracture-affected zone). After exposure for 24-48 h, the arrays were retrieved, and the DGT-Cd and DGT-As fluxes were determined by inductively coupled plasma mass spectrometry (ICP-MS) to quantify the spatiotemporal distribution of bioavailable heavy metals. Simultaneously, a portable micro-area X-ray fluorescence spectrometer (μ-XRF) was used to perform in-situ scanning of the paddy field soil profile, obtaining two-dimensional spatial co-distribution maps of Cd, As, Fe, Mn, and P elements. Principal component analysis was then used to identify the spatial coupling relationships between heavy metals and passivation components (such as Fe / Mn oxides and phosphate phases). Three representative plants were randomly selected from each zone, and samples were collected from different organs (roots, stems, leaves, and brown rice). After microwave digestion, the total Cd and As contents were determined by ICP-MS, and the translocation coefficient (TF = brown rice concentration / root concentration) and enrichment coefficient (BCF = brown rice concentration / soil DGT concentration) were calculated. Synchrotron X-ray absorption fine structure (SR-XAFS) analysis was performed on soil samples from high-risk areas before and after remediation to resolve the coordination environment of Cd (whether Cd3(PO4)2, CdS, or cadmium-oxygen-iron (Cd-O-Fe) bonds were formed) and the valence state and adsorption forms of As (As(V)-FeOOH, As(III)-thiols, etc.), verifying the molecular-scale passivation mechanism of modified biochar. DGT data, μ-XRF images, phytoaccumulation results, and SR-XAFS mechanism evidence were fused from multiple sources to construct a four-dimensional evaluation matrix of effectiveness, mobility, bioaccumulation, and molecular fixation, used to determine remediation achievement and guide material optimization in the next cycle.

[0029] Step S6: Based on the remediation effect assessment data, dynamically update the risk zoning results and modified biochar application strategy for the next cycle, forming a six-step closed-loop graded remediation mechanism of risk assessment, zoning customization, material application, process control, effect monitoring, and strategy optimization.

[0030] Specifically, the DGT-Cd, DGT-As, grain Cd / As content, Eh control records, and modified biochar type and application rate obtained during this remediation cycle were structured and stored to form a standardized remediation case archive. The actual remediation efficiency of each risk zone was calculated: for high-risk zones, grain As ≤ 0.3 mg / kg and Cd ≤ 0.2 mg / kg were used as the criteria for compliance; for medium-risk zones, a reduction rate of available Cd ≥ 40% or an As reduction rate ≥ 35% were used as the criteria for effectiveness; and for low-risk zones, no new leaching flux was used as the criterion for stability. If a zone did not achieve the expected results, its dominant limiting factors were analyzed retrospectively: if As still exceeded the standard, it was determined to be due to insufficient redox regulation or low iron and manganese loading; if Cd did not decrease significantly, it was determined to be due to insufficient phosphate release or Eh window shift; if pollutant downward migration intensified, it was determined to be due to insufficient microporous biochar control capacity. Based on the above attribution results, the remediation parameters for the next cycle are dynamically adjusted: for areas with insufficient redox regulation, the Eh control window is activated earlier and the aeration time is extended; for areas with low loading, the iron-manganese dual oxide loading ratio is increased to Fe:Mn = 4:1; for areas with insufficient resistance control, the alkaline heat treatment temperature is increased from 90℃ to 110℃ to enhance micropore development. The input conditions (Initial Pollution Risk Composite Index (CPRI), fracture density), intervention measures (biochar type, Eh control scheme), output results (remediation efficiency, limiting factors), and optimization suggestions are encoded into triplets and stored in the local remediation knowledge base to support intelligent recommendation of remediation strategies for similar fields in the future.

[0031] Simultaneously, a basic tiered remediation system for cadmium and arsenic contaminated paddy fields using modified biochar with synergistic passivation is also provided, including:

[0032] The intelligent pollution risk assessment module is used to integrate soil DGT available heavy metal data, karst fissure three-dimensional detection results, groundwater level dynamics and vegetation remote sensing indices to construct an eight-dimensional pollution-hydrology-ecology coupled feature vector and generate high / medium / low risk zoning maps.

[0033] The modified biochar customized supply unit has three independent silos, which store iron-manganese double oxide co-supported biochar, phosphate impregnated modified biochar and alkali-heat treated microporous biochar respectively, and are equipped with metering screws and pneumatic conveying devices.

[0034] The variable precision application actuator integrates a GNSS positioning terminal and an electro-hydraulic proportional valve. It receives risk zoning map signals, automatically switches hoppers and adjusts the feeding rate according to geographical coordinates, and realizes differentiated variable application of modified biochar in different risk zones.

[0035] The Eh-pH coordinated regulation irrigation subsystem includes a wireless sensor array deployed in the rhizosphere zone, an intelligent water valve controller, and a water source allocation device. Based on the rice growth stage and measured Eh / pH values, it dynamically adjusts the field surface water depth, drainage frequency, and SO4 levels in the irrigation water.2- / NO3 - The ratio was used to stabilize the soil redox potential within a window of +100mV to -150mV.

[0036] The in-situ assessment terminal for remediation effects integrates a portable μ-XRF spectrometer and a DGT probe recovery workstation to simultaneously acquire the spatial distribution of heavy metals and bioavailable fluxes in soil profiles.

[0037] A closed-loop optimization decision engine is used to compare the heavy metal detection results of grains with the preset compliance threshold, identify and repair failure factors, and generate biochar type adjustment instructions and Eh control parameter optimization schemes for the next cycle.

[0038] Furthermore, the three types of biochar in the customized modified biochar supply unit meet the following physicochemical properties: the Fe / Mn oxide nanoclusters of the iron-manganese dioxide co-supported biochar have a particle size of 5–20 nm, a maximum adsorption capacity for As(V) ≥45 mg / g, and can promote Cd adsorption. 2+ Forms Cd(OH)2 or CdFe2O4 precipitates; phosphate-impregnated modified biochar has a mesoporous content ≥65% and a phosphorus content of 5-15wt%, which can slowly release PO4. 3- Induces Cd3(PO4)2 precipitation and inhibits As(III) desorption; alkaline-thermal treated microporous biochar has a micropore volume ratio ≥70% and a specific surface area of ​​200-300 m². 2 / g, after alkaline heat treatment at 90-110℃, the density of oxygen-containing functional groups on the surface increases by more than 30%, which is used to enhance the stability of soil aggregates and inhibit the vertical migration of Cd / As along karst fissures.

[0039] Based on the above-mentioned basic implementation method, and considering the special geological background of a certain region with high calcium content, strong karst fissure development, and violent hydrological fluctuations, the following preferred embodiment, which includes a more refined mechanism model and intelligent control algorithm, is further proposed.

[0040] like Figure 1 As shown, a method for graded remediation of cadmium and arsenic contaminated paddy fields using modified biochar with synergistic passivation includes the following steps:

[0041] Multi-source heterogeneous environmental data were collected from paddy field areas to construct a coupled feature vector that combines pollution characteristics, hydrogeological parameters, and ecological indicators.

[0042] Specifically, multi-source heterogeneous environmental data refers to comprehensive datasets derived from different sensors, sampling frequencies, and spatial dimensions. For example, it can include soil chemical property data obtained through physical sampling, subsurface structure data obtained through geophysical exploration, and surface vegetation data obtained through remote sensing technology. The process of constructing coupled feature vectors essentially involves spatiotemporal alignment and standardization of the multi-source heterogeneous environmental data, and extracting feature components that play a crucial indicative role in the migration and transformation of cadmium and arsenic. These feature components not only include traditional pollutant concentration information but also incorporate parameters that characterize specific geological backgrounds, providing a complete data foundation for subsequent risk assessment and precise regulation. In practical applications, a gridded monitoring network can be deployed in paddy field areas, combined with drone inspections and fixed-point observation stations, to achieve digital perception of all elements of the paddy field environment.

[0043] The coupled feature vector is input into a pre-built risk assessment model to calculate the composite risk probability of grid cells, and the paddy field area is divided into risk control zones of different levels based on the composite risk probability.

[0044] In this embodiment, the risk assessment model is a computational model capable of establishing a nonlinear mapping relationship between feature vectors and pollution risks. The composite risk probability refers to the comprehensive likelihood of excessive accumulation, activation, or migration of heavy metals such as cadmium and arsenic in paddy field soil under current environmental conditions. Based on the calculated composite risk probability, cluster analysis or threshold determination methods can be used to discretize the continuous paddy field area into several risk management zones with clearly defined spatial boundaries. For example, the risk management zones include high-risk, medium-risk, and low-risk zones. For instance, areas with a composite risk probability higher than a preset high-risk threshold can be designated as high-risk management zones, areas lower than a preset low-risk threshold as low-risk management zones, and the remaining areas as medium-risk management zones. This zoning management strategy can break away from the traditional one-size-fits-all remediation model, achieving differentiated governance for areas with different levels of pollution and improving the utilization efficiency of remediation resources.

[0045] Based on the risk management zone level, modified biochar materials with differentiated functions are deployed. Modified biochar materials include composite-supported biochar for high-risk areas.

[0046] Specifically, differentiated application refers to using different types, formulations, or application rates of remediation materials for different risk management zones. For high-risk areas, the focus is on simultaneously controlling the combined pollution of cadmium and arsenic; therefore, the selected composite-supported biochar needs to be capable of simultaneously passivating both cationic cadmium and anionic arsenic. For example, this composite-supported biochar can achieve arsenic adsorption by loading specific metal oxides onto its surface, while simultaneously using the biochar's alkalinity and porous structure to immobilize cadmium. For medium- and low-risk areas, the focus can be on the targeted regulation of specific pollutants or preventing pollutant leaching into deeper soil layers, respectively. Through a functionally matched application strategy, the physicochemical properties of the remediation materials can be highly compatible with the dominant pollution processes in the area, ensuring remediation effectiveness while reducing treatment costs.

[0047] During the critical growth period of rice, the rhizosphere microenvironment is monitored in real time, and the boundary of the target regulation feasible domain for plots is constructed based on coupled feature vectors.

[0048] In other words, during the critical growth period of rice, the rhizosphere microenvironment is monitored in real time after the differentiated application of modified biochar materials, and the boundary of the target regulation feasible domain for plots is constructed based on the coupled feature vector.

[0049] In this embodiment, the key growth stages of rice typically include the tillering and heading stages, which are highly sensitive to heavy metal absorption. The rhizosphere microenvironment refers to the physicochemical properties of the soil surrounding the rice roots, mainly including key indicators such as redox potential and pH. The site-adaptive target control feasible domain boundary refers to the safe range for inhibiting heavy metal activation, dynamically calculated based on the specific geological, hydrological, and pollution characteristics of each site. Unlike traditional fixed threshold control, this feasible domain boundary is dynamically changing and can adaptively adjust to different site characteristic parameters such as the degree of fissure development and groundwater level changes. For example, for sites with a high degree of fissure development, the control boundary may be tightened accordingly to prevent rapid pollutant migration; while for sites with high carbonate content, the boundary may be modified based on buffering capacity. This adaptive mechanism effectively overcomes the high spatial heterogeneity of the region, ensuring that the control target remains within the optimal safe range.

[0050] Based on the feasible boundary of target regulation, a multivariate collaborative regulation strategy is implemented to maintain the rhizosphere microenvironment state within a safe range that inhibits cadmium and arsenic activation.

[0051] In other words, based on the boundary of the feasible domain of target regulation, a multivariate collaborative regulation strategy is implemented to generate regulation instructions and control the actions of irrigation and drainage actuators, so as to maintain the state of the rhizosphere microenvironment within the safe domain of inhibiting cadmium and arsenic activation.

[0052] Specifically, the multivariate synergistic regulation strategy refers to the comprehensive use of multiple regulatory means (such as water level regulation and the addition of chemical regulators) to synergistically control the rhizosphere microenvironment, rather than relying solely on single irrigation and drainage measures. The goal of implementing this strategy is to ensure that the real-time monitored rhizosphere microenvironment state always falls within the safe zone defined by the feasible boundary of the target regulation domain. Within this safe zone, cadmium mainly exists in a stable precipitate form, while arsenic is primarily adsorbed onto the soil solid surface, minimizing the effective concentration and bioavailability of both in the soil solution. For example, when the rhizosphere state is detected to deviate from the safe zone, the system automatically calculates the optimal combination of regulatory actions, such as adjusting irrigation water volume or adding specific regulators, to quickly bring the rhizosphere environment back to a safe state. Through synergistic regulation based on dynamic boundaries, the contradictory redox requirements of cadmium and arsenic are resolved, achieving precise treatment of complex pollution.

[0053] In a preferred embodiment, the coupled feature vector includes: bioavailable cadmium and bioavailable arsenic flux data; fracture flux index and karst water redox potential characteristic values ​​characterizing the degree of karst development; seasonal variation of groundwater level characterizing groundwater dynamics; carbonate buffer strength index and calcium carbonate saturation index characterizing soil chemical buffering capacity; pH ​​value and arsenic valence ratio characterizing soil physicochemical properties; and normalized vegetation index characterizing vegetation growth status.

[0054] In this embodiment, the different dimensions of the coupled feature vectors together constitute a complete dataset describing the paddy field environment. The bioavailable cadmium (DGT-Cd) and bioavailable arsenic (DGT-As) flux data were obtained through an in-situ diffusion gradient film (DGT) device, directly reflecting the bioavailability of heavy metals, rather than just their total amount. The fracture flux index (KFI), characterizing the degree of karst development, is a comprehensive geological parameter. Its calculation formula can be expressed as fracture density multiplied by the average fracture aperture, used to quantify the potential for rapid migration of water and pollutants through karst fractures. The karst water redox potential characteristic value refers to the statistical characteristic value of the redox potential (Eh) measured in karst fracture water, used to indicate the redox background of the groundwater environment.

[0055] Furthermore, the seasonal variation in groundwater level (ΔGW) refers to the difference between the groundwater level during the dry and wet seasons within a hydrological year. This parameter directly affects the risk of vertical migration of pollutants due to groundwater fluctuations. The carbonate buffer strength index (CBI) is a key parameter specific to certain regions, used to characterize the soil's ability to resist changes in pH. Its calculation can be derived based on the concentrations of calcium ions, bicarbonate ions, and the solubility product constant of calcium carbonate in the soil solution. The calcium carbonate saturation index (SIc) reflects the precipitation or dissolution trend of calcium carbonate in the soil solution. The arsenic valence ratio (As(III) / As(V)) refers to the ratio of trivalent arsenic to pentavalent arsenic concentrations. Since trivalent arsenic is far more toxic than pentavalent arsenic, this ratio is an important indicator for assessing ecological risk. The normalized difference vegetation index (NDVI) is a vegetation growth status indicator extracted from remote sensing imagery, used to help identify crop areas under heavy metal stress. The fusion of these multidimensional features enables the system to capture complex physical, chemical, and biological coupling processes in the environment.

[0056] In one possible implementation, after calculating the composite risk probability of the grid cells, the following is also included:

[0057] Based on the composite risk probability and the preset sampling cost factor, the control failure value index of each grid cell is calculated.

[0058] The control failure value index is compared with the preset trigger threshold.

[0059] When the control failure value index of any grid cell exceeds the trigger threshold, a supplementary sampling instruction is generated for that grid cell, and the coupling feature vector is updated using the supplementary sampling instruction.

[0060] Alternatively, based on the composite risk probability and the sampling cost factor pre-set according to terrain accessibility, the control failure value index of each grid cell is calculated; the control failure value index is compared with the preset trigger threshold; when the control failure value index of any grid cell exceeds the trigger threshold, a supplementary sampling instruction for that grid cell is generated and sent to the sampling execution terminal; the supplementary sampling data returned by the sampling execution terminal is received, and the historical data of the corresponding grid cell in the coupled feature vector is replaced with the supplementary sampling data to complete the update.

[0061] Specifically, a closed-loop active sampling logic based on control failure risk is constructed. Traditional sampling strategies are often based on predictive uncertainty, while the strategy proposed in this embodiment is based on control failure value. The control failure value index (I...) iThe control failure value index (I) is a comprehensive quantitative indicator that measures the potential loss value resulting from the failure of subsequent control strategies if the risk assessment of the current grid cell deviates. The preset sampling cost factor reflects the human, material, or time costs required for supplementary sampling in that grid cell. i When the preset trigger threshold is exceeded, it means that the benefit of supplementary sampling at that point (reducing the risk of control failure) outweighs the cost. The system then generates a supplementary sampling instruction containing the sampling coordinates and sampling type. By backfilling the database with the latest high-precision data obtained from the supplementary sampling and updating the coupled feature vector, the bias of the risk assessment model can be corrected in real time, realizing the online evolution of the model.

[0062] In one exemplary embodiment, the control failure value index I i Calculated using the following formula:

[0063] I i =(P violate，i ·G i ·(1-C i )) / (Cost i +ε);

[0064] Where P violate，i G is the probability of regulation exceeding the limit calculated based on a pre-configured prediction model. i C is a risk boundary gradient index calculated based on composite risk probabilities. i To predict the repair compliance confidence level, Cost i ε is a preset sampling cost factor, and ε is a positive number.

[0065] In this embodiment, a specific mathematical definition of the control failure value index is given. Where P violate，i This represents the probability, based on the current predictive model, that the control target (such as the Eh value) of this grid cell will exceed the safe and feasible region in the future. i The risk boundary gradient index is used to characterize the degree to which a grid cell is located in the risk level transition zone. Specifically, G... i This can be obtained by calculating the absolute value of the difference between the composite risk probability of the grid cell and the mean of the composite risk probabilities of its neighboring grid cells. For example,

[0066] G i = | P risk，i - (1 / N) *∑ j=1 N (P risk，j )|;

[0067] Among them G i P represents the risk boundary gradient index for the i-th grid cell; risk，iLet ∑ be the composite risk probability value of the i-th grid cell; j=1 N (P risk，j ) represents the sum of the composite risk probabilities of all adjacent grid cells j of the i-th grid cell; N is the total number of adjacent grid cells, usually 4 or 8; | represents the absolute value operation. G i The larger the value, the more critical it is to the point being on the edge of risk change; its state determination is crucial for defining the partition boundary. C i This represents the confidence level of the repair meeting the standard, that is, the degree to which the model believes that the repair at this point will achieve the safety standard. (1-C) i This indicates the risk of failing to meet the standards. Cost i This is the sampling cost factor; for example, this value is set higher for mountainous grids with inaccessible terrain, and lower for easily accessible plains grids. ε is a very small positive number (e.g., 10). -6 (), used to prevent the denominator from being zero.

[0068] For example, suppose that for a specific grid cell i, the prediction model calculates the regulation overshoot probability P. violate，i A value of 0.4 indicates a 40% probability that the regulation will fail; the risk boundary gradient index G is calculated based on spatial analysis. i A value of 0.8 indicates that this point is on the verge of a strong risk change; the predicted repair achievement confidence level C i The preset sampling cost factor is 0.6. i A value of 10 represents a medium sampling difficulty. Substituting the above values ​​into the formula, we can obtain I. i =(0.4×0.8×(1-0.6)) / (10+10 -6 )≈0.0128. If the preset trigger threshold I th Since 0.0128 > 0.01, the system will determine that the grid cell has high resampling value and automatically trigger a resampling command. This ensures that sampling resources are always concentrated in critical areas with high risk, ambiguous boundaries, and severe consequences if out of control, thereby improving the overall system's control reliability at minimal cost.

[0069] In a preferred implementation, for grid cells designated as high-risk zones, the composite-loaded biochar is iron-manganese-lanthanum ternary co-loaded modified biochar.

[0070] In this embodiment, high-risk areas typically refer to regions with severe cadmium and arsenic combined pollution and high bioavailability. The iron-manganese-lanthanum ternary co-supported modified biochar tailored for high-risk areas can solve the problem that traditional iron-manganese biochar has a negatively charged surface in high-pH (typically greater than 7.5) soils, resulting in reduced adsorption capacity for anionic arsenate.

[0071] Iron-manganese-lanthanum ternary co-supported modified biochar was prepared by in-situ growth of iron-manganese oxides on the surface of a biochar support and doping with lanthanum, wherein the molar ratio of iron, manganese and lanthanum was 3:1:0.5. It was used to enhance the adsorption of arsenate in high pH soil environment by maintaining the surface positive charge.

[0072] Specifically, the preparation process of iron-manganese-lanthanum ternary co-supported modified biochar can be carried out using a co-precipitation-calcination method. Using common agricultural wastes in a certain region, such as corn stalks or mulberry branches, as raw materials, raw biochar is prepared by slow pyrolysis at 500°C under nitrogen protection. A mixed precursor solution containing ferric nitrate, manganese sulfate, and lanthanum nitrate is prepared, with the molar ratio of metal ions in the solution controlled at ferric:manganese:lanthanum (La) = 3:1:0.5. The raw biochar is impregnated in this mixed solution and continuously stirred under a 60°C water bath. During this period, sodium hydroxide solution is slowly added dropwise to adjust the pH to approximately 10, inducing in-situ precipitation of metal hydroxides in the pores and surface of the biochar. After 12 hours of aging, washing, and drying, a second calcination is performed at 300°C to obtain modified biochar with a nano-cluster distribution. The introduction of lanthanum is particularly crucial because lanthanum oxide has a high isoelectric point of about 10.0, which ensures that the material maintains a positive surface potential in alkaline soil environments, thereby strongly adsorbing arsenate ions through electrostatic attraction and coordination exchange. At the same time, iron and manganese oxides promote the conversion of cadmium ions into stable cadmium ferrite or cadmium hydroxide precipitates through co-precipitation.

[0073] In a further embodiment, the modified biochar material with differentiated delivery function matching further includes:

[0074] For grid cells classified as medium-risk areas, phosphate-impregnated modified biochar was applied to release phosphate ions to induce cadmium phosphate precipitation.

[0075] In this embodiment, the medium-risk area refers to an area where the total amount of heavy metals exceeds the standard but the activity is relatively low, or where only a single pollutant exceeds the standard. The preparation method of phosphate-impregnated modified biochar can be a liquid-phase impregnation method. Specifically, the raw biochar is impregnated in a potassium dihydrogen phosphate solution with a concentration of 1.0 mol / L, controlling the solid-liquid ratio at 1:5, and ultrasonically assisted for 30 minutes to promote the entry of phosphate into the mesoporous structure of the biochar. It is then dried at 80 degrees Celsius and calcined at 350 degrees Celsius in an inert atmosphere, so that the phosphate is immobilized on the surface of the biochar in the form of metaphosphate or pyrophosphate. When this material is applied to the soil, it can slowly release phosphate ions, which react with cadmium ions in the soil solution to form a cadmium phosphate precipitate with extremely low solubility, achieving long-term cadmium fixation. Simultaneously, the phosphate groups on the surface of the biochar can also inhibit the desorption of trivalent arsenic through surface complexation.

[0076] For grid cells classified as low-risk zones, microporous modified biochar treated with alkali heat was applied to control the vertical leaching of heavy metals along karst fissures by utilizing the microporous structure.

[0077] In this embodiment, the primary focus in the low-risk area is preventing pollutants from migrating downwards along karst fissures with rainfall or irrigation water and contaminating groundwater. The preparation process of alkaline-thermal-treated microporous modified biochar emphasizes the regulation of pore structure. Specifically, the raw biochar is placed in a sodium hydroxide solution with a concentration of 1.5 mol / L and subjected to alkaline-thermal treatment at 90°C to 110°C for 3 hours. The strong alkaline etching effect can open up the pores of the biochar, increasing the volume ratio of micropores (pore size less than 2 nanometers). The specific surface area of ​​the prepared material can typically reach 200 to 300 square meters per gram. The microporous structure can adsorb dissolved organic matter-heavy metal complexes in the soil solution through physical interception and capillary coagulation, enhancing the water stability of soil aggregates and effectively blocking the vertical migration channels of pollutants to deep groundwater. In actual application, it can be combined with risk zoning map and a variable spreader with satellite positioning function can be used to automatically switch the hopper for precise application according to the preset formula. For example, the application amount in high-risk areas is 8 to 12 tons per hectare, and in medium-risk areas it is 5 to 8 tons per hectare.

[0078] like Figure 2 As shown, according to one aspect of this application, implementing a multivariate collaborative regulation strategy includes:

[0079] Obtain the carbonate buffer strength index in the coupled feature vector;

[0080] The carbonate buffer strength index is compared with a preset buffer threshold.

[0081] When the carbonate buffer strength index is greater than the buffer threshold, it enters a regulation mode based on Eh-pCO2 synergy.

[0082] In this embodiment, the carbonate buffer strength index (CBI) is the core criterion for determining whether to activate the Eh-pCO2 synergistic regulation mode. Specifically, the formula for calculating CBI can be expressed as:

[0083] CBI=([Ca 2+ ]×[HCO3 - ] / K sp )×(1 / (1+10 pH-pKa2 ));

[0084] Where [Ca 2+ ] and [HCO3 - [ ] represents the molar concentrations of calcium ions and bicarbonate ions in the soil solution, respectively; K sp is the solubility product constant of calcium carbonate, which is approximately 3.36 × 10⁻⁶ at 25°C.-9 pKa2 is the second-order dissociation constant of carbonic acid, approximately 10.33; pH is the measured soil acidity / alkalinity. The CBI calculation formula quantifies the sensitivity of the soil solution to changes in carbon dioxide partial pressure. For example, when the measured soil solution calcium ion concentration is 2 mmol / L, bicarbonate concentration is 3 mmol / L, and pH is 7.8, the calculated CBI value is approximately 2.5. If the preset buffer threshold is set to 2.0, the site is determined to be in a high buffer state of calcium carbonate supersaturation. Traditional single Eh regulation will fail due to pH lock-in, therefore the system automatically switches to a regulation mode based on Eh-pCO2 synergy.

[0085] like Figure 3 As shown, in a preferred embodiment, the Eh-pCO2 synergistic regulation mode includes:

[0086] Based on the real-time monitored rhizosphere redox potential, the preset Eh baseline, and the karst water redox potential characteristic value in the coupled feature vector, the karst water intrusion index is calculated.

[0087] A dual closed-loop control architecture was constructed, including an Eh control loop triggered by the karst water intrusion index and a pCO2 control loop triggered by the carbonate buffer strength index.

[0088] In this embodiment, the Karst Water Intrusion Index (KWII) is used to quantitatively characterize the degree of disturbance of the rhizosphere redox environment by fissure water intrusion. Its calculation formula is as follows:

[0089] KWII=(Eh measured -Eh baseline ) / (Eh karst -Eh baseline )×100%;

[0090] Among them Eh measured For real-time readings from the rhizosphere sensor; Eh baseline The steady-state baseline value during periods without rainfall, for example -100 millivolts; Eh karst The redox potential (Eh) is a characteristic value of karst water, typically in an oxygen-rich state, for example, +300 mV. If Eh is measured to be +100 mV in real time, then KWII = (100 - (-100)) / (300 - (-100)) = 50%, indicating that the rhizosphere environment has been disturbed by fissure water intrusion. The dual closed-loop control architecture refers to the system operating two control loops in parallel: the Eh control loop is mainly responsible for counteracting external oxidation disturbances identified by KWII through water management, thereby inhibiting arsenic activation; the pCO2 control loop is mainly responsible for maintaining the thermodynamic stability of cadmium carbonate precipitation by adjusting the partial pressure of carbon dioxide, thereby fixing cadmium.

[0091] In one possible implementation, the specific execution logic of the pCO2 control loop is as follows:

[0092] Maintain the pCO2 of the rhizosphere soil within the target window of 0.003 atm to 0.03 atm;

[0093] When pCO2 is below the lower limit of the target window, the partial pressure is increased by stimulating rhizosphere microbial respiration with the application of organic acid solution.

[0094] When pCO2 is higher than the upper limit of the target window, the partial pressure is reduced by decreasing the thickness of the water layer on the field surface to promote gas dissipation.

[0095] Specifically, the target window of 0.003 atm to 0.03 atm is a safe range derived from the cadmium carbonate precipitation-dissolution equilibrium. When the rhizosphere pCO2 is detected to be below 0.003 atm, it indicates a risk of cadmium carbonate dissolution. At this point, the system generates an instruction to apply an organic acid solution to the rhizosphere via the drip irrigation system. Preferably, a 0.5% citric acid solution can be used at a rate of 500 liters per hectare. Citric acid, as an readily available carbon source, can rapidly stimulate rhizosphere microbial respiration, releasing a large amount of carbon dioxide in a short time, thereby increasing the partial pressure. Conversely, when the pCO2 is above 0.03 atm, it indicates a risk of soil acidification. The system will then instruct the water level to be lowered to below 2 cm or drained, utilizing atmospheric and soil gas exchange to remove excess carbon dioxide.

[0096] In a further embodiment, the Eh control loop also includes an active prevention mechanism, specifically:

[0097] Obtain meteorological and rainfall forecast data, and calculate the time delay parameters of karst fissure water intrusion into the rhizosphere by combining the fissure flux index in the coupled feature vector;

[0098] Based on the time delay parameter, preventive control measures are initiated within a preset time window before the predicted time of fissure water intrusion. These preventive control measures include increasing the depth of the surface water layer to use hydrostatic pressure to block the lateral infiltration of fissure water.

[0099] In this embodiment, the time lag characteristic of karst hydrological processes is utilized to achieve preventative control. For example, the time lag parameter T... lag It can be estimated using empirical formulas:

[0100] T lag =a×R -0.5 ×(1+b×θ);

[0101] Where R represents the forecast rainfall, θ represents the current soil moisture content, and a and b are geological coefficients related to the fissure flux index (KFI). For example, if a heavy rainstorm is forecast for the next 24 hours, the model calculates that fissure water will reach the rhizosphere 36 hours after the rain, i.e., T.lag =36h, the system will initiate preventative irrigation 24 hours after the rain, or 12 hours before intrusion, raising the water depth on the field surface to 10 cm. The hydrostatic pressure created by the deep water acts like a hydraulic shield, physically suppressing the lateral inflow of groundwater from fissures and protecting the rhizosphere's restorative environment from damage.

[0102] The coupling correction rule is executed to identify and coordinate the control conflict between the Eh control loop and the pCO2 control loop, and the final control command is generated.

[0103] In a preferred implementation, the coupling correction rule is executed, including:

[0104] When the Eh control loop requests flooding measures and the pCO2 control loop detects that the carbon dioxide partial pressure is close to the upper limit of the window, shallow water irrigation is implemented in conjunction with surface mulching of the field water body to isolate oxygen while allowing some carbon dioxide to escape through the edge of the water body.

[0105] Specifically, when the KWII indicator suggests flooding to reduce Eh, but pCO2 is already close to the dangerous upper limit of 0.03 atm, simple deep flooding will prevent CO2 from escaping and cause it to exceed the limit. In this case, the system implements a shallow water + mulching strategy: maintaining a shallow water layer of 2 to 3 cm to allow some gas exchange, while covering the surface of the water body with a biodegradable liquid film or agricultural film. Mulching can physically isolate atmospheric oxygen from entering (meeting the Eh reduction requirement), while the permeability at the edge of the shallow water layer can prevent excessive CO2 accumulation in the rhizosphere (meeting the pCO2 control requirement).

[0106] When the Eh control loop requests aeration and the pCO2 control loop detects that the carbon dioxide partial pressure is close to the lower limit of the window, intermittent short-term aeration is implemented to reduce the amount of carbon dioxide escape while increasing the redox potential.

[0107] Specifically, when Eh is too low, for example, below -200 mV, increasing the risk of arsenic release and requiring aeration, but pCO2 is already low at around 0.003 atm, continuous aeration will blow away all CO2, leading to cadmium carbonate dissolution. In this case, the system implements an intermittent aeration strategy: for example, turning on the microporous aeration device for 15 minutes every 2 hours. Pulsed oxygen injection can increase Eh while allowing time for microbial respiration and CO2 accumulation, maintaining a safe balance between the two indicators.

[0108] Furthermore, to achieve accurate prediction of future states, a system of kinetic differential equations based on physical driving forces is preferred as the core of the prediction model. For example, the rate of change of Eh, d[Eh] / dt, can be expressed as the difference between the oxidation rate term and the reduction rate term, i.e.:

[0109] d[Eh] / dt=k ox [O2]-k red [Reductant];

[0110] The rate of change of pCO2, d[logpCO2] / dt, describes the dynamic balance between respiration production, gas dissipation, and carbonate buffering, i.e.:

[0111] d[logpCO2] / dt=R resp -k evade ×(1-θ water )-k buffer ×CBI;

[0112] Where k ox [O2] is the oxidation rate constant; [O2] is the soil dissolved oxygen concentration; k red R is the reduction rate constant; [Reductant] is the total activity of soil reducing substances; resp k is the rate of gas production by rhizosphere microbial respiration. evade θ is the gas escape coefficient; water The water cover of the field surface (0 to 1); k buffer is the buffer consumption coefficient; CBI is the carbonate buffer strength index. Solving this system of equations through numerical integration can predict the environmental response trajectory under complex regulatory actions more accurately than a purely data-driven model. The Eh-pCO2 synergistic regulation model is suitable for sites with high carbonate content, strong buffering capacity, and significant karst fissure development.

[0113] Otherwise, it will enter a control mode based on model predictive control.

[0114] In this embodiment, the path selection logic is based on the carbonate buffer strength index (CBI). When the CBI is less than a preset buffer threshold, it means that the soil carbonate content is low, the pH value is no longer strongly locked, and the redox potential Eh and pH value regain a strong coupling relationship. Furthermore, the environment is more significantly affected by hydrological fluctuations. Therefore, the system automatically switches to model predictive control mode, utilizing the algorithm's optimization capabilities to cope with the dynamic changes of multivariate coupling.

[0115] It should be noted that the setting of the buffer threshold affects the selection of the control mode. If the threshold is set too high, plots that should use Eh-pCO2 synergistic control may mistakenly enter the MPC mode, thus failing to effectively address the carbonate buffering effect. If the threshold is set too low, some low-buffered plots may require overly complex dual-loop control, increasing unnecessary control costs. In practical applications, it is recommended to determine the threshold based on the statistical distribution of soil calcium carbonate content in the target area, typically ranging from 1.5 to 2.5.

[0116] like Figure 4 As shown, in one embodiment of this application, the control mode based on model predictive control includes:

[0117] Based on the fracture flux index, seasonal variation of groundwater level, and calcium carbonate saturation index in the coupled feature vector, an adaptive Eh feasible domain boundary function is constructed for the land parcel.

[0118] Specifically, the plot-adaptive Eh feasible region boundary function can solve the problem that traditional fixed thresholds cannot adapt to different geological backgrounds. The plot-adaptive Eh feasible region boundary function dynamically outputs the optimal control interval for the current plot by establishing a quantitative mapping relationship between geological parameters and the safe Eh boundary. The lower boundary function of the Eh feasible region (Eh...) low ) and upper boundary function (Eh high Eh can be represented as a linear combination of eigenvector components. For example, Eh low The calculation formula is:

[0119] Eh low =α0+α1×KFI+α2×ΔGW+α3×SIc;

[0120] Where KFI is the fracture flux index; ΔGW is the seasonal variation of groundwater level; SIc is the calcium carbonate saturation index; and α0, α1, α2, and α3 are regression weight coefficients obtained through training with historical data. A higher KFI (indicating well-developed fractures) means a greater risk of vertical leaching of pollutants, and coefficient α1 is usually positive, raising the lower boundary of Eh. This requires a weaker reducing environment to prevent the rapid reduction and dissolution of iron and manganese oxides, which could lead to the release of adsorbed heavy metals. If ΔGW is large (indicating large water level fluctuations), coefficient α2 is used to correct the boundary to reserve an oxidation buffer when the water level falls. Similarly, Eh... high The calculation formula is:

[0121] Eh high =β0+β1×KFI+β2×ΔGW+β3×SIc;

[0122] Where β0, β1, β2, and β3 are regression weight coefficients obtained through training with historical data. By substituting the feature values ​​of the current plot, such as KFI of 0.6, ΔGW of 2.0 meters, and SIc of 0.5, the current dynamic feasible region can be calculated to be [-80 mV, +120 mV]. The dynamic window is more in line with the actual safety requirements of the plot than the fixed window.

[0123] With the goal of minimizing the pre-defined cadmium-arsenic composite risk proxy value and control cost, the optimal control sequence in the future prediction time domain is solved under the condition of satisfying the boundary function constraint of the feasible region Eh of the land parcel.

[0124] In this embodiment, the core lies in constructing an optimization objective function J that includes risk and cost terms. The objective function can be expressed as:

[0125] J=Σ k=1 T [w1×R Cd (k)+w2×R As (k)+w3×||U(k)|| 2 ];

[0126] Where k is the time step size in the prediction time domain, for example, a step size of 1 hour for the next 24 hours; T is the total number of steps in the prediction control time domain; R Cd (k) and R As (k) represents the cadmium migration risk proxy value and the arsenic activation risk proxy value at time k, respectively; U(k) is the control vector at step k; ||U(k)|| 2 The norm square of the control vector represents energy consumption and reagent cost; w1, w2, and w3 are weighting coefficients. The solution process typically employs quadratic programming (QP) or nonlinear programming algorithms, solving the problem on a rolling basis within each control cycle, executing only the first control action in the sequence.

[0127] In one possible implementation, the site-adaptive Eh feasible region boundary function is constructed by establishing a linear or nonlinear combination relationship between the fracture flux index, the seasonal variation of groundwater level, the calcium carbonate saturation index, and the initial pH value. The site-adaptive Eh feasible region boundary function is used to output the lower boundary threshold and the upper boundary threshold of Eh that vary with the geological characteristics of the site, as state constraints in the process of solving the optimal control sequence.

[0128] Specifically, this means that when solving for the minimum value of the objective function J, the constraint condition Eh must be satisfied. low ≤Eh predicted (k)≤Eh high Eh predicted (k) represents the predicted soil redox potential Eh at time k. If the prediction model finds that the future trajectory of Eh may break through this boundary, the optimizer will force an adjustment to the control variable U(k), such as by increasing irrigation or adding chemicals, to pull the state back into the feasible region.

[0129] In one embodiment of this application, the optimal control sequence includes irrigation and drainage duration, sulfate dosage, and nitrate dosage;

[0130] Solving for the optimal control sequence in the future prediction time domain includes:

[0131] We introduce a cadmium migration risk proxy function and an arsenic activation risk proxy function based on the Sigmoid activation function;

[0132] The output values ​​of the cadmium migration risk proxy function and the arsenic activation risk proxy function being less than a preset safety threshold are used as gating constraints.

[0133] Gating constraints are satisfied by adjusting the dosage of sulfate and nitrate.

[0134] Specifically, the optimal control sequence is used to control the irrigation and drainage actuators and chemical dosing devices to perform corresponding actions. To address the nonlinear relationship between heavy metal risk and environmental parameters, a risk surrogate function is constructed. The cadmium migration risk surrogate function R... Cd It can be represented as:

[0135] R Cd =1 / (1+exp(-(c0+c1×ln(DGT Cd )+c2×Eh)));

[0136] Arsenic activation risk proxy function R As It can be represented as:

[0137] R As =1 / (1+exp(-(a0+a1×ln(DGT As )-a2×Eh)));

[0138] Where R Cd and R As These are the cadmium migration risk proxy value and the arsenic activation risk proxy value, respectively, ranging from 0 to 1; exp() represents the exponential function; ln() represents the natural logarithm function; DGT Cd and DGT As R0 represents the most recently measured bioavailable flux; Eh represents the predicted redox potential; c0, c1, c2 and a0, a1, a2 are model fitting parameters. The introduction of the Sigmoid function normalizes the risk value to between 0 and 1, simulating the abrupt change in the probability of risk occurrence. Gating constraints refer to the requirements for R0... Cd <γ and R As <γ; γ is a preset safety threshold, such as 0.2. When simply adjusting the moisture content cannot meet this constraint, the algorithm will automatically search for the optimal dosage of sulfate and nitrate.

[0139] Furthermore, the process of satisfying the gating constraints by adjusting the dosage of sulfate and nitrate follows the following branch logic:

[0140] When the real-time monitored redox potential is lower than the lower boundary threshold of the Eh feasible domain, it is determined to be a state of preferential control for arsenic activation, and the amount of nitrate added is preferentially increased as an alternative electron acceptor.

[0141] When the real-time monitored redox potential is higher than the upper boundary threshold of the feasible Eh domain, it is determined to be a cadmium migration priority control state, and the sulfate dosage is increased to promote cadmium sulfide precipitation.

[0142] When the real-time monitored redox potential is between the lower and upper thresholds of Eh, the sulfate and nitrate dosages are adjusted synchronously according to a preset ratio.

[0143] In this embodiment, the branching logic clearly defines the priority of electron acceptor regulation. When Eh falls below the lower boundary, the system faces a severe risk of arsenic activation, requiring oxidation to raise Eh. This is because nitrate (NO3) - The standard electrode potential of arsenic oxide is relatively high, making it an excellent alternative electron acceptor. Therefore, the system preferentially increases the dosage of nitrate to inhibit the reduction and dissolution of iron and manganese oxides, thereby locking in arsenic. Conversely, when Eh is above the upper boundary, the system faces the risk of cadmium activation. In this case, the dosage of sulfate is preferentially increased to utilize the sulfide ions (S₂O₃) generated by sulfate reduction. 2- This induces the formation of cadmium sulfide precipitates. When Eh is in the intermediate region, the system maintains a low-dose equilibrium of both to preserve the chemical stability of the system. The model-based predictive control mode is suitable for complex sites with strong hydrological heterogeneity, drastic groundwater level fluctuations, and difficulty in control through simple rules.

[0144] As an alternative implementation, for hardware environments with limited computing resources, the aforementioned Model Predictive Control (MPC) optimization process can be replaced by a simplified rule lookup table. For example, the optimal control actions under different combinations of KFI and ΔGW can be pre-calculated offline and stored as a lookup table. During real-time execution, the control command can be retrieved directly from the table based on the current state.

[0145] According to one aspect of this application, it also includes a step of in-situ evaluation of the repair effect, specifically:

[0146] The diffusion flux of bioavailable cadmium and bioavailable arsenic in soil was determined by using a diffusion gradient thin-film probe array pre-deployed in the rhizosphere and non-rhizosphere zones.

[0147] Specifically, the diffusion gradient film (DGT) probe array consists of multiple DGT samplers, which are embedded in the rhizosphere region with dense rice roots and the non-rhizosphere region without roots, respectively, for a duration of 24 to 48 hours. Based on Fick's first law of diffusion, DGT technology calculates the average flux of heavy metals in the soil solution by measuring the mass of heavy metals accumulated on the resin film. This method more accurately reflects the actual absorption potential of heavy metals by plants than traditional methods involving total soil mass or single extraction. In practice, after the DGT device is recovered, the metals in the resin layer are extracted using nitric acid solution, and quantitative analysis is performed using inductively coupled plasma mass spectrometry (ICP-MS).

[0148] In-situ scanning of soil profiles was performed using a portable micro-area X-ray fluorescence spectrometer to obtain two-dimensional spatial co-distribution maps of heavy metal elements and passivation components.

[0149] In this embodiment, micro-area X-ray fluorescence spectrometry (μ-XRF) scanning can reveal the microscopic binding of heavy metals and passivating materials. By attaching a portable spectrometer probe to the undisturbed soil profile, elemental distribution heatmaps with micrometer-level resolution can be generated. By overlaying and analyzing the spatial distribution correlations of cadmium, arsenic, iron, manganese, and lanthanum (characteristic elements of the passivating agent), it is possible to visually determine whether the passivating agent has successfully captured the target pollutants. Furthermore, as a preferred mechanism verification method, synchrotron X-ray absorption fine structure spectroscopy (SR-XAFS) analysis can also be used. By resolving the XAFS spectrum, newly formed cadmium-containing mineral phases in the soil can be identified, such as the detection of Cd-O-Fe or Cd-S coordination bonds, confirming the effectiveness of the iron-manganese oxide or sulfide precipitation mechanism at the molecular level. The diffusion flux and two-dimensional spatial co-distribution spectrum are used as remediation effect assessment data for subsequent closed-loop feedback updates.

[0150] In a further embodiment, a closed-loop feedback update step is also included, specifically:

[0151] The measured diffusion fluxes of bioavailable cadmium and bioavailable arsenic were compared with the preset remediation thresholds.

[0152] Based on the comparison results and the two-dimensional spatial co-distribution map, the weighting coefficients of the adaptive Eh feasible domain boundary function of the land parcel are updated online, as well as the parameters of the risk assessment model used to calculate the composite risk probability, so as to optimize the control strategy for the next remediation cycle.

[0153] Specifically, assuming the cadmium flux measured by DGT is still higher than the preset safety threshold, it indicates that the current control strategy is too weak. The system uses the deviation signal to correct the boundary function coefficients α online using the recursive least squares (RLS) method. i and β iFor example, if the actual cadmium release risk is found to be higher than the model prediction, the algorithm will automatically increase the weight coefficient α1 of KFI, thereby tightening the lower boundary of Eh and forcing a more conservative control strategy in the next cycle, favoring an oxidizing environment. Simultaneously, for the risk assessment model, an incremental learning algorithm can be used to add the feature vector-pollution result data pairs obtained in the current assessment cycle to the training set, fine-tuning the edge weight parameters of the graph convolutional network (GCN) to improve the model's prediction accuracy for that specific site. This achieves a closed-loop process based on actual measurement, bias, and correction. For example, the updated boundary function coefficient vector, i.e., α... i and β i It can be represented as:

[0154] θ new = θ old +K gain * (y actual - y predicted );

[0155] Where θ new The updated boundary function coefficient vector; θ old K is the coefficient vector of the previous period; gain The gain matrix can be calculated based on the covariance matrix; y actual The ideal Eh boundary value is obtained based on the DGT measured inversion; y predicted These are the boundary values ​​predicted by the model. This reflects the dynamic correction process of geological parameter weights using measured deviations.

[0156] This embodiment not only employs chemical probe technology but also combines spectroscopic methods to ensure the accuracy of the evaluation results and the interpretability of the mechanism.

[0157] It should be noted that after the aforementioned closed-loop feedback update, the corrected risk assessment model parameters will be automatically synchronized to the intelligent risk assessment module for use in calculating the composite risk probability in the next remediation cycle. Simultaneously, the updated Eh feasible domain boundary function coefficients will be written to the local storage of the intelligent control terminal to ensure that subsequent control strategies can be optimized based on the latest land parcel characteristics.

[0158] According to another aspect of this application, a modified biochar synergistic passivation and graded remediation system for cadmium and arsenic contaminated paddy fields is provided for performing the modified biochar synergistic passivation and graded remediation method for cadmium and arsenic contaminated paddy fields described in any of the above embodiments, comprising:

[0159] The multi-source sensing module is used to acquire multi-source heterogeneous environmental data of the paddy field area and generate coupled feature vectors.

[0160] Specifically, the multi-source sensing module includes wireless sensor network nodes deployed in the field. Each node is equipped with an industrial-grade pH / Eh electrode for monitoring Eh and pH, a non-dispersive infrared (NDIR) carbon dioxide sensor for monitoring rhizosphere pCO2, and a pressure-type water level gauge for monitoring surface and groundwater levels. In addition, the multi-source sensing module includes a remote sensing interface for periodically receiving multispectral image data transmitted from satellites or drones. All sensor data is transmitted to an edge computing gateway via LoRa or 5G narrowband IoT (NB-IoT), where a built-in data preprocessing algorithm cleans, aligns, and generates feature vectors.

[0161] The risk intelligence assessment module is used to run risk assessment models to calculate composite risk probabilities and generate risk management zoning instructions.

[0162] In this embodiment, the risk intelligent assessment module can be deployed on a cloud server or a high-performance edge computing box. A pre-trained graph convolutional neural network (GCN) model is loaded. After receiving the coupled feature vector, the risk intelligent assessment module quickly infers the composite risk probability of each grid cell and generates a risk management zoning instruction map in GeoJSON format based on a preset grading standard. This map includes the spatial coordinates of each plot and its corresponding risk level (high, medium, low).

[0163] The modified biochar customized supply unit includes multiple independent silos and variable application mechanisms, used to deliver remediation materials, including iron-manganese-lanthanum ternary co-loaded modified biochar, in a differentiated manner according to risk control zoning instructions.

[0164] Specifically, the customized modified biochar supply unit is typically mounted on an unmanned agricultural vehicle or tractor platform. It features three independent hoppers, respectively storing ferromanganese-lanthanum biochar for high-risk areas, phosphate-impregnated biochar for medium-risk areas, and alkaline-thermal microporous biochar for low-risk areas. A variable-rate dispensing mechanism precisely controls the rotation speed of the screw feeder at the bottom of each hopper via a servo motor. When the vehicle travels to different risk zones, the onboard controller reads GPS coordinates and matches zone commands, automatically adjusting the discharge rate of each hopper to achieve millisecond-level formula switching and precise dispensing.

[0165] The intelligent regulation and control terminal is used to connect the rhizosphere monitoring sensor network and irrigation and drainage execution mechanism, and executes multi-variable collaborative regulation strategies based on the boundary of the adaptive target regulation feasible domain of the plot.

[0166] In this embodiment, the intelligent control terminal is the core of the field control, typically employing a programmable logic controller (PLC) or an embedded industrial computer. It is hard-connected to solenoid valves, water pumps, and pesticide dosing devices (for adding organic acids, sulfates, or nitrate solutions) in the field. The terminal internally runs a control algorithm (state machine or MPC solver). It reads sensor network data in real time, calculates the deviation between the current state and the target feasible region boundary, and outputs pulse width modulation (PWM) signals to drive actuators, such as opening solenoid valve 1 for 30 minutes of irrigation, or starting metering pump 2 for 5 minutes to add calcium nitrate solution. The system in this embodiment adopts a modular design, integrating sensing, decision-making, execution, and feedback functions, enabling it to adapt to harsh field operating environments.

[0167] In a preferred embodiment, the intelligent control terminal is also equipped with an anomaly detection and fault tolerance module. When the rhizosphere monitoring sensor detects abnormal data, such as a sudden change in reading exceeding five times the normal range, the system automatically marks the sensor data as invalid and uses the weighted average of neighboring sensors as a substitute input. When the communication link is interrupted for more than a preset timeout threshold, the system switches to local autonomous mode, continues to execute the most recently effective control strategy, and automatically updates upon communication restoration. In the event of an extreme rainstorm warning, the system prioritizes preventative deep-water irrigation to ensure the stability of the rhizosphere environment.

[0168] In summary, the graded remediation method for cadmium and arsenic contaminated paddy fields using modified biochar for synergistic passivation includes: collecting multi-source environmental data to construct feature vectors coupled with hydrogeological parameters, and dividing risk control zones accordingly; applying iron-manganese-lanthanum ternary co-loaded modified biochar in differentiated manner to high-risk areas; constructing adaptive feasible domain boundaries for regulation targets during the critical growth stages of rice, and implementing multivariate synergistic regulation strategies. The system automatically selects the regulation mode based on geological characteristics: implementing Eh-pCO2 dual closed-loop regulation for high-carbonate plots, and implementing model predictive control (MPC) based on electron acceptor ratio for plots with strong hydrological variability.

[0169] In a detailed embodiment, a typical cadmium-arsenic co-contaminated paddy field in a certain region was selected as the experimental site. The soil pH in this area was 7.8, the calcium carbonate content was 45 g / kg, and the underground karst fissures were moderately developed. The experimental field was divided into three blocks: Block A adopted the graded synergistic remediation method described in this application; Block B adopted the existing lime passivation + traditional flooding remediation method; and Block C served as a blank control group. During the implementation in Block A, an adjacency matrix construction strategy based on fissure connectivity was adopted for the construction of the graph convolutional neural network (GCN) in the risk assessment model. Specifically, resistivity tomography (ERT) data was used to identify the underground fissure network. If there was a connected water-conducting fissure between two grid cells, the edge weight in the adjacency matrix was set to 1.0; if it was only connected by soil matrix, it was set to 0.1. This enabled the model to capture the spatial dependence of pollutant rapid migration along fissures.

[0170] After a full rice growing season, the remediation effects of each block were compared as follows: Regarding the compliance rate of heavy metal content in rice, Block A had an average Cd content of 0.12 mg / kg and an average As content of 0.15 mg / kg, achieving a compliance rate of 96% for both indicators. In contrast, although Block B met the Cd standard (0.18 mg / kg), the long-term flooding led to severe As activation, resulting in an average As content as high as 0.45 mg / kg (exceeding the standard), resulting in an overall compliance rate of only 40%. Regarding environmental stability control capabilities, monitoring data showed that during heavy rain, Block A successfully controlled the fluctuation range of rhizosphere Eh within ±30 mV through hydrological lag prediction and preventative regulation; while Block B, affected by fissure water intrusion, experienced Eh fluctuations exceeding ±150 mV, leading to the redissolution of precipitated Cd. Regarding material and water costs, thanks to the active sampling and tiered application strategy, the total amount of modified biochar added in Block A was reduced by 35% compared to the traditional uniform application model; thanks to MPC optimized control, irrigation water consumption was reduced by 45% compared to traditional flood irrigation. After continuous monitoring for three rice growing seasons, the compliance rates of heavy metals in rice in Block A were 96%, 94%, and 95%, respectively, indicating that the remediation effect of this application has long-term stability.

[0171] The experimental data above demonstrate that the synergistic passivation and graded remediation method of modified biochar for cadmium and arsenic contaminated paddy fields proposed in this application can effectively overcome the two major challenges of pH lock-in caused by high calcium and drastic hydrological changes caused by fissures in a certain area, achieving multiple benefits such as simultaneous attainment of cadmium and arsenic standards, reduced remediation costs, and water conservation.

[0172] This invention abandons the traditional single Eh control mode and adopts a dual-variable synergistic regulation strategy of Eh and pCO2. By introducing carbonate buffer strength diagnosis and utilizing precise organic acid adjustment, intermittent aeration, and water surface film, it successfully regulates the partial pressure of carbon dioxide in the rhizosphere under the premise that the pH is locked by the carbonate buffer, maintaining the thermodynamic stability of cadmium carbonate precipitation and effectively preventing cadmium reactivation at the safe potential. An adaptive target regulation feasible domain boundary is constructed for the site, replacing the rigid fixed threshold. Combining model predictive control algorithms and hydrological lag prediction mechanisms, the system can dynamically adjust control constraints according to the fracture flux index and initiate preventive intervention in advance before fracture water invades the rhizosphere, achieving proactive adaptation and stable control to drastic fluctuations in groundwater. Furthermore, the differentiated application of Fe-Mn-La ternary co-loaded modified biochar further solves the problem of decreased arsenic adsorption capacity of conventional materials under high pH conditions, achieving long-term treatment of cadmium and arsenic simultaneously.

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

Claims

1. A method for graded remediation of cadmium and arsenic contaminated paddy fields using modified biochar with synergistic passivation, characterized in that... include: Collect multi-source heterogeneous environmental data from paddy field areas and construct a coupled feature vector that combines pollution characteristics, hydrogeological parameters, and ecological indicators; The coupled feature vector is input into a pre-built risk assessment model to calculate the composite risk probability of the grid cell, and the paddy field area is divided into risk control zones of different levels accordingly. Based on the risk management zone level, differentiated modified biochar materials with matching functions are deployed, including composite supported biochar for high-risk areas; During the critical growth period of rice, the rhizosphere microenvironment was monitored in real time, and the boundary of the target regulation feasible domain for plots was constructed based on coupled feature vectors. Based on the boundary of the feasible domain of target regulation, a multivariate collaborative regulation strategy is implemented to maintain the rhizosphere microenvironment state within the safe domain of inhibiting cadmium and arsenic activation. Implementing multivariate collaborative regulation strategies, including: Obtain the carbonate buffer strength index in the coupled feature vector; The carbonate buffer strength index is compared with a preset buffer threshold. When the carbonate buffer strength index is greater than the buffer threshold, the system enters the Eh-pCO2 synergistic control mode; otherwise, it enters the model predictive control mode. Eh-pCO2 synergistic regulation modes include: Based on the real-time monitored rhizosphere redox potential, the preset Eh baseline, and the karst water redox potential characteristic value in the coupled feature vector, the karst water intrusion index is calculated. A dual closed-loop control architecture was constructed, including an Eh control loop triggered by the karst water intrusion index and a pCO2 control loop triggered by the carbonate buffer strength index. The coupling correction rule is executed to identify and coordinate the control conflict between the Eh control loop and the pCO2 control loop, and the final control command is generated.

2. The method according to claim 1, characterized in that, After calculating the composite risk probability of the grid cells, the following is also included: Based on the composite risk probability and the preset sampling cost factor, the control failure value index of each grid cell is calculated. The control failure value index is compared with the preset trigger threshold. When the control failure value index of any grid cell exceeds the trigger threshold, a supplementary sampling instruction is generated for that grid cell, and the coupling feature vector is updated using the supplementary sampling instruction.

3. The method according to claim 2, characterized in that, Control Failure Value Index I i Calculated using the following formula: I i =(P violate，i ·G i ·(1-C i )) / (Cost i +ε); Where P violate，i G is the probability of regulation exceeding the limit calculated based on a pre-configured prediction model. i C is a risk boundary gradient index calculated based on composite risk probabilities. i To predict the repair achievement confidence level, Cost i ε is a preset sampling cost factor, and ε is a positive number.

4. The method according to claim 1, characterized in that, For grid cells designated as high-risk areas, the composite-loaded biochar deployed was iron-manganese-lanthanum ternary co-loaded modified biochar. Iron-manganese-lanthanum ternary co-supported modified biochar was prepared by in-situ growth of iron-manganese oxide on the surface of a biochar support and doping with lanthanum, wherein the molar ratio of iron, manganese and lanthanum was 3:1:0.

5.

5. The method according to claim 1, characterized in that, The Eh control loop also includes an active prevention mechanism, specifically: Obtain meteorological and rainfall forecast data, and calculate the time delay parameters of karst fissure water intrusion into the rhizosphere by combining the fissure flux index in the coupled feature vector; Based on the time delay parameter, preventive control measures are initiated within a preset time window before the predicted time of fissure water intrusion, including increasing the depth of the surface water layer to use hydrostatic pressure to block the lateral infiltration of fissure water.

6. The method according to claim 1, characterized in that, The specific execution logic of the pCO2 control loop is as follows: Maintain the pCO2 of the rhizosphere soil within the target window of 0.003 atm to 0.03 atm; When pCO2 is below the lower limit of the target window, the partial pressure is increased by stimulating rhizosphere microbial respiration with the application of organic acid solution. When pCO2 is higher than the upper limit of the target window, the partial pressure is reduced by decreasing the thickness of the water layer on the field surface to promote gas dissipation.

7. The method according to claim 1, characterized in that, Model-based predictive control (MMDC) modes of regulation include: Based on the fracture flux index, seasonal variation of groundwater level, and calcium carbonate saturation index in the coupled feature vector, an adaptive Eh feasible domain boundary function for the land parcel is constructed. With the goal of minimizing the pre-defined cadmium-arsenic composite risk proxy value and control cost, the optimal control sequence in the future prediction time domain is solved under the condition of satisfying the boundary function constraint of the feasible region Eh of the land parcel.

8. The method according to claim 7, characterized in that, The optimal control sequence includes irrigation and drainage duration, sulfate dosage, and nitrate dosage; Solving for the optimal control sequence in the future prediction time domain includes: We introduce a cadmium migration risk proxy function and an arsenic activation risk proxy function based on the Sigmoid activation function; The output values ​​of the cadmium migration risk proxy function and the arsenic activation risk proxy function being less than a preset safety threshold are used as gating constraints. Gating constraints are satisfied by adjusting the dosage of sulfate and nitrate.

Images