A method for tracing based on coupling hydrodynamics and pollutant degradation equations

Abstract

The present application belongs to the cross field of environmental engineering and hydrology and hydrodynamics, and specifically is a tracing method based on coupling hydrology and pollutant degradation equation, which comprises the following steps: firstly, obtaining multi-source heterogeneous monitoring data and preprocessing; secondly, selecting one-dimensional Saint-Venant equation or two-dimensional shallow water equation to solve the hydrodynamic field according to the water body form; thirdly, coupling multi-component convection-diffusion-reaction equation to simulate the migration and transformation of pollutants; fourthly, introducing LSTM module to identify suspected pollution events; fifthly, through the double-channel framework of PDE constraint optimization and Bayesian inference, the pollution source parameters are inversed and the uncertainty is quantified; sixthly, combining data assimilation to update the model parameters online; finally, quantifying the uncertainty and verifying it. The present application takes into account the tracing accuracy, efficiency and compliance, supports multiple water bodies and multiple data sources, and is suitable for complex water body pollution tracing.

Description

Technical Field

[0001] This invention belongs to the interdisciplinary field of environmental engineering and hydrology and hydrodynamics, specifically a source tracing method based on the coupling of hydrodynamics and pollutant degradation equations. Background Technology

[0002] In the interdisciplinary field of environmental engineering and hydrology and hydrodynamics, source tracing of water pollution is a core prerequisite for the precise treatment of industrial wastewater, urban sewage, and agricultural non-point source pollution. Its goal is to clarify the source location, emission sequence, intensity, and time-varying contribution rate of pollutants, providing a basis for pollution accountability and the formulation of treatment plans. As the pollution scenarios faced by various types of water bodies, such as rivers, lakes, river networks-lake groups, and estuaries, become increasingly complex (e.g., dam scheduling, tidal alternation, and non-steady-state inflows), higher demands are placed on the accuracy, stability, and interpretability of source tracing technologies. Traditional source tracing methods are no longer adequate for the practical needs of complex hydraulic conditions and the coexistence of multiple pollution sources, necessitating an integrated technical solution that combines mechanistic models and intelligent computing.

[0003] Existing pollution source tracing technologies have significant limitations: Traditional methods rely on tracer administration or empirical back-inference from a few monitoring sections, which cannot accurately depict the migration and diffusion patterns of pollutants under complex hydrodynamic fields. When multiple pollution sources coexist, equivalent source confusion is likely to occur, leading to deviations in pollution source location and misjudgment of timing. Purely data-driven tracing solutions (such as end-to-end neural networks) can process time-series data, but they lack physical mechanism support, have poor interpretability, are insufficient in generalization to rare operating conditions (such as sudden emissions and sudden changes in boundary conditions), and are difficult to meet the audit compliance requirements required by regulators. Pure mechanistic models that rely solely on the Saint-Venant equation or shallow water equation and convection-diffusion-reaction equation face "illness" and high-dimensional uncertainty problems in the source term inversion process. Furthermore, complex factors such as wet-dry boundary, wind stress, and tides significantly increase the difficulty and time consumption of numerical calculations, making it impossible to achieve efficient online source tracing.

[0004] The aforementioned technical shortcomings make it difficult for existing solutions to balance traceability accuracy, efficiency, and compliance: traditional methods lack accuracy, purely data-driven methods lack interpretability, and purely mechanistic models are inefficient. Furthermore, existing technologies generally do not consider alternative mechanisms in case of intelligent module failure, making them prone to traceability interruptions in scenarios with limited algorithm permissions or data drift.

[0005] Therefore, this invention provides a source tracing method based on coupled hydrodynamics and pollutant degradation equations. Summary of the Invention

[0006] In order to overcome the shortcomings of the prior art, at least one technical problem raised in the background art is solved.

[0007] The technical solution adopted by this invention to solve its technical problem is: a source tracing method based on coupled hydrodynamics and pollutant degradation equations, comprising the following steps: S1. Data acquisition and preprocessing: This is used to acquire and integrate multi-source heterogeneous monitoring data and perform standardization processing. The processing includes: filling in missing observation data using interpolation algorithms, and identifying and removing outliers using the interquartile range method combined with physical constraints. S2. Hydrodynamic solution is used to accurately simulate the water flow state of the target water body. The simulation adopts an adaptive model selection strategy, and solves the equations using the one-dimensional Saint-Venant equation or the two-dimensional shallow water equation according to the water body morphology. The finite volume method is used to discretize the governing equations, and the Runge-Kutta method with high-order total variation reduction is used for time advancement. S3, Reaction-Migration Coupling, is used to simulate the migration and transformation of pollutants in water bodies based on the solved hydrodynamic field, and to solve multi-component convection-diffusion-reaction equations through process coupling. S4 and LSTM diffusion event recognition are used to introduce a diffusion event recognition module based on long short-term memory network to identify emission events of suspected external pollution sources from multi-index time series of multiple monitoring sections. The module is trained to perform joint classification and regression tasks and output emission time window, pollutant type and preliminary location clues. S5, Source Inversion and Contribution Rate Assessment, is used to estimate the location, start-up and shutdown time and emission intensity of pollutant source terms based on the dual-channel inversion framework, quantify the uncertainty of the inversion results, and calculate the real-time time-varying contribution rate of each candidate pollution source to the target section or region through the contribution rate assessment module. S6, Assimilation and Parameter Calibration, is used to update key hydrodynamic parameters, pollutant degradation rates, and boundary condition error terms online using data assimilation technology, thereby achieving online model updates and improving prediction accuracy; S7, Uncertainty Quantification and Verification, is used to propagate perturbations to uncertainties in model parameters, boundary conditions, and model structure, thereby outputting the confidence interval and robustness score of the source tracing results, and combining them with historical pollution events for posterior verification.

[0008] As a preferred approach, S1 data acquisition and preprocessing specifically includes: Collect hydrological data, hydraulic data, meteorological data, geographic information data, sediment characteristic parameters, a list of known discharge outlets and their locations, as well as various water quality data; All collected data underwent unified calibration of coordinate system and time reference. For missing observation data, methods based on Kalman filtering, expectation-maximization algorithm or neighborhood regression are used to fill in the gaps; For outliers, the interquartile range method is used in combination with physical constraints of maximum flow velocity and minimum water depth for identification and removal. After alignment, the multi-source data is quality-labeled to form a standardized input dataset and a set of prior parameters containing geographic, hydrological, and water quality information.

[0009] As a preferred method, the S2 hydrodynamic solution specifically includes: For water bodies with significant one-dimensional characteristics, the one-dimensional Saint-Venant equation is used for simulation. The one-dimensional Saint-Venant equation includes the mass continuity equation and the momentum equation. For water bodies with significant two-dimensional characteristics, the solution is obtained through the depth-averaged two-dimensional shallow water equation, which includes the mass conservation equation, the momentum conservation equation in the x-direction, and the momentum conservation equation in the y-direction. In terms of numerical implementation, the finite volume method is used to discretize the governing equations, and the flux is calculated using an HLL or HLLC Riemann solver. The time progression employs a second- or third-order Runge-Kutta method to reduce higher-order total variation; The source term balance, wet-dry boundary treatment, and dry bed conservation treatment were optimized. Boundary conditions are set including upstream given flow rate or water level, downstream given water level or free outflow, and the effects of wind stress, rainfall, evaporation and tidal harmonics or measured tidal sequences are taken into account.

[0010] As a preferred approach, S3 reaction-migration coupling specifically includes: Flux reconstruction of the convection term employs a MUSCL limiter to suppress numerical oscillations; The dispersion term considers longitudinal and lateral dispersion as well as turbulent diffusion effects; The reaction term is multi-component, including at least first-order decay kinetics and biochemical degradation processes, and further includes Arrhenius temperature correction, dissolved oxygen limitation, Monod or multi-substrate kinetic models, adsorption-desorption and sedimentation-resuspension processes, depending on the pollutant characteristics and environmental conditions. To address the strong nonlinear coupling between convection, diffusion, and reaction terms, an operator splitting method is employed.

