Methods for Simulating and Predicting Heavy Metal Concentration in Water

Abstract

This invention discloses a method for simulating and predicting heavy metal concentrations in water, comprising: integrating raw monitoring data, hydrodynamic data, total suspended solids, image remote sensing data, and human activity data to generate a multi-source cleaning sequence data package; using the image remote sensing data in the data package to perform cross-modal adsorption capacity estimation, inferring particle chemical composition and adsorption isotherm parameters from image texture, and generating a capacity feature package containing an upper bound of adsorption capacity; calculating time-varying travel times based on hydrodynamic data and human activity data, causally aligning the capacity feature package and upstream signals to generate a travel time aligned feature package; and combining metal fingerprint parameters, applying the upper bound of adsorption capacity as a physical constraint in the morphology allocation constraint header, explicitly decoupling and predicting its morphology, and generating a prediction result package. This invention deeply couples hydrodynamic physical mechanisms with particle adsorption chemical mechanisms, improving the prediction accuracy and physical consistency of the model under unsteady-state conditions.

Description

Technical Field

[0001] This invention belongs to the field of water environment monitoring technology, and in particular, it is a method for simulating and predicting the concentration of heavy metals in water. Background Technology

[0002] Heavy metal pollution is one of the key environmental problems threatening the safety of aquatic ecosystems. Due to their high toxicity, recalcitrant nature, and bioaccumulation, the migration and transformation of heavy metals in water bodies are subject to complex influences from multiple factors, including hydrological, meteorological, and biochemical factors. Therefore, developing technologies that can dynamically simulate and accurately predict heavy metal concentrations in water bodies can enable timely early warning of pollution incidents, assessment of water environmental capacity, and watershed risk management.

[0003] Currently, various technical approaches have been developed for simulating heavy metal concentrations in water. One type is mechanistic models based on physicochemical processes, such as water migration and diffusion models. These models characterize the spatiotemporal evolution of pollutants by solving convection-diffusion equations, typically requiring detailed river network topology, cross-sectional parameters, and accurate boundary conditions. Another type is data-driven statistical or machine learning models. For example, autoregressive integrated moving average (ARIMA) models or long short-term memory networks (LSTM) are used to predict the historical concentration series of a single monitoring section over time; or spatial interpolation techniques (such as Kriging) are used to estimate the spatial distribution of concentrations at different sections at the same time. In addition, some studies have begun to explore the integration of multi-source data such as water quality, water quantity, and meteorological data, using feature stitching to train models and improve prediction accuracy.

[0004] However, existing technologies still suffer from limitations in simulating unsteady hydrodynamic and complex biochemical coupling processes, including a lack of real-time chemical mechanism constraints and difficulty in achieving causal alignment of upstream and downstream signals under dynamic flow conditions. Therefore, further research and innovation are needed to address these issues in existing technologies. Summary of the Invention

[0005] Purpose of the invention: In view of the above-mentioned problems of the prior art, this application provides a method for simulating and predicting the concentration of heavy metals in water.

[0006] Technical solution: According to one aspect of this application, a method for simulating and predicting heavy metal concentrations in water includes:

[0007] Integrate raw monitoring data, hydrodynamic data, total suspended solids, remote sensing image data, and human activity data to generate a multi-source cleaning sequence data package;

[0008] Multi-source cleaning sequence data packets are used to perform cross-modal adsorption capacity estimation, infer particle chemical composition and adsorption parameters, and generate capacity feature packets.

[0009] Based on the multi-source cleaned sequence data packets, the travel time is calculated, and the upstream signals in the capacity feature packets and multi-source cleaned sequence data packets are causally aligned to generate travel time aligned feature packets.

[0010] By combining travel time alignment feature packets, capacity feature packets, and metal fingerprint parameters, physical constraints are applied to the morphological allocation constraint header to decouple and generate fractional concentration and total concentration prediction result packets.

[0011] Beneficial effects: This invention deeply couples the hydrodynamic physical mechanism with the particle adsorption chemical mechanism, improving the prediction accuracy and physical consistency of the model under unsteady-state conditions. The related technical effects will be described in detail below with reference to specific embodiments. Attached Figure Description

[0012] Figure 1 A flowchart illustrating a method for simulating and predicting heavy metal concentrations in water provided in this application embodiment.

[0013] Figure 2 This is a flowchart illustrating the generation of a capacity feature packet, as provided in an embodiment of this application.

[0014] Figure 3 This is a flowchart illustrating how constrained regression is applied to image texture morphology feature vectors to infer particle composition vectors representing chemical composition, as provided in an embodiment of this application.

[0015] Figure 4 This is a flowchart illustrating how to parse image remote sensing data from a multi-source cleaning sequence data packet and obtain image texture morphology feature vectors, as provided in an embodiment of this application.

[0016] Figure 5 This is a flowchart illustrating the decoupling and generation of fractional concentration and total concentration prediction result packages provided in the embodiments of this application. Detailed Implementation

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

[0018] It should be noted that the terms "first," "second," etc., in the specification and accompanying drawings of this invention are used to distinguish similar objects and are not necessarily used to describe a specific order or sequence. It should be understood that such data can be interchanged where appropriate so that embodiments of the invention described herein can be implemented in orders other than those illustrated or described herein. Furthermore, the terms "including" and "having," and any variations thereof, are intended to cover non-exclusive inclusion; for example, a process, method, system, product, or apparatus that includes a series of steps or units is not necessarily limited to those steps or units explicitly listed, but may include other steps or units not explicitly listed or inherent to such processes, methods, products, or apparatus.

[0019] To address the aforementioned issues, the applicant conducted in-depth searches and analyses, and discovered:

[0020] On the one hand, existing models typically simplify the adsorption-desorption process of heavy metals or use static, offline-calibrated partition coefficients. It is impossible to dynamically infer the upper bound of the total adsorption capacity of water bodies based on the real-time chemical composition (e.g., organic matter, clay ratio) and physical form (which can be characterized by image texture) of suspended matter and total suspended solids (TSS) in the water.

[0021] On the other hand, when dealing with the spatiotemporal dependence of upstream and downstream sections, existing methods mostly employ fixed time delays or pure statistical correlations. When river flow velocity changes significantly due to gate regulation or rainfall, fixed time delay windows cannot accurately capture the time-varying travel time of pollutants as they propagate with the water flow, leading to a causal misalignment between input features and the prediction target, affecting the model's response accuracy to sudden pollution events. The lack of chemical mechanisms and the misalignment between these processes and hydrodynamic physical processes make the model unreliable in predicting the distribution ratios of key forms such as dissolved, colloidal, and particulate states, especially under extreme scenarios of high concentrations or high TSS, where the prediction results are prone to violating the physical upper bound of adsorption saturation.

[0022] To solve these problems, combined with Figures 1 to 5 The present invention will be specifically described through the following embodiments.

[0023] Example 1: This invention relates to the field of water environment monitoring and data analysis, and in particular to a method for simulating and predicting the concentration of heavy metals in water.

[0024] In the context of this invention, to maintain consistency in description, technical terms and data items in some embodiments are defined as follows:

[0025] in, # _hat represents the symbol. Я represents the partial derivative.

[0026] C _d_m_fallback C _c_m_fallback、 C_p_m_fallback These represent the standby / alternative concentrations for the three states (dissolved, colloidal, and particulate). C t obs This indicates the actual observed concentration.

[0027] q _max_m q represents the maximum adsorption capacity per unit mass of particles. _max_clay q _max_org To represent the maximum adsorption capacity of a single adsorbent component (such as pure clay or pure organic matter) for the target substance (i.e., the maximum amount that a unit mass of the component can adsorb when it exists alone), w _clay w _org This represents the corresponding ratio.

[0028] The value of b is the affinity coefficient.

[0029] Multi-source data can comprehensively capture the complex factors influencing the migration and transformation of heavy metals; it may include:

[0030] Raw monitoring data, such as time series data of historical heavy metal concentrations (e.g., hourly or daily concentrations of lead, chromium, arsenic, etc.).

[0031] Hydrodynamic data, such as variables reflecting the physical movement of water bodies, such as flow velocity and water level at monitoring sections.

[0032] Total suspended solids (TSS) is the concentration of total suspended solids in water bodies and one of the main carriers for the adsorption of heavy metals.

[0033] Image remote sensing data, such as water images acquired through satellite remote sensing or on-site camera equipment. This image data can be used to analyze the spatial distribution of suspended matter, the image texture of particulate matter, etc.

[0034] Human activity data, such as variables that have a significant impact on hydrological conditions, such as records of the opening and closing status of upstream reservoirs or sluice gates, and records of sewage discharge from coastal industries.

[0035] Furthermore, traditional water heavy metal monitoring mainly relies on on-site sampling and laboratory chemical analysis. This method has inherent drawbacks, such as high cost, long cycle time, and difficulty in conducting high-frequency, large-scale monitoring. At the same time, the monitoring data is discrete in time and space, failing to capture the dynamic changes in pollution processes. Although predictions based on traditional mechanistic models (such as water migration and diffusion models) exist, these models typically require accurate fluid dynamic parameters and complex boundary conditions, making modeling difficult and resulting in slow response to sudden pollution events.

[0036] The method proposed in this invention can run on various computing devices. For example, it can be deployed in one or more server, workstation, personal computer, embedded system, or cloud platform environments. Exemplarily, the implementation environment may include a data acquisition module, a data storage module, and a prediction module. The data acquisition module is responsible for acquiring data from various sensors, remote sensing satellites, or third-party data platform interfaces. The data storage module (such as a database) is used to store massive amounts of historical data and real-time incoming data. The prediction module loads and executes the method of this invention to simulate and predict future heavy metal concentrations.

[0037] Optionally, an online dynamic update mechanism is provided. Following conventional methods in the art, the latest multi-source data is continuously acquired over time for preprocessing and feature extraction. Model parameters are updated based on the new data, for example, by retraining a single-source model weekly or monthly using new data, or updating the weight calculations for Bayesian fusion. During the online update process, both the old and new models are retained, and A / B testing is used to compare the performance of the two models in real-world applications (e.g., comparing the deviation rate between predicted and measured values), ensuring a smooth transition and optimizing prediction accuracy. It should be understood that parameters in the aquatic environment and external influences are dynamically changing; the online update mechanism enables the model to adapt to environmental changes and new data characteristics, maintaining long-term prediction accuracy and stability.

