Agricultural non-point source pollution inversion method and system based on multi-level topological constraint

Abstract

Aiming at the problem that the existing agricultural non-point source pollution inversion method is unstable and has unclear physical meaning under the condition of limited monitoring section, this paper proposes a field scale agricultural non-point source pollution inversion method and system based on multi-level topological constraints. Firstly, the study area is divided into three levels of field, catchment unit and river section by using digital elevation model and spatial data, and the pollution path relationship of "field → catchment unit → section" is established. Then, by integrating the data of terrain slope, flow distance, soil properties, vegetation coverage and rainfall, a matrix is constructed to express the pollution transmission process from field to section. Taking the measured flux of section as the main condition, three rules of non-negative, mass conservation and spatial smoothing are introduced: non-negative ensures the physical reasonableness of load value, mass conservation ensures the conservation of matter, and spatial smoothing avoids the sharp jump of load in adjacent fields. The three rules form a compound objective function, which is used to invert the pollution load of each field. Finally, according to the inversion results, the pollution contribution rate calculation and spatial distribution mapping are carried out. Under the condition of limited section data, this method can stably and accurately analyze the pollution source at field scale, and support precise management decision-making of the watershed.

Description

Technical Field

[0001] This invention belongs to the field of agricultural environmental protection and watershed water pollution control technology. Specifically, it relates to a technical method and system for constructing a transmission matrix using multi-scale spatial topological relationships and solving the field-scale agricultural non-point source pollution load into the river under multiple constraints such as river cross-section flux. Background Technology

[0002] With increasing agricultural intensification, nutrients such as nitrogen and phosphorus accumulate in farmland. Driven by rainfall and runoff, these substances migrate and diffuse into river systems, forming a non-point source pollution pattern with complex sources and dispersed spatial distribution. How to reverse-analyze the pollution contribution of different fields from limited river monitoring section data, and accurately characterize the migration relationship of pollution between spatial units, is a technical challenge currently facing the refined management of watershed water environment. Existing methods for estimating agricultural non-point source pollution loads mainly fall into three categories: output coefficient models, hydrological process models, and statistical regression models. While these methods can estimate pollution loads to some extent, they share several common problems. First, most methods use watersheds or sub-watersheds as the basic calculation unit, resulting in coarse spatial granularity that fails to meet the requirements of field-level responsibility allocation. Second, these models are mostly forward simulations of "source → sink," lacking a structured mathematical expression of the multi-level spatial topological relationships of "field-catchment-section," leading to the inability to effectively feed back monitoring section information into the analysis of source loads. Third, there is currently no inversion framework that can simultaneously integrate conditions such as cross-sectional flux, non-negativity, mass conservation, and spatial smoothness; therefore, the inversion results are often unstable or lack clear physical meaning. To address the aforementioned problems, this paper proposes a novel approach: first, the spatial topological and hydraulic connections between "field-catchment unit-cross-section" are encoded into a computable transport matrix; then, based on this matrix, a composite solution framework integrating cross-sectional flux, mass conservation, and smoothing conditions is constructed; finally, within this framework, the pollution load entering the river at the field scale is stably inverted. This method provides a feasible technical path to overcome the problems of low resolution, weak topological representation, and lack of conditional mechanisms in traditional methods. Summary of the Invention

