Agricultural non-point source pollution prediction and tracing method

Abstract

The application discloses a kind of agricultural non-point source pollution prediction and tracing method, it is related to agricultural monitoring technical field, to gravity standard triangulation network as skeleton, time-lapse gravity residual is converted to grid virtual deformation, and determines edge level priority channel and convergence / divergence;In the same edge-node space coupling water / soil / air and cross-phase parameter.Adopt programmable expert sequence aggregation: gravity deformation, three-phase flow field, interannual variation, microbial community memory index, groundwater lag and observation consistency, and by forecast hook, rainfall, wind, temperature is mapped as parameter increment, and rolling update.Method is carried out mass conservation and threshold rejection at edge level, output main diffusion ridge line, sensitive section p10 / p50 / p90 arrival time and source area contribution rate, realize the explainable prediction, verification and accurate tracing of non-point source pollution.

Description

Technical Field

[0001] This invention relates to the field of agricultural monitoring technology, specifically a method for predicting and tracing agricultural non-point source pollution. Background Technology

[0002] Agricultural non-point source pollution is characterized by spatial dispersion, sudden temporal occurrence, and concealed transport pathways. It is often influenced by rainfall / irrigation, surface roughness, soil permeability, near-surface wind, and surface-to-subsurface coupling. Pollutants are easily diluted, adsorbed, and exchanged during migration, resulting in weakened chemical signals downstream, making prediction and source tracing difficult. Current engineering practices largely rely on distributed hydrological models based on DEM flow direction and empirical parameters, or use chemical / particle tracing and microbial source tracking (MST / eDNA) for validation. Separate groundwater numerical models are also used for long-tail effect assessment. However, these approaches still have shortcomings in the unified representation of multi-source heterogeneous data, the closure of evidence chains, and rapid updates.

[0003] Analysis of existing technologies reveals the following: First, flow direction determination based primarily on topographic slope is prone to path ambiguity in flat areas and areas with dense artificial ditches, lacking independent constraints that can reveal hidden priority channels; second, model parameters are mostly static or annual averages, making it difficult to respond to interannual differences and short-term weather processes, resulting in delayed prediction updates; third, the line-area structure of GIS raster fields and river / road / field plots makes it difficult to uniformly solve flux, dilution, and arrival time at the "edge" scale, making it difficult to stably guarantee mass conservation and non-negativity; fourth, microbial evidence is mostly used for identifying pollution types or host sources, and has not yet been effectively coupled with the consistency of physical paths and timelines ("functional memory" and decay rhythms); fifth, the time-lag linkage between surface runoff and groundwater is insufficiently characterized, leading to an underestimation of the risk of late arrival in distant areas and riverbank zones; sixth, multi-source observations and model results lack weight gating and threshold removal within the same state space, resulting in insufficient interpretability and stability of the output main diffusion paths, arrival time windows, and source area contribution rates.

[0004] Existing technologies predict and validate runoff based on DEM-based runoff confluence analysis, empirical / physical hybrid distributed hydrological models, chemical or particle tracer experiments, microbial source tracing, and independent groundwater simulation. However, they still have certain limitations (data and models are not coupled on the same spatial object, there is a lack of independent physical constraints that can indicate convergence / divergence and main direction, there is a lack of explicit mapping of interannual variation and short-term forecasts, and there is no consistent gating between the evidence layer and the dynamic layer), such as unstable path orientation, large arrival time error, high source area contribution rate deviation, difficulty in closing the mass conservation loop, insufficient identification of groundwater wake, and difficulty in forming an auditable chain of evidence.

[0005] Therefore, there is an urgent need for a method that can use time-shifted gravity information for priority channel identification in a "gravity standard triangulation network" within a unified edge-node state space, establish a "edge-level physical parameter coupling system", use "programmable expert sequential aggregation" to fuse three-phase flow fields, interannual variations and "forecasting linkage", and introduce "microbial community memory index", groundwater lag and observation consistency for calibration. Finally, under the mass conservation and threshold elimination rules, the method outputs the main diffusion ridge, arrival time window and source area contribution rate for agricultural non-point source pollution prediction and tracing. Summary of the Invention

[0006] The purpose of this invention is to overcome the shortcomings of the prior art and propose a method for predicting and tracing agricultural non-point source pollution to solve the above-mentioned problems.

[0007] The objective of this invention is achieved through the following technical solution: a method for predicting and tracing agricultural non-point source pollution, comprising the following steps:

[0008] S1. Constructing a standard gravity triangulation network: Within the agricultural area, topographic information is acquired based on a high-resolution digital elevation model. A constrained Delaunay triangulation network is established, using canals, fields, roads, watersheds, and riverbanks as boundary constraints. The network must meet quality control requirements, with a minimum interior angle of 30 degrees and an aspect ratio not exceeding 3. Time-lapsed gravity observation data is loaded onto the nodes and edges of the triangulation network. Gravity residuals are calculated and converted into virtual deformation information for the grid, determining the priority channel weights and convergence / divergence characteristics of each edge.

[0009] S2. Establish a border-level physical parameter coupling system: Establish a unified physical parameter package for each edge of the triangular network. The parameter package includes water phase transport parameters, soil phase transport parameters, gas phase transport parameters, and interphase exchange parameters. Among them, water phase transport parameters include roughness coefficient, channel capacity, head loss coefficient, inflow bias weight, sediment carrying capacity, and dilution mixing length; soil phase transport parameters include boundary permeability, permeability adjustment factor, water storage coupling strength, and adsorption label; gas phase transport parameters include effective roughness length, ventilation corridor coefficient, shear shielding coefficient, and settlement bias; and interphase exchange parameters include interphase exchange weight and reactivity label.

[0010] S3. Execute programmable expert sequential aggregation: On the same side node state space of the triangular network, six types of programmable expert modules are called sequentially in a predetermined order to update the edge state, including:

[0011] Gravity Deformation Expert: Generates gravity weights, inflow biases, and priority channel sequences based on gravity residuals and virtual deformation information; Three-Phase Flow Field Expert: Calculates water, soil, and air flux distribution at each edge based on shallow water runoff, soil seepage, and near-surface wind, and provides preliminary arrival times; Interannual Variation Expert: Corrects edge parameters year by year based on multi-year historical meteorological data and multi-temporal satellite remote sensing images; Microbial Functional Memory Expert: Implements path consistency gating and timeline verification for candidate paths based on the Microbial Community Memory Index (SMI); Groundwater Lag Expert: Establishes a delay queue for surface-subsurface composite transport at the boundary between the riparian zone and the groundwater area; Observational Consistency Expert: Calibrates edge-level confidence based on measured flow velocity, water level, pollutant concentration, and groundwater level data;

[0012] S4. Integrated meteorological forecast-driven prediction: Short-term meteorological forecast data is projected onto side parameters through a forecast-linking mechanism. The short-term meteorological forecast data includes hourly or three-hourly rainfall, wind direction and speed, and temperature. Rainfall is mapped to the increment of inflow bias, wind direction and speed are mapped to the adjustment of ventilation corridor coefficient, and temperature is mapped to the correction of interphase exchange or reaction activity, so as to advance the prediction of pollution diffusion path and arrival time in future periods.

[0013] S5. Multi-expert aggregation and result output: On each edge, the outputs of six types of experts are aggregated and adjudicated according to weight gating and threshold elimination rules to generate a single edge flux, pollutant mass distribution, arrival time probability distribution and confidence score; at the edge level, mass conservation and non-negativity constraints are implemented, and flux is preferentially deducted for edges with low confidence; at the node level, the edge-level results are summarized and the main diffusion ridge of agricultural non-point source pollution, the arrival time window of sensitive sections and the source area contribution rate are output.

[0014] In the process of constructing the gravity standard triangulation network in step S1, a triangle quality control strategy is adopted to ensure that the minimum interior angle of the triangular unit is not less than 30 degrees and the aspect ratio does not exceed 3. Watersheds, canals, field boundaries, roads, and riverbanks are used as constraint lines. At the same time, the edges are semantically classified and labeled as ditch channels, slope confluence, field boundaries, road drainage, riverbank exchange, or ventilation corridors, which are used to drive the subsequent differential setting of three-phase parameters and weights.

[0015] In the gravity residual transformation process of step S1, by calculating the changes in the area and angle of the triangular unit, the main elongation direction is determined as the dominant direction of pollution diffusion, the area with reduced area is marked as the convergence zone, and the area with increased area is marked as the divergence zone, forming a priority channel sequence at the edge level.

[0016] In the interannual variation expert module of step S3, by analyzing historical meteorological records of more than ten years and multi-temporal satellite remote sensing images, interannual variation characteristics such as rainfall pattern spectrum, vegetation cover index NDVI, bare land rate, soil moisture and surface water frequency are extracted, and the boundary parameters are corrected year by year.

[0017] In step S3, the microbial functional memory expert module collects soil and water microbial samples along the main diffusion ridge, bifurcation point, riparian zone and key cross-section. The microbial community memory index (SMI) is constructed by 16S gene sequencing and metagenomic analysis, and the path consistency is verified by the spatial distribution gradient of functional genes and enzyme activities.

[0018] In the groundwater lag expert module in step S3, the coupling relationship between the surface and groundwater is established at the riverbank zone, alluvial fan boundary and shallow groundwater observation point. The groundwater lag transport process is simulated through a delayed queue mechanism to supplement the long-term pollution risk assessment.

[0019] In the observation consistency expert module of step S3, flow velocity monitoring, water level observation, pollutant concentration detection and groundwater level data are collected in real time and compared with the edge-level prediction results. The confidence weight of edges that are inconsistent with the observation data is reduced to ensure the quality conservation constraint.

[0020] In the meteorological forecast fusion in step S4, the forecast hooking mechanism includes explicit mapping rules: mapping the forecast rainfall to a weighted increment of the edge inflow bias, mapping the forecast wind direction and speed to an adjustment of the ventilation corridor coefficient and shear shielding coefficient, mapping the forecast temperature to an adjustment of the interphase exchange weight or reactive label, and limiting the above mappings to maintain edge mass conservation within the time step.

[0021] In the multi-expert aggregation process in step S5, the outputs of the six types of experts are superimposed and aggregated through weight gating, threshold elimination and normalization to generate a comprehensive confidence score of the edge, and the probability arrival time window of p10, p50 and p90 is output at the sensitive section.

[0022] The source region contribution rate decomposition is calculated by aggregating, normalizing, and aligning the edge fluxes from different fields or management units at the target cross section, and simultaneously outputs the main diffusion ridge map, high-risk convergence zone identification, and composite path uncertainty assessment; the source region contribution rate is used for governance priority ranking and source tracing decision support.

[0023] The beneficial effects of this invention are:

[0024] This invention constructs a "gravity standard triangulation network," converting time-shifted gravity observation data into grid virtual deformation information within a constrained Delaunay triangulation network. This allows for the determination of priority channel weights and convergence / divergence characteristics, enabling quantitative identification of hidden confluence channels and dominant diffusion directions, overcoming the ambiguity caused by relying solely on topographic slope aspect. Furthermore, by establishing a "border-level physical parameter coupling system," water phase transport parameters, soil phase transport parameters, gas phase transport parameters, and interphase exchange parameters are uniformly configured within the same edge-node state space. This allows flux, dilution, and arrival time solutions to be completed directly at the border level, achieving path-level physical consistency and controllable closure in conjunction with mass conservation and non-negativity constraints.