[0038] Example 2: A method for simulating and predicting the concentration of heavy metals in water is provided. In this example, the method can be executed by the prediction module described above.

[0039] Step S2.1: Integrate raw monitoring data, hydrodynamic data, total suspended solids, remote sensing image data, and human activity data to generate a multi-source cleaning sequence data package. This provides a foundation for subsequent analysis.

[0040] Specifically, data from multiple sources, different modalities, and different sampling frequencies are collected and preprocessed. Preprocessing may include time alignment (e.g., standardizing to hourly or daily granularity), quality control (e.g., removing obvious outliers), and imputation of missing data. This results in a unified, spatiotemporally aligned training sample stream, i.e., a multi-source cleaned sequence data package.

[0041] Step S2.2: Perform cross-modal adsorption capacity estimation using multi-source cleaning sequence data packets, infer particle chemical composition and adsorption parameters, and generate capacity feature packets.

[0042] In this embodiment, image remote sensing data (visual features) from the data package is used to estimate the physicochemical properties of particulate matter, such as the particle type composition (e.g., the proportion of clay, organic matter, and silt), through a cross-modal inference model. Based on this, the adsorption isotherm parameters of this particle composition under a specific metal are further retrieved or corrected from a knowledge base. Then, combined with the total suspended solids (TSS) concentration, the theoretical maximum adsorption capacity of the water body for that metal (i.e., the upper bound of the total adsorption capacity Q) is calculated. _max_m The inferred chemical parameters are organized into a capacity feature package, providing chemical mechanism constraints for subsequent steps.

[0043] Step S2.3: Calculate the travel time based on the multi-source cleaning sequence data packets, and perform causal alignment on the capacity feature packet and the upstream signals in the multi-source cleaning sequence data packets to generate a travel time aligned feature packet. This step introduces physical constraints from hydrodynamics.

[0044] Specifically, using hydrodynamic data (flow velocity, water level) and human activity data (gate status) from the multi-source cleaning sequence data packet, the travel time (T) from the upstream section (e.g., section i) to the downstream section (e.g., section j) at the current time (t) is calculated. _ij，t Furthermore, it is time-varying.

[0045] Next, based on this travel time, causal alignment is performed on the upstream signals. For example, the downstream concentration at time t should be influenced by sources such as upstream pollution discharge or upstream adsorption capacity state occurring within the range of time t. _ij，t This step involves constructing a physically reachable time window and weighting and aggregating upstream signals (including capacity feature packets, sewage discharge records, etc.) to prevent the model from incorrectly using future information or temporally misaligned correlations. The aggregated result is the travel time aligned feature packet.

[0046] Step S2.4: Combining the travel time alignment feature package, capacity feature package, and metal fingerprint parameters, apply physical constraints to the morphology allocation constraint header to decouple and generate the partial state concentration and total concentration prediction result package. This step is the final prediction and decoupling step. The metal fingerprint parameters refer to parameters characterizing the specific chemical behavior of each metal (such as valence state, complexing ability, partition coefficient, etc.). The morphology allocation constraint header is a special prediction model structure that does not directly predict the total concentration but explicitly models the three main states of heavy metals: dissolved state, colloidal state, and particulate state. When performing morphology allocation, this header applies physical constraints from the capacity feature package; for example, the upper limit of particulate state concentration must never exceed the upper limit of the total adsorption capacity of the water body. Simultaneously, it uses the alignment features from the travel time alignment feature package as input. Based on this, the header outputs the partial state concentrations of the three states, which are then synthesized into the total concentration, forming the partial state concentration and total concentration prediction result package.

[0047] It is understood that, in this embodiment, the capacity feature package generated by the cross-modal adsorption capacity estimation (specifically, the upper bound of the adsorption capacity) will serve as the upper limit of the physical constraint on the particulate concentration in the morphology allocation constraint head, thereby achieving dual constraint collaborative prediction of physical mechanism (travel time alignment) and chemical mechanism (capacity constraint).

[0048] Example 3 describes an optional implementation of data integration and uncertainty labeling. As an optional implementation, raw monitoring data, hydrodynamic data, total suspended solids, remote sensing image data, and human activity data are integrated to generate a multi-source cleaned sequence data package, which also includes spatiotemporal alignment and cleansing of multi-source heterogeneous data.

[0049] Accordingly, a unified time axis and a baseline time granularity (e.g., 1 hour) are established, and data with different sampling frequencies (e.g., daily average hydrodynamic data, weekly image data) are resampled or windowed. Simultaneously, a unique spatial index ID is created for each monitoring section.

[0050] Furthermore, the alignment process can also handle data conflicts and delays. For example, when there are multiple sources recording duplicates at the same cross-section and at the same time (such as simultaneous sensor direct readings and manual recordings), a primary / secondary source arbitration rule can be adopted (for example, automatic readings from the monitoring instrument take precedence over manual recordings, and manual recordings take precedence over interpolated estimates). In addition, the inherent acquisition delays and clock drift of different sensors can be corrected to ensure phase consistency of cross-source events (such as rainfall and flow velocity changes).

[0051] Based on this, this step performs missing measurement imputation and uncertainty labeling on the cleaned data, which may specifically include:

[0052] Step S3.1: For missing values ​​in the original monitoring data, use a combination of local K-nearest neighbors and Bayesian autoregression to generate imputed mean values.

[0053] In this embodiment, the local K-nearest neighbor algorithm is used to capture the spatial correlation between cross sections, that is, to estimate the missing values ​​using contemporaneous data from K geographically close cross sections. A Bayesian autoregressive model is used to capture the time-series dependence of the cross section itself. Using both methods in combination can generate relatively robust imputed means within a strict time roll-forward window (even if it does not use information from future times).

[0054] Step S3.2: Estimate the interpolation variance of the interpolation mean using multiple interpolation. Wherein, the interpolation variance (σ) imp 2 The interpolation mean is a quantitative indicator of the uncertainty or reliability of the interpolation. Multiple interpolation is an optional method to achieve this estimation.

[0055] Based on this, for the same missing data segment, the imputation process is independently repeated M times (e.g., M=5 or 10 times). Each imputation process will produce slightly different imputed values ​​due to the randomness of the model (such as Bayesian posterior sampling). Then, by calculating the variance of the M imputation results, the true imputation variance for the missing value is estimated. If the variance is large, it indicates that the imputation reliability at that point is very low.

[0056] Step S3.3: Incorporate the interpolation variance into the multi-source cleaning sequence data packet as the basis for calculating the data quality weight.

[0057] Based on this, the interpolated mean, binary labels indicating whether a value is interpolated, and the estimated interpolation variance, along with all other cleaned, aligned, and standardized data (such as hydrodynamic data and image features), are merged and packaged to form the final multi-source cleaned sequence data package. The interpolation variance field in this data package will be used in subsequent alignment steps to calculate data quality weights, automatically reducing the influence of low-confidence (high interpolation variance) data when the model aggregates upstream signals.

[0058] Example 4 describes the specific implementation process of cross-modal adsorption capacity estimation (VCAP), used to establish a bridge between the visual characteristics (images) of water bodies and their chemisorption capacity. Exemplarily, this example can be implemented as follows:

[0059] Step S4.1: Parse the image remote sensing data in the multi-source cleaning sequence data package to obtain the image texture morphology feature vector. Accordingly, read the image remote sensing data (e.g., particulate matter images acquired by on-site camera equipment or remote sensing satellite imagery) from the multi-source cleaning sequence data package.

[0060] Optionally, denoising, flat-field correction, and scale normalization are performed on the remote sensing image data. In some embodiments, preprocessing is performed to accurately extract effective features from the original image; specifically, this includes: applying mean-mode filtering or bilateral filtering to denoise the image and eliminate random noise; performing flat-field correction or histogram equalization to eliminate the effects of uneven illumination; and identifying and segmenting granular regions through threshold adaptation and morphological opening and closing operations.

[0061] On the preprocessed image, multi-scale texture and morphological features are extracted to form an image texture and morphological feature vector. Features may optionally include: gray-level co-occurrence matrix statistics (e.g., contrast, homogeneity, energy, entropy, used to describe the coarseness and regularity of the texture), local binary mode histogram (e.g., rotation-invariant mode, used to describe local texture), morphological grain size spectrum or equivalent diameter distribution (as a proxy index for particle size distribution), and fractal dimension (used to describe the complexity of particle edges).

[0062] Step S4.2: Apply constrained regression to the image texture morphology feature vector to infer the particle composition vector representing the chemical composition.

[0063] In this embodiment, the particle composition vector is used to characterize the relative proportions of different chemical components in the total suspended solids of water. For example, this vector can be represented as {w _clay w _org w _silt w _sand}, where w _clay For clay ratio, w _org The proportion of organic matter, w _silt The ratio of silt to sand, w _sand The proportion of sand particles is [not specified]. This type of component has a much higher adsorption capacity for heavy metals (especially clay and organic matter) than other components.

[0064] Constrained regression is used to establish a mapping from image texture morphology feature vectors (visual domain) to particle composition vectors (chemical domain). Its constrained nature is reflected in the following steps: applying a regressor to the image texture morphology feature vectors yields the original composition vectors; then, a simplex projection is performed on these (original composition vectors) to ensure that the particle composition vectors satisfy physical constraints, specifically, the constraints of non-negativity and summation to one (e.g., w). _clay +w _org +w _silt +w _sand =1, and all terms are greater than or equal to 0).

[0065] As an alternative implementation, the constrained regressor can be implemented as a cross-modal mapping network, such as a deep neural network. This network may include: one or more feature extraction layers (such as convolutional or fully connected layers) for receiving image texture morphology feature vectors; one or more cross-modal mapping layers (such as fully connected hidden layers) for learning a nonlinear mapping from visual features to chemical composition; and constrained projection layers (e.g., using a Softmax activation function or an explicit simplex projection operator) that satisfy the physical constraints of nonnegativity and summation to one for the output particle composition vector.

[0066] Step S4.3: Based on the particle composition vector and metal fingerprint parameters, retrieve and correct the adsorption isotherm parameters from the knowledge base.