[0011] As a preferred method, S4LSTM diffusion event identification specifically includes: The LSTM diffusion event identification module takes multi-index interval time series from multiple monitoring sections as input. The sliding window length of the time series is adaptively adjusted between 30 and 180 steps, and the step size is selected between 1 and 5 steps. The network structure uses two or more layers of LSTM units, followed by a fully connected layer; The LSTM diffusion event recognition module is trained to perform joint classification and regression tasks. The classification task identifies the type of pollution event, and the regression task predicts the start and end time, peak arrival time, and peak amplitude of the event. The training dataset consists of historically labeled pollution event samples and samples synthesized by the mechanistic model; The loss function is designed as a weighted combination of classification cross-entropy loss, Huber regression loss, and physical consistency regularization term. The physical consistency regularization term is used to constrain the event recognition results to conform to the principle of mass conservation and the physical laws of pollutant propagation speed.

[0012] As a preferred approach, the S5 source tracing and contribution rate assessment employs a dual-channel inversion framework, which includes: The PDE-constrained optimization channel is used to construct an objective functional that includes a spatiotemporal weighted error term and apply L1 norm or total variation sparse regularization to the pollutant source term. The optimization problem is solved iteratively by calculating the gradient of the objective functional with respect to the source term through the adjoint equation method and combined with the optimization algorithm to jointly estimate the location, start-up and shutdown time and emission intensity of the pollution source. The Bayesian inference channel combines geographical prior information with the observation likelihood function. It samples the posterior probability distribution of pollution source parameters using the Markov chain Monte Carlo method or the sequential Monte Carlo algorithm, thereby outputting the posterior mean, variance, and high-density interval of the pollution source parameters.

[0013] As a preferred approach, the objective functional in the PDE-constrained optimization channel Defined as: ; in, These represent the simulated and observed concentrations at spatial location i and time t, respectively. Spatiotemporal weights; Denotes the L1 norm of the source term s. Indicates source item The total variation is used to apply sparse regularization; This indicates a penalty term based on prior physical knowledge; Gradient calculation is achieved through the adjoint equation, and the optimization algorithm uses L-BFGS or Adam.

[0014] As a preferred approach, the observation likelihood function in the Bayesian inference channel is characterized by a Gaussian or Student's t-distribution to represent the observation error; Geographic prior information includes known outlet locations, emission characteristics of specific industries, and pollution accessibility analysis results.

[0015] As a preferred method, S6 assimilation and parameter calibration specifically include: Data assimilation techniques such as ensemble Kalman filtering, ensemble Kalman smoothing, or particle filtering are employed. The assimilation process involves the simultaneous assimilation of state quantities such as water depth, flow velocity, and pollutant concentration. The assimilation process supports multi-frequency assimilation within a range of 5 to 60 minutes and can integrate multiple data sources from online monitoring, manual sampling, and remote sensing monitoring.

[0016] As a preferred approach, mechanistic equivalent substitution mechanisms are also included, specifically: When the long short-term memory network module or surrogate model is unavailable or untrustworthy due to data extrapolation, model drift, or algorithm permission not being met, the system automatically switches to a pure physical mechanism solution. The alternative is: event recognition is achieved by combining threshold judgment with cumulative sum or mutation test, and adaptive mesh refinement and time step splitting strategy is enabled after the proxy model is turned off; The entire process records data versions, parameter snapshots, model code verification, and key decision logs, forming an audit evidence chain.

[0017] The beneficial effects of this invention are as follows: By tightly coupling 1D / 2D hydrodynamics with multi-component pollutant degradation equations, the "black box" problem of traditional data-driven models is avoided, providing clear physical basis and logical chain. The inversion process incorporates prior physical knowledge and sparse regularization techniques, effectively mitigating the ill-conditioned nature of inverse problem solving, improving the accuracy and stability of pollution source identification, and making the source tracing results highly interpretable. By introducing an LSTM-based diffusion event identification and physical consistency proxy model, the number of computations of expensive mechanistic models required for inversion iterations is significantly reduced. Combined with a data assimilation link, online updates of system parameters (such as roughness and degradation rate) and state variables (such as water level, flow velocity, and concentration) are achieved, enabling hourly-level rapid updates in water body simulations with a scale of 100,000 to 1 million computational units, meeting the needs of real-time or near-real-time source tracing. By using sparse observation data and complex boundary conditions, this invention effectively identifies and distinguishes multiple coexisting pollution sources through a joint constraint mechanism of PDE constraint optimization and Bayesian posterior inference, making the estimation of source location, start-up and shutdown sequence and emission intensity more robust and reducing the risk of misjudgment caused by equivalent source effects. By introducing a mechanism-equivalent substitution strategy for AI modules, a seamless switch to a purely mechanism-driven backup solution can be achieved when the AI ​​model is unavailable or poses licensing / infringement risks. Simultaneously, the system retains data versions, parameter snapshots, model code checksums, and decision logs throughout the entire process, forming a complete and auditable chain of evidence, greatly facilitating regulatory review and evidence presentation for intellectual property compliance. The method of this invention supports the fusion of multiple pollutants, multiple water body forms (such as rivers, lakes, and estuaries), and multiple data sources (online, manual, and remote sensing), exhibiting good versatility and scalability. Its modular architecture makes it easy to integrate with existing SCADA (Supervisory and Data Acquisition) systems and other regulatory platforms, enhancing the flexibility of practical applications. Attached Figure Description

[0018] The invention will now be further described with reference to the accompanying drawings.

[0019] Figure 1 This is a schematic diagram of the overall process of a source tracing method based on coupled hydrodynamics and pollutant degradation equations according to the present invention.

[0020] Figure 2 This is a schematic diagram of the system architecture of a source tracing method based on coupled hydrodynamics and pollutant degradation equations according to the present invention. Detailed Implementation

[0021] To make the technical means, creative features, objectives and effects of this invention easier to understand, the invention will be further described below in conjunction with specific embodiments.

[0022] like Figure 1 and Figure 2 As shown in the figure, an embodiment of the present invention provides a source tracing method based on coupled hydrodynamics and pollutant degradation equations. S1: Engineering implementation of the data acquisition and preprocessing module: The data acquisition module aims to integrate multi-source heterogeneous monitoring data, and its implementation involves the integration of physical sensor networks, geographic information systems (GIS) and remote sensing platforms.

[0023] Specifically, the system collects data via distributed sensor nodes (e.g., integrated flow meters, water level gauges, velocity meters, and online water quality analyzers for DO, pH, conductivity, turbidity, COD, ammonia nitrogen, and total phosphorus) at a sampling frequency of 5 to 30 minutes, and uploads it in real time to a cloud data platform via LoRaWAN, NB-IoT, or 4G / 5G cellular networks. Geographic information data (including elevation model DEM, river cross-section data, and land use types) is integrated through a GIS platform and subjected to unified WGS84 or CGCS2000 coordinate system projection conversion. Meteorological data is obtained through an interface with the meteorological bureau, including hourly wind speed, wind direction, rainfall, temperature, relative humidity, and evaporation. Known discharge outlet lists, their geographic coordinates, design emission standards, and historical emission data are also included in the prior database.

[0024] In the data preprocessing stage, all collected data undergoes rigorous time calibration and coordinate system alignment. For missing data, this invention implements an adaptive interpolation strategy: For short-term missing data (<2 hours): linear interpolation or cubic spline interpolation should be used first. If the data is highly volatile, a local regression model based on signal decomposition (such as Empirical Mode Decomposition, EMD) should be used to fill the gaps in order to preserve the nonlinear and non-stationary characteristics of the signal.

[0025] For medium- to long-term missing data (2 hours to 24 hours): The system automatically selects either Kalman filtering or the Expectation-Maximization (EM) algorithm based on the data type. For data with significant dynamic processes, such as flow rate and water level, Kalman filtering is configured as a linear or nonlinear state-space model with physical constraints (the rate of change of water level must not exceed the historical maximum instantaneous rate of change) for real-time or offline prediction and imputation. Its process noise and observation noise covariance matrices are adaptively estimated based on historical data. For data influenced by multiple factors, such as pollutant concentration, the EM algorithm estimates missing values ​​by iteratively maximizing the likelihood function and can simultaneously estimate the model parameters, making it particularly suitable for scenarios with multivariate correlations.

