Method for simulating and predicting concentration of heavy metals in water body

Abstract

The invention discloses a water heavy metal concentration simulation and prediction method, which comprises the following steps: integrating original monitoring data, hydrodynamic data, total suspended solids, image remote sensing data and human activity data, and generating a multi-source cleaning sequence data packet; executing cross-modal adsorption capacity estimation by using image remote sensing data in the data packet, inferring particle chemical composition and adsorption isotherm parameters from image textures, and generating a capacity feature packet containing an adsorption capacity upper bound; time-varying travel time is calculated based on the hydrodynamic data and the human activity data, causal alignment is performed on the capacity feature packet and the upstream signal, and a travel time alignment feature packet is generated; and in combination with metal fingerprint parameters, applying an adsorption capacity upper bound as a physical constraint on a form distribution constraint head, explicitly decoupling and predicting the form, and generating a prediction result packet. According to the method, the hydrodynamic physical mechanism and the particle adsorption chemical mechanism are deeply coupled, and the prediction precision and the physical consistency of the model under the unsteady state condition are improved.

Description

TECHNICAL FIELD

[0001] The application belongs to the technical field of water environment monitoring, and particularly relates to a water body heavy metal concentration simulation and prediction method. BACKGROUND

[0002] Heavy metal pollution is one of the key environmental problems threatening the safety of water ecological systems. Due to the high toxicity, non-degradability and bioaccumulation of heavy metals, their migration and transformation processes in water bodies are complexly affected by multiple factors such as hydrology, meteorology and biochemistry. Therefore, constructing a technical means that can dynamically simulate and accurately predict the concentration of heavy metals in water bodies can realize timely early warning of pollution events, water environment capacity assessment and watershed risk control.

[0003] At present, for the simulation of water body heavy metal concentration, various technical paths have been developed. One is the mechanism model based on physical and chemical processes, such as the water quality migration and diffusion model, which describes the spatio-temporal evolution of pollutants by solving the convection-diffusion equation set, usually requiring detailed river network topology, cross-section parameters and accurate boundary conditions. Another is the data-driven statistical or machine learning model, such as using the autoregressive integrated moving average (ARIMA) model or long short-term memory network (LSTM) to make time prediction on the historical concentration sequence of a single monitoring section; or using spatial interpolation technology (such as the Kriging method) to estimate the concentration of different sections at the same time.

[0004] However, the existing technology still has the problem that the model lacks real-time chemical mechanism constraints when simulating non-steady-state hydrodynamic and complex biochemical coupling processes, and it is difficult to achieve causal alignment of upstream and downstream signals under dynamic water flow conditions. Therefore, further research and innovation are needed to solve the above problems existing in the prior art. SUMMARY

[0005] The application provides a water body heavy metal concentration simulation and prediction method.

[0006] Technical scheme: According to one aspect of the application, a water body heavy metal concentration simulation and prediction method comprises:

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

[0008] Perform cross-modal adsorption capacity estimation using the multi-source cleaned sequence data package to infer particle chemical composition and adsorption parameters, and generate a capacity feature package;

[0009] Calculate the travel time based on the multi-source cleaning sequence data packet, perform causal alignment on the upstream signal in the capacity feature packet and the multi-source cleaning sequence data packet, and generate a travel time alignment feature packet;

[0010] In combination with the travel time alignment feature packet, the capacity feature packet and the metal fingerprint parameter, a physical constraint is applied in the morphology distribution constraint head, decoupling and generating a state concentration and total concentration prediction result packet.

[0011] Beneficial effects: the present application deeply couples the hydrodynamic physical mechanism and the particle adsorption chemical mechanism, improves the prediction accuracy and physical consistency of the model under non-steady state conditions. Related technical effects will be described in detail below in conjunction with specific embodiments. BRIEF DESCRIPTION OF DRAWINGS

[0012] Figure 1 A flowchart of a water body heavy metal concentration simulation and prediction method provided by the embodiment of the present application.

[0013] Figure 2 A flowchart of generating a capacity feature packet provided by the embodiment of the present application.

[0014] Figure 3 A flowchart of applying a constrained regression to an image texture morphology feature vector and inferring a particle composition vector representing chemical composition provided by the embodiment of the present application.

[0015] Figure 4 A flowchart of analyzing image remote sensing data in the multi-source cleaning sequence data packet and obtaining an image texture morphology feature vector provided by the embodiment of the present application.

[0016] Figure 5 A flowchart of decoupling and generating a state concentration and total concentration prediction result packet provided by the embodiment of the present application. DETAILED DESCRIPTION

[0017] In order to enable the personnel in the technical field to better understand the present application scheme, the technical solutions in the embodiments of the present application will be described clearly and completely in conjunction with the drawings in the embodiments of the present application. Obviously, the described embodiments are only a part of the embodiments of the present application, not all. Based on the embodiments in the present application, all other embodiments obtained by the person skilled in the art without creative labor should belong to the scope of protection of the present application.

[0018] It is to be understood that the terminology used herein is for the purpose of describing particular embodiments only and is not intended to be limiting, unless the context dictates otherwise. It will be further understood that the use of relational terms such as first and second, and the like are used solely to distinguish one from another entity without necessarily implying a relationship or order. It is to be understood that where the disclosure is specifically directed to embodiments, examples, or the like that include, but do not require, concurrent, serial, or sequential steps or procedures that can have been performed chronologically, in some cases, the steps or procedures can be performed in different order, concurrently, or with other steps or procedures intervening. Also, the terms "comprises," "comprising," "includes," "including," "contains," "containing," "has," "having," "united," "unites," "coupled," "coupling," and the like can be used synonymously, unless the context dictates otherwise.

[0019] To solve the above problems, the applicant has conducted in-depth retrieval and analysis, and found that:

[0020] On the one hand, the existing model usually simplifies the adsorption-desorption process of heavy metals, or uses static and offline calibrated distribution coefficients. It is impossible to dynamically infer the upper bound of the total adsorption capacity of the water body according to the real-time chemical composition (such as the proportion of organic matter and clay) and physical form (which can be represented by image texture) of the suspended solids and total suspended solids concentration (TSS) in the water body.

[0021] On the other hand, when dealing with the spatiotemporal dependence of upstream and downstream sections, the existing method usually uses fixed time lag or pure statistical correlation. When the flow velocity of the river changes greatly due to gate regulation or rainfall, the fixed time lag window cannot accurately capture the time-varying travel time of the pollutants along with the water flow, resulting in a causal misplacement between the input characteristics and the prediction target, and affecting the response accuracy of the model to sudden pollution events. The lack of chemical mechanism and the misplacement of hydrodynamic physical process lead to unreliable prediction of the distribution proportion of key forms such as dissolved state, colloidal state and particulate state, especially in extreme situations of high concentration or high TSS, and the prediction results easily violate the physical upper bound of adsorption saturation.

[0022] To solve these problems, in combination with Figures 1 to 5 The present application is specifically illustrated by the following embodiments.

[0023] Embodiment one, the present application relates to the field of water environment monitoring and data analysis, especially a water body heavy metal concentration simulation and prediction method.

[0024] In the context of the present application, in order to keep the description consistent, some technical terms and data items in some embodiments are defined:

[0025] Among them, # represents the hat symbol. Я represents the partial derivative.

[0026] C _d_m_fallback , C _c_m_fallback、 C_p_m_fallback Reserve / alternative concentration representing the concentration of the three states (dissolved state, colloidal state and particulate state) respectively. C t obs Actual observed concentration.

[0027] q _max_m Maximum adsorption capacity per unit mass of particles, q _max_clay , q _max_org Maximum adsorption capacity of the target substance corresponding to a single adsorption component (such as pure clay, pure organic matter) (i.e. the maximum amount of adsorption per unit mass when the component exists alone), w _clay , w _org Corresponding proportion.

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