[0025] This invention employs "programmable expert sequential aggregation," where a gravity deformation expert provides the framework and inflow bias, a three-phase flow field expert generates initial values ​​for water, soil, and air fluxes and arrival times, an interannual variation expert corrects the boundary parameters year by year, a microbial functional memory expert performs path consistency gating and timeline verification using the microbial community memory index, a groundwater lag expert characterizes the wake effect through the delay queue of surface-subsurface composite transmission, and an observation consistency expert calibrates the boundary-level confidence based on measured flow velocity, water level, pollutant concentration, and groundwater level. The above sequential aggregation generates a single-sided state under weighted gating and threshold elimination rules, effectively improving the robustness of orientation and the reliability of arrival times under conditions of chemical signal dilution, observation sparsity, and multiple solutions.

[0026] This invention projects short-term weather forecast data onto side parameters through a "forecast-linked mechanism," explicitly mapping rainfall to inflow bias increments, wind direction and speed to ventilation corridor coefficient adjustments, and temperature to interphase exchange or reactivity corrections. This enables rolling forecasts and rapid updates before and during events, significantly shortening the response chain from weather changes to pollution pathway rearrangement. Finally, at the node level, side-level results are aggregated, outputting the arrival time windows and source area contribution rates of the main diffusion ridge and sensitive sections. It provides channel-level interpretable evidence decomposition and governance priority ranking, offering an integrated and implementable technical solution for the prediction, source tracing, and precise intervention of agricultural non-point source pollution. Attached Figure Description

[0027] Figure 1 The process of this invention Figure 1 ;

[0028] Figure 2 The process of this invention Figure 2 ;

[0029] Figure 3 The process of this invention Figure 3 . Detailed Implementation

[0030] The technical solution of the present invention will be clearly and completely described below with reference to the embodiments. 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 are within the scope of protection of the present invention.

[0031] It should be noted that the directional concepts of "left", "right", "up", "down", "front", "back", "inner", and "outer" in the following scheme are all relative directions, and will not be listed one by one here.

[0032] Example 1:

[0033] like Figure 1 As shown, this embodiment targets agricultural areas with typical complex landforms of ditches, fields, roads, and riverbanks. It utilizes a high-resolution digital elevation model, vector data of watersheds and riverbanks from ditches, fields, roads, and two periods of time-shifted gravity observation data to establish a constrained Delaunay triangulation as a unified spatial framework. Within the state space of nodes on the same side of this triangulation, edge-level physical parameter packages are configured. Edge fluxes, pollutant mass distributions, arrival time probability distributions, and confidence scores are generated through sequential aggregation by gravity deformation experts and three-phase flow field experts. Subsequently, mass conservation and non-negativity constraints are applied at the edge level and summarized at the node level to obtain the arrival time windows and source region contribution rates of the main diffusion ridge and sensitive sections. The steps, data structures, and key computational relationships disclosed in this embodiment are sufficient for those skilled in the art to implement.

[0034] Before implementation, a high-resolution digital elevation model (DEM) of 1 to 20 meters was acquired. Rivers, fields, roads, watersheds, and riverbanks were extracted as constraint lines. A constrained Delaunay triangulation was used to divide the region into triangulation networks. Triangle quality control was implemented to ensure that the minimum interior angle was not less than 30 degrees and the aspect ratio did not exceed 3. Semantic classification was performed on the generated edge set, labeling them as ditch channels, slope confluences, field boundaries, road drainage, riverbank exchange, or ventilation corridors, facilitating subsequent differentiated parameter configuration. Subsequently, two periods of time-shifted gravity surveys were conducted 48 to 72 hours before and after heavy rain or agricultural irrigation events. Combined with base station corrections, solid tide corrections, and terrain corrections, a regional gravity residual grid was obtained, which was then interpolated or extrapolated to the grid nodes. Within this node domain, to convert gravity residual information into grid virtual deformation information, node displacements were constructed according to the following formula, and unit deformation and convergence / divergence characteristics were calculated accordingly: ;

[0035] in, Let be the virtual displacement vector of node p; is the gravity residual scalar field at node p processed by the smoothing operator (obtained from the observation residual field through spatial smoothing); This represents the spatial gradient operator over the node domain; This is the displacement scaling factor, used to map the dimensions of the gravity residual to the deformation scale. After updating the three nodes of the triangle using the above displacement, the deformation gradient is used to measure the scaling and orientation information of the unit:

[0036] ;

[0037] in, Let T be the deformation gradient tensor of the triangular unit T; The initial node coordinates of the element; These are the node coordinates after virtual deformation; The unit area ratio (Jacobi). Representatives gathered, This represents divergence. Further, the principal elongation direction is obtained using the right Cauchy-Green tensor to indicate the dominant direction of potential diffusion: ;

[0038] in, Let T be the right Cauchy–Green tensor of cell T; The unit eigenvector corresponding to the largest eigenvalue; The largest eigenvalue. Project along the outer edge of the element, and combine Assigning direction and convergence degree to adjacent edges enables the determination of priority channel weights.

[0039] In the initialization phase of the edge-level physical parameter coupling system, a unified parameter package is established for each edge. Water phase transport parameters include roughness coefficient, channel capacity, head loss coefficient, inflow bias weight, sediment carrying capacity, and dilution-mixing length; soil phase transport parameters include boundary permeability, permeability adjustment factor, water storage coupling strength, and adsorption label; gas phase transport parameters include effective roughness length, ventilation corridor coefficient, shear shielding coefficient, and settlement bias; and interphase exchange parameters include interphase exchange weight and reactivity label. Initial parameter values ​​are retrieved from databases of land types, soils, channel cross-sections, road materials, vegetation, and remote sensing textures, and differentiated according to edge semantic classification. Subsequently, gravity deformation experts provide a "skeleton orientation / priority order" on the edge domain, synthesizing deformation measurement, slope aspect consistency, and semantic bias into edge-priority channel weights.

[0040] ;

[0041] in, The priority channel weight for edge e (between 0 and 1). It is an S-shaped compression function; The angle between the side direction and the main extension direction of the adjacent unit; The Jacobian of the unit adjacent to edge e; The category bias assigned based on the semantic classification of the edge (ditch passage, slope confluence, field boundary, road drainage, riverbank exchange or ventilation corridor); This is a weighted composite coefficient used to balance the contributions of directional consistency, convergence, and semantic bias. This weight is used for subsequent three-phase flux allocation and preliminary time-of-arrival estimation.