[0003] The purpose of this invention is to address the problems of inaccurate pollution source inversion and unstable solutions caused by the limited number of monitoring sections in existing technologies. To this end, an inversion method and system based on multi-level spatial topological conditions is proposed. The core idea is to matrixify the spatial relationships and solve them under multiple physical and mathematical constraints. The specific technical solution is as follows: (1) Multi-level spatial unit division and topological relationship construction First, the target watershed is divided into three nested levels: field units (the smallest decision-making units), catchment units (extracted from DEM hydrological analysis, serving as an intermediate scale), and river cross-section units (water quality and flow monitoring points). Then, using the flow direction and cumulative runoff extracted from the DEM, a directed connection is established from the field, through the catchment units, and finally into a specific cross-section. This forms a structured set of pollution transport paths. This multi-level structure is the foundation for transferring cross-section condition information to the field scale, and compared to single-grid inversion, it better balances field-level accuracy with the consistency of cross-section observations. (2) Construction of multi-level topological constraint transfer matrix Based on the aforementioned spatial topology and transmission paths, this paper encodes the pollution migration relationship into a quasi-linear transmission matrix A. The specific approach involves three steps: First, based on the DEM and hydrological flow direction, determine the primary (or unique) transport path for each field to reach the target section. Second, for each path, a transmission (or attenuation) coefficient is calculated by comprehensively considering topographic slope, path length, underlying surface type (such as crop cover and soil texture), and rainfall driving factors. This coefficient represents the proportion of a unit load from a field that ultimately reaches the cross-section. Third, the transmission coefficients from all fields to all cross sections are summarized to form a transmission matrix A, whose elements are... (1) here, It is the path length. It's the slope. is the underlying surface factor, and k is the attenuation coefficient. The A matrix transforms the constraints of multi-level topology on pollution migration into a mathematically solvable form. (3) Construction and solution of inversion model based on cross-sectional flux composite constraint The core of this step is to construct an inversion model that uses river cross-section monitoring data as a hard condition and integrates various physical and mathematical conditions to stably solve for the field pollution load vector L. First, establish a basic linear relationship model. Let the field load vector be L (n×1) and the cross-sectional observation flux vector be F (m×1), then there exists... (2) Where ε is the error term. The goal is to solve for L. Next, the objective function under the combined conditions is constructed. To ensure the uniqueness, stability, and physical meaning of the solution, the following multiple conditions are introduced and integrated into a single function for solution: a. Cross-sectional flux constraints: The core is to make the simulated flux A·L as close as possible to the measured value F, i.e. (3) b. Non-negativity constraint: requires L ≥ 0, which is a basic condition to ensure the physical rationality of the load. c. Mass Conservation Constraint: At the catchment unit scale, the output load (after attenuation) of all fields within the unit must be substantially consistent with the estimated flux at the unit's outlet. This is achieved through a penalty term. The form is added to the model. d. Spatial smoothness constraint: Introducing a regularization term (4) Q is the Laplace matrix that characterizes the adjacency relationship between fields. Its function is to suppress drastic changes in the load of adjacent fields, making the results more continuous and smooth, and reducing noise interference. The final objective function is integrated as follows: (5) Constraints in, and This is a regularization parameter used to balance the weights of different constraint terms in the objective function; The penalty term for the mass conservation constraint is expressed as follows: (6) Finally, an optimization algorithm with nonnegativity constraints (such as nonnegative least squares or its variants) is used to iteratively solve the above objective function until convergence, yielding the optimal estimate of L. This framework integrates monitoring data, spatial structure, physical laws, and mathematical priors, ensuring the stability of the inversion. (4) Analysis and spatial allocation of pollution contribution rate After obtaining the field load vector L, and combining it with the transfer matrix A, the pollution contribution rate of any field to any cross section can be quantitatively calculated: (7) At the same time, it can aggregate the contributions of different scales, such as the catchment units, to form a multi-level spatial distribution structure of pollution load. (5) Result output After analyzing the results, a pollution load distribution map at the field scale and a spatial distribution map of the contribution rate into rivers are generated. Spatial analysis tools are then used to delineate areas with significant contributions for reference in governance decisions. (6) System composition and data flow organization methods This method is also accompanied by a system, which consists of six modules: data acquisition, topology construction, transmission matrix construction, flux condition modeling, inversion solution, and result analysis. The data flow follows the sequence of "spatial description → structure construction → mathematical expression → condition solution → result analysis," and key intermediate steps can be adjusted to allow observational data to correct the inversion process. Attached Figure Description

[0004] Figure 1 A general technical flow diagram of an embodiment of the present invention is shown; Figure 2 A schematic diagram of the multi-level (field-catchment unit-section) spatial unit topology structure in an embodiment of the present invention is shown; Figure 3 A schematic diagram illustrating the construction process of encoding pollution migration relationships as a transmission matrix A in an embodiment of the present invention is shown; Figure 4 The diagram illustrates the inversion model structure and solution logic diagram that integrates multiple constraints in an embodiment of the present invention. Figure 5 An example diagram showing the results of field-scale pollution load inversion and spatial allocation of contribution rate in an embodiment of the present invention is provided. Figure 6 A step-by-step flowchart of an embodiment of the present invention is shown. Detailed Implementation