[0026] Extreme missing data or sensor failure (>24 hours): In this case, the system will combine geospatial neighborhood (historical data from upstream or downstream adjacent monitoring stations) and machine learning methods (regression models based on random forest or XGBoost, using historical data as the training set) to fill in the missing data, while imposing strict physical constraints on the filling results.

[0027] Outlier identification and removal employs an improved interquartile range (IQR) method, typically marking data points outside the range [Q1−2.5×IQR,Q3+2.5×IQR] as outliers, where Q1Q1 and Q3Q3 are the first and third quartiles, respectively.

[0028] Furthermore, outlier identification also incorporates physical constraints: instantaneous flow velocity must not exceed the maximum empirical flow velocity of the river channel (3-5 m / s in mountainous rivers, 0.5-2 m / s in plain rivers); pollutant concentration must not be negative; and DO concentration must not exceed the saturated dissolved oxygen concentration. Identified outliers will be removed, replaced with nearby normal values, or marked as low-confidence data based on their characteristics and impact.

[0029] All preprocessed data are accompanied by detailed quality annotations, including metadata such as data source, sampling frequency, processing method, and confidence score. These standardized input datasets and prior parameter sets (topography, channel roughness, sediment chemistry, known pollution source locations, and historical emission ranges) are ultimately stored in a data lake in a unified structured format, providing a reliable and high-quality data foundation for subsequent hydrodynamic, reaction-migration, and source tracing inversion modules.

[0030] S2: In-depth analysis of the hydrodynamic solution module. The hydrodynamic solution module is the physical foundation of the entire source tracing method, and its accuracy directly affects the reliability of subsequent pollutant migration simulations. The hydrodynamic model of this invention has adaptive selection capabilities, and can intelligently switch according to the geometric characteristics and hydraulic properties of the target water body.

[0031] Implementation details of the one-dimensional Saint-Venant equation: For water bodies with significant one-dimensional flow characteristics, such as typical rivers and irrigation canals, this invention uses the one-dimensional Saint-Venant equation for simulation. In terms of numerical discretization, this invention employs the finite volume method (FVM), dividing the computational domain into a series of non-overlapping control volumes (river segments or computational units). Mass continuity equation. and momentum equation Integration is performed on each control volume. Flux calculations employ the HLL (Harten-Lax-van Leer) or HLLC (Harten-Lax-van Leer Contact) Riemann solver to effectively capture discontinuities such as flood waves and hydraulic jumps, and to ensure flux conservation and stability at steep slopes or variable cross sections. The HLLC solver handles water interfaces more accurately by introducing contact waves, making it particularly suitable for cases with complex bottom topography and wet-dry interfaces. Time propagation uses a second- or third-order Runge-Kutta (RK) method with higher-order total variation reduction (TVD), such as TVD-RK2, whose expression is: in, It represents a spatial discrete operator, which uses a limiter to ensure that the total variation does not increase, thereby effectively suppressing numerical oscillations.

[0032] For optimization in complex terrain and dynamic water levels: the source term balancing employs hydrostatic reconstruction technology to accurately balance the bed slope force and pressure terms, ensuring accurate solutions under still water surfaces or gentle slopes and avoiding the generation of spurious velocities. This technology is particularly critical because traditional methods may lead to the accumulation of numerical errors and instability under complex terrain.

[0033] Wet-Dry Boundary Treatment: A very small threshold water depth hmin (e.g., 1 × 10⁻³ m) is introduced to determine the wet / dry state of a cell. When the water depth of a cell is below hmin, it is considered a dry bed, and cross-boundary flux is restricted. When a wet cell is re-submerged, a local mass conservation algorithm is used to avoid instantaneous negative water depths.

[0034] Dry bed conservation treatment: After the water body has completely receded, a near-zero static water depth is maintained on the dry bed to ensure the conservation of mass and momentum in this area, preventing computational instability or spurious fluctuations. Boundary condition setting: The upstream boundary can be given by flow rate (e.g., reservoir discharge curve) or water level (e.g., lake or reservoir control water level); the downstream boundary can be given by water level (e.g., tidal curve, downstream dam control water level) or free outflow (e.g., processed by the characteristic line method). For artificially controlled facilities (e.g., dams, pumping stations), their operating curves (e.g., gate opening, pumping head and flow rate relationship) are integrated in real time as dynamic boundary conditions, coupled through hydraulic structure equations or equivalent resistance models. The Coriolis effect can be used as a source term in the momentum equation according to the river width and latitude. (especially for large-scale water bodies).

[0035] Refinement of the two-dimensional shallow water equation: For water bodies with significant two-dimensional flow characteristics, such as lakes, estuaries, floodplains, or wide river sections, this invention uses a depth-averaged two-dimensional shallow water equation for simulation. The FVM discretization, flux calculation (HLLC), and time progression (TVD-RK2 / RK3) strategies for the equation are consistent with those of the one-dimensional model.

[0036] Momentum source term Composition and calculation of bottom friction: The Manning formula or Chezy formula is usually used. The Manning formula is as follows: Where n is the Manning roughness coefficient, spatially distributed according to the type of sediment (e.g., sandy riverbed n≈0.025, gravelly riverbed n≈0.03−0.04, lush vegetation area n≈0.06−0.08) and water depth. The Xie Cai formula is... .

[0037] Bed slope force: directly derived from water depth gradient Components in the x, y directions.

[0038] Wind stress: in air density, drag coefficient , This is the wind speed vector.

[0039] Coriolis force: In a two-dimensional plane, it is expressed as fv in the x-direction and −fu in the y-direction.

[0040] Controlled boundaries: such as dams and pumping stations, whose effects are integrated into the grid cells through local momentum balance equations or equivalent resistance models.

[0041] Challenges and solutions for wetting-drying treatment in 2D scenarios: In 2D, the wetting-drying boundary is more complex, potentially involving irregular receding and flooded areas. This invention employs a cell-based mass conservation method: when the water depth in a cell is below a threshold, it is marked as a dry bed, and the flux between dry bed cells and wet cells is ensured to be zero. During flooding, a "water pushing" mechanism is used to inject water from neighboring wet cells into dry bed cells in a mass-conserving manner, preventing numerical collapse caused by excessively low local water depths.

[0042] S3: Implementation strategy of the reaction-migration coupling module: This module solves the hydrodynamic field (including water depth h, flow velocity) in module S2. Based on this, the spatiotemporal evolution of pollutants in water bodies is simulated. This invention couples and solves the multi-component convection-diffusion-reaction (ADR) equations. .

[0043] Numerical processing of the convection term: Convection term Left and right states Obtained through MUSCL reconstruction, such as Where ϕ is the limiter function. This represents the local gradient ratio.

[0044] Consideration of diffusion terms: diffusion terms This describes the diffusion and mixing processes of pollutants in water bodies. The dispersion coefficient tensor D considers longitudinal dispersion. Horizontal diffusion And turbulent diffusion. Among them , Typically, these parameters are estimated using empirical formulas based on hydraulic parameters (such as flow velocity, water depth, and Manning roughness), such as the Taylor-Aris dispersion theory or Elder's formula. For example, in rivers, It can be estimated as Where Q is the flow rate, W is the river width, H is the water depth, R is the hydraulic radius, U* is the shear velocity, and k is the KAMAN constant. The turbulent diffusion coefficient is calculated using the local turbulent kinetic energy (TKE) or shear velocity model.

[0045] Multicomponent reaction kinetic model: The reaction terms R(C,T,DO,…) are designed as a configurable multicomponent model, containing at least: First-order decay: −kC, where k is the first-order decay rate constant, which typically ranges from 0 to 1. The exact number depends on the type of pollutant (e.g., biodegradable organic matter, bacteria) and environmental conditions.

[0046] Arrhenius temperature correction: This is used to describe the change in reaction rate with water temperature t. The rate constant at 20°C This is a temperature correction factor, typically ranging from 1.02 to 1.09.

[0047] Dissolved oxygen (DO) limitation: For aerobic biochemical reactions, when the DO concentration is below a certain critical threshold... At concentrations of 1.0–2.0 mg / L, the reaction rate will be inhibited or even stopped to simulate biochemical processes under hypoxic conditions. This is typically achieved through a monod function or a step function.