[0029] Multi-source data can comprehensively capture complex factors affecting heavy metal migration and transformation; can include:

[0030] Raw monitoring data, for example, historical time series data of heavy metal concentration (such as hourly or daily lead, chromium, arsenic, etc. concentration values).

[0031] Hydrodynamic data, for example, variables reflecting the physical movement of the water body, such as flow rate, water level, etc. at the monitoring section.

[0032] Total suspended solids (TSS), i.e. the concentration of total suspended solids in the water body, one of the main carriers of heavy metal adsorption.

[0033] Image remote sensing data, for example, water body images collected through satellite remote sensing or on-site camera equipment. These image data can be used to analyze the spatial distribution of suspended solids, image texture of particulate matter, etc.

[0034] Human activity data, for example, variables that have a significant impact on hydrological regime, such as records of the opening and closing status of upstream reservoirs or gates, records of industrial pollution along the coast, etc.

[0035] Further, traditional water body heavy metal monitoring mainly relies on on-site sampling and laboratory chemical analysis, which has inherent defects such as high cost, long period, and difficulty in high-frequency and large-scale monitoring. At the same time, the monitoring data is discrete in time and space, and cannot capture the dynamic change process of pollution. Although there are predictions based on traditional mechanism models (such as water quality migration and diffusion models), such models usually require accurate fluid dynamics parameters and complex boundary conditions, making modeling difficult, and the response to sudden pollution events is slow.

[0036] The method can be implemented on various computing devices. For example, the method can be deployed in one or more servers, workstations, personal computers, embedded systems, or cloud platform environments. Illustratively, the implementation environment can include a data acquisition module, a data storage module, and a prediction module. The data acquisition module is responsible for obtaining data from various sensors, remote sensing satellites, or third-party data platform interfaces. The data storage module (e.g., a database) is used to store massive amounts of historical data and real-time streaming data. The prediction module loads and executes the method of the present application to simulate and predict future heavy metal concentrations.

[0037] Optionally, an online dynamic updating mechanism is provided. According to conventional methods in the art, the latest multi-source data is continuously obtained over time, and pre-processing and feature extraction processing are performed. The model parameters are updated according to the new data, for example, the single-source model is retrained with new data every week or month, or the weight calculation of Bayesian fusion is updated. During the online updating process, the old model and the new model are retained at the same time, and the performance of the two models in actual application is compared (such as comparing the deviation rate of the predicted value and the measured value) through A / B testing, so that the transition is smooth, and the prediction accuracy is optimized. It should be understood that the parameters in the water environment and the external influences are dynamically changing, and the online updating mechanism enables the model to adapt to environmental changes and new data characteristics, maintaining the accuracy and stability of long-term prediction.

[0038] Embodiment two, a water body heavy metal concentration simulation and prediction method is provided, in this embodiment, the method can be executed by the prediction module described above.

[0039] Step S2.1, integrate the 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. This provides a basis for subsequent analysis.

[0040] Specifically, data from various sources, different modalities, and different sampling frequencies are collected and pre-processed. The pre-processing process can include time alignment (e.g., unified to hourly or daily granularity), quality control (e.g., removing obvious outliers), and interpolation of missing data. A unified, spatio-temporally aligned training sample stream, i.e., a multi-source cleaning sequence data package, is formed.

[0041] Step S2.2, perform cross-modality adsorption capacity estimation using the multi-source cleaning sequence data package, infer the chemical composition of particles and adsorption parameters, and generate a capacity feature package.

[0042] In this embodiment, the physico-chemical properties of the particulate matter, such as the type composition of the particles (e.g. the proportion of clay, organic matter, silt) are estimated from the image remote sensing data (visual features) in the data packet, through a cross-modal inference model. On this basis, the adsorption isotherm parameters of the particulate composition under a specific metal are further retrieved or corrected from the knowledge base. On this basis, combined with the total suspended solids (TSS) concentration, the theoretical maximum adsorption capacity (i.e. the upper bound of total adsorption capacity Q _max_m ) of the water body for the metal and the current saturation degree are calculated. The inferred chemical parameters are organized into a capacity feature packet, providing chemical mechanism constraints for subsequent steps.

[0043] Step S2.3, calculate the travel time based on the multi-source cleaning sequence data packet, and causally align the capacity feature packet and the upstream signals in the multi-source cleaning sequence data packet to generate a travel time aligned feature packet. This step introduces the physical constraints of hydrodynamics.

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

[0045] Next, based on the travel time, the upstream signals are causally aligned. For example, the concentration at time t downstream, whose influence source (such as upstream pollution, upstream adsorption capacity state) should occur around (t-T _ij，t ). This step constructs a physically reachable time window to weight and aggregate the upstream signals (including the capacity feature packet, pollution records, etc.), so that the model does not mistakenly use future information or time-displaced correlations. The aggregated result is the travel time aligned feature packet.

[0046] Step S2.4, combine the travel time aligned feature packet, the capacity feature packet and the metal fingerprint parameters, apply physical constraints in the morphological distribution constraint head, decouple and generate the state concentration and total concentration prediction result packet. This step is the final prediction and decoupling step. The metal fingerprint parameters refer to parameters that characterize the specific chemical behavior of each metal (such as valence state, complexing ability, distribution coefficient, etc.). The morphological distribution constraint head is a special prediction model structure that does not directly predict the total concentration, but explicitly models the three main forms of heavy metals: dissolved state, colloidal state and particulate state. When performing morphological distribution, the head will apply the physical constraints in the capacity feature packet, for example, the upper limit of the particulate concentration cannot exceed the upper bound of the total adsorption capacity of the water body. At the same time, the alignment features in the travel time aligned feature packet are used as input. On this basis, the head outputs the three-state state concentration, which is combined into the total concentration to form the state concentration and total concentration prediction result packet.

[0047] It can be understood that in the embodiment, the capacity feature package generated by the cross-modal adsorption capacity estimation (specifically, the adsorption capacity upper bound) will serve as the upper limit of the physical constraint of the particle state concentration in the modal allocation constraint header, realizing the dual constraint of the physical mechanism (travel time alignment) and the chemical mechanism (capacity constraint).

[0048] Embodiment three describes an optional implementation of data integration and uncertainty annotation. As an optional implementation, the original monitoring data, hydrodynamic data, total suspended solids, image remote sensing data, and human activity data are integrated to generate a multi-source cleaning sequence data package, which also includes the spatio-temporal alignment and cleaning of multi-source heterogeneous data.

[0049] Correspondingly, a unified time axis and reference time granularity (e.g., 1 hour) are established, and data of different sampling frequencies (e.g., daily average hydrodynamic data, image data collected once a week) are resampled or window aggregated. At the same time, a unique spatial index ID is established for each monitoring section.

[0050] Further, data conflicts and delays can also be handled during the alignment process. For example, when there are multiple source records at the same section and the same time (such as simultaneous sensor direct reading and manual supplementary recording), a primary and secondary source arbitration rule can be used (for example, monitoring instrument automatic reading is prior to manual recording, and manual recording is prior to interpolation estimated value). In addition, the inherent collection delay and clock drift of different sensors can also be corrected to ensure the phase consistency of cross-source events (such as rainfall and flow rate change).

[0051] On this basis, this step performs missing value interpolation and uncertainty annotation on the cleaned data, which can specifically include:

[0052] Step S3.1, for the missing values in the original monitoring data, a local K-nearest neighbor and Bayesian autoregressive joint generation interpolation mean is used.

[0053] In the embodiment, the local K-nearest neighbor algorithm is used to capture the spatial correlation between sections, that is, the synchronous data of K sections with similar geographical positions is used to estimate the missing value. The Bayesian autoregressive model is used to capture the time series dependence of the section itself. The joint use of the two methods can generate a relatively robust interpolation mean within a strict time rolling window (even if it does not use future time information).