[0042] In the three-phase flow field expert stage, flux allocation and time-scale estimation are performed with the edges as first-class citizens. For the aqueous flux, a joint characterization of edge channel capacity, roughness coefficient, and inflow bias weight is used, and the edge aqueous volumetric flux is calculated using the following formula: ;

[0043] in, Let e ​​be the water phase volume flux. This is the equivalent coefficient for the side channel capacity; This is an estimate of the water surface potential difference between the two nodes at the two ends of the edge; Indicates the non-negative part; This is the equivalent frictional resistance synthesized from the roughness coefficient and the local head loss coefficient. For both the soil and gas phases, boundary permeability, permeability adjustment factor, effective roughness length, and ventilation corridor coefficient are used as the main control parameters to calculate the soil phase volume flux. Gas phase volume flux And synthesize them according to the same weight-resistance mechanism based on information completeness. The three-phase volumetric flux is converted into the side-average phase velocity and the initial arrival time is given:

[0044] ;

[0045] in, Estimate the median arrival time of edge e; The length of the side; For edge equivalent phase velocity; The equivalent cross-sectional area used for conversion (determined by edge channel capacity and semantic category). When cross-phase exchange of pollutants is significant, the mass distribution between phases is performed within the edge according to the cross-phase exchange parameters to ensure consistency in subsequent concentration distribution calculations.

[0046] To ensure that the edge flux and contaminant mass distribution satisfy physical constraints, mass conservation and nonnegativity constraints are enforced at the edge and node levels, respectively. For each node n, the following constraint is enforced within one time step:

[0047] ;

[0048] in, and Sum the sets of edges flowing into and out of node n, respectively; This represents the volumetric flux of the corresponding phase; This represents the net source or sink term for the node (or zero if unknown). When imbalances occur, flux deduction is prioritized for edges with low confidence to restore conservation. Pollutant mass dilution-transport within edges is updated based on edge mixing length and flux. ;

[0049] in, and These represent the inflow and outflow pollutant quality of edge e within the time step; The equivalent dilution factor is determined by the dilution mixing length, edge volume flux, and adsorption label. By each side By performing additive superposition along the path, the median arrival time estimate from the source node to any target section is obtained, and a probability arrival time distribution is constructed based on the edge arrival time uncertainty and parameter interval propagation.

[0050] During the aggregation and output phase, the multi-source computational load of each edge is weighted and thresholded to form a single edge state, including edge flux, pollutant mass distribution, arrival time probability distribution, and confidence score. Confidence increases with the intensity of gravity residuals, deformation consistency, and the completeness of the three-phase parameters, and decreases with the increase of semantic conflicts and conservation corrections. At the node level, the edge results are aggregated, and the main diffusion ridge is extracted along the path set with the highest confidence, continuous flux, and monotonic arrival time. At sensitive cross-sections, the arrival time window and the edge flux aggregation ratio of each upstream field or management unit at that cross-section are output as the source area contribution rate. This embodiment, without enabling interannual variation, weather forecast, microbial functional memory, and groundwater hysteresis modules, can already provide physically consistent and interpretable main directions and arrival times of pollution diffusion based on virtual deformation derived from gravity residuals and edge flux driven by three-phase parameters, laying the framework for the subsequent superposition of evidence layers and time-varying layers.

[0051] The working process of this embodiment is as follows: First, in the constrained Delaunay triangulation, the time-shifted gravity residual is mapped to the nodes. Based on the displacement-deformation relationship, virtual deformation of the elements is generated, and the main elongation direction and Jacobian are extracted. Then, priority channel weights are assigned to the edges. Subsequently, the edge-level physical parameter package is loaded according to the edge semantic classification. The gravity deformation expert is called to confirm the skeleton direction and priority order. Then, the three-phase flow field expert allocates the volume flux of water phase, soil phase, and gas phase on the edges and estimates the median arrival time of the edges. Next, mass conservation and non-negativity constraints are performed at the edge level to eliminate the imbalance caused by incomplete parameters or local instability. At the node level, the constrained edge states are summarized, and the output includes a set containing a single main diffusion ridge, multiple suboptimal paths, and corresponding arrival time windows. The contribution rate of the source region of sensitive sections is calculated. Through two phases of gravity surveys and several basic hydrological points for a single event, a stable macroscopic diffusion pattern can be obtained without relying on chemical concentration and biological evidence.

[0052] First, the gravity residual, after being quantified by virtual deformation and deformation gradient, directly serves to determine the weights of edge-priority channels, significantly improving the ability to identify hidden channels compared to relying solely on topographic slope aspect to determine flow direction. Second, the edge-level physical parameter package integrates flux, dilution, and arrival time calculations in the spatial topology, avoiding the secondary assumptions from nodal fields to line fluxes and improving the controllability of mass conservation. Third, the sequentially aggregated expert chain first determines the skeleton and then allocates fluxes, reducing path oscillations caused by parameter indiscernibility, resulting in more stable and interpretable main diffusion ridges. Fourth, the conservation strategy of deducting fluxes from low-confidence edge-priority channels maintains cross-sectional flow closure even when local observations are scarce or parameters are incomplete, ensuring the engineering usability of the results. Overall, this embodiment achieves closed-loop modeling from skeleton to dynamics using the minimum available link of "gravity virtual deformation edge-priority channels" and "edge-level three-phase fluxes," providing a reusable unified edge node state space foundation for further overlaying of interannual variation, weather forecasting, microbial functional memory, and groundwater hysteresis modules.

[0053] Example 2:

[0054] like Figure 1 and Figure 2 As shown, based on the same experimental area and mesh framework, the constrained Delaunay triangulation, gravity residual mapping, and edge semantic classification, consistent with the previous embodiment, are used to continue time-varying modeling and prediction with edges as first-class citizens. This embodiment adds two new parts: an interannual variation network and a forecast linkage mechanism. The goal is to adaptively update the edge-level physical parameter package with the year and future weather fields while maintaining edge-level mass conservation and non-negativity constraints, thereby obtaining a time-adaptive main pollution diffusion path and arrival time window, and maintaining consistency with the priority channel weights driven by gravity virtual deformation. To avoid repetition, the mesh construction, virtual deformation, and basic three-phase flux calculation disclosed in the previous embodiment are omitted here; only the specific implementation method, working process, and core calculation relationships of the interannual variation expert and forecast linkage mechanism are fully disclosed.

[0055] First, an interannual variation network is constructed, selecting historical meteorological and multi-temporal satellite remote sensing data of at least ten years. For each year, interannual variation characteristics such as rainfall pattern, vegetation cover index, bare land rate, soil moisture, and surface water frequency are generated. The sampling domain is the geometric range of the edges or the buffer zone at the midpoint of the edges. Edge-level interannual feature vectors are obtained through spatial resampling and temporal aggregation. Based on this, the interannually varying components in the edge-level physical parameter package (permeability adjustment factor, roughness coefficient, interphase exchange weight, ventilation corridor coefficient, and shear shading coefficient, etc.) are corrected year by year. To avoid abnormal parameter drift, a gated vectorized correction model is used for each edge. ;

[0056] in, Let e ​​be the parameter vector of year y (components include permeability adjustment factor, roughness coefficient, phase exchange weight, ventilation corridor coefficient, shear shielding coefficient, etc.). The parameter vector for the base year (or the static initialization of the previous embodiment); This is a mapping matrix obtained from offline learning of historical events; Let e ​​be the interannual feature vector of year y (components include standardized rainfall pattern index, NDVI percentile, bare land rate, soil moisture and surface water frequency, etc.). This is a trimming operator that limits parameters to the minimum / maximum physical range allowed by the operating conditions on a component basis. The meanings of the above symbols are as described above, ensuring the interpretability and controllability of parameter correction.

[0057] To further suppress erroneous corrections in years with weak signals, an adaptive gating system is superimposed on the parameter vector layer, allowing features from high-variance or extreme years to gain greater influence.

[0058] ;

[0059] in, These are the parameters that take effect after gating; The gating factor; It is an S-shaped compression function; Represents the variance of the interannual eigenvectors (after scale normalization); This is a gated overparameter. This gate ensures that the parameter returns to the baseline in stable years and appropriately increases the correction amount in abnormal years.

[0060] In the forecast-linked mechanism, hourly or three-hourly rainfall, wind direction and speed, and temperature fields are projected onto the boundary parameter package, forming a time-step driven dynamic increment. The inflow bias weights are mapped using direct rainfall increment mapping. ;

[0061] in, The inflow bias increment for edge e at time step t; The rainfall-confluence sensitivity coefficient (different values ​​can be set according to edge semantic classification, such as confluence of ditch channels above slopes). The predicted rainfall intensity is projected onto edge e. A coupled mapping of wind direction alignment factor and dimensionless wind speed is applied to the ventilation corridor coefficient and shear shielding coefficient:

[0062] ;

[0063] in, and These are the time step increments for the ventilation corridor coefficient and the shear shielding coefficient, respectively. For wind-induced mapping coefficients; This is the prevailing wind direction at that time step; Let e ​​be the direction angle of side e; This refers to the wind speed at that time step. Normalized reference wind speed; m, n are empirical indices. Temperature linear correction is applied to the interphase exchange weights and reactivity labels. ;

[0064] in, The weight increment for the interphase exchange of edge e at time step t; This refers to the temperature sensitivity coefficient. The temperature at that time step; The reference temperature is used. The above three mappings satisfy the time step constraint of edge-level mass conservation. Specifically, the estimated volume flux of the edge is first calculated using the updated parameters, and then conservation calibration is performed at the nodes.

[0065] At each time step, based on the three-phase flux formula of the previous embodiment, the edge-level volumetric flux and equivalent phase velocity are updated with interannually corrected parameters and forecast increments, and then the median time of arrival estimate of the path is calculated. For any set of directed paths pp from the upstream source region to the sensitive section, the time of arrival estimate at time step tt is: ;

[0066] in, Let p be the midpoint arrival time of path p at time step t; Let e ​​be the length of edge e; Let e ​​be the equivalent phase velocity of the edge, synthesized from the fluxes of the aqueous, soil, and gas phases. To ensure that the edge-level mass closes within the time step, a local minimum intervention calibration optimization is performed on each node:

[0067] ;

[0068] in, This is the calibrated edge volume flux; The estimated flux is obtained directly from the updated parameters; Adjust the weights for flux (take smaller values ​​for low-confidence edges so they are adjusted first); Let n be the set of edges connected to node n. The set of edges flowing into / out of the node; For external source and sink terms of the nodes. This convex optimization problem is solved independently at each node, ensuring the conservation and non-negativity of edge mass within the time step, while approximating the physical flux obtained by forecasting as closely as possible.

[0069] During operation, the system first performs batch corrections and gating of the edge-level physical parameter packages for each year using an interannual variation network, ensuring that the permeability adjustment factor, roughness coefficient, interphase exchange weight, ventilation corridor coefficient, and shear shielding coefficient reflect the long-term background differences of that year. Then, it enters the rolling forecasting phase. At each forecast time step, a forecasting hook mechanism maps the increments in rainfall, wind field, and temperature to inflow bias increments, ventilation corridor / shear shielding adjustments, and interphase exchange weight corrections, respectively, obtaining the edge-level dynamic parameters within the time step. Based on this, the edge volume flux and equivalent phase velocity are calculated, the median arrival time of the path is given, and mass conservation and non-negativity constraints are enforced through node local optimization. The main diffusion ridge is updated based on edge confidence, flux continuity, and arrival time monotonicity, and new arrival time windows are generated at sensitive sections. The entire process is completed within the same edge node state space, without disrupting the gravity virtual deformation and priority channel weight system established in the previous embodiment.