[0048] Monod or multi-substrate kinetics: used to describe the dependence of microbial growth or pollutant degradation on the concentration of a single or multiple limiting substrates, such as... ,in For the maximum specific degradation rate, It is the half-saturation constant.

[0049] Adsorption-desorption: describes the distribution of pollutants between the aqueous phase and sediment or suspended particulate matter. This is typically achieved using linear partition coefficients Kd or more complex Freundlich and Langmuir adsorption isotherm models. For example, particulate concentration... in This represents the concentration in the aqueous phase.

[0050] Sedimentation-resuspension: describes the vertical migration process of particulate pollutants between the bottom of a water body and the water surface, quantified by sedimentation flux. and resuspension flux Composition, in which For the settling velocity, The resuspension rate is the resuspension rate.

[0051] Coupled solution strategy: To handle the strong nonlinear coupling between convection, diffusion, and reaction terms, this invention employs two strategies to ensure computational efficiency and stability: Operator splitting methods (such as Strang splitting) decompose the complex ADR equations into independent convection, diffusion, and reaction subproblems, which are then solved alternately. For example, within a time step, it can be decomposed as: convection → diffusion → reaction → diffusion → convection. This method uses the most suitable numerical algorithm for each substep and is easily parallelized, but it may introduce splitting errors.

[0052] Implicit Newton-Krylov methods solve all terms simultaneously within a single time step. This method linearizes the nonlinear system through Newton iterations and then uses Krylov subspace methods (such as GMRES or BiCGSTAB) to solve large sparse linear systems. Its advantages include high stability, especially for systems with "hard" reactive terms (i.e., characteristic timescales much smaller than the time step), but its computational cost is relatively high. This invention prioritizes this method when reactive terms are strong or time steps are large.

[0053] S4: In-depth analysis of the LSTM diffusion event recognition module To efficiently identify potential pollutant emission events, this invention introduces a diffusion event identification module based on Long Short-Term Memory (LSTM) networks. This module aims to automatically extract feature information of pollution events from multi-channel, high-frequency monitoring data.

[0054] Input data: This module takes as input a multi-index time series from at least 3 (usually 5-8) monitoring sections at equal intervals. These indicators typically include COD, ammonia nitrogen, total phosphorus, DO, pH, water temperature, flow rate, and water level, with a sampling frequency of 10 minutes. The input series usually uses a sliding window method, with the window length W adaptively adjusted between 30 and 180 time steps (corresponding to 5 hours to 30 hours), and the step size s selected between 1 and 5 time steps.

[0055] Network Architecture: The network structure employs two or more layers (e.g., three layers) of LSTM units, each containing 64 to 256 memory units, and introduces Dropout layers (e.g., Dropout rate 0.3) to prevent overfitting. The LSTM layers capture long-term dependencies in the time series. Following the LSTM layers are one or more fully connected layers (e.g., a ReLU activation layer containing 128 neurons), and the final output layer branches into two branches: Classification branch: The Softmax activation function is used to output the probability distribution of pollution event types (e.g., point source emissions, non-point source runoff, endogenous release, background fluctuations, etc.), which is a multi-classification task.

[0056] Regression branch: Employing the Sigmoid activation function (mapping the output to the [0,1] interval and then obtaining the actual value through a linear transformation) or direct output, this branch is used to predict the start and end times of events. Peak arrival time Peak amplitude Start and end times and peak times are usually expressed as offsets relative to the sliding window start time, while the peak amplitude is the standardized concentration increment.

[0057] Construction of the training dataset: The training dataset consists of two parts: Historically labeled pollution event samples: Real pollution events identified by experts and confirmed by laboratory analysis, including precise labels such as event type, start and end time, and emission intensity.

[0058] Mechanism Model Synthesis Samples: To address the problem of sparse real-world pollution event samples that are difficult to comprehensively cover all scenarios, this invention utilizes the forward mechanism model of modules S2 and S3 to synthesize a large number of pollution diffusion scenarios under different hydrodynamic conditions, pollutant types, emission locations, emission intensities, and time series. Corresponding monitoring section time series and precise pollution event labels are then generated for these scenarios. This approach significantly improves the model's generalization ability and its ability to identify rare events.

[0059] Loss function design: The loss function is designed as the classification cross-entropy loss. Huber regression loss and physical consistency regularization terms Weighted combination: .

[0060] Classification cross-entropy loss: used for classification tasks, its form is... .

[0061] Huber regression loss: Used for regression tasks, it is insensitive to outliers and takes the following form: ; Where e is the prediction error The threshold (e.g., 1.01.0).

[0062] Physical Consistency Regularization: This term is crucial for ensuring that the LSTM model output conforms to physical laws. It constrains the identified events to be consistent with the physical processes of the water body in terms of mass increment and propagation velocity. Specifically, it penalizes: a) significant deviations between the total mass of the identified events at the downstream monitoring section (obtained by integrating the concentration curve) and the total mass of the events identified at the upstream section; b) deviations of the predicted pollutant propagation velocity from the physically reasonable range of the mainstream average flow velocity calculated by the S2 module. For example, a penalty is imposed if the propagation velocity predicted by the LSTM exceeds the maximum flow velocity in the water body or falls below a certain lower dispersion limit.

[0063] Output: The output of the LSTM module includes: emission time windows of suspected external pollution sources (e.g., from...). The most likely pollutant type label (e.g., heavy metals, organic matter, nutrients) and preliminary location clues calculated based on the peak concentration time difference at different monitoring sections (e.g., the pollution source may be located between section A and section B, approximately XX kilometers from section A). These preliminary location clues will serve as prior information input to the S5 module to narrow down the inversion search space.

[0064] S5: The core algorithm and engineering implementation of the source tracing and contribution rate assessment module. This step is the core of the invention to solve the pathological nature of pollution source inversion. It adopts a dual-channel inversion framework to improve the accuracy, robustness and uncertainty quantification ability of the source tracing results.

[0065] PDE-constrained optimization channel: This channel treats pollution source terms (including location, time series, and intensity) as optimization variables and estimates source terms by minimizing the difference between simulation results and observation data.

[0066] Target functional Defined as: ; Error term: The simulated concentration at spatial location i and time t was quantified. With observed concentration The L2 norm difference between them. Spatiotemporal weights. The value can be dynamically adjusted based on the quality of the observation data, the uncertainty of the sensor (e.g., the weight of online monitoring data is higher than that of manually sampled data), and the level of attention to specific pollution events or areas (e.g., the weight of peak concentration periods or sensitive water areas can be set higher). The typical value is the reciprocal of the observation error variance.

[0067] Sparse regularization term: Used to apply sparsity regularization to suppress the occurrence of spurious sources, prompting the model to identify a limited number of primary pollution sources.

[0068] L1 norm Regularization coefficient: Encourages a majority of elements in the source term s to be zero, thus identifying a small number of active pollution sources. The typical range of values ​​is .

[0069] Total variation Encouraging source terms s to be piecewise constants or piecewise smooth in time and space effectively suppresses noise and makes the start-stop timing and spatial distribution of source terms clearer.

[0070] Physical Prior Penalty Term (Φ(Physical Prior)): This term incorporates known physical prior knowledge into the optimization process, further constraining the inversion results. For example, if the source term obtained from the inversion has a high intensity at a location where it is known that no discharge outlet could exist, or if the emission intensity exceeds the design emission capacity of a known discharge outlet, a penalty is applied. A specific implementation could be: in and These represent the maximum and minimum reasonable strengths of source j, respectively. , which is the penalty weight.