[0054] Step S3.2, the interpolation variance of the interpolation mean is estimated in a multiple interpolation manner. The interpolation variance (σ imp 2 ) is a quantitative indicator of the uncertainty or reliability of the interpolation mean. Multiple interpolation is an optional way to achieve this estimation.

[0055] On this basis, for the same segment of missing data, the imputation process is independently repeated M times (for example, M = 5 or 10 times), and each time the imputation process produces slightly different imputation values due to the randomness of the model (such as Bayesian posterior sampling). On this basis, by calculating the variance of the M times of imputation results, the true imputation variance of the missing value is estimated. If the variance is large, it means that the imputation credibility of this point is very low.

[0056] Step S3.3, incorporate the imputation variance into the multi-source cleaned sequence data package as the basis for calculating the data quality weight.

[0057] On this basis, the imputation mean, the binary label indicating whether it is an imputed value, and the estimated imputation variance are combined with all other cleaned, aligned and standardized data (such as hydrodynamic data, image features, etc.) to form the final multi-source cleaned sequence data package. The imputation variance field in the data package will be used to calculate the data quality weight in the subsequent alignment step, to automatically reduce the influence of those low credibility (high imputation variance) data when aggregating upstream signals in the model.

[0058] Embodiment Four, describes the specific implementation process of cross-modal adsorption capacity estimation (VCAP), which is used to establish a bridge between water body visual features (images) and chemical adsorption capacity. Exemplarily, the present embodiment can be implemented as follows:

[0059] Step S4.1, parse the image remote sensing data in the multi-source cleaned sequence data package to obtain the image texture and morphology feature vector. Correspondingly, read the image remote sensing data (for example, particulate matter images collected by on-site camera equipment or remote sensing satellite images) from the multi-source cleaned sequence data package.

[0060] Optionally, denoising, flat field correction and scale standardization are performed on the image remote sensing data. In some embodiments, preprocessing is performed to accurately extract effective features from the original image; specifically including: applying median filtering or bilateral filtering for image denoising to eliminate random noise; performing flat field correction or histogram equalization to eliminate the influence of uneven illumination; performing threshold self-adaptation and morphological opening and closing operation for particle region identification and segmentation.

[0061] On the preprocessed image, multi-scale texture and morphology features are extracted to form the image texture and morphology feature vector. The features can optionally include: gray level co-occurrence matrix statistics (such as contrast, homogeneity, energy, entropy, used to describe the fineness and regularity of the texture), local binary pattern histogram (such as rotation invariant pattern, used to describe local texture), morphological granularity spectrum or equivalent diameter distribution (as a proxy indicator of particle size distribution), and fractal dimension (used to describe the complexity of the particle edge).

[0062] Step S4.2, apply constrained regression to it (image texture morphological feature vector) to infer a particle composition vector that characterizes the chemical composition.

[0063] In this embodiment, the particle composition vector is used to represent the relative proportions of different chemical components in the total suspended solids of a water body. Exemplarily, the vector can be represented as {w _clay , w _org , w _silt , w _sand}, where w _clay is the proportion of clay, w _org is the proportion of organic matter, w _silt is the proportion of silt, and w _sand is the proportion of sand. The adsorption capacity of such components (especially clay and organic matter) for heavy metals is much higher than that of other components.

[0064] Constrained regression is used to establish a mapping relationship from the image texture morphological feature vector (visual domain) to the particle composition vector (chemical domain). Its constrained nature is reflected in that: applying a regressor to the image texture morphological feature vector to obtain an original composition vector; next, performing a simplex projection on it (original composition vector) to make the particle composition vector satisfy the physical constraints, specifically, the constraints of non-negativity and sum to one (for example, w _clay +w _org +w _silt +w _sand =1, and each item is greater than or equal to 0).

[0065] As an optional implementation, the constrained regressor can be implemented as a cross-modal mapping network, such as a deep neural network. The network can include: one or more feature extraction layers (such as convolutional layers or fully connected layers) for receiving the image texture morphological feature vector; one or more cross-modal mapping layers (such as fully connected hidden layers) for learning the non-linear mapping of visual features to chemical composition; and a constraint projection layer (for example, using a Softmax activation function or an explicit simplex projection operator) to make the output particle composition vector satisfy the physical constraints of non-negativity and sum to one.

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