[0005] To make the objectives, solutions, and advantages of the present invention clearer, a detailed description is provided below in conjunction with specific embodiments. Figure 1 The diagram shows the overall technical flow of an embodiment of the present invention, and this embodiment is implemented step by step according to this flow. This example selects a typical hilly agricultural watershed in Area A. The watershed area is approximately 150 km², with terrain mainly consisting of gently sloping farmland, and agricultural non-point source pollution characteristics are quite obvious. Conventional water quality monitoring sections have already been set up in the area, providing the basic conditions for conducting inversion analysis. (1) Basic data acquisition and preprocessing The data used included: a 12.5 m resolution DEM (for hydrological analysis), Sentinel-2 multispectral images (for extracting NDVI and crop cover), agricultural statistics (fertilization intensity and planting structure), and continuous monitoring data from three hydrological and water quality monitoring sections (TN, TP concentration and flow rate, with a time resolution of monthly scale). Data processing procedure: First, depression filling, flow direction, and runoff accumulation calculations were performed on the DEM. The river network was extracted using the threshold method, and 37 catchment units were delineated. Then, the plot boundary vector data and land use classification results were overlaid to subdivide the study area into 1256 plot units. Figure 2A schematic diagram of the multi-level (field-catchment unit-section) spatial unit topology in an embodiment of the present invention is shown. The hierarchical relationship between the above-mentioned field, catchment unit and monitoring section is shown in the figure. Next, the water quality concentration and flow data of each section are converted into pollution flux time series, and then the adult-scale flux data are summarized to construct the target vector F for inversion. (2) Construction of transmission matrix A Based on the water flow direction extracted from the DEM, the water flow path of each field is traced to determine the catchment units it passes through and the final target cross-section (cross-section A, B, or C) into which it flows. For the path from each field j to its target cross-section i, a transport coefficient between 0 and 1 is calculated by comprehensively considering path length, average slope, crop cover along the path, and soil erodibility factor (K factor), using the following formula: in For path length, For average soil erodibility, The average slope Let be the average crop coverage factor, and k be a rate constant. In this way, a 3 × 1256 dimensional transmission matrix A is constructed. Figure 3 The diagram illustrates the process of constructing a transfer matrix A to encode the pollution migration relationship in an embodiment of the present invention. (3) Solving the composite constraint inversion model Based on the TN concentration C and flow rate Q at three monitoring sections, the annual total nitrogen flux was calculated: (Unit: tons / year) Then construct the inversion objective function: The constraint condition is L ≥ 0. Here, L is the annual TN output load vector for the 1256 fields to be determined. This is the mass conservation penalty term, calculated as the sum of squares of (the sum of all field loads after decay within the unit - the unit outlet flux estimated from the field loads) at the scale of 37 catchment units. Matrix Q is a Laplace matrix constructed based on the shared edge or concurrent adjacency relationships of the fields. λ1 = 0.5 and λ2 = 0.1 are determined using the L-curve method. Figure 4 The diagram illustrates the inversion model structure and solution logic diagram of the embodiment of the present invention, which integrates multiple constraints. The objective function, constraint terms and solution process are all reflected in the diagram. Finally, an iterative shrinkage-threshold algorithm with non-negativity constraints is used to solve the objective function. It converges after approximately 200 iterations, yielding the optimal estimate of the field load vector L. (4) Results analysis and comparative verification The inversion yielded the annual TN output load for each of the 1256 fields. Figure 5 An example diagram showing the spatial distribution results of pollution load inversion and contribution rate at the field scale in an embodiment of the present invention is provided, which visually illustrates the spatial distribution of the load and the differences in contribution rate among different fields. To verify the effectiveness, the inversion results are aggregated by catchment unit and compared with the estimation results of the traditional output coefficient method. Some results are shown in the table below: <![CDATA[SU 05 ]]> 12.5 9.8 high <![CDATA[SU 12 ]]> 28.3 36.1 high <![CDATA[SU 21 ]]> 15.7 14.2 high Total of the entire basin 45.2 tons 41.5 tons Completely consistent with the total cross-sectional flux (41.5 tons). Analysis shows that: 1) The total load of the entire basin retrieved by this invention is in perfect agreement with the sum of the measured total fluxes at the three cross sections, indicating that the model strictly adheres to the mass conservation condition. 2) Compared with the traditional output coefficient method, this method significantly adjusts the load distribution at the catchment unit scale. This is because the spatial distribution of the load is recalibrated through cross-sectional conditions and the transmission matrix, making it more consistent with actual hydrological paths and monitoring results. 3) SU was identified. 12 SU 08 Five catchment units, accounting for approximately 45% of the total pollution load in the watershed, were identified as priority control areas. This embodiment verifies the effectiveness of this method in the refined inversion of agricultural non-point source pollution in small watersheds. Table 1. Core Input, Parameters, and Output Datasets for Section Constraint Inversion in Area A 1 Cross-sectional water quality concentration (TN / TP) mg / L The pollutant concentration at the river monitoring section is used to calculate the flux observation vector F. 2 Cross-sectional flow (Flow) m³ / s The cross-sectional water flow rate monitored simultaneously is used to calculate the flux observation vector F. 3 Cross-sectional pollution flux (Flux) kg / d or t / a Concentration × Flow Rate constitutes the target vector F of the inversion model. 4 Transfer Matrix (A) Dimensionless Pollution migration coefficient matrix from field to cross section, dimension m×n, core parameters