[0071] Gradient Calculation and Optimization Algorithm: The gradient of the objective functional with respect to the source term is calculated using the adjoint equation method. This method efficiently calculates the gradient of the objective function with respect to a large number of parameters by inversely integrating the adjoint equations of the original hydrodynamics and ADR equations, without needing to simulate perturbations for each parameter, thus significantly reducing computational costs. After obtaining the gradient, this invention combines optimization algorithms such as L-BFGS (Limited-memory Broyden–Fletcher–Goldfarb–Shano) or Adam for iterative solution. L-BFGS is a quasi-Newton algorithm suitable for large-scale nonlinear optimization problems; it accelerates convergence by approximating the Hesian matrix. Adam is an adaptive learning rate optimizer with lower requirements for hyperparameter tuning and faster convergence speed. The iterative convergence condition is typically set to a gradient norm less than... Or the maximum number of iterations (e.g., 200) can be reached. This method supports hybrid representations of discrete point sources (by applying the Diracdelta function to grid points) and distributed surface sources (by defining piecewise constant sources over multiple grid cells).

[0072] Bayesian Inference Channel: The Bayesian inference channel quantifies the uncertainty of pollution source parameters using a probabilistic framework. It combines geographical prior information with the observational likelihood function and samples the posterior probability distribution of pollution source parameters using the Markov Chain Monte Carlo (MCMC) method or the Sequential Monte Carlo (SMC) algorithm.

[0073] Likelihood Function: The observation likelihood function P(D|M) uses either a Gaussian distribution or a Student's t-distribution to characterize the observation error. The Gaussian distribution assumes that the error follows a normal distribution, while the Student's t-distribution is more robust to outliers in the observation data, and its degrees of freedom parameter can be adjusted according to the data characteristics. For example, .

[0074] Prior information: Geographic prior information (such as preliminary location clues provided by the S4 module, probability distributions of known discharge outlet locations, emission characteristic ranges of specific industries, and pollution accessibility analysis results—i.e., whether pollutants can reach the observation point from a certain location via water flow) is integrated with historical emission records to form a prior probability distribution P(M) of the pollution source parameters. For example, for source location, a Gaussian prior centered on known discharge outlets can be set, or a binary prior based on river connectivity can be set. For emission intensity, a log-normal prior based on industry statistics or historical events can be set.

[0075] Sampling Algorithm: Delayed Rejection Adaptive MCMC (DRAM): This algorithm can adaptively adjust the covariance matrix of the proposal distribution to improve sampling efficiency, and after the first rejection of a proposal, it will try again with a proposal with a smaller perturbation, thereby reducing unnecessary rejections and improving the acceptance rate.

[0076] Preconditional CoherentNoiseMCMC (pCN-MCMC): Applicable to high-dimensional parameter spaces, it accelerates convergence by introducing a precondition matrix.

[0077] The Sequential Monte Carlo (SMC) algorithm extracts samples from complex posterior distributions through a series of importance sampling and resampling steps, making it particularly suitable for non-stationary or multi-peak posterior distributions. Using these methods, the system generates posterior probability distributions of pollution source parameters (location, start / stop time, emission intensity) and outputs the posterior mean, variance, and high-density interval (HPD), thereby quantifying the uncertainty of the inversion results and providing a confidence assessment of the source tracing results.

[0078] Contribution rate assessment: After obtaining pollution source information, the contribution rate assessment module is used to calculate the real-time time-varying contribution rate of each candidate pollution source to the target section or area. This invention provides three optional schemes: Sensitivity superposition method: Based on the assumption of a linearized model, this method estimates the contribution rate by calculating the sensitivity matrix of emission changes from each pollution source at different times to the target section concentration, and then superimposing the effects of different sources. This method is computationally efficient, but requires the system response to be approximately linear.

[0079] Single-source perturbation experiment: By activating or deactivating a pollution source one by one (or reducing its emission intensity to zero), and running a forward model, the marginal contribution of that source to the total pollution load is evaluated. This method is more accurate in nonlinear systems, but it is computationally expensive (requiring N+1 forward model runs, where n is the number of pollution sources).

[0080] The Shapley value attribution method, derived from cooperative game theory, assigns a "fair" contribution to the total pollution load for each pollution source. The Shapley value considers the marginal contributions under all possible combinations of pollution sources, handles nonlinear interactions between sources, and ensures that the sum of the contribution rates equals the total load. Although it has the highest computational cost (theoretically requiring 2N²N simulations, but in practice, the Monte Carlo approximation is usually used), it provides the most impartial attribution results. Ultimately, the system outputs detailed source attribution conclusions, including the location coordinates (longitude, latitude) of each pollution source, emission start and stop times, emission intensity curves (e.g., pollutant mass flow rate per hour, in kg / h), and real-time time-varying contribution rate curves and contribution decomposition pie charts for peak pollution periods for each source to a specific monitoring section in the receiving water body.

[0081] S6: Online update mechanism for assimilation and parameter calibration module: In order to realize online updating of the model and improve prediction accuracy, this invention adopts advanced data assimilation technology to correct key hydrodynamic parameters, pollutant degradation rate and boundary condition error terms online.

[0082] Assimilation technique selection: This invention supports data assimilation techniques such as ensemble Kalman filtering (EnKF), ensemble Kalman smoothing (EnKS), or particle filtering.

[0083] Ensemble Kalman Filtering (EnKF): This method propagates uncertainty and assimilates states and parameters by generating a "set" of model states and updating each member of the set using observation data. EnKF does not require computation of adjoint equations and is suitable for large nonlinear systems.

[0084] Ensemble Kalman Smoothing (EnKS): An extension of EnKF, it can use observational data from a future period to correct historical states and parameters, thereby providing a more accurate estimate of the "past" model state.

[0085] Particle filtering: It represents the posterior probability distribution using a set of "particles" (random samples of model states and parameters) and updates it through importance sampling and resampling processes. Particle filtering is highly applicable to non-Gaussian and nonlinear systems, but it has high computational costs, especially in high-dimensional state spaces.

[0086] Assimilation parameters and state variables: The assimilation process can simultaneously assimilate state variables such as water depth (h), flow velocity (u / v), and pollutant concentration (c), as well as Manning roughness n (typically ranging from 0.015 to 0.08), longitudinal and transverse dispersion coefficients, etc. (Typical range) Pollutant degradation rate kk (typical range) Key parameters, etc.

[0087] Multi-frequency and multi-source data fusion: This invention supports multi-frequency assimilation within a range of 5 to 60 minutes to adapt to sampling intervals from different data sources. For example, online monitoring data can be assimilated every 5 minutes, while manually sampled data is batch-processed and assimilated immediately after acquisition. The system can fuse multiple data sources, including online monitoring (high frequency, low precision), manual sampling analysis (low frequency, high precision), and remote sensing monitoring (large scale, medium precision, such as chlorophyll a and suspended matter concentration), achieving optimal fusion by setting different observation error covariance matrices for different data sources.

[0088] S7: Implementation details of the uncertainty quantification and verification module: In order to comprehensively evaluate the reliability of the traceability results, this invention quantifies and propagates the inherent uncertainty in the model.

[0089] Uncertainty propagation method: This invention uses Latin hypercube sampling (LHS) or Monte Carlo perturbation propagation to address uncertainties in key parameters and boundary conditions.

[0090] Parameter perturbations: Manning roughness n, longitudinal and transverse dispersion coefficients The values ​​of key parameters such as degradation rate k are often uncertain. This invention addresses this by setting a probability distribution for these parameters (e.g., uniform distribution, normal distribution, log-normal distribution) and performing extensive sampling (e.g., 1000 times). For example, the Manning roughness n can be perturbed by ±20% from a baseline value determined based on riverbed sediment surveys.

[0091] Boundary condition errors: Observational errors or prediction uncertainties also exist in boundary conditions such as upstream flow rate and downstream water level. This invention simulates this uncertainty by superimposing random noise (e.g., following a Gaussian distribution with zero mean and variance determined by historical errors) onto the given boundary conditions.

[0092] Model structure uncertainty: For simplified models or physical processes that are not fully considered, this invention can represent the uncertainty of the model structure by introducing an error term or parameterizing a subgrid effect and perturbing it.

[0093] Output: By propagating these uncertainties, the system runs multiple sets of forward simulations and source inversions, ultimately generating the probability distribution of the source tracing results (source location, intensity, contribution rate), and outputting the confidence interval (e.g., 95% confidence interval) and robustness score (e.g., the ratio of confidence interval width to posterior mean).