[0067] Specifically, the baseline isotherm parameters are retrieved from the knowledge base based on the particle composition vector. Alternatively, based on the particle composition vector (e.g., {w... _clay w _orgThe baseline isotherm parameters are retrieved or obtained from a pre-established parameter knowledge base. This knowledge base stores the adsorption parameters of different single components (such as pure clay, pure organic matter) for different metals (such as lead, chromium) (e.g., the maximum adsorption capacity q of the Langmuir model). _max_base And affinity coefficient b _base , or the two feature parameters K of the Freundlich model _base and n _base By weighting the particle composition vector (e.g., q...) _max_m =w _clay ×q _max_clay +w _org ×q _max_org +...), to obtain the reference isotherm parameter set of the mixed particles.

[0068] Based on this, a competition correction coefficient based on metal fingerprint parameters is introduced to reduce the reference isotherm parameters (especially the affinity coefficient) to obtain adsorption isotherm parameters that take into account the multi-metal competition effect.

[0069] As an optional implementation method, the competing correction coefficient γ _m，m＇ (The competitive strength of metal m' against metal m) can be constructed in the following way:

[0070] ;

[0071] Among them, b _m b is the reference affinity coefficient for the target metal m (derived from reference isotherm parameters); _m＇ θ is the baseline affinity coefficient for the competing metal m'; θ and φ are dimensionless exponential parameters calibrated through a small number of competitive adsorption experiments; a _m＇ The effective activity of the competing metal m' (which can be calculated based on its concentration, pH, ionic strength, and other hydrochemical conditions); a _ref The reference activity is (e.g., 1 mol / L). This competition correction factor is then used to reduce the affinity coefficient, for example, the corrected effective affinity coefficient b. _eff_m It can be calculated as: b _eff_m =b _m / (1+Σ _m＇≠m γ _m，m＇ ×a _m＇ ); where Σ represents the summation over all coexisting competing metals m'.

[0072] Furthermore, to improve the accuracy of the parameters in the field, the corrected isotherm parameter set can be normalized for pH / ionic strength / organic matter conditions. That is, based on the water chemistry conditions (such as pH, temperature, and dissolved organic carbon) in the multi-source cleaning sequence data package, the parameters are corrected to the actual chemical conditions of the current water body to obtain the adsorption isotherm parameters under field conditions.

[0073] Step S4.4: Combine the total suspended matter and adsorption isotherm parameters in the multi-source cleaning sequence data package to calculate the adsorption capacity and generate a capacity feature package.

[0074] Accordingly, the theoretical upper bound Q for calculating the capacity of a unit volume of water for metal m is calculated. _max_m (Unit: mg / L). This upper limit is determined by the maximum adsorption capacity per unit mass of particles (i.e., q). _max_m Q is determined by multiplying the concentration of total suspended solids (TSS, g / L, from the isotherm parameter set) with the total suspended solids concentration (TSS). Correspondingly, Q... _max_m =q _max_m ×TSS.

[0075] At the same time, determine the capacity response function q _m (C _d_m This function is used to calculate the solubility concentration C. _d_m Mapped to the actual adsorption amount q per unit mass of particles _m If the Langmuir model is used, then the function is: q _m (C _d_m )=(q _max_m ×b _eff_m ×C _d_m ) / (1+b _eff_m ×C _d_m );

[0076] If the Freundlich model is used, then the function is: q _m (C _d_m )=K _m ×(C _d_m ) n_m Among them, b _eff_m K _m n _m All parameters are for on-site conditions.

[0077] In this embodiment, a branching strategy is adopted to obtain the saturation anchor signal during model training and avoid the logical contradiction of using the value to be predicted (particle concentration). When there is measured fractional phase data, if the historical data includes the measured particle phase concentration C... _p，m_observed Then the determined anchor point S is calculated. _m_anchor =min(C _p，m_observed / Q _max_m,1).

[0078] If no fractional measured data is available, the observed saturation is not calculated. Alternatively, the total prior C, constructed based on upstream information in subsequent steps, can be used. _proxy_m Generate soft anchor point range (e.g., S) _m_anchor Belongs to the interval [0, min(C)] _proxy_m / Q _max_m ,1)]), or skip this item.

[0079] Based on this, the calculated theoretical upper bound of capacity and capacity response function (and its specific parameter q) will be used to... _max_m b _eff_m K _m n _m ), and optional anchor point S _m_anchor The confidence level and its associated value are encapsulated together as a capacity feature packet. This feature packet will be aligned as an upstream signal and can also serve as a constraint for morphological assignment.

[0080] Example 5 describes an exemplary scheme for hydrodynamic-based travel time alignment that introduces physical constraints from hydrodynamics, such that downstream concentration prediction is based on its physically reachable upstream sources of influence. Specifically:

[0081] Step S5.1: Using the hydrodynamic data and human activity data in the multi-source cleaning sequence data packet, determine the accessibility of the cross-section and calculate the travel time. In other words, based on the hydrodynamic data and human activity data, determine whether the cross-section is accessible and further calculate the travel time.

[0082] In this embodiment, travel time refers to the time required for a pollutant (or its carrier) to be transported from upstream section i to downstream section j. It should be understood that travel time is not a fixed constant, but a time-varying travel time that is dynamically updated with flow velocity, water level, and gate status.

[0083] Optionally, the accessibility of the upstream section i to the downstream section j at time t can be determined based on the river network topology, real-time water level difference (from hydrodynamic data), and gate opening / closing status (from human activity data). For example, if the gate is closed or reverse pressure prevents the water flow from reaching the downstream section, it is determined to be inaccessible.

[0084] When reachability is determined, the convective travel time is calculated using the convection method. One possible calculation method is to use a piecewise constant velocity approximation of the travel time T. _ij，t ≈Σ _e∈path(i→j) (Δx _e / u _t (e));

[0085] Where path(i→j) is the set of all river segments e along the path from i to j; Δx _eIt is the length of river segment e; u _t (e) is the average flow velocity of the river segment at time t (from hydrodynamic data).

[0086] Furthermore, monitoring shows that the convective travel time exceeds a preset physical limit (e.g., T) due to extremely low flow velocities. _max This upper limit can be estimated based on the river length and historical p95 velocity (the case where this upper limit is used to solve numerical stability problems, i.e., u). _t (e) The problem that travel time tends to infinity in extremely low temperatures.

[0087] When the trigger condition (i.e., T) is met _ij，t >T _max The diffusion nucleus rollback mechanism is then enabled. At this point, convection time is no longer used; instead, the effective diffusion coefficient κ is used. _eff (This can be obtained by matching historical data or using an experience base) Estimate the diffusion backtracking travel time T _diff For example, T _diff ≈L _ij 2 / (2×κ _eff ), where L _ij Let be the total path length between i and j. Based on this, the diffusion backtracking travel time is used as the travel time. Furthermore, to suppress anomalous spikes in the travel time series, temporal smoothing, such as median smoothing or exponential smoothing, can be performed on the calculated travel time series.

[0088] Step S5.2: Construct a physically reachable time window based on travel time and generate a causal mask. The causal mask is a binary or weighted tensor that masks upstream information that is physically impossible to affect downstream predictions (i.e., time-displaced).

[0089] Optionally, based on T _ij，t To construct a physically reachable time window [tT] for the downstream section j at time t. _ij，t -δ, tT _ij，t +δ].

[0090] The window width δ is used to cover the uncertainties of the hydrodynamic model and the diffusion effect of the water body. Optionally, the window width can be determined using an adaptive rule, for example, δ=κ×IQR(T) _ij_history ), where κ is the adjustment coefficient (e.g., 1.5), IQR(T) _ij_history κ represents the interquartile range of the travel time from upstream section i to downstream section j within the historical window. When diffusion nucleus regression is triggered, κ can be appropriately amplified (e.g., κj). _diff >κ), covering greater diffusion uncertainty.