[0070] By employing vectorized correction and gating mechanisms in the interannual variation network, the edge-level physical parameter package reflects the systematic differences in meteorological and surface conditions over multiple years, avoiding systematic biases caused by extrapolating empirical parameters from a single year to other years, thereby improving the robustness of arrival time and flux allocation. On the other hand, through explicit mapping rules in the forecast-linked mechanism, rainfall, wind direction, wind speed, and temperature are projected in real-time as inflow bias increments, ventilation corridor / shearing adjustments, and interphase exchange weight corrections. This enables real-time updates of edge-level flux and arrival time estimates before and during events, significantly shortening the link from forecast to warning. Regarding mass conservation, local minimum intervention optimization at nodes forces closed budgets within the time step, while edge confidence weights prioritize adjustments to low-confidence edges, ensuring the consistency and interpretability of the overall output. Compared to approaches that rely solely on static parameters, this embodiment maintains the continuity of path morphology and time windows across different years and weather events. It can rapidly rearrange the main diffusion ridge and arrival time sequence in scenarios such as strong convection or sudden wind changes, and provides a reliable temporal basis for subsequent superposition of microbial functional memory, groundwater lag, and observation consistency.

[0071] Example 3:

[0072] like Figures 1 to 3As shown, based on the constrained Delaunay triangulation, gravity residual-generated grid virtual deformation information, edge semantic classification, and edge-level physical parameter package shared with the aforementioned embodiments, microbial functional memory experts, groundwater hysteresis experts, and observation consistency experts are introduced to perform evidence fusion and time-delay expansion on edge-level flux, arrival time probability distribution, and confidence scores within the same edge node state space. The grid construction, gravity deformation-driven priority channel weight determination, and basic calculations of the three-phase flow field in this embodiment follow the existing process and will not be repeated here. The key points disclosed are the specific implementation methods and core calculation relationships for the construction of the microbial community memory index and path consistency gating, the delay queue mechanism for surface-subsurface composite transport, the calibration and update of edge-level confidence scores based on observation data, the calculation of arrival time quantiles, and the output of source region contribution rates. Furthermore, mass conservation and non-negativity constraints are continuously executed at the edge level, and flux is preferentially deducted from low-confidence edges to obtain the arrival time windows (p10 / p50 / p90) of the main diffusion ridge and sensitive sections, the source region contribution rate, and the composite path uncertainty.

[0073] In the implementation of microbial functional memory, sampling points were set up along the main diffusion ridge, bifurcation points, riverbank zones, and key cross-sections. 16S gene sequencing and metagenomic functional annotation were performed on water and soil samples to form a standardized composition vector with functional pathways as the dimension. Using edges as first-class citizens, the functional pathway distribution of each edge within the sampling period was represented as follows: The distribution of functional pathways in regions without fertilization references or upstream areas was selected. The boundary-level microbial community memory index is defined as...

[0074] ;

[0075] in, The microbial community memory index for time step t; The weights for functional pathway k are set based on the pathway's sensitivity to the target pollutant, satisfying the following conditions: ; A symmetric divergence metric between 0 and 1, used to measure the difference between the current and reference functional compositions. For path consistency gating, the unit pointing vector of the edges is used. (From upstream node to downstream node), define the difference in the microbial community memory index between the two ends of the edge. And construct a bio-evidence consistency score: ;

[0076] in, Score the path consistency of edge e at time step t; It is an S-shaped compression function; These are the gating parameters. If... This indicates that there exists a monotonic gradient along the edge direction from "high memory" to "low memory". Increases, and vice versa. To utilize metabolic "clock" information, functional markers or enzyme activity levels sensitive to time decay are selected. The timescale of pollution input is estimated by the activity ratio of upstream and downstream sides; ;

[0077] in, Estimating the metabolic clock of edge e at time step t; This is the attenuation rate parameter for this function in a borderland environment; and These represent the functional activity intensities of the upstream and downstream sides, respectively. and The biological gain mapped to the edge confidence is updated incrementally using log odds.

[0078] ;

[0079] in, The edge confidence level after introducing biological evidence; Prior confidence level before biological evidence; ; This is the gain coefficient; To map the metabolic clock as a monotonic function of confidence gain (e.g., taking positive values ​​within a reasonable time window and negative values ​​outside the window).

[0080] When groundwater lag experts implement the procedure, a delay queue for surface-subsurface composite transport is established for edges labeled as associated with riverbank exchange or alluvial fans. Let the surface volume flux time series of the edge be... The lag flux of groundwater quantity is obtained by convolving it with a normalized kernel function.

[0081]

[0082] in, Let e ​​be the volume component of the groundwater at edge e; The delay kernel for edge features is taken as a gamma distribution. For shape parameters; Let be the scale parameter; Γ(⋅) be the gamma function. Let the groundwater sharing coefficient of the edge be... Then the mean and variance of the equivalent travel time for that side can be written as:

[0083] ;

[0084] in, Let e ​​be the mean of the equivalent travel time at time step t; This represents the equivalent travel time variance. The length of the side; Let be the equivalent phase velocity of the edge at time step t (composed of the three-phase flux). The travel time of the path is approximately an additive variable of the travel times of each edge. For any directed path p from the upstream source region to the sensitive section, we have: ;

[0085] in, and Let be the mean and variance of the travel time of path p at time step t, respectively. The quantile arrival times are obtained using a normal approximation: ;

[0086] in, For quantiles Arrival time; These are the quantiles of a standard normal distribution. Therefore, p10 / p50 / p90 are directly output at the sensitive cross-section within the time window.

[0087] During the implementation of the observation consistency expert program, flow velocity monitoring, water level observation, pollutant concentration detection, and groundwater level sequences during the event period are collected. Observation operators are established according to instrument and cross-sectional geometry, mapping the boundary-level states to the predicted values ​​of each observation point. Let the set of observation points be denoted as... The observed value is The predicted value is Sensitivity weights for edge pairs of observations Constructing the edge likelihood:

[0088] ;

[0089] in, Let e ​​be the consistency likelihood of edge e at time step t; Let r be the noise standard deviation at observation point r. The likelihood is incorporated into the confidence score using a log-odds summation method.

[0090] ;

[0091] in, The edge-level confidence level after observation consistency fusion; This refers to the observation gain coefficient. To ensure the overall map quality is conserved and non-negativity, the imbalance between income and expenditure is considered at each node. Prioritize deductions and allocate them according to the "deductible weight" ratio of each edge: ;

[0092] in, This is the edge volume flux after deduction; The edge volume flux before deduction; This represents the imbalance between income and expenditure at node n at time step t (a positive value indicates that it needs to be deducted from the outflow side). Let the set of edges flowing out of the node be denoted as 'node'. Non-negative weights associated with edge capacity or resistance are used to avoid forcibly deducting excessive throughput to edges with insufficient capacity. The aforementioned fusion confidence level. This formula prioritizes low-confidence edges for flux adjustment, avoiding excessive perturbation to high-confidence paths.

[0093] At the output layer, edge flux, travel time distribution, and confidence scores, calibrated for consistency between biological evidence and observations within the same time step, are aggregated and adjudicated. Edges with confidence scores below the threshold and path segments with non-monotonic arrival times are removed. The main diffusion ridge is extracted along paths with high confidence and continuous flux. At sensitive cross-sections, the p10 / p50 / p90 arrival time windows are estimated using quantiles. The contribution rate of the source region is calculated based on the pollutant mass flux passing through the edge at the cross-section. Let the set of management units (fields) be M, and the cross-section contribution of a certain unit m is defined as follows: ;

[0094] in, The contribution rate of the time step t management unit m at the sensitive section S; Let m be the set of edges originating from element m and reaching section S; Let e ​​be the pollutant mass flux at section S. The identification of high-risk convergence zones integrates convergence indices and mass density based on virtual deformation at the edge unit level, using: ;

[0095] in, Risk score for unit T; For unit Jacobi (less than 1 indicates convergence); Estimate the mass density of pollutants within the cell; This is a weighting factor. For Areas exceeding a set threshold are marked as high-risk and used for governance priority ranking.

[0096] The working process of this embodiment is as follows: Based on the interannual variation network and forecasting linkage formed in the previous embodiment, edge-level dynamic parameters are periodically generated and edge flux and initial values ​​of arrival time are calculated on a rolling basis; microbial functional pathway data are acquired within the sampling period to form edge-level microbial community memory index, path consistency score and metabolic clock estimation, and biological gain is applied to the edge-level confidence; surface-subsurface delayed convolution is performed on the relevant edges of riverbanks and alluvial fans to obtain the mean and variance supplements of groundwater content and path travel time; at each time step, observational data such as flow velocity / water level / concentration / well water level are collected, consistency likelihood is calculated and edge-level confidence is updated, and then mass conservation calibration with priority deduction is performed at the nodes; finally, the main diffusion ridge, p10 / p50 / p90 arrival time window and source area contribution rate are output with aggregation adjudication, while high-risk convergence zones are identified and composite path uncertainty is given. Compared to scenarios without biological and groundwater evidence, this embodiment can improve the reliability of main diffusion path direction determination by gating the path with the microbial community memory index under conditions of significant dilution of chemical concentration or sparse observation; improve the estimation of time windows for distant and delayed arrivals by explicitly characterizing the surface-subsurface coupling wake through delayed queues; and maintain overall map quality conservation by preferentially deducting flux from low-confidence edges through observation consistency, ensuring that the output main diffusion ridge and quantile arrival times remain stable and interpretable throughout the entire event evolution process.

[0097] The above description is merely a preferred embodiment of the present invention. It should be understood that the present invention is not limited to the forms disclosed herein and should not be construed as excluding other embodiments. It can be used in various other combinations, modifications, and environments, and can be modified within the scope of the concept described herein by means of the above teachings or the technology or knowledge in related fields.

Claims