Claims

1. A method for inverting agricultural non-point source pollution load at the field scale based on multi-level topological conditions, characterized in that, Includes the following steps: S1: Obtain basic field-scale data within the study area, including plot spatial boundaries, digital elevation models, soil properties, crop types, hydrological flow direction, and measured values ​​of cross-sectional pollution fluxes; S2: Based on the spatial adjacency relationship between fields and the hydrological connectivity rules, construct a spatial topology network with fields as nodes and directed edges representing the confluence relationship, denoted as G=(V, E), and define node attributes and edge weight functions; S3: Based on the aforementioned spatial topology network, establish a topological transport matrix A to describe the migration process of pollutants between fields. Matrix elements... The calculation formula is in, Indicates path length. Indicates the slope factor. Let represent the underlying surface factor, and k be the attenuation coefficient; S4: Based on the cross-sectional pollution flux observation data, construct the coupling constraint equation between the topological transport matrix and the cross-sectional flux, and establish a basic linear relationship model. Define the field pollution load vector L (dimension n × 1) and the cross-sectional pollution flux observation vector F (dimension m × 1), then the following relationship exists: Where ε is the error term, the goal of the model is to solve L, and further construct a multi-constraint coupled system, including: (1) Cross-sectional fitting constraint: Minimizing the residual between the simulated flux A·L and the observed flux F is taken as the core objective, which is usually expressed as minimizing the sum of squared residuals, i.e.: (2) Non-negativity constraint: It is stipulated that the pollution load of all field units must be non-negative, that is, L ≥ 0; (3) Mass conservation condition: For each catchment unit, the output load of all fields within it, after attenuation along the flow path, should be basically consistent with the estimated flux at the outlet of the unit; (4) Spatial smoothing condition: Introduce a smoothing regularization term R(L) between adjacent fields, for example, based on the quadratic form of the Laplace matrix, i.e.: Where Q is a positive semidefinite matrix representing the spatial adjacency relationship of fields. The purpose of this constraint term is to suppress drastic changes in the load of adjacent fields, making the inversion results more continuous and smooth in space, and effectively reducing outliers caused by observation noise or model errors. S5: Combine the above conditions into a unified objective function for solution, the function form of which is: The constraints are defined by the regularization parameter, which is used to balance the weights of different constraint terms in the objective function. The penalty term for the mass conservation constraint is expressed as follows: S6: Based on the pollution load results obtained from the inversion, calculate the pollution contribution rate of the field to the cross section: And perform spatial allocation and expression.

2. The method according to claim 1, characterized in that, In step S2, only directed edges that satisfy the hydrological flow direction constraints are retained:

3. The method according to claim 1, characterized in that... In step S3, the transfer matrix satisfies the normalization condition:

4. The method according to claim 1, characterized in that, The edge weight function adopts an exponential decay form:

5. The method according to claim 1, characterized in that, The spatial smoothing matrix Q is defined as: Q=D−W Where D is the degree matrix and W is the adjacency weight matrix.

6. The method according to claim 1, characterized in that, The mass conservation constraint is embedded in the objective function in the form of a penalty term.

7. The method according to claim 1, characterized in that, The regularization parameter is determined by the L-curve method or cross-validation.

8. The method according to claim 1, characterized in that, Further, this includes identifying high-contribution fields: L i >t 9. The method according to claim 1, characterized in that, Output multi-scale pollution allocation results, including field scale, catchment unit scale, and watershed scale.

10. A field-scale agricultural non-point source pollution inversion system based on multi-level topological constraints, characterized in that, include: The module comprises a data acquisition module, a topology construction module, a transfer matrix construction module, a flux constraint modeling module (corresponding to S4), an inversion solution module (corresponding to S5), and a result parsing module. The flux constraint modeling module is used to construct the following constraint relationships: F=AL The inversion solution module is used to solve the following optimization problem:

11. The system according to claim 10, characterized in that: The modules are connected in series through a data flow-driven approach, forming a closed-loop computation process of "data → topology → matrix → constraint → optimization → output".

Images