[0094] Model validation: Cross-validation and posterior validation using historical pollution events are employed to ensure the model's predictive performance under different scenarios. Commonly used validation metrics include root mean square error (RMSE), Nash efficiency coefficient (NSE), and coefficient of determination (R²²).

[0095] One-dimensional to two-dimensional adaptive coupling interface: This invention employs advanced adaptive coupling interface technology at the boundaries of different dimensional model regions such as river networks, lake areas, and estuaries to ensure the continuity of mass and momentum. This interface can dynamically switch between strong and weak coupling based on the actual physical characteristics of the water body connection.

[0096] Strongly coupled methods: When significant hydraulic interactions exist at the interface (e.g., narrow openings where rivers flow into lakes, or floodplains during floods), an interface solution based on the Riemannian problem is employed. Specifically, the cells at the coupling boundary are treated as virtual control volumes, and their flux exchanges are calculated using a local Riemannian problem solver, ensuring precise conservation of mass and momentum on both sides of the interface. This requires consistent time steps between models of different dimensions, typically achieved by finding the greatest common divisor time step or employing asynchronous parallel computation and data interpolation strategies.

[0097] Weak coupling method: When hydraulic interactions are relatively weak or higher computational efficiency is required (e.g., two large water bodies connected by a sluice gate), a weak coupling method based on water level-discharge characteristic curves is adopted. In this method, the outflow of one model serves as the inflow boundary condition for another model, and vice versa. Characteristic curves (e.g., QH relationship) can be generated independently from simulation or measured data before coupling. This method allows models of different dimensions to use their own independent time steps, but requires additional iteration or filtering steps to ensure long-term stability. This coupling strategy ensures uniform simulation in complex water systems, effectively avoiding computational errors and instabilities caused by dimensional mismatch.

[0098] This invention cleverly integrates neural network models and introduces a "mechanism-equivalent substitution" strategy to ensure compliance and reliability.

[0099] Detailed Explanation of LSTM Parameters and Output: The adaptive adjustment of the sliding window length (WW) in the LSTM event recognition module is typically based on the statistical distribution of historical event durations. A smaller W is chosen if most events have short durations, and a larger W is chosen otherwise. The choice of step size (s) balances computational efficiency and the timeliness of event capture. The Adam optimizer is used during network training, and the learning rate is typically set to... This is combined with a learning rate decay strategy (e.g., decreasing by 0.1 every 50 epochs). The training batch size is typically 32-128. The output parameters precisely include the start and end times of the events. Peak arrival time Peak amplitude And the corresponding event type label. These time parameters are expressed in seconds relative to the start of the sliding window, and the peak amplitude is the normalized increment of pollutant concentration.

[0100] Construction and application of proxy models: This invention uses convolutional long short-term memory networks (ConvLSTM), U-Net, or physical information neural networks (PIn) to construct proxy models.

[0101] ConvLSTM introduces convolutional operations into LSTM units, enabling them to process both spatial and temporal features simultaneously, making them particularly suitable for spatiotemporal sequence data such as pollutant concentration fields. The ConvLSTM model learns the local mapping relationship of pollutant concentration (cc) under given water depth and flow velocity (u,v,hu,v,h). Its architecture consists of multiple ConvLSTM layers followed by convolutional layers for output prediction.

[0102] U-Net: A classic encoder-decoder architecture, adept at image-to-image conversion tasks. Here, it is trained to learn the mapping from hydrodynamic inputs (e.g., flow field, water depth field) to pollutant concentration fields (as "images"), making it particularly suitable for concentration field prediction in two-dimensional water bodies.

[0103] Physical Information Neural Network (PIn): By incorporating the residual terms of the ADR equation into the loss function of the neural network, PIn forces its predictions to conform to physical laws during training. This allows PIn to generate predictions with high physical consistency even in sparse data situations and effectively avoids non-physical phenomena such as negative concentrations.

[0104] These surrogate models are primarily used for short-term (e.g., less than or equal to 6 hours) concentration field predictions, or to replace computationally expensive ADR solutions in the inner loop of the inversion algorithm, thereby significantly improving computational efficiency. For example, in the S5 inversion iteration, the surrogate model can be used for rapid prediction in the first 10-20 iterations. When the iteration approaches convergence or the surrogate model prediction deviates from the preset physical constraint threshold (e.g., negative concentrations occur, pollutant mass non-conservation error exceeds 5%), the system will automatically revert to the accurate mechanistic model solver to avoid error accumulation and result drift.

[0105] Mechanism Equivalent Substitution and Compliance Strategy: This invention proposes a "mechanism equivalent substitution" mechanism to address the limitations or compliance requirements of AI models in specific contexts. When an AI module (such as an LSTM event recognition or proxy model) becomes unavailable or untrusted due to reasons such as data extrapolation (i.e., encountering extreme scenarios outside the training data range), model drift (model performance deteriorates over time), or unmet third-party algorithm licensing requirements, the system will automatically and seamlessly switch to a purely physical mechanism-driven equivalent substitution solution.

[0106] Alternatives for event recognition: When the LSTM event recognition module fails, the event recognition function will be implemented using multi-threshold criteria (e.g., a sudden increase in concentration exceeding three standard deviations from the historical mean and lasting for more than 30 minutes), combined with a cumulative sum (CUSUM) algorithm (for detecting small but persistent changes in time series) or abrupt change detection algorithms (e.g., sliding t-test or Lepage test). These methods are based on statistical principles and physical thresholds, offering high interpretability and auditability.

[0107] Replacement of surrogate models: Surrogate models will be disabled, and the system will automatically enable adaptive mesh refinement and time-step splitting strategies. For example, when a steep concentration gradient or drastic flow velocity change is detected, the local mesh will be dynamically refined, and the time step will be dynamically reduced according to the CFL number to complete the simulation without sacrificing accuracy. At the same time, more robust but computationally more expensive numerical schemes (e.g., implicit solvers) will be enabled.

[0108] Construction of an auditable chain of evidence: To meet stringent regulatory and auditing requirements, the system of this invention records data versions (through data checksums and version numbers), parameter snapshots (MD5 hashes of model parameter files), model code checksums (Git version control and code repository hashes), and key decision logs (including the enabled / disabled status of AI modules, triggering conditions and switching times for alternative solutions, and user intervention records) throughout the entire traceability analysis chain. These records constitute a complete auditable chain of evidence, stored in a distributed ledger or timestamp service, which can be used to support regulatory audits and provide legal basis for potential algorithm licensing or infringement disputes.

[0109] Example 1: Pollutant Source Tracing in a One-Dimensional River Scenario This embodiment focuses on a narrow river channel located in the mountainous region of southern China, approximately 85 kilometers long, with an average width of 30 meters and an average depth of 3.5 meters. The river is regulated by a reservoir upstream and has an ecological control gate downstream. Six online monitoring sections (M1 to M6) were established along the river channel, with M1 located approximately 5 kilometers downstream of the reservoir and M6 approximately 3 kilometers upstream of the ecological control gate. Each section was equipped with a high-precision flow meter, water level gauge, and multi-parameter water quality analyzer (COD, ammonia nitrogen, total phosphorus, DO), with a sampling frequency of 15 minutes.

[0110] Scenario Description: From 14:00 to 18:00 on July 15, 2023, an abnormal increase in COD concentration was observed at monitoring section M3 (approximately 30 km from M1), rising from a baseline of approximately 15 mg / L to a peak of 120 mg / L, before gradually decreasing. A similar increase in COD was observed at monitoring section M4 (approximately 45 km from M1) from 17:00 to 21:00 on July 15, 2023, with a peak of approximately 95 mg / L.

[0111] (1) Data Acquisition and Preprocessing (S1): Collect water level, flow rate, COD, and water temperature data from sections M1 to M6. Obtain rainfall and evaporation data for the area through an interface with the meteorological bureau. Extract river channel DEM data and roughness distribution maps from the geographic information database. Perform timestamp calibration and outlier detection on all data. For example, the flow rate data of section M2 at 16:30 on July 15 was missing due to a momentary sensor malfunction. The system used Kalman filtering to predict and fill in the missing data, with the filled value being... The confidence level is marked as high.