1. A method for predicting and tracing agricultural non-point source pollution, characterized in that, Includes the following steps: S1. Constructing a standard gravity triangulation network: Within the agricultural area, topographic information is acquired based on a high-resolution digital elevation model. A constrained Delaunay triangulation network is established, using canals, fields, roads, watersheds, and riverbanks as boundary constraints. This network meets quality control requirements, with a minimum interior angle of 30 degrees and an aspect ratio not exceeding 3. Time-lapse gravity observation data is loaded onto the nodes and edges of the triangulation network. Gravity residuals are calculated and converted into virtual deformation information for the grid, determining the priority channel weights and convergence / divergence characteristics of each edge. The nodal displacements are constructed according to the following formula, and the deformation and convergence / divergence characteristics of the mesh elements are calculated accordingly: ； in, Let be the virtual displacement vector of node p; The gravity residual scalar field at node p is processed by the smoothing operator and is obtained by spatial smoothing of the observation residual field. This represents the spatial gradient operator over the node domain; The displacement scale coefficient is used to map the dimensions of the gravity residual to the deformation scale. After updating the three nodes of the triangle using the virtual displacement vector, the deformation gradient tensor is used to measure the scaling and orientation information of the mesh cells. ； in, Let T be the deformation gradient tensor of the triangular mesh element T; The initial node coordinates of the element; These are the node coordinates after virtual deformation; This refers to the ratio of unit area, i.e., the Jacobian. Representatives gathered, Represents divergence; The principal elongation direction is obtained using the right Cauchy-Green tensor to indicate the dominant direction of potential diffusion: ； in, Let T be the right Cauchy–Green tensor of cell T; The unit eigenvector corresponding to the largest eigenvalue; The largest eigenvalue; Project along the outer edge of the element, and combine Assign direction and degree of convergence to adjacent edges; The edge-priority channel weights are calculated using the following formula: ； in, Let e ​​be the priority channel weight; It is an S-shaped compression function; The angle between the direction of edge e and the main extension direction of the adjacent element; The Jacobian of the unit adjacent to edge e; The category bias assigned to the edge semantic classification based on ditch channels, slope confluence, field boundaries, road drainage, riverbank exchange, or ventilation corridors; These are the weighted composite coefficients; S2. Establish a border-level physical parameter coupling system: Establish a unified physical parameter package for each edge of the triangular network. The parameter package includes water phase transport parameters, soil phase transport parameters, gas phase transport parameters, and interphase exchange parameters. Among them, water phase transport parameters include roughness coefficient, channel capacity, head loss coefficient, inflow bias weight, sediment carrying capacity, and dilution mixing length; soil phase transport parameters include boundary permeability, permeability adjustment factor, water storage coupling strength, and adsorption label; gas phase transport parameters include effective roughness length, ventilation corridor coefficient, shear shielding coefficient, and settlement bias; and interphase exchange parameters include interphase exchange weight and reactivity label. S3. Execute programmable expert sequential aggregation: On the same side node state space of the triangular network, six types of programmable expert modules are sequentially called in a predetermined order to update the edge state, including: Gravity Deformation Expert: Generates gravity weights, inflow biases, and priority channel sequences based on gravity residuals and virtual deformation information; Three-Phase Flow Field Expert: Calculates water, soil, and air flux distribution at each edge based on shallow water runoff, soil seepage, and near-surface wind, and provides preliminary arrival times; Interannual Variation Expert: Corrects edge parameters year by year based on multi-year historical meteorological data and multi-temporal satellite remote sensing images; Microbial Functional Memory Expert: Implements path consistency gating and timeline verification for candidate paths based on the Microbial Community Memory Index (SMI); Groundwater Lag Expert: Establishes a delay queue for surface-subsurface composite transport at the boundary between the riparian zone and the groundwater area; Observational Consistency Expert: Calibrates edge-level confidence based on measured flow velocity, water level, pollutant concentration, and groundwater level data; S4. Integrated meteorological forecast-driven prediction: Short-term meteorological forecast data is projected onto side parameters through a forecast-linking mechanism. The short-term meteorological forecast data includes hourly or three-hourly rainfall, wind direction and speed, and temperature. Rainfall is mapped to an increment of inflow bias, wind direction and speed are mapped to an adjustment of ventilation corridor coefficient, and temperature is mapped to a correction of interphase exchange or reaction activity, so as to advance the prediction of pollution diffusion paths and arrival times in future periods. S5. Multi-expert aggregation and result output: On each edge, the outputs of six types of experts are aggregated and adjudicated according to weight gating and threshold elimination rules to generate a single edge flux, pollutant mass distribution, arrival time probability distribution and confidence score; at the edge level, mass conservation and non-negativity constraints are implemented, and flux is preferentially deducted for edges with low confidence; at the node level, the edge-level results are summarized and the main diffusion ridge of agricultural non-point source pollution, the arrival time window of sensitive sections and the source area contribution rate are output.

2. The method according to claim 1, characterized in that, In the process of constructing the gravity standard triangulation network in step S1, the edges are semantically classified and labeled as ditch channels, slope confluence, field boundaries, road drainage, riverbank exchange, or ventilation corridors, which are used to drive the subsequent differential setting of three-phase parameters and weights.

3. The method according to claim 1, characterized in that, In the gravity residual transformation process of step S1, by calculating the changes in the area and angle of the triangular unit, the main elongation direction is determined as the dominant direction of pollution diffusion, the area with reduced area is marked as the convergence zone, and the area with increased area is marked as the divergence zone, forming a priority channel sequence at the edge level.

4. The method according to claim 1, characterized in that, In the interannual variation expert module of step S3, the interannual variation characteristics of rainfall pattern, vegetation cover index NDVI, bare land rate, soil moisture and surface water frequency are extracted by analyzing historical meteorological records of more than ten years and multi-temporal satellite remote sensing images, and the edge parameters are corrected year by year.

5. The method according to claim 1, characterized in that, In the microbial functional memory expert module of step S3, soil and water microbial samples are collected along the main diffusion ridge, bifurcation point, riparian zone and key cross-section. The microbial community memory index (SMI) is constructed by 16S gene sequencing and metagenomic analysis, and the path consistency is verified by the spatial distribution gradient of functional genes and enzyme activities.

6. The method according to claim 1, characterized in that, In the groundwater lag expert module of step S3, the coupling relationship between the surface and groundwater is established at the riverbank zone, alluvial fan boundary and shallow groundwater observation point. The groundwater lag transport process is simulated through a delayed queue mechanism to supplement the long-term pollution risk assessment.

7. The method according to claim 1, characterized in that, In the observation consistency expert module of step S3, flow velocity monitoring, water level observation, pollutant concentration detection and groundwater level data are collected in real time and compared with the edge-level prediction results. The confidence weight of edges that are inconsistent with the observation data is reduced to ensure the quality conservation constraint.

8. The method according to claim 1, characterized in that, In the meteorological forecast fusion in step S4, the forecast hooking mechanism includes explicit mapping rules: mapping the forecast rainfall to a weighted increment of the edge inflow bias, mapping the forecast wind direction and speed to an adjustment of the ventilation corridor coefficient and shear shielding coefficient, mapping the forecast temperature to an adjustment of the interphase exchange weight or reactive label, and limiting the above mapping to maintain edge mass conservation within the time step.

9. The method according to claim 1, characterized in that, In the multi-expert aggregation process in step S5, the outputs of the six types of experts are superimposed and aggregated through weight gating, threshold elimination and normalization to generate a comprehensive confidence score for the edges.

10. The method according to claim 1, characterized in that, The source region contribution rate decomposition is calculated by aggregating, normalizing, and aligning the edge fluxes from different fields or management units at the target cross section, and simultaneously outputs the main diffusion ridge map, high-risk convergence zone identification, and composite path uncertainty assessment; the source region contribution rate is used for governance priority ranking and source tracing decision support.

Images