[0067] Specifically, according to the particle composition vector, the reference isotherm parameters are retrieved from the knowledge base. Or, according to the particle composition vector (such as {w _clay , w _org,...} ) are retrieved or regressed from a pre-established parameter knowledge base. The knowledge base stores adsorption parameters (e.g., maximum adsorption capacity q _max_base and affinity coefficient b _base for different single components (e.g., pure clay, pure organic matter) on different metals (e.g., lead, chromium) (or two characteristic parameters K _base and n _base for Freundlich model). The reference isotherm parameter set for mixed particles is obtained by weighted average of the particle composition vector (e.g., q _max_m = w _clay × q _max_clay + w _org × q _max_org +... ).

[0068] On this basis, a competition correction coefficient is introduced, which is constructed based on the metal fingerprint parameters, to reduce the reference isotherm parameters (especially the affinity coefficient) and obtain the adsorption isotherm parameters considering the multi-metal competition effect.

[0069] As an optional implementation, the competition correction coefficient γ _m，m＇ (indicating the competition intensity of metal m' on metal m) can be constructed as follows:

[0070] ;

[0071] where b _m is the reference affinity coefficient of the target metal m (from the reference isotherm parameters); b _m＇ is the reference affinity coefficient of the competing metal m'; θ and φ are dimensionless index parameters calibrated by a small amount of competition adsorption experiments; a _m＇ is the effective activity of the competing metal m' (which can be calculated according to its concentration, pH, ion strength, and other hydrochemical conditions); a _ref is the reference activity (e.g., 1 mol / L). The competition correction coefficient is then used to reduce the affinity coefficient, for example, the effective affinity coefficient b _eff_m after correction can be calculated as: b _eff_m = b _m / (1 +∑ _m＇≠m γ _m，m＇ × a _m＇ ); where ∑ represents the summation of 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 there is no observed data, then the observed saturation is not calculated. Instead, the total mass prior C _proxy_m is constructed based on upstream information in the following steps, and used to generate soft anchor ranges (e.g. S _m_anchor ) that belong to the interval [0, min(C _proxy_m / Q _max_m , 1)] or skip this item.

[0079] On this basis, the calculated theoretical capacity upper bound, capacity response function (and its specific parameters q _max_m , b _eff_m , K _m , n _m ), and optional anchor S _m_anchor and its confidence are packaged together as a capacity feature package. This feature package will be aligned as an upstream signal, and can also be used as a constraint for morphological allocation.

[0080] Example Five, describes an exemplary scheme based on water-driven travel time alignment, which introduces the physical constraints of hydrodynamics, so that the prediction of downstream concentration is based on its physically reachable upstream influence source. Specifically as follows:

[0081] Step S5.1, use the water dynamics data and human activity data in the multi-source cleaning sequence data package to determine the reachability of the section and calculate the travel time. In other words, based on the water dynamics data and human activity data, determine whether the section is reachable, and further calculate the travel time.

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

[0083] Optionally, according to the topological relationship of the river network, the real-time water level difference (from the water dynamics data) and the gate opening and closing state (from the human activity data), the reachability of the upstream section i to the downstream section j at time t is determined. For example, when the gate is closed or the reverse pressure occurs, it is determined that the water flow cannot reach.

[0084] When it is determined to be reachable, the convective travel time is calculated in a convective manner. One optional calculation method is to use the piecewise constant speed approximation 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 sections e on the path from i to j; Δx _eis the length of the river reach e; u _t (e) is the average flow velocity of the river reach at time t (from hydrodynamic data).

[0086] Further, the case where the advection travel time exceeds a pre-set physical upper limit (e.g. T _max , which can be estimated from the length of the river reach and the historical p95 velocity) due to extremely low flow velocity is monitored. _t (e) the problem that the travel time tends to infinity when the flow velocity is extremely low,

[0087] When the case is triggered (i.e. T _ij，t > T _max ), a diffusion kernel fallback mechanism is enabled. At this time, the advection travel time is no longer used, but the diffusion fallback travel time T _eff is estimated based on the effective diffusion coefficient K _diff , which can be given by historical data matching or an empirical library, e.g. T _diff ≈ L _ij 2 / (2 x K _eff ), where L _ij is the total path length between i and j. On this basis, the diffusion fallback travel time is used as the travel time. Further, in order to suppress abnormal spikes in the travel time series, time smoothing, such as median smoothing or exponential smoothing, can also be performed on the calculated travel time series.

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

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

[0090] Wherein, the window width δ is used to cover the uncertainty of the hydrodynamic model and the dispersion effect of the water body. Optionally, the window width can be determined by an adaptive rule, e.g. δ = κ x IQR(T _ij_history ), where κ is an adjustment coefficient (e.g. 1.5) and IQR(T _ij_history ) is the interquartile range of the travel time from upstream section i to downstream section j within the historical window. When the diffusion kernel fallback is triggered, κ can be appropriately enlarged (e.g. κ _diff > κ) to cover greater diffusion uncertainty.

[0091] The causal mask M(i, j, t, t') is located at the upstream time t' in the time window [t-T _ij，t- δ, t - T _ij，t = 1 (or a weight value based on the position within the window, such as a Gaussian weight) if it is within the time window, and strictly 0 if it is outside the time window.

[0092] Further, if it is determined that the current period is not reachable (travel time is invalid and no diffusion fallback), a full zero mask is generated for this (i, j, t); if it is intermittently reachable (e.g. recently reachable but not currently reachable due to gate closure), the last window in which it was reachable can be kept as a weak reference, and the mask weight is downgraded.

[0093] Step S5.3, applying the causal mask, a weighted aggregation including freshness decay and data quality is applied to the capacity feature package and upstream signals (e.g. upstream pollution records, upstream metal concentration, upstream TSS) derived from the multi-source cleaning sequence data package. In some embodiments, the causal mask, the freshness decay weight and the data quality weight are combined to form a comprehensive weight; the comprehensive weight is applied to aggregate the upstream signals. All relevant upstream signals (possibly from multiple upstream sections i and multiple time instants t' within the physical time window) are aggregated into a feature vector x ~ _j (t) for the downstream section j at time t. The weighted aggregation employs the comprehensive weight w. The comprehensive weight w is optionally obtained by multiplying three components:

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

[0095] where M(i, j, t, t') is the causal mask (0 or 1). w _fresh (Δ) is the freshness decay weight, the farther in time from the travel time, the lower the information value of the signal; specifically, the freshness decay weight w _ij，t = exp(-λ x Δ) can be derived based on the travel time to calculate the signal lag Δ = |(t - t') - T _fresh |, where λ is the decay parameter. _quality (σ imp 2 ) is the data quality weight. It is used to reduce the influence of interpolated data or abnormal data. Specifically, the interpolation variance (σ imp 2 ) in the multi-source cleaning sequence data package resulting from the missing data interpolation step is extracted, based on which the data quality weight is derived, for example w _quality = 1 / (1 + σ imp 2 ).

[0096] Optionally, the freshness decay parameter initial value λ is set_0 approximately equal to 1 / median(T _ij，t ) (i.e. scaled by the inverse of the median travel time); next, the optimal value is selected among multiple candidate values (e.g. {0.25λ _0 , 0.5λ _0 , λ _0 , 2λ _0 , 4λ _0} ) via a time-series rolling validation.

[0097] In another alternative implementation, a set of feature-specific freshness parameters can be employed. For example, a larger λ is used for fast-varying variables (e.g. rainfall, gate position change) so that they decay faster, and a smaller λ is used for slow-varying variables (e.g. upstream concentration related to sediment release) so that they decay slower.

[0098] On this basis, the aligned downstream feature x ~ _j (t) can be obtained via 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 upstream i section at time t' (e.g. adsorption capacity upper bound from the capacity feature pack, or discharge record from the multi-source flushing pack).

[0101] Step S5.4, merge the aligned downstream feature stream generated from the weighted aggregation with the travel time and the causal mask to generate a travel time aligned feature pack for the morphological allocation constraint header. This step is an encapsulation step. The aligned downstream feature stream, the travel time and the causal mask are merged to form a structured data pack, i.e. the travel time aligned feature pack. This feature pack will be directly input to the morphological allocation constraint header.

[0102] According to an aspect of the present application, a case of computing time-varying travel time is provided, as follows:

[0103] Suppose a pollutant needs to flow from upstream section i to downstream section j. The path is composed of two consecutive river reaches e _1 and e _2 . Reach e _1 has a length Δx _e_1 = 18,000 meters (18 km); reach e _2 has a length Δx _e_2= 12,000 meters (12 km); total path L _ij = 30,000 meters (30 km). According to historical hydrological data, the default physical upper limit T _max = 172,800 seconds (48 hours). When triggered, the effective spread factor K _eff = 60 m 2 / s.

[0104] When t = 10:00 (normal flow rate), at this time, the water flow is smooth, and the gate is open. According to the hydrodynamic data, e _1 segment flow rate u _t=10 (e _1 ) = 1.5 m / s; e _2 segment flow rate u _t=10 (e _2 ) = 1.2 m / s.

[0105] Calculated in a convection manner. e _1 segment time T _e_1 = Ax _e_1 / u _t=10 (e _1 ) = 18,000 / 1.5 = 12,000 seconds. e _2 segment time T _e_2 = Ax _e_2 / u _t=10 (e _2 ) = 12,000 / 1.2 \text{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 (about 6.1 hours). T _ij，t=10 (22,000 s) < T _max (172,800 s), not exceeding the upper limit. The convection travel time T _ij，t=10 = 22,000 seconds is adopted.

[0106] When t = 14:00 (upstream gate opening), at this time, the upstream gate opening causes the flow rate to increase. According to the hydrodynamic data, e _1 segment flow rate u_{t=14}(e _1 ) = 3.0 m / s; e _2 segment flow rate u _t=14 (e _2 ) = 2.5 m / s.

[0107] e _1 segment time T _e_1 = 18,000 / 3.0 = 6,000 seconds. e _2 segment 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,800 s) < T _max (172,800 s), not exceeding the upper limit.

[0110] Adopting convective travel time T _ij，t=14 = 10,800 seconds. It can be seen that T _ij，t varies with the time-varying characteristics of the flow velocity u _t .

[0111] When t = 18:00 (dry season and downstream gate closed), at this time, the water flow is very slow, close to stagnation. According to the 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.05 m / s.

[0112] Calculated in the convective way. e _1 section time T _e_1 = 18,000 / 0.1 = 180,000 seconds. e _2 section 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,000 s) > T _max (172,800 s).

[0115] Since the convective travel time exceeds the preset physical upper limit, the diffusion kernel rollback mechanism is triggered. The rollback calculation is enabled: T _diff ≈ L ij 2 / (2 x K _eff ) = (30,000) 2 / (2 x 60) ≈ 900,000,000 / 120 = 7,500,000 seconds (about 86.8 days).

[0116] Diffusive backtracking travel time T _ij，t=18 =7,500,000 seconds, the transmission of this period is dominated by diffusion (e.g., using a wider time window or lower weight) in subsequent weighted aggregation.

[0117] According to one aspect of the present application, the metal fingerprint parameter of each metal refers to a set of parameters used to characterize the intrinsic or slowly changing physicochemical specific properties of a specific heavy metal (such as lead, chromium, arsenic, etc.) in the water environment. The parameter is the basis for distinguishing the migration and transformation behavior (such as adsorption, complexation, precipitation) of different metals.

[0118] In optional embodiments of the present application, the metal fingerprint parameter includes, but is not limited to, one or a combination of the following:

[0119] Valence state distribution ratio, for example, for arsenic (As), the ratio of As(III) to As(V); for chromium (Cr), the ratio of Cr(III) to Cr(VI). The toxicity, solubility and adsorption of different valence states differ greatly.

[0120] Complexation stability constant family, which characterizes the stability constant of the metal ion forming a soluble or insoluble complex with common ligands in water (such as carbonate, chloride, sulfate, humic acid or dissolved organic carbon).

[0121] Carrier distribution coefficient family, which characterizes the distribution coefficient (Kd) or adsorption affinity sequence of the metal on different adsorption carriers (such as organic matter, clay minerals, iron and manganese oxides, etc.).

[0122] Precipitation / hydrolysis related parameters: for example, the solubility product (K _sp ) interval of the metal to form hydroxide, carbonate or sulfide precipitates, and the hydrolysis constant of the metal ion.

[0123] Redox sensitivity, which characterizes the influence factor of the change of the oxidation-reduction potential (Eh) of the water body on the valence state distribution or solubility of the metal.

[0124] Colloid affinity parameter, which characterizes the tendency of the metal to bind to colloidal particles (usually referring to particles with a particle size of less than 1 micron).

[0125] Sedimentation tendency parameter, which characterizes the empirical coefficient related to the sedimentation rate of particles (especially particulate) attached with the metal.

[0126] In some embodiments, the acquisition path of the metal fingerprint parameters can be diverse. Optionally, the typical values are acquired from professional literature, chemical manuals or public databases as priori intervals; if the on-site conditions allow, small-scale laboratory isothermal adsorption or competitive adsorption experiments can be carried out by collecting water samples and sediment samples to calibrate the key parameters (such as the maximum adsorption capacity of the unit mass of particles of the local TSS or the b value); for the parameters that cannot be directly measured, the inverse adjustment or Bayesian update within the priori interval can also be carried out in the model training through the historical data.

[0127] Further, since some parameters (such as the complexation constant, the solubility product) can slowly change with the water temperature, the pH and other environmental factors, or differ in different river sections (such as industrial and agricultural areas), these parameters can be spatially and temporally interpolated to form a time-series and spatialized metal fingerprint parameter set for real-time calling by the model. The metal fingerprint parameter set can be used to construct the competition correction coefficient; it can also be used to establish the initial proportion of the morphological distribution for determining the colloidal state proportion coefficient.

[0128] Embodiment six, describes an optional implementation scheme of the morphological distribution constraint head (MFCN_head), i.e. the metal fingerprint coupled morphological distribution constraint head. This embodiment integrates physical and chemical constraints to achieve robust decoupling of heavy metal forms. Specifically, it includes:

[0129] Step S6.1, in combination with the travel time alignment feature package, the capacity feature package and the metal fingerprint parameters, applying physical constraints in the morphological distribution constraint head to decouple and generate the state concentration and total concentration prediction result package.

[0130] Correspondingly, the morphological distribution constraint head is a special model structure (for example, a deep neural network module or a coupled calculation equation), the input of which includes the travel time alignment feature package (including the aligned upstream signal, TSS, water chemical interpretation quantity, etc.); the capacity feature package (including the adsorption capacity upper bound and the capacity response function); the metal fingerprint parameter set. The head is used to explicitly model and decouple the dissolved state concentration (C _d_m ), the colloidal state concentration (C _c_m ) and the particulate state concentration (C _p_m ).

[0131] Step S6.2, establishing the initial morphological distribution proportion.

[0132] Before applying the constraints, as an optional implementation, the head establishes the initial morphological distribution proportion based on the metal fingerprint parameter set and the alignment features (such as TSS, pH, dissolved organic carbon, etc.) in the travel time alignment feature package. _p_m_init The initial proportion can be learned by a regression sub-network, which represents the particulate state distribution tendency determined only by the chemical affinity and the water body conditions without considering the adsorption capacity saturation limit.

[0133] Step S6.3, constraint shaping on the initial morph fraction.

[0134] In some embodiments, the initial morph fraction is constraint shaped by applying the upper bound on adsorption capacity in the capacity feature package and the alignment signal in the travel time alignment feature package.

[0135] Optionally, the upstream contribution is aggregated with the local feature quantity from the travel time alignment feature package to construct the total quantity prior C _proxy_m Here, the total quantity prior is an estimated total concentration based on the upstream alignment signal and local features (e.g. rainfall, sewage discharge), which is only used as an internal variable for constraint calculation, not the final output.

[0136] Next, the total quantity prior C _proxy_m is used to calculate the upper bound on the particulate fraction π _max_m with the upper bound on adsorption capacity Q _p_m_cap in the capacity feature package.

[0137] Illustratively, the upper bound can be calculated as:

[0138] ;

[0139] where ε is a numerical stability constant to prevent the denominator from being zero. This avoids using the dissolved concentration. It should be understood that the circular dependency is caused by using the dissolved concentration to be solved in the solving process.

[0140] Further, to reflect the inhibitory effect when the adsorption capacity is close to saturation, the saturation constraint is introduced in this embodiment.

[0141] Correspondingly, when the capacity saturation S (S _m_anchor or the capacity saturation estimated by the total quantity prior and q _m function) approaches 1, the particulate fraction is dynamically adjusted downward by the saturation inhibition function.

[0142] Illustratively, the function can be expressed as:

[0143] ;

[0144] where π _p_m_init is the initial fraction obtained in the previous step, η (e.g. in the interval [0, 1]) and ρ (e.g. greater than or equal to 1) are learnable parameters for controlling the inhibition intensity and nonlinearity.

[0145] At the same time, the dissolved and colloidal fractions are adjusted compensatorily so that the sum of the three fractions is still 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 a non-negative monotonically increasing function of C _d_m , the right side of the mass conservation equation is monotonically increasing with respect to the dissolved concentration. This guarantees the uniqueness of the numerical root-finding.

[0157] On this basis, the dissolved concentration solution is obtained by solving the mass conservation equation by numerical root-finding method. In other words, the equation can be solved by numerical root-finding method (for example, bisection method, Newton method or Broyden method, with certain convergence guarantee), to obtain a unique dissolved concentration solution.

[0158] Further, based on the dissolved concentration solution, the colloidal and particulate concentrations are calculated, and the multi-state concentration and total concentration prediction result package is synthesized. That is, based on the obtained C _d_m , the following is calculated:

[0159] ;

[0160] ;

[0161] ;

[0162] The three-state concentrations and their sum (or the sum of the three) are collectively packaged as a multi-state concentration and total concentration prediction result package.

[0163] Step S6.5, failure fallback mechanism. In order to enable the system to still provide a solution when the input is extreme or the numerical calculation is unstable, the embodiment includes a fallback mechanism.

[0164] Correspondingly, the convergence state of the numerical root-finding method is monitored. When it is determined that it is not convergent, the failure fallback mechanism is enabled. The triggering conditions of non-convergence can optionally include that the number of iterations exceeds a preset upper limit N _max (for example, N _max = 100 iterations), or the relative error of the residual is greater than the convergence threshold ε (for example, ε = 1e-5).

[0165] When the fallback is triggered, a proportional locking linear approximation method is used to decouple the generation of the multi-state concentration and total concentration prediction result package based on the unconstrained estimate of the total concentration and the capacity-constrained morphological distribution proportion.

[0166] Correspondingly, the capacity-constrained morphological distribution proportion π _p_m and the unconstrained estimate of the total concentration C _m_tilde , the concentration is distributed 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, after preprocessing, feature extraction and standardization of data from different data sources (such as environmental indicators, social indicators, water quantity indicators, water quality indicators, image features), the data is spliced in the feature dimension. A unified, high-dimensional feature matrix is formed. Next, the unified feature matrix is taken as input to train a single prediction model (for example, 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 in an end-to-end manner. This method has a simple structure and can directly utilize information from all data sources; however, when there is a large amount of missing data in any data source, additional imputation processing is required, and the data standardization requirement is higher.

[0183] As another optional implementation, a decision-level fusion (Bayesian fusion) can be used. This method is suitable for handling cases where there are differences in sampling frequency, granularity and partial data missing between different data sources.

[0184] Correspondingly, for each independent data source (for example, one group is water power and water quality data, and the other group is image and remote sensing data), an independent single-source prediction model is constructed (for example, partial least squares PLS, LASSO or random forest can be used). For each single-source model, the mean square error (MSE) is calculated according to its historical prediction performance in the most recent time window (for example, the past i time points). The mean square error (MSE) is regarded as a measure of the prediction uncertainty of the model, and based on this, the contribution weight of each single-source model is dynamically allocated through the Bayesian weight calculation formula.

[0185] Exemplarily, the weight w _i (j) of the jth model at time i can be set as the inverse of its recent mean square error MSE _i|i-1 (j) for normalization w _i (j)=(1 / MSE _i|i-1 (j)) / ∑ _k (1 / MSE _i|i-1 (k)); where ∑ _k represents the sum of all k single-source models. The weighted prediction results of each single-source model are taken as the final fusion prediction output.

[0186] This decision-level fusion method can dynamically adjust the weight according to the real-time performance of different data source models, and has strong robustness. Even if part of the data source (i.e., a single-source model) is temporarily disabled or the data quality decreases (manifested as an increase in MSE), its weight will be automatically reduced, ensuring the stability and accuracy of the overall prediction.

[0187] Example 9. An exemplary scheme based on differentiable physical propagation and sparse source inversion is provided, which not only enables forward prediction, 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 to numerically discretize (e.g., using finite difference method or finite element method) the one-dimensional or two-dimensional convection-diffusion equation describing the migration of pollutants, so as to make it a differentiable computational operator.

[0189] Exemplarily, the operator can be represented in the form of matrix time step iteration as follows:

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

[0191] where C _t is a vector representing the heavy metal concentration at all river grid points at time t; C _t+1 is the concentration vector at t+1. A(u, K) is the state transition matrix, which is determined by the flow field u and the diffusion coefficient K of the river grid points (u and K can come from hydrodynamic calculation or data package); the A matrix encapsulates the physical processes of convection and diffusion. S _t is a vector representing the pollution source release intensity at each grid point at time t. B is the source injection matrix, which is used to distribute the source term S _t to the corresponding river grid points.

[0192] Since all matrix operations (multiplication, addition) in this equation are differentiable, this PDE layer can be embedded into a neural network to realize end-to-end training. In forward prediction, given S _t (e.g., known pollution records), the future C _t can be predicted by forward propagation at each time step.

[0193] Further, this embodiment can utilize the characteristics of the differentiable PDE layer to realize sparse multi-hypothesis source inversion. After observing the abnormal concentration of the downstream monitoring point, the unknown pollution source S _t (position, start and end time, and intensity) can be inferred backward. It can be constructed as an optimization problem, and its objective function can optionally include:

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

[0195] where S ≥ 0 is a physical constraint that the pollution source strength cannot be negative. Σ _t || C # _t (S) - C t obs || 2 2 is an L2 loss term that represents the Euclidean distance (i.e., prediction error) between the predicted concentration C # _t (S) and the actual observed concentration C t obs . λ _1 × || S || 1 _1 is an L1 regularization term (|| S || 1 _1 is the L1 norm of S _t ) that is used to impose a sparsity prior. Assuming that a sudden pollution source usually only occurs at a few locations and for a short time, the L1 regularization tends to produce a sparse solution (i.e., most elements in S _t are zero). λ _2 × Σ _t TV(S _t ) is a total variation (TV) regularization term that is used to impose a temporal smoothness prior. Assuming that the release intensity of a pollution source does not jump drastically between adjacent time instants. λ _1 and λ _2 are hyperparameters that are used to balance the weights of the terms.

[0196] Since the entire objective function is differentiable (or can be solved by proximal gradient descent), the optimization problem can be solved by a backpropagation algorithm that backpropagates from the observed C t obs to infer the most likely pollution source S _t .

[0197] Optionally, to improve the robustness of the inversion, a multiple hypothesis tracking can be used, i.e., the optimization problem is solved in parallel for multiple possible initial pollution locations or reporting time instants, and the Top-k source hypotheses are ranked and confidence output based on the final loss value and information criteria (e.g., AIC / BIC).

[0198] In another alternative implementation, the differentiable PDE layer can be combined with a machine learning model (e.g., a residual CNN), where the PDE layer is responsible for modeling the interpretable advection-diffusion dominated process, while the ML layer is used to learn and compensate for non-ideal effects (e.g., local backwater, temporary gate effects, etc.) that are not modeled by the PDE layer.

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

[0200] Suppose at a time t, for a metal m (e.g. lead) at a downstream j cross section, the model has completed the preceding steps and obtained the following input parameters (the numerical values are only examples and do not limit the present application):

[0201] Concentration unconstrained estimate C _m_tilde = 50.0 µg / L. Upper bound on adsorption capacity Q _max_m = 25.0 µg / L. Total suspended solids concentration TSS = 100.0 mg / L. Colloidal fraction coefficient 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 ). Suppose parameters are q _max_m = 0.3 mg / g, b _eff_m = 0.2 L / µg.

[0202] Further, calculate the particulate concentration C _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); 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 _p_m = 25.69 µg / L is slightly larger than the capacity upper bound Q _max_m = 25.0 µg / L. In the actual solver, due to the min() function, the solution will converge to the point where C _p_m = Q _max_m .

[0212] Correspondingly, the modified solution (the actual convergence point of the solver) is calculated. C _p_m = 25.0 µg / L (reached the capacity upper bound 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. The total concentration C _m = 50.0 µg / L.

[0214] Suppose the numerical root-finding method determines that it has not converged (e.g., the number of iterations exceeds N _max = 100).

[0215] The failure fallback mechanism is enabled. The system enables the linear approximation method of proportion locking. Suppose that the (constraint reshaping) has calculated the proportion of morphological distribution π _p_m = 0.45 (45%) under the capacity constraint. (This proportion itself has taken into account the constraints of Q _max_m and C _proxy_m ).

[0216] Using linear approximation decoupling, the particulate state C _p_m = π _p_m × C _m_tilde = 0.45 × 50.0 = 22.5 µg / L. The dissolved state 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. The colloidal state 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 the present application, the metal fingerprint parameters can also be (i.e., inverse tuning) include:

[0223] Accordingly, the typical value of the target parameter P (e.g., the stability constant K of a certain complex reaction or the solubility product Ksp of a certain precipitation _sp ) is obtained from professional literature, chemical manuals or databases, and a reasonable prior interval [P _min , P _max ] in the physical and chemical sense is set.

[0224] Optionally, the parameter P is regarded as a trainable weight of the model (e.g., the morphology distribution constraint head). In order to keep P within the prior interval during the training process, a constraint transformation can be used. For example, an unconstrained trainable variable P _trainable is defined inside the model, which is mapped to the target interval by a nonlinear function (such as a scaled Sigmoid function or a Tanh function), i.e.: P = P _min + (P _max - P _min ) x Σ (P _trainable ), where Σ (x) is a Sigmoid activation function with a value range of [0, 1].

[0225] During the model training process, the internal variable P _trainable is regarded as a trainable parameter together with other weights (W, b) of the neural network. When calculating the loss function (e.g., the threshold consistency loss function), the gradient of the total loss Loss with respect to P, ЯLoss / ЯP, is back-propagated to P _trainable by the chain rule, obtaining ЯLoss / ЯP _trainable . The optimizer (e.g., Adam or SGD) updates P _trainable using the gradient ЯLoss / ЯP _trainable . Based on this, the parameter P is automatically inverse tuned within its pre-set physical prior interval [P _min , P _max ] under the constraint of historical observation data, towards the direction of minimizing the total loss of the model, to obtain a parameter value that conforms to the characteristics of the local water environment.

[0226] According to still another aspect of the present application, the total quantity prior C _proxy_m can also be constructed as follows. Accordingly, the total quantity prior is a predicted value generated by a dedicated total quantity prior regression sub-model. The regression sub-model can optionally be a machine learning model, such as a gradient boosting machine (GBDT), a random forest (RF), or a small, lightweight multi-layer perceptron (MLP) network. The input feature vector X _proxyis a feature set specially constructed to circumvent the circular dependency, which does not contain the current time concentration to be solved.

[0227] Optionally, the feature set includes: all upstream water quality and hydrodynamic variables in the aligned downstream feature flow (for example, upstream metal concentration, upstream TSS concentration, upstream flow rate, etc. after travel time alignment); local environmental features (for example, local TSS concentration, local flow rate, local water level, local rainfall, etc.) of the current downstream section j at time t from the multi-source cleaning sequence data packet; corresponding human activity data (for example, local pollution record, upstream gate alignment signal, etc.).

[0228] On this basis, the training target of the regression sub-model is the true total concentration observation C _observed in the original monitoring data. The model can be jointly trained with the main parameters of the morphological allocation constraint head (MFCN_head) in the time series forward rolling training process, or two-stage training is adopted; that is, first train the regression sub-model independently, freeze it and use it as a feature extractor to train MFCN_head.

[0229] Further, the regression sub-model output ;

[0230] is used as input to calculate the physical upper bound of the particulate allocation ratio. A reasonable estimate of the total concentration is provided, which does not depend on the dissolved concentration to be solved at the current step, solving the circular dependency problem.

[0231] The technical scheme of the present application provides a method for cross-modal adsorption capacity estimation (VCAP) to solve the problems of lack of real-time chemical mechanism constraint and inability to dynamically infer the upper bound of total adsorption capacity from images. The method analyzes the image texture morphological feature vector, uses a constrained regression (such as a cross-modal mapping network) to infer the particle composition vector representing the chemical composition, and then combines the competition correction to calculate the dynamic adsorption capacity upper bound. The method solves the limitations of traditional methods that rely on static allocation coefficients.

[0232] It should be understood that the present scheme provides a hydrodynamic alignment strategy. The strategy calculates the time-varying travel time according to the real-time flow rate and gate state, introduces a diffusion kernel rollback mechanism to deal with non-steady flow. A physically reachable time window and a causal mask are constructed to make the upstream and downstream signals physically consistent, solving the problem of causal misplacement.

[0233] On this basis, in view of the consequence that the predicted morphology is unreliable and easy to violate the physical upper bound of adsorption saturation, the scheme adopts a morphology allocation constraint head (MFCN_head). The head applies the adsorption capacity upper bound calculated in the foregoing steps as a physical constraint when decoupling the calculation of the three-state concentration, introduces a saturation suppression function, makes the prediction results comply with the physical saturation upper bound, and improves the prediction reliability in extreme scenarios.

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

[0235] Optionally, read the image and remote sensing data of the multi-source cleaning sequence data packet _ particulate matter image / texture and the corresponding total suspended matter, perform denoising (median or bilateral filtering), flat field correction, and scale standardization (unify the spatial resolution to the reference pixel size) on the image, perform particle region segmentation (threshold self-adaptation + morphological opening and closing operation), extract multi-scale texture and morphological features, form an image texture morphological feature vector, and use it to regress the particle composition vector. The features include: gray level co-occurrence matrix statistics (contrast, homogeneity, energy, entropy, etc.), local binary pattern histogram (rotation invariant / multiple radius), edge direction histogram, particle size proxy (morphological particle size spectrum / equivalent diameter distribution), shape spectrum (circularity, aspect ratio, convexity), fractal dimension, etc.

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

[0237] Optionally, read the particle composition vector and the metal fingerprint parameters (metal fingerprint parameters) of each metal (valence, complexing ability, distribution coefficient, settling tendency, etc.), retrieve or regress the isotherm parameters of each metal m on the given particle type from the knowledge base, and obtain the effective parameter set by weighting the composition. For the Langmuir model, retrieve q _max_m (mg / g) and b _m (L / mg); for the Freundlich model, retrieve K _m and n _m . For competitive adsorption of multiple metals coexisting, introduce a competition correction coefficient γ _m，m( indicates the competitive strength of metal m' over metal m, varying with ion strength and pH), constructing effective affinity coefficient b _eff_m = b _m _m＇≠m _m，m＇ _m＇ ; where a _m＇ is the activity of metal m' or the relative effective concentration derived from per-metal metal fingerprint parameter and particle composition vector. The upper bound for surface site saturation is constrained by total site density N _s (given by particle type library, unit mmol / g) for q _max_m not exceeding N _s × M _m (M _m is the molar mass of metal m). The output is isotherm parameter set_per metal (q _max_m , b _eff_m or K _m , n _m ) and competitive correction parameter set, for capacity and upper bound calculation, morphological distribution constraints.

[0238] Optionally, read isotherm parameter set_per metal and total suspended solids, calculate the theoretical capacity upper bound of unit volume of water body to metal m (unit mg / L): Q _max_m = q _max_m × TSS. Define capacity response function q _m (C _d_m ), which maps dissolved concentration C _d_m to unit mass particle adsorption amount.

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

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

[0241] Particle phase concentration upper bound C _p_m_cap = min(q _m (C _d_m ) × TSS, Q _max_m ​​​) is given by the training period. If there exists the same-epoch raw monitoring data_ metal concentration sequence, the assigned prior estimate is simplified to calculate the provisional saturation S _m_obs = C _p_m _max_m (not as a hard constraint). The output contains Q _max_m , b _eff_m , or K _m , n _m , the capacity response function, and the optional S _m_obs , the capacity feature package, which can be used as an upstream signal for alignment, and as a core input for morphology assignment and upper bound constraint.

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

[0243] Correspondingly, read the raw hydrodynamic data_ flow rate and water level and human activity data_ gate status in the multi-source cleaning sequence data package, establish a directed graph according to the latitude and longitude of the section and the connectivity relationship of the river section, and according to the water level difference and the gate opening and closing state, determine the directionality and continuity for each time t; when the gate is closed or the reverse pressure occurs, the corresponding directed edge is set to zero. Save the length Δx and the current effective flow rate u _t for each directed edge. Output the directed adjacency and reachability matrix, which can be used to calculate the travel time and also to construct the causal mask.

[0244] Correspondingly, read the directed adjacency and reachability matrix and the raw hydrodynamic data_ flow rate and water level, and perform travel time integration on all reachable paths from the upstream section i to the downstream section j according to the piecewise constant speed:

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

[0246] When there are multiple reachable paths, select the minimum travel time; set a lower limit for the extremely low flow rate section to avoid numerical divergence, and smooth the mutation. The travel time of the unreachable pair is recorded as infinity or marked as missing. The travel time can be used to construct a physical window, and also to align the weight and freshness decay.

[0247] Correspondingly, read the travel time, and construct the physical reachable time window [t-T _ij，t -δ, t-T _ij，t +δ] for the downstream section j at time t. δ is used to cover the hydrodynamic uncertainty, and δ = κ × IQR(T _ij ​or by flow rate variance setting; for un-reachable or missing travel time cases, set the window to empty. Generate a causal mask (binary mask) M(i,j,t,t') for each upstream time index t', 1 within the window, 0 outside the window. Used for feature alignment and weighting under physical constraints.

[0248] Correspondingly, read the time-sequenced capacity feature package, upstream explanatory variables (upstream metal concentration, pollution record, rainfall, etc.) in the multi-source cleaning sequence data package, and the causal mask and travel time, and perform alignment and aggregation within the physically reachable time window. Apply freshness decay weight w _fresh =exp(-λ×Δ) to samples with upstream time offset Δ=t-t' _quality imp 2 , where σ imp 2 is the interpolation variance; integrate the weight w=M×w _fresh ×w _quality . For each feature component 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'); obtain the aligned downstream feature flow (including upstream concentration prior, aligned capacity feature, aligned pollution / weather explanatory variables). Output the aligned downstream feature flow, which can be packaged as a travel time alignment feature package.

[0249] The above describes optional embodiments of the present application in detail, but the present application is not limited to the specific details in the above embodiments, and various equivalent transformations can be made to the technical solutions of the present application within the technical concept of the present application, and these equivalent transformations all belong to the protection scope of the present application.​

Claims

1. A method for simulating and predicting heavy metal concentrations in water, characterized in that, include: 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; 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. 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. 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.

2. The method according to claim 1, characterized in that, 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, including: Analyze the image remote sensing data in the multi-source cleaning sequence data package to obtain the image texture morphology feature vector; By applying a constrained regression, the particle composition vector characterizing the chemical composition is inferred; Based on the particle composition vector and metal fingerprint parameters, the adsorption isotherm parameters are retrieved from the knowledge base and corrected. By combining the total suspended matter and adsorption isotherm parameters in the multi-source cleaning sequence data package, the adsorption capacity is calculated and a capacity feature package is generated.

3. The method according to claim 2, characterized in that, By applying constrained regression to the image texture morphology feature vector, the particle composition vector representing the chemical composition is inferred, including: A regressor is applied to the image texture morphology feature vector to obtain the original composition vector; Perform a simplex projection on it to make the particle composition vector satisfy the physical constraints of non-negativity and summation to one.

4. The method according to claim 2, characterized in that, Based on the particle composition vector and metal fingerprint parameters, adsorption isotherm parameters are retrieved from the knowledge base and corrected, including: Based on the particle composition vector, retrieve the baseline isotherm parameters from the knowledge base; Introducing a competing correction coefficient γ based on metal fingerprint parameters _m，m＇ By subtracting the reference isotherm parameters, adsorption isotherm parameters that take into account the multi-metal competition effect are obtained. The competition correction factor is then used to reduce the affinity factor, resulting in the effective affinity factor b. _eff_m ; in, ; b _eff_m =b _m / (1+Σ _m＇≠m c _m，m＇ ×a _m＇ ); Among them, b _m b is the reference affinity coefficient for the target metal m; _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 competing metal m'; a _ref For reference activity.

5. The method according to claim 2, characterized in that, Parse the image remote sensing data in the multi-source cleaning sequence data package to obtain image texture and morphological feature vectors, including: Denoising, flattening correction, and scale normalization are performed on remote sensing image data. Multi-scale texture and morphological features are extracted to form an image texture and morphological feature vector. The features include: gray-level co-occurrence matrix statistics, local binary pattern histogram, morphological granularity spectrum, and fractal dimension.

6. The method according to claim 1, characterized in that, 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 partial concentration and total concentration prediction result packets, including: An initial morphology allocation ratio is established based on the feature packet aligned with metal fingerprint parameters and travel time. By applying the upper bound of adsorption capacity contained in the capacity feature package and the alignment signal in the travel time alignment feature package, the initial morphology allocation ratio is constrained and shaped to obtain a capacity-constrained morphology allocation ratio. 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 packet of predicted concentrations of different states and total concentrations.

7. The method according to claim 6, characterized in that, Constraint shaping of the initial shape allocation ratio includes: By aggregating upstream contributions and local feature quantities from the travel time-aligned feature package, a total prior C is constructed. _proxy_m ; Using total prior and the upper bound Q of adsorption capacity in the capacity feature bag _max_m Calculate the physical upper bound π of the particle state distribution ratio. _p_m_cap To obtain the morphological allocation ratio under capacity constraints ; in, ; ； Where min() is the objective function to be minimized, ε is the numerical stability constant, and π _p_m_init The initial ratio is obtained from the aforementioned steps, η and ρ are learnable parameters used to control the suppression strength and nonlinearity, and S is the capacity saturation.

8. The method according to claim 6, characterized in that, The concentrations of dissolved, colloidal, and particulate states are decoupled and calculated, and then synthesized into a package of partial-state concentration and total concentration prediction results, including: Unconstrained estimates of total concentration are obtained by aligning feature packets with travel time. Based on unconstrained estimation of total concentration C _m_tilde Based on the capacity response function and upper bound of adsorption capacity in the capacity characteristic package, as well as the colloidal state proportion coefficient determined by the metal fingerprint parameters, a mass conservation equation for the concentration of dissolved state is established. Solve the mass conservation equation to obtain the dissolved concentration C. _d_m ; The concentrations C of colloidal and particulate matter were calculated based on the solution of the dissolved state concentration. _p_m The results of the synthesis fractional concentration and total concentration prediction package; Among them, 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 Let Q be the capacity response function, TSS be the total suspended solids concentration, and Q be the total suspended solids concentration. _max_m The upper bound of the particulate capacity and k _c_m C is the colloid proportion coefficient. _c_m Corresponding colloidal concentration.

9. The method according to claim 1, characterized in that, Based on multi-source cleaned sequence data packets, travel time is calculated, and causal alignment is performed on the capacity feature packet and upstream signals to generate a travel time aligned feature packet, including: By using hydrodynamic data and human activity data from multi-source cleaning sequence data packets, the accessibility of cross sections can be determined and travel time can be calculated; Construct a physically reachable time window based on travel time and generate a causal mask; By applying causal masking, a weighted aggregation including freshness decay and data quality is applied to the capacity feature packet and the upstream signal from the multi-source clean sequence data packet to generate a travel time aligned feature packet. Weighted aggregation uses a comprehensive weight w(i,j,t,t′). w(i, j, t, t') = M(i, j, t, t') × w _fresh (Δ) × w _quality (σ imp 2 ); where w _quality = 1 / (1 + σ imp 2 ); w _fresh = exp(-λ × Δ); Δ = |(t - t') - T _ij，t |; Where M(i, j, t, t') is the causal mask (0 or 1), w _fresh (Δ) represents the freshness decay weight, i and j are the index numbers of the upstream and downstream sections, t is the target time for predicting or constructing features for the downstream section j, t' is the upstream time, and T _ij,t Let λ be the time-varying travel time, λ be the freshness decay parameter, Δ be the signal hysteresis, and σ be the signal delay. imp 2 To interpolate the variance, w _quality (σ imp 2 ) represents the data quality weight.

10. The method according to claim 9, characterized in that, Using hydrodynamic data and human activity data from multi-source cleaning sequence data packets, the accessibility of cross sections is determined and travel time is calculated, including: The convection travel time is calculated based on the convection patterns in the hydrodynamic data. Monitor situations where convective travel time exceeds a preset physical limit due to extremely low flow rates; When the above conditions are triggered, the diffusion kernel backoff mechanism is activated, the diffusion backoff travel time is estimated based on the effective diffusion coefficient, and the diffusion backoff travel time is used as the travel time.

Images