[0112] (2) Hydrodynamic Solution (S2): Since the river width is much smaller than its length, this invention automatically selects the one-dimensional Saint-Venant equation for hydrodynamic simulation. The model mesh is divided into 250 elements, each approximately 340 meters long. The upstream M1 section is given the reservoir discharge curve (e.g., 30 m³ / s), and the downstream M6 section is given the control water level of the ecological gate (e.g., 50.2 m). The Manning roughness n is a piecewise constant, such as n=0.035 for the upstream section and n=0.030 for the mid-to-lower reaches. The time step is dynamically adjusted to ensure that the CFL number is around 0.75.

[0113] (3) Reaction-migration coupling (S3): Based on the solved hydrodynamic field, the migration and degradation of COD are simulated. A first-order decay kinetic model is used for COD, and the decay rate is... The temperature correction factor θ = 1.05. The longitudinal dispersion factor Dl is estimated to be 50 m² / s. The convection term is reconstructed using a VanLeer limiter, with the time step consistent with the hydrodynamic model.

[0114] (4) LSTM Dispersion Event Identification (S4): The LSTM module uses a 15-minute COD time series from sections M1 to M6 as input, with a sliding window length W=60 (15 hours) and a step size s=2. The network identifies a potential pollution event between sections M3 and M4. The LSTM outputs the start time of the event as 13:50±10min on July 15, 2023, the stop time as 18:10±15min on July 15, 2023, the pollution type as organic matter, and the preliminary location clues pointing to a location approximately 10-15 km downstream of section M3.

[0115] (5) Source tracing and contribution rate assessment (S5): Using the preliminary clues provided by S4, the source tracing search area is limited to the range of 10 kilometers downstream of M3 to 5 kilometers upstream of M4.

[0116] PDE-constrained optimization channel: Construct an objective functional, where the error term weights wi,t are set as the reciprocal of the variance of COD observation errors at each cross-section (e.g., online monitoring error). .

[0117] TV regularization coefficient The physical prior penalty term Φ imposes a high penalty on sources located more than 100 meters from the riverbank. The optimization algorithm uses L-BFGS with a maximum of 150 iterations.

[0118] Bayesian inference channel: The likelihood function adopts a Gaussian distribution. Prior information is fused from: a) the location clues identified by LSTM as Gaussian prior centers; b) the known discharge characteristics (COD emission intensity range of 50–500 kg / h) of two small dyeing and printing factories in this river section as intensity priors. DRAM-MCMC is used for sampling to generate the posterior distribution of source parameters. Through dual-channel inversion, the system identifies a single major pollution source located approximately 12.5 km downstream of M3, with its emission start-stop time from 14:15 to 17:45 on July 15, 2023, and a peak emission intensity of 180 kg / h. This source contributes 92% to the increase in COD at section M4.

[0119] (6) Assimilation and Parameter Calibration (S6): EnKF data was used for assimilation, with flow rate, water level, and COD concentration at sections M1 to M6 assimilated hourly. The assimilation also updated the Manning roughness n (within ±10% of the baseline value) and the COD decay rate k (within ±20% of the baseline value). Through assimilation, the model's prediction error for COD concentration over the next 24 hours was reduced by 15%.

[0120] (7) Uncertainty Quantification and Validation (S7): Latin hypercube sampling was performed on the Manning roughness n and COD decay rate k (1000 times each), with a perturbation range of ±15%. By propagating these uncertainties, the 95% confidence interval for identifying the source location was calculated to be downstream of M3 [11.8 km, 13.2 km], and the 95% confidence interval for the peak emission intensity was [165 kg / h, 195 kg / h]. Cross-validation results of historical pollution events showed that the model's source tracing accuracy RMSE for similar events was 0.7 km (location) and 8% (intensity).

[0121] Example 2: Pollutant Source Tracing and Comparison in a Two-Dimensional Lake Scene This embodiment is applied to a shallow lake with an area of ​​approximately 100 km³. The lake has multiple rivers flowing into it, a buoy station is deployed in the center of the lake, and there are several aquaculture farms and industrial wastewater outlets around the lake. The monitoring data includes 30-minute online data of DO, COD, total nitrogen (TN), and total phosphorus (TP) from the buoy station, as well as weekly sampling data near several wastewater outlets.

[0122] Scenario description: From 06:00 to 18:00 on August 20, 2023, the COD concentration was observed to rise abnormally from 20 mg / L to a peak of 85 mg / L at the lake center buoy station.

[0123] (1) Data acquisition and preprocessing (S1): Collect data on lake topography, sediment type, wind field, and inflow river flow. Process missing and abnormal data.

[0124] (2) Hydrodynamic Solution (S2): Since the lake is a two-dimensional water body, this invention automatically selects a depth-averaged two-dimensional shallow water equation. The computational grid uses an unstructured triangular grid with an average cell area of ​​approximately 200 m², totaling 500,000 cells. Boundary conditions include a given flow rate of the inflowing river and a given water level at the lake outlet. Wind stress is introduced into the momentum equation.

[0125] (3) Reaction-migration coupling (S3): COD adopts first-order decay, decay rate .

[0126] (4) LSTM diffusion event identification (S4): The LSTM module uses 30 minutes of water quality data from the lake center buoy station as input to identify the COD rise event in the lake center. The initial clues point to the southwest direction of the lake, and the duration is about 10 hours.

[0127] (5) Source tracing and contribution rate assessment (S5): The search area is set on the southwest side of the lake. There are three known potential pollution sources in this area (P1: discharge outlet of large-scale farm, P2: direct discharge outlet of small-scale industrial wastewater, P3: overflow outlet of nearby urban sewage treatment plant).

[0128] Invention scheme (dual-channel inversion): PDE constraint optimization parameter settings: Physical prior penalties apply to sources located more than 200 meters from the existing discharge outlet location.

[0129] Bayesian prior inference: The emission intensities of P1, P2, and P3 each follow a log-normal distribution, with the mean and variance determined based on historical monitoring data and industry emission standards. For example, the prior mean of COD intensity for P1 is 150 kg / h, for P2 it is 300 kg / h, and for P3 it is 200 kg / h.

[0130] Source tracing results: P2 was identified as the main pollution source, with emissions starting and stopping between 06:30 and 17:00 on August 20, 2023, and a peak emission intensity of 380 kg / h. P1 and P3 had normal emission intensities during this event, with extremely low contribution rates. P2 contributed 85% to the peak COD level at the lake's central monitoring point.

[0131] The comparative scheme (single data-driven inversion) employs a backpropagation method based on a traditional artificial neural network (An), directly mapping monitoring data to source parameters without introducing mechanistic models or physical constraints, nor quantifying uncertainties. The network training data primarily comes from historical events and a small amount of mechanistic model simulation data.

[0132] Source tracing results: Three potential sources were identified, but there is significant uncertainty regarding the emission timing and intensity, and it is impossible to distinguish between primary and secondary sources. The An inversion yielded average contribution rates of 30%, 40%, and 30% for P1, P2, and P3, respectively, which deviates from the actual situation (P2 is the primary source), and it cannot provide emission timing information.

[0133] Table 1: Comparison of Source Tracing Results for Two-Dimensional Lake Scenes Identification of major pollution sources P2 (Small-scale industrial wastewater direct discharge outlet) Unclear, 3 potential sources identified Source position accuracy 25m from actual location 95% CI: 15m-35m 150m from actual location Emission start-stop time 2023 / 08 / 2006:30-17:00(95%CI:±30min) Specific start-stop times were not provided; only average emission periods were estimated. Peak emission intensity 380kg / h 95% CI: 360kg / h-405kg / h 350±120 kg / h (high uncertainty) Peak COD contribution rate of the lake center P2:85% P1: 30%, P2: 40%, P3: 30% (This does not reflect the actual situation) Quantification of uncertainty Provides 95% confidence intervals (CI) and posterior distributions. Not provided Explainability Based on the physical mechanism model, the results have clear physical meaning. The "black box" model makes the results difficult to interpret. computational efficiency The inversion process took approximately 3 hours (accelerated using a surrogate model). The inversion process takes approximately 1.5 hours (no mechanistic model required, but accuracy is low). As can be seen from the comparison results in Table 1, the proposed method significantly outperforms traditional single-data-driven inversion methods in terms of pollution source identification accuracy, emission characteristic estimation, uncertainty quantification, and interpretability. Although the computational efficiency is slightly lower than that of the pure data-driven model, the efficiency has been greatly improved by introducing a surrogate model, and the reliability and credibility of the results are higher.