[0091] The causal mask M(i, j, t, t') is located in the time window [tT] at the upstream time t'. _ij，t-δ, tT _ij，t When within the window, its value is set to 1 (or a weight value based on the position within the window, such as Gaussian weights); while outside the window, its value is strictly set to 0.

[0092] Furthermore, if it is determined that the current period is unreachable (the travel time is invalid and there is no diffusion backoff), then an all-zero mask is generated for (i, j, t); if it is intermittently reachable (e.g., reachable in the near future but unreachable in the current period due to gate closure), the window that is most recently reachable can be retained as a weak reference, and the mask weight is reduced.

[0093] Step S5.3: Apply a causal mask to the capacity feature packet and upstream signals (e.g., upstream sewage records, upstream metal concentration, upstream TSS) originating from multi-source cleaning sequence data packets, and apply a weighted aggregation that includes freshness decay and data quality. In some embodiments, the causal mask, freshness decay weight, and data quality weight are combined to form a comprehensive weight; the comprehensive weight is then applied to aggregate the upstream signals. All relevant upstream signals within the physical time window (potentially from multiple upstream sections i and multiple time points t' within the window) are aggregated into a feature vector x for the downstream section j at time t. ~ _j (t). Weighted aggregation uses a composite weight w. This composite weight w can optionally be obtained by multiplying three components:

[0094] w(i,j,t,t')=M(i,j,t,t')×w _fresh (Δ)×w _quality (σ imp 2 );

[0095] Where M(i, j, t, t') is the causal mask (0 or 1). _fresh (Δ) represents the freshness decay weight; the further the signal deviates from the travel time in time, the lower its information value. Specifically, the signal lag Δ = |(tt') - T can be calculated based on the travel time. _ij，t |, Export the freshness decay weight w _fresh =exp(-λ×Δ). λ is the attenuation parameter. _quality (σ imp 2 ) represents the data quality weight. It is used to reduce the influence of imputed or outlier data. Specifically, it extracts the imputation variance (σ) generated by the missing data imputation step from the multi-source cleaned sequence data packets. imp 2 Based on this, data quality weights are derived, such as w. _quality =1 / (1+σ imp 2 ).

[0096] Optionally, an initial value λ is set for the freshness decay parameter._0 Approximately equal to 1 / median(T) _ij，t (i.e., measured by the reciprocal of the median travel time); next, through time-series roll-forward verification, among multiple candidate values ​​(e.g., {0.25λ})... _0 0.5λ _0 , λ _0 ,2λ _0 ,4λ _0 Choose the optimal value.

[0097] In another alternative implementation, a set of characteristic freshness parameters can be used. For example, a larger λ can be used for fast variables (such as rainfall, gate position changes) (to make them decay faster), and a smaller λ can be used for slow variables (such as upstream concentrations related to sediment release) (to make them decay more slowly).

[0098] Based on this, the aligned downstream feature x ~ _j (t) can be obtained using the following weighted average formula:

[0099] x ~ _j (t)=Σ _i，t＇ (w(i, j, t, t') × x _i (t�)) / Σ _i，t＇ w(i, j, t, t');

[0100] Where, x _i (t') is the original signal at time t' of the upstream section i (e.g., the upper limit of the adsorption capacity from the capacity feature pack, or the discharge record from the multi-source cleaning pack).

[0101] Step S5.4: Merge the aligned downstream feature stream generated by weighted aggregation with the travel time and causal mask to generate a travel time aligned feature packet for use in the morphology assignment constraint header. This step is an encapsulation step. The aligned downstream feature stream, travel time, and causal mask are merged to form a structured data packet, namely the travel time aligned feature packet. This feature packet will serve as the direct input to the morphology assignment constraint header.

[0102] According to one aspect of this application, a case study for calculating time-varying travel time is provided, as follows:

[0103] Assume the pollutants need to flow from upstream section i to downstream section j. This path consists of two consecutive river segments e. _1 and e _2 Composition. River section e _1 Length Δx _e_1 =18,000 meters (18 km); river section e _2 Length Δx _e_2=12,000 meters (12km); Total path L _ij =30,000 meters (30 km). The default physical upper limit T is set based on historical hydrological data. _max =172,800 seconds (48 hours). The effective diffusion coefficient κ used when triggering rollback. _eff =60m 2 / s.

[0104] At t=10:00 (normal flow rate), the water flow is stable, and the gate opens. According to hydrodynamic data, e _1 Section velocity u _t=10 (e _1 ) = 1.5 m / s; e _2 Section velocity u _t=10 (e _2 =1.2m / s.

[0105] Calculated based on convection. _1 For a period of time T _e_1 =Δx _e_1 / u _t=10 (e _1 ) = 18,000 / 1.5 = 12,000 seconds. e _2 For a period of time T _e_2 =Δx _e_2 / u _t=10 (e _2 =12,000 / 1.2 m / s = 10,000 seconds. Total travel time T_{ij, t=10} = T _e_1 +T _e_2 =12,000 + 10,000 = 22,000 seconds (approximately 6.1 hours). T _ij，t=10 (22,000s) <T _max (172,800s), not exceeding the upper limit. Convection travel time T is used. _ij，t=10 =22,000 seconds.

[0106] At t=14:00 (when the upstream sluice gate is opened to release water), the upstream sluice gate opening causes the flow velocity to increase. According to hydrodynamic data, e _1 Section flow velocity u_{t=14}(e _1 ) = 3.0 m / s; e _2 Section flow velocity u _t=14 (e _2 =2.5m / s.

[0107] e _1 For a period of time T _e_1 =18,000 / 3.0 = 6,000 seconds. e _2 For a period of time T _e_2=12,000 / 2.5=4,800 seconds.

[0108] Total travel time T _ij，t=14 =T _e_1 +T _e_2 =6,000 + 4,800 = 10,800 seconds (3.0 hours).

[0109] T _ij，t=14 (10,800s) <T _max (172,800s), not exceeding the upper limit.

[0110] Using convection travel time T _ij，t=14 =10,800 seconds. It can be seen that T _ij，t With flow velocity u _t The changing time-varying characteristics.

[0111] At t=18:00 (dry season and downstream gate closed), the water flow is extremely slow, almost stagnant. According to hydrodynamic data, e _1 Section flow velocity u _t=18 (e _1 ) = 0.1 m / s; e _2 Section flow velocity u _t=18 (e _2 =0.05m / s.

[0112] Calculated based on convection. _1 For a period of time T _e_1 =18,000 / 0.1=180,000 seconds. e _2 For a period of time T _e_2 =12,000 / 0.05=240,000 seconds.

[0113] Total travel time T _ij，t=18 =T _e_1 +T _e_2 =180,000 + 240,000 = 420,000 seconds.

[0114] T _ij，t=18 (420,000s)>T _max (172,800s).

[0115] The diffusion kernel rollback mechanism is triggered because the convection travel time exceeds the preset physical limit. Rollback calculation enabled: T _diff ≈L ij 2 / (2×κ _eff ) = (30,000) 2 / (2×60)≈900,000,000 / 120=7,500,000 seconds (approximately 86.8 days).

[0116] Using diffusion backtracking travel time T _ij，t=18 =7,500,000 seconds, in which the transmission during this period is marked as diffusion-dominated (e.g., using a wider time window or lower weights) in subsequent weighted aggregation.

[0117] According to one aspect of this application, the metal fingerprint parameters (hereinafter referred to as metal fingerprint parameters) refer to a set of parameters used to characterize the inherent or slowly changing physicochemical specific properties of a particular heavy metal (e.g., lead, chromium, arsenic, etc.) in an aquatic environment. These parameters serve as the basis for distinguishing different metal migration and transformation behaviors (e.g., adsorption, complexation, precipitation).

[0118] In an optional embodiment of the present invention, the metal fingerprint parameters include, but are not limited to, one or more of the following combinations:

[0119] Valence state distribution ratios, for example, for arsenic (As), the ratio of As(III) to As(V); for chromium (Cr), the ratio of Cr(III) to Cr(VI). Different valence states exhibit vastly different toxicities, solubilities, and adsorption properties.

[0120] The complexation stability constant family characterizes the stability constants of soluble or insoluble complexes formed by metal ions with common ligands in water (such as carbonate, chloride, sulfate, humic acid, or dissolved organic carbon).

[0121] The partition coefficient family of carriers characterizes the partition coefficient (Kd) or adsorption affinity order of metals on different adsorbent carriers (such as organic matter, clay minerals, iron and manganese oxides, etc.).

[0122] Precipitation / hydrolysis related parameters: For example, the solubility product (K) of the metal in which it forms hydroxide, carbonate, or sulfide precipitates. _sp The range of ) and the hydrolysis constant of metal ions.

[0123] Redox sensitivity is a factor that characterizes the influence of changes in the redox potential (Eh) of water on the distribution of metal valence states or solubility.

[0124] Colloidal affinity parameters characterize the tendency of a metal to bind to colloidal particles (usually those with a particle size of less than 1 micrometer).

[0125] The settling tendency parameter is an empirical coefficient characterizing the settling rate of particles (especially particulates) with attached metal.

[0126] In some embodiments, the acquisition path for metal fingerprint parameters can be diverse. Optionally, typical values ​​can be obtained from professional literature, chemical handbooks, or public databases as prior intervals; if field conditions permit, small-scale laboratory isothermal adsorption or competitive adsorption experiments can be conducted by collecting water samples and sediment samples to calibrate key parameters (such as the maximum adsorption capacity per unit mass of particles or the b-value of the local TSS); for parameters that cannot be directly measured, inversion fine-tuning or Bayesian updates can also be performed within the prior interval using historical data during model training.

[0127] Furthermore, since certain parameters (such as complexation constant and solubility product) may change slowly with environmental factors such as water temperature and pH, or differ in different river sections (such as industrial and agricultural areas), these parameters can be spatiotemporally interpolated to form a temporally sequenced and spatially allocated set of metal fingerprint parameters for real-time model use. This set of metal fingerprint parameters can be used to construct competing correction coefficients; it can also be used to establish the initial proportion of morphological allocation to determine the colloidal proportion coefficient.

[0128] Example 6 describes an optional implementation scheme for the morphology allocation constraint head (MFCN_head), namely, a morphology allocation constraint head with metal fingerprint coupling. This example integrates physical and chemical constraints to achieve robust decoupling of heavy metal morphology. Specifically, it includes:

[0129] Step S6.1: Combining the travel time alignment feature package, capacity feature package, and metal fingerprint parameters, apply physical constraints to the morphological allocation constraint header to decouple and generate the fractional concentration and total concentration prediction result packages.

[0130] Accordingly, the morphology assignment constraint header is a specialized model structure (e.g., a deep neural network module or coupled computational equations) whose inputs include a travel time alignment feature package (containing aligned upstream signals, TSS, water chemistry interpretations, etc.); a capacity feature package (containing the upper bound of adsorption capacity and the capacity response function); and a set of metal fingerprint parameters. This header is used to explicitly model and decouple dissolved state concentration (C0). _d_m ), colloidal concentration (C _c_m ) and particulate concentration (C _p_m ).

[0131] Step S6.2: Establish the initial morphological allocation ratio.

[0132] Before applying constraints, as an optional implementation, the header establishes an initial morphological allocation ratio π based on alignment features (such as TSS, pH, dissolved organic carbon, etc.) in the metal fingerprint parameter set and travel time alignment feature package. _p_m_init This initial ratio can be learned through a regression subnetwork and represents the particle distribution tendency determined solely by chemical affinity and aquatic conditions, without considering adsorption capacity saturation limitations.

[0133] Step S6.3: Constrain and reshape the initial shape allocation ratio.

[0134] In some embodiments, the adsorption capacity upper bound contained in the capacity feature package and the alignment signal in the travel time alignment feature package are used to constrain and shape the initial morphology allocation ratio to obtain a capacity-constrained morphology allocation ratio.

[0135] Optionally, upstream contributions and local feature quantities are aggregated from the travel time aligned feature package to construct a total prior C. _proxy_m Here, the total prior is the total concentration estimate obtained by regression of upstream alignment signals and local features (such as rainfall and sewage discharge), and is used only as an internal variable for constraint calculation, not as the final output.

[0136] Next, using the total prior C _proxy_m The upper bound of adsorption capacity Q in the capacity feature package _max_m Calculate the physical upper bound π of the particle state distribution ratio. _p_m_cap Used to avoid circular dependencies.

[0137] For example, this upper bound can be calculated as:

[0138] ;

[0139] Here, ε is a numerical stability constant used to prevent the denominator from reaching zero. The use of dissolved concentrations is avoided. It should be understood that cyclic dependence arises from the use of the dissolved concentrations to be determined during the solution process.

[0140] Furthermore, to reflect the inhibitory effect when the adsorption capacity is close to saturation, this embodiment introduces a saturation constraint.

[0141] Accordingly, when the capacity saturation S (S _m_anchor Or through total prior and q _m When the capacity saturation obtained by the function estimation approaches 1, the particle state allocation ratio is dynamically reduced by the saturation suppression function.

[0142] For example, the function can be represented as:

[0143] ;

[0144] Where, π _p_m_init The initial ratio is obtained from the aforementioned steps, and η (e.g., in the [0, 1] interval) and ρ (e.g., greater than or equal to 1) are learnable parameters used to control the suppression strength and nonlinearity.

[0145] At the same time, the ratio of dissolved and colloidal states is adjusted to compensate for the loss of 1, so that the sum of the ratios of the three states remains 1 after adjustment.

[0146] Based on this, the physical upper bound and the ratio after saturation suppression are combined (for example, taking the minimum of the two) to obtain the final capacity-constrained morphological allocation ratio.

[0147] Step S6.4: Decouple and calculate the fractional concentrations. In some embodiments, based on the capacity-constrained morphological distribution ratio, the capacity response function in the capacity feature packet, and the travel time alignment feature packet, the concentrations of dissolved, colloidal, and particulate states are decoupled and calculated, and synthesized into a fractional concentration and total concentration prediction result packet.

[0148] Furthermore, by aligning feature packets with travel time and regressing subheadings, an unconstrained estimate of the total concentration C is obtained. _m_tilde The C _m_tilde It is an internal intermediate quantity, representing the model's best fit to the total concentration at the current moment.

[0149] Next, based on the unconstrained estimation of total concentration, the capacity response function and upper bound of adsorption capacity in the capacity characteristic packet, and the colloidal state proportion coefficient determined by the metal fingerprint parameters, a mass conservation equation for the dissolved state concentration is established. This can be described as the following formula: C _m_tilde =C _d_m +C _c_m +C _p_m .

[0150] Among them, the colloidal concentration C _c_m Approximated as proportional to the concentration of dissolved substances: C _c_m =k _c_m ×C _d_m .

[0151] It should be noted that the colloidal proportion coefficient k _c_m It is not a fixed constant, but is optionally dynamically determined by the set of metal fingerprint parameters and the water chemistry interpretations (such as dissolved organic carbon and ionic strength) in the travel time aligned feature package.

[0152] Particulate concentration C _p_m Then, based on the capacity response function q _m (C _d_m ) and the upper limit of adsorption capacity (Q) _max_m )Decide:

[0153] C _p_m =min(q _m (C _d_m )×TSS, Q _max_m ).

[0154] Substituting the above relationship into the mass conservation equation, we obtain the result only regarding the unknown C. _d_m Nonlinear equations:

[0155] ;

[0156] It can be proven that, due to k _c_m ≥0, and q _m (C _d_m ) is about C _d_m A nonnegative monotonically increasing function (such as the Langmuir or Freundlich function) ensures that the global right-hand side of the mass conservation equation is monotonically increasing with respect to the dissolved concentration. This guarantees the uniqueness of the numerical root.

[0157] Based on this, the mass conservation equation is solved using numerical root-finding methods to obtain the solution for the dissolved concentration. In other words, the equation can be solved using numerical root-finding methods (e.g., the bisection method, Newton's method, or Broyden's method, which have a definite convergence guarantee) to obtain a unique solution for the dissolved concentration.

[0158] Furthermore, based on the solution of the dissolved state concentration, the concentrations of colloidal and particulate states are calculated, and a package of predicted partial and total concentrations is synthesized. That is, based on the obtained C... _d_m Calculate in reverse:

[0159] ;

[0160] ;

[0161] ;

[0162] The three-state concentrations and their sum (or the sum of the three) are packaged together into a package of predicted results for both state concentration and total concentration.

[0163] Step S6.5, Failure Rollback Mechanism. To ensure that the system can still provide a solution under extreme input conditions or when numerical calculations are unstable, this embodiment includes a rollback mechanism.

[0164] Accordingly, the convergence status of the numerical root-finding method is monitored. When convergence is determined to be non-convergent, a failure fallback mechanism is activated. The triggering conditions for non-convergence may optionally include: the number of iterations exceeding a preset upper limit N. _max (e.g. N) _max =100 iterations), or the relative error of the residual is greater than the convergence threshold ε (e.g., ε=1e-5).

[0165] When a rollback is triggered, a proportionally locked linear approximation method is used to decouple the fractional concentration and total concentration prediction result packages based on the unconstrained total concentration estimation and the capacity-constrained morphological allocation ratio.

[0166] Correspondingly, the morphological allocation ratio π is constrained by capacity. _p_m Unconstrained estimation of total concentration C _m_tilde The concentration is allocated according to the following linear approximation rule: C _p_m_fallback=π _p_m ×C _m_tilde C _d_m_fallback =(1-π _p_m ) / (1+k _c_m )×C _m_tilde ;

[0167] C _c_m_fallback =k _c_m ×C _d_m_fallback ;

[0168] Optionally, C can be further modified. _p_m_fallback Perform an upper bound check, i.e., C _p_m_final =min(C _p_m_fallback Q _max_m The truncated difference will be proportionally replenished to C. _d_m and C _c_m .

[0169] This rollback method keeps the output value stable, and its distribution ratio π _p_m It is still constrained by physical capacity, which ensures the rationality of the results.

[0170] Example 7 provides optional implementations for model training, calibration, and export, describing the training and deployment process of a prediction model. Specifically, it includes:

[0171] Step S7.1: Using the fractional concentration and total concentration prediction results package and the original monitoring data (training labels), construct a threshold consistency loss function that is sensitive to the regulatory threshold.

[0172] In this embodiment, model training is used not only to minimize the root mean square error (RMSE) but also to ensure the accuracy of predictions for high-concentration events (i.e., events approaching or exceeding regulatory thresholds). To this end, an optional threshold-consistent loss function is employed.

[0173] For example, this loss function can be composed of a quantile loss and an extremum correction term. The quantile loss allows the model to directly learn to predict a high quantile (e.g., the 0.90 or 0.95 quantile), penalizing upper tail errors much more than lower tail errors. The extremum correction term is used to better model fat-tailed distributions; for example, it can be achieved by fitting the upper tail of the predicted residuals to a generalized Pareto distribution (GPD) and incorporating its negative log-likelihood as part of the loss function, allowing the model to learn the statistical properties of extreme events.

[0174] Step S7.2: Apply the threshold consistency loss function to perform temporal roll-forward training and probability calibration.

[0175] Furthermore, a loss function is applied to optimize the model (including all trainable parameters such as the VCAP network, the initial scaling network, and the unconstrained estimation network).

[0176] Optionally, a time-order roll-forward training method can be used. Specifically, the chronological order of the data is strictly adhered to to avoid future information leakage. For example, the dataset is divided into K folds (Fold1, Fold2, ..., Fold...) according to time. K The model was trained on Fold1 and validated on Fold2; on Fold... 1+2 Train on Fold3, validate on Fold3, and so on. This method simulates the online learning and prediction performance of the model in a real-world deployment environment.

[0177] Furthermore, after the model parameters have been trained, probability calibration is performed on the model's probability output. This is because the confidence or probability output of many machine learning models (especially neural networks) is often uncalibrated (e.g., the model predicts a 70% probability of overshooting), meaning that the predicted probability does not match the actual frequency of the event.

[0178] Based on this, on an independent validation set, techniques such as temperature scaling (adjusting the sharpness of the Softmax output through a scalar T) or ordinal-preserving regression (a non-parametric method) can be used to post-process the model's original probability output so that the calibrated probability (e.g., the overthreshold probability) can reflect the actual frequency of the event.

[0179] Step S7.3: Export the calibrated model and inference pipeline. Solidify and package all necessary components for easy deployment to the production environment.

[0180] The inference pipeline is an executable, end-to-end processing flow. Optionally, it encapsulates data preprocessing and interpolation logic; the VCAP model and its learned parameters; travel time calculation and causal alignment logic; the MFCN_head (including a numerical solver or backoff mechanism) and its learned parameters; and probability calibration parameters (e.g., the T-value for temperature scaling). The derived calibrated model and inference pipeline constitute the final delivery of this invention, which can be used in online prediction systems to receive real-time multi-source data input and output calibrated fractional concentrations, total concentrations, and overthreshold probabilities.

[0181] Example 8 provides an exemplary scheme based on a multi-source data fusion strategy, which integrates multi-source data using different machine learning fusion strategies. This example can be based on data preparation and feature engineering, and its prediction model can adopt a fusion method.

[0182] In this embodiment, one possible implementation is a feature-level fusion method. Specifically, data from different data sources (such as environmental indicators, social indicators, water quantity indicators, water quality indicators, and image features) are preprocessed, feature extracted, and standardized, and then concatenated along the feature dimension to form a unified, high-dimensional feature matrix. Next, this unified feature matrix is ​​used as input to train a single prediction model (e.g., random forest, gradient boosting tree, or multilayer perceptron MLP), which directly learns the mapping relationship from the fused features to the target heavy metal concentration end-to-end. This method has a simple structure and can directly utilize information from all data sources; however, when any data source has large-scale missing data, additional imputation processing is required, and the data standardization requirements are high.

[0183] As an alternative implementation, decision-level fusion (Bayesian fusion) can be used. This method is suitable for handling situations where there are differences in sampling frequency, granularity, and partial data loss among the data sources.

[0184] Accordingly, for each independent data source (e.g., one set of hydrodynamic and water quality data, and another set of image and remote sensing data), an independent single-source prediction model is constructed (e.g., using Partial Least Squares (PLS), LASSO, or Random Forest, etc.). For each single-source model, its mean squared error (MSE) is calculated based on its historical prediction performance within the most recent time window (e.g., the past i time points). The MSE is considered as a measure of the model's prediction uncertainty, and based on this, the contribution weights of each single-source model are dynamically allocated using a Bayesian weighting formula.

[0185] For example, the weight w of the j-th model at time i _i (j) can be set to its recent mean square error MSE _i|i-1 The reciprocal of (j) is normalized to w. _i (j)=(1 / MSE _i|i-1 (j)) / Σ _k (1 / MSE _i|i-1 (k)); where Σ _k This represents the summation over all k single-source models. The weighted prediction results of each single-source model are used as the final fused prediction output.

[0186] This decision-level fusion method can dynamically adjust weights based on the real-time performance of different data source models. It is highly robust. Even if some data sources (i.e., a single-source model) temporarily fail or the data quality deteriorates (manifested as an increase in MSE), their weights will be automatically reduced, ensuring the stability and accuracy of the overall prediction.

[0187] Example 9 provides an exemplary scheme based on differentiable physics propagation and sparse source inversion. This method can not only make positive predictions, but also has the ability to perform source inversion for sudden pollution events.

[0188] In this embodiment, the model can be constructed as a differentiable physical propagation layer (PDE-Layer). Specifically, this layer is a differentiable computational operator that numerically discretizes the one-dimensional or two-dimensional convection-diffusion equations describing pollutant migration (e.g., using the finite difference method or the finite element method).

[0189] For example, the operator can be represented in the following matrix time-step iteration form:

[0190] C _t+1 =A(u, κ)×C _t +B×S _t ;

[0191] Among them, C _t It is a vector representing the heavy metal concentration at all river grid points at time t; C _t+1 S is the concentration vector at time t+1. A(u, κ) is the state transition matrix, determined by the velocity field u and diffusion coefficient κ at the river grid points (u and κ can be obtained from hydrodynamic calculations or data packages); the A matrix encapsulates the physical processes of convection and diffusion. _t S is a vector representing the pollution source release intensity at each grid point at time t. B is the source injection matrix, used to represent the source term S. _t Assign them to the corresponding river grid points.

[0192] Since all matrix operations (multiplication and addition) in this equation are differentiable, this PDE layer can be embedded into a neural network to achieve end-to-end training. During forward prediction, given S... _t (For example, known sewage discharge records) can be forward-propagated hourly to predict future C. _t .

[0193] Furthermore, this embodiment can utilize the characteristics of the differentiable PDE layer to achieve sparse multi-hypothesis source inversion. This is used to infer the unknown pollution source S after observing concentration anomalies at downstream monitoring points. _t (Its location, start and end times, and intensity). It can be constructed as an optimization problem, with its objective function optionally including:

[0194] min _S≥0 Σ _t |C # _t (S)-C t obs ||2 2 +λ _1 ×∣∣S∣∣ _1+λ _2 ×Σ _t TV(S _t );

[0195] Here, S≥0 is a physical constraint, indicating that the intensity of the pollution source cannot be negative. Σ _t |C # _t (S)-C t obs ||2 2 This is the L2 loss term, representing the predicted concentration C obtained through forward propagation through a differentiable PDE layer. # _t (S) and actual observed concentration C t obs The Euclidean distance between them (i.e., the prediction error). λ _1 ×∣∣S∣∣ _1 It is an L1 regularization term (∣∣S∣∣) _1 For S _t The L1 norm is used to impose sparse priors. Assuming that sudden pollution sources typically occur only at a few locations and for a short period, L1 regularization tends to produce sparse solutions (i.e., S1). _t (Most elements in λ are zero). _2 ×Σ _t TV(S _t ) is the total variation (TV) regularization term, used to apply a time-smoothing prior. It is assumed that the release intensity of the pollution source does not change drastically between adjacent time points. λ _1 and λ _2 It is a hyperparameter used to balance the weights of various components.

[0196] Since the entire objective function is differentiable (or can be solved by proximal gradient descent), the optimization problem can be solved using the backpropagation algorithm, starting from the observed C. t obs Inversely deduce the most likely source of pollution, S _t .

[0197] Optionally, to improve the robustness of the inversion, multiple hypothesis tracking can be adopted, that is, the optimization problem is solved in parallel for multiple possible initial pollution locations or reporting times, and the Top-k source hypotheses are sorted and their confidence scores are output based on the final loss value and information criteria (such as AIC / BIC).

[0198] In another alternative implementation, the differentiable PDE layer can be combined with a machine learning model (such as a residual CNN), where the PDE layer is responsible for modeling interpretable convection-diffusion-dominant processes, while the ML layer is used to learn and compensate for non-ideal effects not modeled by the PDE layer (such as local backflow, temporary gate effects, etc.).

[0199] Example 10 provides another specific numerical calculation case, especially the morphological allocation constraint head MFCN_head.

[0200] Suppose that at a certain time t, at the downstream section j, for metal m (e.g., lead), the model has completed the preliminary steps and obtained the following input parameters (the values ​​are only examples and are not intended to limit the invention):

[0201] Unconstrained concentration estimation C _m_tilde =50.0µg / L. Upper limit of adsorption capacity Q _max_m =25.0µg / L. Total suspended solids concentration (TSS) = 100.0mg / L. Colloidal proportioning factor k _c_m =0.1 (dimensionless). Capacity response function (Langmuir model): q _m (C _d )=(q _max_m ×b _eff_m ×C _d ) / (1+b _eff_m ×C _d Assume the parameter is q. _max_m =0.3mg / g, b _eff_m =0.2L / µg.

[0202] Furthermore, the particulate concentration C was calculated using TSS. _p_m (C _d )=min(q _m (C _d )×TSS, Q _max_m )=min(((0.3×0.2×C _d ) / (1+0.2×C _d ))×100.0,25.0)=min((6.0×C _d ) / (1+0.2×C _d ), 25.0); the mass conservation equation is f(C _d )=C _d +C _c +C _p -C _m_tilde

[0203] =C _d +(0.1×C _d )+min((6.0×C _d ) / (1+0.2×C _d ), 25.0)-50.0

[0204] =1.1×C _d +min((6.0×C _d ) / (1+0.2×C _d), 25.0)-50.0=0;

[0205] Furthermore, the equation f(C) is solved using numerical root-finding methods. _d Regarding C _d Monotonically increasing. A bisection method is used to solve this problem, with the search range set to [0, 0, C]. _upper ], where C _upper Option C is acceptable. _m_tilde That is, 50.0.

[0206] Try C _d =(0.0+50.0) / 2=25.0; Calculate f(25.0)=1.1×25.0+min((6.0×25.0) / (1+0.2×25.0),25.0)-50.0=27.5+min(150.0 / 6.0,25.0)-50.0=27.5+min(25.0,25.0)-50.0=27.5+25.0-50.0=2.5

[0207] Since f(25.0)>0, the upper bound C is set. _upper Adjusted to 25.0. The new bracket is [0.0, 25.0].

[0208] Try C _d =(0.0+25.0) / 2=12.5; Calculate f(12.5)=1.1×12.5+min((6.0×12.5) / (1+0.2×12.5), 25.0)-50.0=13.75+min(75.0 / 3.5, 25.0)-50.0=13.75+min(21.43, 25.0)-50.0=13.75+21.43-50.0

[0209] =-14.82; Since f(12.5)<0, the lower bound is adjusted to 12.5. The new brackets are [12.5, 25.0]. ... (After several subsequent iterations, the brackets continue to shrink)...

[0210] Assuming the iteration converges (e.g., the residual is less than ε = 1e-5), the solution concentration C is obtained. _d_m ≈22.1µg / L.

[0211] Further, the concentrations of the three states were calculated and synthesized. Dissolved C _d_m =22.1µg / L. Colloidal C _c_m =k _c_m ×C _d_m =0.1×22.1=2.21µg / L. Particulate C _p_m =C _m_tilde -C _d_m -C_c_m =50.0 - 22.1 - 2.21 = 25.69 µg / L. The calculated C at this point... _p_m =25.69µg / L, slightly higher than the upper capacity limit Q. _max_m =25.0µg / L. In the actual solver, due to the existence of the min() function, the solution will converge to C. _p_m Exactly equal to Q _max_m point.

[0212] Accordingly, the corrected solution (the actual convergence point of the solver) is calculated. C _p_m =25.0µg / L (reaching the upper capacity limit Q) _max_m );

[0213] At this time, C _d_m +C _c_m =C _m_tilde -C _p_m =50.0-25.0=25.0; C _d_m +0.1×C _d_m =25.0; 1.1×C _d_m =25.0; C _d_m =22.73µg / L. C _c_m =0.1 × 22.73 = 2.27 µg / L. Total concentration C _m =50.0µg / L.

[0214] Suppose the numerical root-finding method is determined to have not converged (e.g., the number of iterations exceeds N). _max =100 times).

[0215] Enable failure rollback mechanism. The system uses a proportionally locked linear approximation method. It is assumed (constrained shaping) that the capacity-constrained shape allocation ratio π has been calculated. _p_m =0.45 (45%). (This ratio already takes Q into account.) _max_m and C _proxy_m (constraints).

[0216] Using linear approximate decoupling, the particulate C _p_m =π _p_m ×C _m_tilde =0.45 × 50.0 = 22.5 µg / L. Dissolved C _d_m =(1-π _p_m ) / (1+k _c_m )×C _m_tilde =(1-0.45) / (1+0.1)×50.0=0.55 / 1.1×50.0=0.5×50.0=25.0µg / L. Colloidal C _c_m =k _c_m ×C _d_m=0.1×25.0=2.5µg / L. Optionally, the calculated C... _p_m =22.5µg / L, not exceeding Q _max_m =25.0µg / L, the regression result is valid. Therefore, the total concentration C _m =50.0µg / L.

[0217] Through the above calculation process, the model in this embodiment strictly adheres to the physical upper bound during normal convergence and reverts to the linear approximation when the numerical solution fails, thus ensuring the robustness and physical rationality of the method.

[0218] According to one aspect of this application, the following method can also be used for inversion of multiple hypothetical sources:

[0219] Discrete potential pollution source grid points are predefined within the target river network. These grid points can be optionally located at all known industrial discharge outlets, key monitoring areas, tributary confluences, or river hydrodynamic singularities (such as bends or backwaters). Simultaneously, discrete potential reporting times are defined, for example, within 48 hours prior to the observation of an anomaly in concentration at downstream monitoring points, at 1-hour or 3-hour intervals. By combining these spatial grid points and temporal starting points, an initial set of hypotheses is generated, where each hypothesis represents a single-source leakage event occurring at location x and time t.

[0220] For each hypothesis in the initial hypothesis set, it is used as the source term S of the differentiable PDE layer. _t Given the initial conditions (i.e., non-zero at (x, t) and zero elsewhere), perform a sparse source inversion optimization process to solve for the optimal source intensity curve S under this assumption. _t And its corresponding L2 loss value (i.e., the fitting residual between the predicted value and the actual observed value). To improve computational efficiency, an early pruning strategy can be adopted, that is, in the early stage of optimization iteration, the calculation is terminated prematurely if the L2 loss value is much larger than the current optimal solution (or the preset baseline).

[0221] After obtaining (e.g., Top-K) candidate solutions (i.e., S... _t After determining the curve and its L2 loss value, an information criterion (such as Bayesian Information Criterion (BIC) or Akaike Information Criterion (AIC)) is used for final ranking. This criterion is used to balance the L2 loss term (representing goodness of fit) with model complexity (which can be determined by λ). _1 ×∣∣S∣∣ _1 +λ _2 ×ΣTV(S _t The value of the regularization term or S _t The number of non-zero spatiotemporal points is used to characterize the solution. The solution with a lower L2 loss value (good fit) and lower model complexity (sparser and smoother source) will be assigned the highest confidence and used as the most likely pollution source inversion result.

[0222] According to another aspect of this application, the parameters of the metallic fingerprint can also be (i.e., inverse fine-tuning) including:

[0223] Accordingly, the target parameter P (e.g., the stability constant K of a complexation reaction or the solubility product K of a precipitate) is obtained from professional literature, chemical handbooks, or databases. _sp The typical value of ) is used to set a reasonable prior interval in a physicochemical sense [P]. _min P _max ].

[0224] Alternatively, the parameter P can be treated as a trainable weight of the model (e.g., a morphological assignment constraint head). To ensure that P remains within the prior interval during training, constraint transformations can be employed. For example, an unconstrained trainable variable P can be defined within the model. _trainable It is mapped to the target interval using a nonlinear function (such as a scaled sigmoid function or a Tanh function), i.e.: P = P _min +(P _max -P _min )×Σ(P _trainable ), where Σ(×) is the Sigmoid activation function with a range of [0, 1].

[0225] During model training, this internal variable P _trainable Along with the other weights (W, b) of the neural network, these are considered trainable parameters. When computing the loss function (e.g., a threshold-consistent loss function), the gradient ЯLoss / ЯP of the total loss Loss with respect to P is backpropagated to P via the chain rule. _trainable Result in ЯLoss / ЯP _trainable The optimizer (e.g., Adam or SGD) utilizes this gradient to calculate Loss / P. _trainable For P _trainable Update accordingly. Based on this, under the constraints of historical observation data, parameter P is adjusted in the direction that minimizes the total loss of the model within its preset physical prior interval [P]. _min P _max The parameters are automatically inverted and fine-tuned to obtain parameter values ​​that conform to the characteristics of the local water environment.

[0226] According to another aspect of this application, the total prior C _proxy_m Alternatively, it can be constructed using the following method. Accordingly, the aggregate prior is a predicted value generated by a dedicated aggregate prior regression sub-model. This regression sub-model can optionally be a machine learning model, such as a gradient boosting machine (GBDT), random forest (RF), or a small, lightweight multilayer perceptron (MLP) network. The input feature vector X of this regression sub-model... _proxyIt is a feature set specifically constructed to avoid circular dependencies and does not contain the concentration to be solved at the current moment.

[0227] Optionally, the feature set includes: all upstream water quality and hydrodynamic variables in the aligned downstream feature flow (e.g., upstream metal concentration, upstream TSS concentration, upstream flow velocity, etc. after travel time alignment); local environmental characteristics of the current downstream section j at time t from multi-source cleaning sequence data packets (e.g., local TSS concentration, local flow velocity, local water level, local rainfall, etc.); and corresponding human activity data (e.g., local sewage discharge records, upstream gate alignment signals, etc.).

[0228] Based on this, the training target of this regression sub-model is the true total concentration observation value C in the original monitoring data. _observed The model can be jointly trained with the main parameters of the morphological assignment constraint head (MFCN_head) during the time-series roll-forward training process, or it can be trained in two stages: first, the regression sub-model is trained independently, then frozen and used as a feature extractor to train the MFCN_head.

[0229] Furthermore, the output of this regression sub-model ;

[0230] It is used as input to calculate the physical upper bound of the particle state distribution ratio. It provides a reasonable estimate of the total concentration that is independent of the dissolved state concentration to be solved in the current step, thus resolving the cyclic dependency problem.

[0231] This invention addresses the lack of real-time chemical mechanism constraints and the inability to dynamically infer the upper bound of total adsorption capacity from images by providing a method for cross-modal adsorption capacity estimation (VCAP). This method analyzes image texture morphology feature vectors and uses constrained regression (e.g., cross-modal mapping networks) to infer particle composition vectors characterizing chemical composition. Based on this, and combined with competitive correction, a dynamic upper bound of the adsorption capacity is calculated. This overcomes the limitation of traditional methods that rely on static allocation coefficients.

[0232] It should be understood that this scheme provides a hydrodynamic alignment strategy. This strategy calculates the time-varying travel time based on real-time flow velocity and gate status, and introduces a diffusion kernel backoff mechanism to address unsteady flow. A physically reachable time window and a causal mask are constructed to maintain physical consistency between upstream and downstream signals, thus resolving the causal misalignment problem.

[0233] Building upon this, and addressing the consequences of unreliable predicted morphologies and the tendency to violate the physical upper bound of adsorption saturation, this scheme employs a morphology allocation constraint head (MFCN_head). This head, when decoupling the calculation of the three-state concentrations, uses the upper bound of adsorption capacity calculated in the aforementioned steps as a physical constraint, introducing a saturation suppression function to ensure that the prediction results adhere to the physical saturation upper bound, thereby improving the reliability of predictions under extreme scenarios.

[0234] In some embodiments, cross-modal adsorption capacity estimation also includes:

[0235] Optionally, images and remote sensing data—particulate matter images / textures and associated total suspended matter—are read from the multi-source cleaning sequence data package. The images are then subjected to denoising (median or bilateral filtering), flat-field correction, and scale normalization (unifying spatial resolution to a baseline pixel size). Particle region segmentation is performed (threshold adaptive + morphological opening / closing operations), and multi-scale texture and morphological features are extracted to form an image texture and morphological feature vector, which is used to regress the particle composition vector. The included features include: gray-level co-occurrence matrix statistics (contrast, homogeneity, energy, entropy, etc.), local binary mode histograms (rotation invariant / multi-radius), edge orientation histograms, particle size proxy (morphological particle size spectrum / equivalent diameter distribution), shape spectrum (roundness, aspect ratio, convexity), fractal dimension, etc.

[0236] Optionally, the image texture morphology feature vector and total suspended matter are read, and a constrained regressor trained based on experimental calibration (such as X-ray diffraction / burn-off) outputs a particle composition vector (w) that satisfies non-negativity and sums to 1. _clay w _org w _silt w _sand …, the sum of w is 1). To enhance robustness, features are grouped and redundancy is removed, and a simplex projection is applied to the regression output to make it physically feasible; at the same time, the composition uncertainty (such as the confidence score based on the regression residual) is estimated and passed as a weight to the subsequent parameter retrieval process.

[0237] Optionally, the particle composition vector and the metal fingerprint parameters (or simply metal fingerprint parameters) of each metal (valence state, complexing ability, partition coefficient, sedimentation tendency, etc.) are read. The isotherm parameters of each metal m on a given particle type are retrieved from the knowledge base or obtained through regression analysis. The effective parameter set is then obtained by weighting the parameters according to the composition. For the Langmuir model, q is retrieved. _max_m (Unit: mg / g) and b _m (Unit: L / mg); For the Freundlich model, K was retrieved. _m With n _m To address the competitive adsorption of multiple metals in coexistence, a competition correction coefficient γ is introduced. _m，m' (representing the competitive strength of metal m' against metal m, varying with ionic strength and pH), construct the effective affinity coefficient b. _eff_m =b _m / (1+Σ _m＇≠m γ _m，m＇ ×a _m＇ ); where a _m＇ This represents the activity of metal m' or the relative effective concentration derived from the metal fingerprint parameters (or simply metal fingerprint parameters) of each metal and the particle composition vector. The upper bound for surface site saturation is expressed as the total site density N. _s (Given from particle type library, unit mmol / g) Constraint q _max_m No more than N _s ×M _m (M) _m (where m is the molar mass of the metal). Output isotherm parameter set_per metal (q) _max_m b _eff_m or K _m n _m ) and competing correction parameter sets, used for capacity and upper bound calculations, and morphological allocation constraints.

[0238] Optionally, read the isotherm parameter set_per metal and total suspended matter, and calculate the theoretical upper bound (in mg / L) of the capacity of water per unit volume for metal m: Q _max_m =q _max_m ×TSS. Define the capacity response function q. _m (C _d_m ), the dissolved concentration C _d_m This is mapped to the amount of particles adsorbed per unit mass.

[0239] q _m (C _d_m )=(q _max_m ×b _eff_m ×C _d_m ) / (1+b _eff_m ×C _d_m (Solving the Langmuir model);

[0240] q _m (C _d_m )=K _m ×(C _d_m ) n_m (Solving the Freundlich model).

[0241] Upper limit of particulate phase concentration C _p_m_cap =min(q _m (C _d_m )×TSS, Q _max_mTherefore, given this, if the original monitoring data (metal concentration sequence) exists during the training period, simplify the allocation of prior estimates of dissolved state concentrations to calculate the temporary saturation S. _m_obs =C _p_m / Q _max_m (Used only as a monitoring signal, not as a hard constraint). The output includes Q. _max_m b _eff_m or K _m n _m Capacity response function and optional S _m_obs The capacity feature packet can be used as an upstream signal for alignment, or as the core input for morphological assignment and upper bound constraints.

[0242] In some embodiments, hydrodynamic-based travel time alignment also includes:

[0243] Accordingly, the raw hydrodynamic data (flow velocity, water level, human activity data, and gate status) from the multi-source cleaning sequence data packet are read. A directed graph is constructed based on the cross-sectional latitude and longitude and the river segment connectivity. For each time t, the directionality and connectivity are determined based on the water level difference and the gate opening / closing status. When the gate is closed or reverse pressure causes inaccessibility, the corresponding directed edge is set to zero. The length Δx and the current effective flow velocity u of each directed edge are stored. _t The output is a directed adjacency and reachability matrix, which can be used to calculate travel time and to construct causal masks.

[0244] Accordingly, the directed adjacency and reachability matrix and the original hydrodynamic data—velocity and water level—are read, and the travel time is integrated for all reachable paths from upstream section i to downstream section j at piecewise constant velocity:

[0245] T _ij，t ≈∑ _e∈path(i→j) Δx _e u _t (e);

[0246] When multiple reachable paths exist, the minimum travel time is selected; a lower limit is set for extremely low flow rate segments to avoid numerical divergence, and abrupt changes are smoothed out. The travel time of unreachable pairs is recorded as infinity or marked as missing. This travel time can be used to construct physical windows or to align weights with freshness decay.

[0247] Accordingly, the travel time is read to construct a physically reachable time window [tT] for the downstream section j at time t. _ij，t -δ, tT _ij，t +δ]. δ is used to cover hydrodynamic uncertainties, and is taken as δ = κ × IQR(T) _ijThe window is set to _{tW:t+W} or according to the flow rate variance; for cases where it is unreachable or the travel time is missing, the window is set to empty. A causal mask (binary mask) M(i,j,t,t') is generated for each upstream time index t', with 1 inside the window and 0 outside the window. It is used for feature alignment and weighting under physical constraints.

[0248] Accordingly, upstream explanatory variables (upstream metal concentration, sewage discharge records, rainfall, etc.) from the temporally sequenced capacity feature packet and the multi-source cleaning sequence data packet, as well as causal masks and travel times, are read and aligned and aggregated within a physically reachable time window. A freshness decay weight w is applied to samples with an upstream time offset Δ=tt'. _fresh =exp(-λ×Δ), introducing data quality weight w _quality =1 / (1+σ imp 2 ), where σ imp 2 For interpolation variance; the overall weight w = M × w _fresh ×w _quality For each eigencomponent x, construct the aligned value x. ~ _j (t)=Σ _i，t＇ (w(i, j, t, t') × x _i (t�)) / Σ _i，t＇ w(i, j, t, t'); obtains the aligned downstream feature flow (including upstream concentration priors, aligned capacity features, aligned sewage / meteorological explanatory quantities, etc.). Outputs the aligned downstream feature flow, which can be packaged into a travel-time aligned feature package.

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

Claims

1. A method for simulation and prediction of heavy metal concentration in water bodies, characterized in that, The method comprises the following steps: Integrating original monitoring data, hydrodynamic data, total suspended solids, image remote sensing data and human activity data to generate a multi-source cleaning sequence data package; Using the multi-source cleaning sequence data package to perform cross-modal adsorption capacity estimation to infer the chemical composition of particles and adsorption parameters, and generate a capacity feature package; Based on the multi-source cleaning sequence data package, calculate the travel time, and perform causal alignment on the capacity feature package and the upstream signal in the multi-source cleaning sequence data package to generate a travel time alignment feature package; Combine the travel time alignment feature package, the capacity feature package and the metal fingerprint parameters, apply physical constraints at the head of the morphological distribution constraint, decouple and generate a state concentration and total concentration prediction result package; Wherein, using the multi-source cleaning sequence data package to perform cross-modal adsorption capacity estimation to infer the chemical composition of particles and adsorption parameters, and generate a capacity feature package, comprising: analyzing the image remote sensing data in the multi-source cleaning sequence data package to obtain an image texture morphological feature vector; constrained regression is applied to it to infer the particle composition vector representing the chemical composition; according to the particle composition vector and the metal fingerprint parameters, the adsorption isotherm parameters are retrieved and corrected from the knowledge base; combined with the total suspended solids in the multi-source cleaning sequence data package and the adsorption isotherm parameters, the adsorption capacity is calculated to generate the capacity feature package; Based on the multi-source cleaning sequence data package, calculate the travel time, and perform causal alignment on the capacity feature package and the upstream signal in the multi-source cleaning sequence data package to generate a travel time alignment feature package, comprising: using the hydrodynamic data and human activity data in the multi-source cleaning sequence data package to determine the cross-section accessibility and calculate the travel time; based on the travel time, construct a physically accessible time window and generate a causal mask; apply the causal mask to the capacity feature package and the upstream signal from the multi-source cleaning sequence data package to apply weighted aggregation including freshness decay and data quality, and generate a travel time alignment feature package; Wherein, combining the travel time alignment feature package, the capacity feature package and the metal fingerprint parameters, applying physical constraints at the head of the morphological distribution constraint, decoupling and generating a state concentration and total concentration prediction result package, comprising: based on the metal fingerprint parameters and the travel time alignment feature package, establishing an initial morphological distribution ratio; applying the adsorption capacity upper bound contained in the capacity feature package and the alignment signal in the travel time alignment feature package to the initial morphological distribution ratio to obtain the capacity-constrained morphological distribution ratio; according to the capacity-constrained morphological distribution ratio, the capacity response function in the capacity feature package and the travel time alignment feature package, decoupling calculation of dissolved state, colloidal state and particle state concentration, and integration into a state concentration and total concentration prediction result package.

2. The method of claim 1, wherein, Applying constrained regression to the image texture morphological feature vector to infer the particle composition vector representing the chemical composition, comprising: Applying a regressor to the image texture morphological feature vector to obtain an original composition vector; Performing a simplex projection on it to make the particle composition vector satisfy the physical constraints of non-negativity and sum of one.

3. The method of claim 1, wherein, According to the particle composition vector and the metal fingerprint parameters, retrieve and correct the adsorption isotherm parameters from the knowledge base, comprising: Retrieving the reference isotherm parameters from the knowledge base according to the particle composition vector; A competition correction coefficient γ is introduced, which is constructed based on the metal fingerprint parameters _m，m＇ The reference isotherm parameters are reduced to obtain the adsorption isotherm parameters considering the competition effect of multiple metals The competition correction factor is then used to reduce the affinity coefficient to obtain the effective affinity coefficient b _eff_m ; wherein ; b _eff_m =b _m / (1+Σ _m＇≠m γ _m，m＇ ×a _m＇ ); where b _m is the reference affinity coefficient for the target metal m; b _m＇ is the reference affinity coefficient for the competing metal m'; θ and φ are dimensionless exponent parameters calibrated by a small amount of competitive adsorption experiments; a _m＇ is the effective activity of the competing metal m'; a _ref is the reference activity.

4. The method of claim 1, wherein, Analyzing image remote sensing data in a multi-source cleaning sequence data package to obtain an image texture morphological feature vector, comprising: Performing denoising, flat field correction and scale standardization on the image remote sensing data; Extracting multi-scale texture and morphological features to form the image texture morphological feature vector, the features including: gray level co-occurrence matrix statistics, local binary pattern histogram, morphological granularity spectrum and fractal dimension.

5. The method of claim 1, wherein, Constrained reshaping of the initial morphological distribution ratio, including: Aggregating upstream contributions with local characteristic quantities from the travel time alignment feature package to construct the total quantity prior C _proxy_m ; Utilizing the total prior and the upper bound of adsorption capacity Q in the capacity feature bag _max_m , calculating the physical upper bound of the particle distribution ratio π _p_m_cap , obtaining the morphological distribution ratio under the capacity constraint ; wherein ; ； where min() is a minimization objective function, ε is a numerical stability constant, π _p_m_init is an initial proportionality, η and ρ are learnable parameters for controlling the strength of suppression and nonlinearity, and S is a capacity saturation.

6. The method of claim 1, wherein, Decoupling calculation of dissolved, colloidal and particulate concentrations to synthesize a sub-state concentration and total concentration prediction result package, including: Using travel time alignment feature package to obtain total concentration unconstrained estimation; C is estimated based on the total concentration without constraints _m_tilde a mass balance equation for the dissolved concentration is established in terms of the capacity response function in the capacity signature package, the upper bound on adsorption capacity, and the colloidal fraction coefficient determined from the metal fingerprint parameters; Solving the mass conservation equation, the dissolved concentration C is obtained _d_m ; Based on the dissolved concentration solution to calculate the colloidal and particulate concentrations C _p_m , synthetic partial concentration and total concentration prediction results package; wherein C _m_tilde = C _d_m + C _c_m + C _p_m ; C _c_m =k _c_m ×C _d_m ; C _p_m = min(q _m (C _d_m ) × TSS, Q _max_m ); where q _m is the capacity response function, TSS is the total suspended solids concentration, Q _max_m is the particle capacity upper bound and k _c_m is the colloidal proportionality coefficient, C _c_m corresponds to the colloidal concentration.

7. The method of claim 1, wherein, The weighted aggregation employs a comprehensive weight w(i,j,t,t'), w(i,j,t,t') = M(i,j,t,t') x w _fresh (Δ) x w _quality (σ imp 2 ); wherein w _quality = 1 / (1 + σ imp 2 ); w _fresh = exp(-λ x Δ); Δ = |(t-t') - T _ij，t |; where M(i, j, t, t') is a causality mask, taking values of 0 or 1, w _fresh (Δ) is a freshness decay weight, i, j are index numbers of upstream and downstream sections, t is a target time instant for predicting or constructing features for downstream section j, t' is an upstream time instant, T _ij,t is a time-varying travel time, λ is a freshness decay parameter, Δ is a signal lag, σ imp 2 is an interpolation variance, w _quality (σ imp 2 is a data quality weight.

8. The method of claim 1, wherein, Using hydrodynamic data and human activity data in the multi-source cleaning sequence data package to determine cross-section accessibility and calculate travel time, including: Calculating convective travel time according to the convection mode in the hydrodynamic data; Monitoring the situation where the convective travel time exceeds the preset physical upper limit due to extremely low flow rate; When the above situation is triggered, enable the diffusion kernel fallback mechanism, estimate the diffusion fallback travel time based on the effective diffusion coefficient, and use the diffusion fallback travel time as the travel time.

Images