[0134] In practical deployment, this invention will be designed with a microservice architecture, with each module (data preprocessing, hydrodynamics, reaction-transfer, LSTM identification, source tracing and inversion, assimilation, and verification) deployed as an independent service in a containerized environment (such as Docker / Kubernetes). This not only improves the modularity and scalability of the system but also facilitates independent upgrades and maintenance of each module. The system front-end provides a web-based user interface for visually displaying monitoring data, simulation results, source tracing conclusions, and contribution rate decomposition. The back-end API supports data exchange and functional integration with external SCADA (Supervisory and Data Acquisition) systems, environmental management platforms, and early warning systems. Data storage utilizes distributed file systems (such as Ceph / HDFS) and time-series databases (such as InfluxDB / TimescaleDB) to efficiently process and store massive amounts of spatiotemporal monitoring data and model outputs. The entire system architecture ensures high performance, high availability, and easy scalability, meeting the needs of water environment pollution source tracing of varying scales and complexities.

[0135] In summary, this invention provides a comprehensive, accurate, efficient, and auditable solution for tracing the source of water pollutants by tightly coupling a mechanistic model with advanced data intelligence technology and supplementing it with rigorous engineering implementation and compliance assurance mechanisms.

Claims

1. A source tracing method based on coupled hydrodynamics and pollutant degradation equations, comprising the following steps: S1. Data acquisition and preprocessing: This is used to acquire and integrate multi-source heterogeneous monitoring data and perform standardization processing. The processing includes: filling in missing observation data using interpolation algorithms, and identifying and removing outliers using the interquartile range method combined with physical constraints. S2. Hydrodynamic solution is used to accurately simulate the water flow state of the target water body. The simulation adopts an adaptive model selection strategy, and solves the equations using the one-dimensional Saint-Venant equation or the two-dimensional shallow water equation according to the water body morphology. The finite volume method is used to discretize the governing equations, and the Runge-Kutta method with high-order total variation reduction is used for time advancement. S3, Reaction-Migration Coupling, is used to simulate the migration and transformation of pollutants in water bodies based on the solved hydrodynamic field. It is achieved by solving multi-component convection-diffusion-reaction equations through coupling. S4, LSTM diffusion event identification, is used to introduce a diffusion event identification module based on long short-term memory network to identify emission events of suspected external pollution sources from multi-index time series of multiple monitoring sections. The module is trained to perform joint classification and regression tasks and output emission time window, pollutant type and preliminary location clues; the preliminary location clues include the possible location range of pollution sources calculated based on the peak concentration time difference of different monitoring sections. S5, Source Inversion and Contribution Rate Assessment, is used to estimate the location, start-up and shutdown time and emission intensity of pollutant source terms based on the dual-channel inversion framework, quantify the uncertainty of the inversion results, and calculate the real-time time-varying contribution rate of each candidate pollution source to the target section or region through the contribution rate assessment module. S6, Assimilation and Parameter Calibration, is used to update key hydrodynamic parameters, pollutant degradation rates, and boundary condition error terms online using data assimilation technology, thereby achieving online model updates and improving prediction accuracy; S7. Uncertainty Quantification and Verification: This function is used to propagate perturbations to uncertainties in model parameters, boundary conditions, and model structure, thereby outputting the confidence interval and robustness score of the source tracing results, and combining them with historical pollution events for posterior verification. The S5 source tracing and contribution rate assessment adopts a dual-channel inversion framework, which includes: The PDE-constrained optimization channel is used to construct an objective functional that includes a spatiotemporal weighted error term and apply L1 norm or total variation sparse regularization to the pollutant source term. The optimization problem is solved iteratively by calculating the gradient of the objective functional with respect to the source term through the adjoint equation method and combined with the optimization algorithm to jointly estimate the location, start-up and shutdown time and emission intensity of the pollution source. The Bayesian inference channel combines geographical prior information with the observation likelihood function. It samples the posterior probability distribution of pollution source parameters using the Markov chain Monte Carlo method or the sequential Monte Carlo algorithm, thereby outputting the posterior mean, variance, and high-density interval of the pollution source parameters.

2. The source tracing method based on coupled hydrodynamics and pollutant degradation equations according to claim 1, characterized in that, S1 data acquisition and preprocessing specifically include: Collect hydrological data, hydraulic data, meteorological data, geographic information data, sediment characteristic parameters, a list of known discharge outlets and their locations, as well as various water quality data; All collected data underwent unified calibration of coordinate system and time reference. For missing observation data, methods based on Kalman filtering, expectation-maximization algorithm or neighborhood regression are used to fill in the gaps; For outliers, the interquartile range method is used in combination with physical constraints of maximum flow velocity and minimum water depth for identification and removal. After alignment, the multi-source data is quality-labeled to form a standardized input dataset and a set of prior parameters containing geographic, hydrological, and water quality information.

3. A source tracing method based on coupled hydrodynamics and pollutant degradation equations according to claim 1, characterized in that, The S2 hydrodynamic solution specifically includes: For water bodies with significant one-dimensional characteristics, the one-dimensional Saint-Venant equation is used for simulation. The one-dimensional Saint-Venant equation includes the mass continuity equation and the momentum equation. For water bodies with significant two-dimensional characteristics, the solution is obtained through the depth-averaged two-dimensional shallow water equation, which includes the mass conservation equation, the momentum conservation equation in the x-direction, and the momentum conservation equation in the y-direction. In terms of numerical implementation, the finite volume method is used to discretize the governing equations, and the flux is calculated using an HLL or HLLC Riemann solver. The time progression employs a second- or third-order Runge-Kutta method to reduce higher-order total variation; The source term balance, wet-dry boundary treatment, and dry bed conservation treatment were optimized. Boundary conditions are set including upstream given flow rate or water level, downstream given water level or free outflow, and the effects of wind stress, rainfall, evaporation and tidal harmonics or measured tidal sequences are taken into account.

4. A source tracing method based on coupled hydrodynamics and pollutant degradation equations according to claim 1, characterized in that, The objective functional J in the PDE-constrained optimization channel is defined as follows: ； in, These represent the simulated and observed concentrations at spatial location i and time t, respectively. Spatiotemporal weights; Denotes the L1 norm of the source term s. represents the total variation of the source term s, used to apply sparse regularization; This indicates a penalty term based on prior physical knowledge; Gradient calculation is achieved through the adjoint equation.

5. A source tracing method based on coupled hydrodynamics and pollutant degradation equations according to claim 4, characterized in that, The observation likelihood function in the Bayesian inference channel is characterized by a Gaussian or Student's t-distribution to represent the observation error; Geographic prior information includes known outlet locations, emission characteristics of specific industries, and pollution accessibility analysis results.

6. A source tracing method based on coupled hydrodynamics and pollutant degradation equations according to claim 1, characterized in that, S6 assimilation and parameter calibration specifically include: Data assimilation techniques such as ensemble Kalman filtering, ensemble Kalman smoothing, or particle filtering are employed. The assimilation process involves the simultaneous assimilation of state quantities such as water depth, flow velocity, and pollutant concentration. The assimilation process supports multi-frequency assimilation within a range of 5 to 60 minutes and can integrate multiple data sources from online monitoring, manual sampling, and remote sensing monitoring.

7. A source tracing method based on coupled hydrodynamics and pollutant degradation equations according to claim 1, characterized in that, It also includes mechanistic equivalent substitution mechanisms, specifically including: When the long short-term memory network module or surrogate model is unavailable or untrustworthy due to data extrapolation, model drift, or algorithm permission not being met, the system automatically switches to a pure physical mechanism solution. The alternative is: event recognition is achieved by combining threshold judgment with cumulative sum or mutation test, and adaptive mesh refinement and time step splitting strategy is enabled after the proxy model is turned off; The entire process records data versions, parameter snapshots, model code verification, and key decision logs, forming an audit evidence chain.

Images