Native sewage water quality prediction method and system based on physical constraint dynamic graph state space network

Abstract

The application belongs to the technical field of intelligent water affairs, and discloses a raw sewage water quality prediction method and system based on a physical constraint dynamic graph state space network, which comprises the following steps: constructing a heterogeneous node classification perception dynamic raw sewage pipe network graph, modeling a specific learnable interaction matrix for dynamic water force correlation by using a node type; embedding a solid-liquid two-phase coupling convection-diffusion equation into a neural network as a physical prior constraint; adopting a double-branch mixed pre-training framework for comparison between noise resistance generation and event perception; using an event-triggered selective state space model to capture long-range time sequence dependence and sudden pollution sensitive characteristics; and outputting water quality prediction values and uncertainty intervals, so that the problems of physical consistency loss, graph topology static assumption limitation and low long-range dependence modeling efficiency of existing methods are solved, high-precision real-time prediction of raw sewage water quality is realized, and real-time early warning requirements of a heat pump heating system are met.

Description

Technical Field

[0001] This invention belongs to the field of intelligent water management and environmental intelligent computing technology, specifically involving a spatiotemporal prediction method for native wastewater quality that deeply integrates prior knowledge of hydraulic physics, dynamic graph topology evolution learning, self-supervised comparison-generation hybrid pre-training, and state space sequence modeling. Background Technology

[0002] The evolution of water quality in urban primary sewage networks is a complex spatiotemporal process influenced by hydraulic laws, biochemical reaction kinetics, and human activities. Traditional methods face a triple dilemma: First, computational cost and accuracy are mutually exclusive. While physical driving models based on computational fluid dynamics (CFD) and SWMM are interpretable, their computation time can take hours to days, which cannot meet the needs of real-time early warning.

[0003] Second, data-driven models lack physical consistency. Existing deep learning methods (LSTM, GRU, Transformer) treat the native sewage pipe network as a purely data-driven "black box," ignoring the physical conservation laws of water mass transport, leading to prediction failures in off-distribution scenarios (such as sudden rainstorms and industrial leaks).

[0004] Third, the scarcity of labeled data restricts the generalization ability of the model. Raw wastewater quality labels (such as pollution source type and cause of concentration peak) need to be labeled by experts, and the scarcity of samples leads to the model's weak generalization ability in cross-pipeline and cross-seasonal scenarios.

[0005] In recent years, spatiotemporal graphical neural networks (ST-GNN) have become the mainstream framework for spatiotemporal prediction, but the following scientific problems remain unresolved: 1. Limitations of the static assumption of graph topology: Existing ST-GNN assumes that the adjacency matrix is ​​fixed or simply interpolated over time, which cannot capture the dynamic topology reconstruction caused by "convergence-diversion-backflow" in the original sewage pipe network.

[0006] 2. Limitations of self-supervised learning paradigms: Existing methods only use a single paradigm of mask reconstruction (MAE) or contrastive learning (CL), ignoring the complementarity between the two - MAE is good at capturing the smoothness of local features, while CL is good at capturing the consistency of global structure.

[0007] 3. The contradiction between long-range dependency and computational efficiency: The self-attention mechanism of the Transformer architecture has a computational complexity of O(N²) on long sequences, which is not suitable for real-time prediction of large-scale native sewage pipe networks. Summary of the Invention

[0008] The present invention proposes a method, device, and storage medium for predicting the quality of raw wastewater based on a physical constraint dynamic graph state space network, which can at least solve one of the technical problems in the background art.

[0009] To achieve the above objectives, the present invention adopts the following technical solution: a method for predicting the quality of raw wastewater based on a physically constrained dynamic graph state-space network, characterized in that the following steps are performed by a computer device: S1: Collect time-series data from multi-source heterogeneous sensors in the urban raw sewage pipe network. The sensor time-series data includes water quality parameters and hydraulic parameters. Perform missing value filling, outlier detection, and feature standardization on the collected data. Label the pipe network nodes according to the sensor installation location. The node types include: pumping station, inspection well, confluence point, and water intake. S2. Based on the node attributes and topology obtained in step S1, construct a heterogeneous graph of the original sewage pipe network; then adopt a classification-aware time-varying adjacency matrix learning mechanism, use the learnable node types to embed specific interaction weight matrices and nodes, combine Gaussian kernel distance weights, calculate the elements of the time-varying adjacency matrix, and use an adaptive update mechanism for time-series smoothing. S3: The solid-liquid two-phase coupled convection-diffusion equation is embedded as a physical prior constraint into the neural network training process; S4: Domain-adaptive dual-branch hybrid pre-training; A collaborative pre-training framework combining a noise-resistant generation branch and an event-aware comparison branch is adopted. S5: Event-triggered selective state-space modeling; Event-triggered selective state-space model is used to capture long-term temporal dependencies and sensitivity to sudden pollution. When the rate of change of flow exceeds the threshold, the event enhancement path is activated. A multi-scale time modeling strategy is adopted to extract features of different time scales with different expansion coefficients. After adaptive fusion with learnable weights, the features are integrated with the event-triggered state through a gating mechanism to obtain the final hidden state vector.

[0010] S6: Output the prediction results; Output the predicted values ​​of water quality parameters and the uncertainty range for future periods; when the predicted value exceeds the preset threshold or the uncertainty range is too large, an early warning signal is triggered.

[0011] A native wastewater quality prediction system based on a physically constrained dynamic graph state-space network includes: Data acquisition module: used to collect time-series data from multi-source heterogeneous sensors in the urban raw sewage pipe network, and to preprocess the data and label the node types; Dynamic graph construction module: used to construct a dynamic graph for classifying heterogeneous nodes. It uses 16 type-specific interaction matrices to model the hydraulic relationships between different node type pairs, and realizes adaptive learning of time-varying adjacency matrices. Physical constraint encoding module: used to transform the solid-liquid two-phase coupled convection-diffusion equation into a differentiable physical constraint loss function, constraining the physical consistency of the neural network output; Hybrid pre-training module: includes a noise-resistant generation branch and an event-aware contrastive branch, used to perform self-supervised pre-training on unlabeled data to learn noise-robust and structurally consistent node representations; State-space modeling module: used to capture long-range temporal dependencies of water quality parameters through event-triggered selective state-space models and multi-scale time fusion mechanisms; Prediction output module: Used to output predicted values ​​and uncertainty ranges of water quality parameters, and generate early warning signals based on the prediction results.

[0012] In another aspect, the present invention also discloses a computer-readable storage medium storing a computer program, which, when executed by a processor, causes the processor to perform the steps of the method described above.

[0013] In another aspect, the present invention also discloses a computer device, including a memory and a processor, wherein the memory stores a computer program, and when the computer program is executed by the processor, the processor performs the steps of the method described above.

[0014] As can be seen from the above technical solution, the original wastewater quality prediction method based on a physically constrained dynamic graph state-space network of the present invention aims to construct a heterogeneous node classification and perception dynamic original wastewater network graph, and to model dynamic hydraulic correlations using node types for specific learnable interaction matrices; it embeds the solid-liquid two-phase coupled convection-diffusion equation as a physical prior constraint into the neural network; it adopts a dual-branch hybrid pre-training framework of noise-resistant generation and event perception comparison; it utilizes an event-triggered selective state-space model to capture long-range temporal dependencies and sudden pollution sensitivity characteristics; and it outputs water quality prediction values ​​and uncertainty intervals. This solves the problems of lack of physical consistency, limitations of static assumptions in graph topology, and low efficiency in long-range dependency modeling in existing methods, and achieves high-precision real-time prediction of original wastewater quality. The prediction MAE is reduced by 52.8% compared with the existing SOTA methods, and the inference time is only 38ms, meeting the real-time early warning requirements of heat pump heating systems. Attached Figure Description

[0015] Figure 1 This is a diagram of the overall spatial network architecture of the present invention; Figure 2 This is a schematic diagram of the dual-branch comparison-generation hybrid pre-training framework of the present invention; Figure 3 This is a comparison chart showing the implementation effect of the present invention for predicting water quality impact during sports events; Figure 4 This is an architecture diagram of the event-triggered selective state-space model of the present invention; Figure 5 This is an architecture diagram of the multi-scale time modeling of the present invention; Figure 6 This is a detailed architecture diagram of the physical constraint encoding of the present invention; Figure 7 This is a diagram showing the overall deployment architecture of the system of this invention; Figure 8 This is a flowchart of the method of the present invention; Figure 9 This is a structural block diagram of the system of the present invention. Detailed Implementation

[0016] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are some embodiments of the present invention, but not all embodiments.

[0017] like Figure 1-9 As shown, the original wastewater quality prediction method based on a physically constrained dynamic graph state-space network described in this embodiment performs the following steps using a computer device: S1: Data acquisition and preprocessing; Time-series data from multi-source heterogeneous sensors in urban raw sewage pipe networks are collected. This sensor time-series data includes water quality parameters (suspended solids concentration SS, chemical oxygen demand COD, ammonia nitrogen NH3-N, pH value, and temperature) and hydraulic parameters (flow rate, velocity, and liquid level). Missing value imputation, outlier detection, and feature standardization are performed on the collected data. Pipeline nodes are labeled according to their sensor installation locations, with node types including: pump stations, manholes, confluence points, and water intakes. S2. Construct a dynamic graph for classifying heterogeneous nodes; S21. Based on the node attributes and topology obtained in step S1, construct a heterogeneous graph of the original sewage pipe network. S22. A classification-aware time-varying adjacency matrix learning mechanism is adopted. The learnable node types are used to embed specific interaction weight matrices and nodes. Combined with Gaussian kernel distance weights, the elements of the time-varying adjacency matrix are calculated, and an adaptive update mechanism is used for temporal smoothing. S3: Physically Constrained Reinforcement Learning; The solid-liquid two-phase coupled convection-diffusion equation is embedded as a physical prior constraint into the neural network training process. S4: Domain-adaptive dual-branch hybrid pre-training; A collaborative pre-training framework combining a noise-resistant generation branch and an event-aware comparison branch is adopted. S5: Event-triggered selective state-space modeling; Event-triggered selective state-space model is used to capture long-term temporal dependencies and sensitivity to sudden pollution. When the rate of change of flow exceeds the threshold, the event enhancement path is activated. A multi-scale time modeling strategy is adopted to extract features of different time scales with different expansion coefficients. After adaptive fusion with learnable weights, the features are integrated with the event-triggered state through a gating mechanism to obtain the final hidden state vector.

[0018] S6: Output the prediction results; Output the predicted values ​​of water quality parameters and the uncertainty range for future periods; when the predicted value exceeds the preset threshold or the uncertainty range is too large, an early warning signal is triggered.

[0019] The specific implementation steps are as follows: Heterogeneous node classification perception dynamic native sewage pipe network map construction Node definition and type encoding Pumping stations, inspection wells, confluence points, and energy station intakes are designated as nodes on the diagram. Each node is associated with a type label.

[0020] Sensor installation location as a node attribute Classification-aware dynamic adjacency matrix learning A learnable time-varying adjacency matrix is ​​used to model the dynamic hydraulic connections between nodes using the following formula:

[0021] In the formula: For a moment node With nodes The time-varying adjacency matrix elements between them take values ​​in the range [0,1]. The Sigmoid activation function maps the output to the [0,1] interval; For a moment node The first set of time-varying node embedding vectors has a dimension of ; superscript This is the matrix transpose operator; For node type pairs The corresponding learnable type interaction weight matrix has the following dimensions: ,common Each matrix corresponds to a hydraulic interaction mode between four node types: pumping station, inspection well, confluence point, and water intake. For a moment node The second set of time-varying node embedding vectors has a dimension of ; This is the scalar multiplication operator; For exponential function operators; For nodes With nodes Geographical distance or pipe segment length between them, in meters; The Gaussian kernel width parameter controls the sensitivity of distance decay. The unit is meters, and the value range is [500, 2000] meters, with a preferred value of 1000 meters.

[0022] in This is a learnable type interaction matrix between node type pairs to capture the differences in hydraulic conduction between different node types.

[0023] Adjacency matrix update frequency: This invention adopts an adaptive update mechanism with a basic update cycle of 15 minutes. When the detected flow change rate exceeds the threshold (|ΔQ / Q|>0.3), an immediate update is triggered.

[0024] Graph signal construction: Node characteristics: flow rate, velocity, temperature, turbidity, COD, SS, pH Side characteristics: pipe segment length, slope, design flow rate Global features: meteorological data (rainfall), calendar features (holidays) Physical constraint coding Embedded solid-liquid two-phase coupled convection-diffusion equations: (e.g.) Figure 6 (As shown) The physical laws governing the transport of suspended solids in the original sewage pipe network are encoded into hard constraints of a neural network. This invention targets solid concentration. Due to the unique characteristics of the transmission of suspended solids concentration (SS) measured by the sensor, the finite difference method is used to calculate the partial derivative terms along the pipe section direction, and a deposition-resuspension source term is introduced. :

[0025] In the formula: This is an approximation of the partial derivative of the solid concentration with respect to time for model prediction, in mg / (L·h). The numerical approximation is equal; The time predicted by the model Solid concentration value, in mg / L; It is the subtraction operator; The time predicted by the model Solid concentration value, in mg / L; The time sampling interval is in hours (h), and in this invention, it is set to 0.25h (i.e., 15 minutes).

[0026] Convection term along pipe section direction Discretization:

[0027] In the formula: This is an approximation of the spatial divergence of the convection term, describing the transport flux of fluid carrying solid particles; The numerical approximation is equal; Along the direction of the pipe section The fluid velocity component, in m / s; This is the concentration of solids, in mg / L. The product of flow rate and concentration, i.e., solid mass flux, is expressed in mg·m / (L·s); subscript Indicates the location of the pipe section Shift downstream The value at the position; the subscript Indicates the location of the pipe section Shifting upstream The value at; The spatial discretization step size is the distance between two adjacent discrete points, expressed in meters (m).

[0028] Diffusion term along pipe section direction Discretization:

[0029] In the formula: This is an approximation of the spatial divergence of the diffusion term, describing the molecular diffusion and turbulent diffusion processes of solid particles; The numerical approximation is equal; The solid diffusion coefficient is expressed in m² / s and is adaptively adjusted based on the solid particle size distribution and the slope of the pipe section. Pipeline section location Solid concentration at the location, in mg / L; Pipeline section location The solid concentration value at the current discrete point, expressed in mg / L. Pipeline section location Solid concentration at the location, in mg / L; The step size for spatial discretization is in meters (m). It is the square of the spatial step size.

[0030] The sedimentation-resuspension source term is modeled as follows:

[0031] In the formula: This is a deposition-resuspension source term function, with output units of mg / (L·h), describing the net deposition or resuspension rate of solid particles at the bottom of the pipe; The independent variable of the function is the current solid concentration, expressed in mg / L. This is the deposition rate coefficient, in h.- ¹ represents the deposition rate constant of solid particles at the bottom of the pipe; This is the deposition term, which describes the deposition process of solid particles at the bottom of the pipe due to gravity, and is proportional to the concentration. This is a subtraction operator, representing the effect of resuspension on the offsetting effect of deposition; This is the resuspension rate coefficient, in h. - ¹ represents the rate constant for the resuspension of deposited particles into the water flow; The term represents the square of the solid concentration, reflecting the nonlinear enhancement effect of resuspension at high concentrations; The maximum solid concentration, expressed in mg / L, represents the upper limit of the concentration that may occur in the pipeline network.

[0032] in The deposition rate coefficient is... The resuspension rate coefficient, This represents the maximum solid concentration. The above parameters are calibrated using a joint strategy of physical prior initialization and gradient backpropagation. Physical prior initialization The physical meaning of is the deposition rate of solid particles at the bottom of the pipe, according to Einstein's deposition formula. With the slope of the pipe section A negative exponential relationship is observed—the deposition rate is low in steep slope sections under gravity. This invention uses an empirical formula for initialization: ; In the formula: This is the deposition rate coefficient, in h. -1 Output quantity to be determined; Reference deposition rate coefficient, in h. -1 In this invention, the value is taken as 0.15h. -1 ; This is the multiplication operator; The natural exponential function operator; The negative operator; The slope attenuation coefficient is dimensionless, and its value is taken as 3.5 in this invention. This is the sine function operator; The slope angle of the pipe section is expressed in rad (radians), and its value range is [value range missing]. .

[0033] in h -1 The sedimentation rate coefficient is a reference value (based on typical values ​​of municipal primary sewage pipe networks). The slope attenuation coefficient (calibrated using publicly available literature [A]). The slope of the pipe section (in radians) is the basis for this formula. The physical principle is that as the slope increases, the shear stress of the water flow increases, making it more difficult for solid particles to settle, thus reducing the sedimentation rate.

[0034] The physical meaning of is the resuspension rate of deposited particles. According to Shields deposition kinetics, when the flow velocity exceeds the critical deposition flow velocity... At this point, the deposited particles begin to resuspend. This invention uses the flow rate ratio as the driving force for resuspension: ; In the formula: This is the resuspension rate coefficient, in h. -1 Output quantity to be determined; For reference resuspension rate coefficient, the unit is h. -1 In this invention, the value is taken as 0.08 hh -1 ; This is the multiplication operator; The ratio of flow rates is dimensionless. The actual flow velocity in the current pipe section, in m / s, is obtained from real-time sensor measurement or hydraulic model calculation. The design flow velocity for the pipe section is expressed in m / s, which is the standard flow velocity specified in the design of the pipe network. The velocity sensitivity index is dimensionless, and its value is 2.0 in this invention.

[0035] in h - ¹ is the reference resuspension rate coefficient. It is a flow velocity sensitivity index (calibrated based on experimental data from native sewage pipelines). For the current flow rate, The design flow velocity for the pipe section is given. The physical basis of this formula is that the resuspension rate is positively correlated with the shear stress of the water flow, and the shear stress is approximately proportional to the square of the flow velocity.

[0036] The initial value is taken as the 99th percentile of historical data to ensure coverage of extreme high concentration events.

[0037] Gradient backpropagation calibration: During training, the three parameters are jointly optimized as learnable variables through gradient descent, while non-negativity constraints are applied. To ensure physical meaning, the optimizer uses Adam with a learning rate of 1e-4 (lower than the main network's learning rate of 1e-3) to ensure that the physical parameters do not lose their physical realism due to overtraining.

[0038] Parameter calibration range

[0039] The physical constraint loss function is defined as: ; In the formula: This represents the expected form of the physical constraint loss function; For time Velocity field and solid concentration The joint mathematical expectation operator represents the average loss over all samples; The time predicted by the model Solid concentration value, in mg / L; The time predicted by the model Solid concentration value, in mg / L; The time sampling interval is in hours (h), and in this invention it is set to 0.25h (i.e., 15 minutes). The fluid velocity vector field predicted or input to the model, in m / s; This represents the solid concentration field predicted by the model, in mg / L. This is the solid-state diffusion coefficient tensor, in m² / s. The sedimentation-resuspension source term function is used to predict concentrations using the model. Calculated as an independent variable; The formula in this case and the formula in claim 3 are essentially the same physical constraint, only expressed differently; The statement in claim 3 is a theoretical definition used to describe physical principles; This formula is a discretized version used in actual calculations, which is the discretized implementation of the physical constraint loss function in claim 3.

[0040] The diffusion coefficient tensor The adjustment is adaptive based on the solid particle size distribution and the slope of the pipe section.

[0041] Physically Constrained Gradient Backpropagation and Joint Parameter Optimization Mechanism Physical constraint loss The neural network is trained by embedding it through the following mechanism: 1. Differentiable Discretization: All finite difference operations in equations (14)-(16) are implemented using the PyTorch automatic differentiation framework, ensuring that the physical constraint loss is converted into the neural network prediction. The gradient chain is fully differentiable. Specifically, The parameters are automatically backpropagated to the GCN encoder and Mamba module through the computation graph.

[0042] 2. Physical parameter learning strategy: physical parameters A joint calibration strategy of "physical prior initialization + gradient backpropagation fine-tuning" is adopted: Pre-training phase: Physical parameters are initialized with prior physical values ​​( h -1 , h -1 , (mg / L), and freezes it to not participate in gradient updates, only the neural network weights. renew; Fine-tuning phase: Unfreeze physical parameters at a small learning rate ( Gradient backpropagation fine-tuning is performed using 1 / 10 of the main network's learning rate, while applying non-negativity constraints ( Ensure physical feasibility.

[0043] 3. Multi-task gradient balancing: The GradNorm algorithm is used to dynamically adjust the weights of the physical constraint loss. This ensures that the gradient norms of each loss term remain on the same order of magnitude, preventing physical constraints from dominating or disappearing.

[0044] Error analysis and convergence proof of physical constraints 1. Error bounds of the finite difference approximation: The time partial derivatives are obtained using the central difference method. The truncation error is

[0045] The space convection term uses central difference, with a truncation error of .

[0046] The diffusion term uses a second-order central difference, with a truncation error of...

[0047] when minute, When m, the total approximation error is approximately The magnitude meets the requirements for real-time prediction accuracy. 2. Gradient stability analysis of physical constraint loss: Due to the use of nonnegativity constraints ( During gradient backpropagation, the gradient of the physical parameters will not cause the parameters to jump out of the physical feasible region. A separation strategy is adopted between the physical parameter learning rate (1e-4) and the main network learning rate (1e-3) to ensure that the convergence speed of the physical parameters matches that of the main network. It has been verified that the physical parameters tend to converge after 50 epochs of training, and the parameter values ​​after convergence are still within the range of the physical prior. 3. Convergence Theorem: Assume the parameters of the neural network... and physical parameters The joint optimization objective is In the formula: To jointly optimize the total loss function; The prediction loss function measures the difference between the model's predicted values ​​and the actual values. It is the addition operator; This is the weighting coefficient for the convection-diffusion physical constraint loss, with a value of 0.5. The convection-diffusion physical constraint loss function; The mass conservation constraint loss weighting coefficient is set to 0.3. Let be the loss function for the mass conservation constraint. If: (a) about and It is a Lipschitz continuum (b) The learning rate satisfies (L is the Lipschitz constant) (c) Physical parameter constraint set It is a convex set The stochastic gradient descent algorithm guarantees convergence to a stable critical point.

[0048] 4. Advantages compared to standard PINN:

[0049] mass conservation constraint Introducing quality conservation loss at the node level, and defining cumulative traffic as a time window. Internal integral traffic: In the formula: For nodes At any moment Cumulative inflow rate, in m³; For definite integral operators; lower limit of integration The start time of the time window; the upper limit of integration. The current moment; For nodes At any moment The inflow rate function, in m³ / h, represents the flow rate from the upstream pipe section into the node. The relationship between flow rate and time; Let be the integral variable, representing any moment within the time window; It is a time infinitesimal element.

[0050] The mass conservation loss is defined as:

[0051] In the formula: The loss function is defined by the mass conservation constraint. It is the set of all nodes in the original sewage pipe network; Use node index variables to iterate through the collection. All nodes in; This is the summation operator; For nodes The set of neighboring nodes, that is, the set of nodes with respect to the node All directly connected neighboring nodes; Use the neighbor node index variable to traverse All nodes in; For a moment From node Flow to Node (or from node) Flow to Node The flow rate of the pipe section, in m³ / h; For nodes In the time window The cumulative inflow rate, in m³, is calculated using formula (20).

[0052] in For pipe section At any moment Traffic, For nodes In the time window The cumulative inflow within.

[0053] Native wastewater domain adaptive dual-branch hybrid pre-training (BCGP) Noise-Robust MAE (Noise-Robust MAE): (e.g.) Figure 2 (As shown) This invention employs a 3-layer stacked Graph Convolutional Network (GCN) as a shared encoder. The propagation formula for layer graph convolution is:

[0054] In the formula: For the first The node feature matrix of the layer has a dimension of ; The ReLU activation function is used to perform element-wise operations on the input. ; The degree matrix is ​​the negative 1 / 2 power, used for symmetric normalization, with dimension 1. ; The time-varying adjacency matrix for adding self-loops has dimensions of [dimensional value missing]. ,in It is the identity matrix; For the first The node feature matrix of the layer, The input feature matrix has a dimension of . ; For the first The learnable weight matrix of the layer has a dimension of The hidden dimensions of the three-layer GCN are 128→128→256 respectively, and the final output node embedding... .

[0055] 1. Automatic identification of abnormal sensor regions: Detects node-time locations with continuous missing values ​​(>3 time steps) or abnormal jumps (>|Δx|>3σ). 2. Segmented masking strategy: Abnormal regions are preferentially masked at 30%, while normal regions are randomly masked at 15%. 3. Node representation after mask extraction by the graph convolutional encoder 4. MLP decoder reconstructs the masked original features 5. Reconstruction Loss:

[0056] In the formula: For noise-resistant mask autoencoder reconstruction loss function; The set of masked nodes The number of elements, i.e., the total number of masked nodes-time positions; This is the summation operator; The index variable is the position of the masked location; iterate through the collection. All elements in; For position The original feature vector at that node contains all sensor measurements of that node at the corresponding time. Location for model reconstruction The feature vector at that location is the predicted value output by the decoder; This is the L2 norm square operator, which calculates the squared Euclidean distance between two vectors.

[0057] Event-Aware Comparison Branch (Event-Aware CL) 1. Event Trigger Detection: When Time-triggered enhanced view generation 2. Extract global graph representations from the original graph and the augmented graph using a graph convolutional encoder.

[0058] 3. Graph augmentation strategies include: (a) Edge Dropout: using probability... (a) Randomly discard some edges; (b) Node feature Dropout: based on probability (c) Add noise to node features; (d) Randomly prune time slices: randomly select continuous subsequences. 4. Employ contrastive loss to maximize the similarity of source views:

[0059] In the formula: The loss function for event-aware contrastive learning is based on the InfoNCE loss form. The operator for negative natural logarithms; For exponential function operators; The similarity calculation function uses cosine similarity. ; The global graph representation vector extracted from the original graph by the graph convolutional encoder has dimensions of . ; To enhance the view (with) The global graph representation vector extracted by the graph convolutional encoder (from the same source) has a dimension of ; The temperature parameter controls the sharpness of the softmax distribution; in this invention, the value is set to 0.5. This refers to the batch size, i.e., the number of samples in a batch. In this invention, the value is 32. This is a sample index variable within the batch, with a value range of [value range missing]. ; The first in the batch The global representation vector of the enhanced view of each sample.

[0060] Among them temperature parameters Batch size .

[0061] 5. Event consistency loss: Ensuring consistency of representation across nodes at the time of event triggering:

[0062] In the formula: Let the event consistency loss function be used. Set of event moments The number of elements, i.e., the total number of time steps in which a sudden pollution event was detected; This is the summation operator; Use time index variables to iterate through the collection. All event moments; For set membership operators; This is a collection of event moments, containing indexes of all moments identified as sudden events by traffic mutation detection. For a moment The global graph representation vector, with dimension . Extracted by a graph convolutional encoder; For a moment The average graph representation vector of the neighboring nodes, i.e. ; This is the L2 norm square operator.

[0063] in For the collection of event moments, For a moment Average representation of neighboring nodes.

[0064] Structural Restructuring Branch Reconstruct the complete adjacency matrix to capture global topological information:

[0065] In the formula: Reconstruct the loss function for the adjacency matrix structure; The actual time-varying adjacency matrix has dimensions of . ,in The total number of nodes; The adjacency matrix for model prediction or reconstruction, with dimension . ; This is the matrix subtraction operator; This is the Frobenius norm operator, which calculates the square root of the sum of the squares of all elements in a matrix.

[0066] Total pre-training loss:

[0067] In the formula: The total loss function for BCGP pre-training is the weighted sum of four sub-loss functions; To mitigate the noise-resistant mask autoencoder reconstruction loss, the reconstruction accuracy of mask features is measured. To compare the loss weighting coefficient, a value of 1.0 was used; For event-aware contrastive learning loss; This is the event consistency loss weighting coefficient, with a value of 0.5; Let the event consistency loss function be used. The structural reconstruction loss weighting coefficient has a value of 0.3. The loss function is reconstructed for the adjacency matrix structure. The type-sharing regularization loss weight coefficient is set to 0.2, which controls the strength of common constraint on the interaction matrix of nodes of the same type. The type-sharing regularization loss function forces all interaction matrices of nodes of the same type to converge to a shared basis vector.

[0068] Among them, weight parameters , , , =0.2; This weighting parameter, selected based on ablation experiments and validated, is designed to control overfitting in the 16 types of interaction matrices. It ensures commonality in interaction patterns within the same type while allowing for differences between types by forcing all interaction matrices of nodes of the same type to converge towards a shared basis vector. The selection criteria and ablation validation for this weighting parameter table are as follows:

[0070] BCGP pre-training and downstream task connection mechanism After pre-training, the model parameters are transferred to the downstream water quality prediction task according to the following strategy: 1. Parameter Preservation: 3-layer graph convolution weights of the GCN encoder All are retained for initialization of downstream tasks; node embedding of the dynamic graph construction module. Interaction matrices of 16 types Keep everything.

[0071] 2. Parameter Reset: The MAE decoder (MLP reconstruction head) and contrastive learning projection head are discarded after pre-training and do not participate in downstream tasks.

[0072] 3. Initialization of new modules: The Mamba state space module (step S4), the multi-scale fusion module and the prediction output head (step S5) are trained from scratch using Xavier uniform initialization.

[0073] 4. Two-stage fine-tuning: - Phase 1 (Freeze Fine-tuning): Freeze the parameters of the GCN encoder and dynamic graph module, and train only the Mamba module and prediction head (about 20 epochs) to adapt the new module to the pre-trained representation; - Second stage (full parameter fine-tuning): Unfreeze all parameters with a small learning rate ( End-to-end fine-tuning is performed (approximately 80 epochs) while preserving the physical constraint loss.

[0074] 5. Physical constraint connection: Physical constraint loss during fine-tuning phase and Continue applying physical parameters In the second phase of thawing, participants will engage in fine-tuning.

[0075] Selective state-space sequence modeling triggered by raw wastewater events Time-series tokenization The time-series data of the original sewage pipe network system is divided into non-overlapping time slices (patch), and each patch contains 32 time steps.

[0076] Event-triggered selective state-space modeling: (e.g.) Figure 4 (As shown) By triggering Mamba modules through events, long-range temporal dependencies of water quality parameters can be captured, while maintaining sensitivity to sudden pollution events. 1. Event detection: Based on sudden changes in traffic. Trigger event indicator

[0077] 2. Event-triggered selective state update:

[0078] In the formula: For a moment The hidden state vector, with dimension . Encoding the cutoff time Temporal dependency information; For a moment A selective state-space gated scalar, with a value range of [0,1], is determined by the input. Adaptive generation; This is Hadamard's element-wise multiplication operator; The state transition weight matrix has dimensions of . This maps the input to the state space; For a moment The input feature vector contains water quality sensor data and environmental variables, with dimensions of [missing information]. ; This is the vector addition operator; For a moment The hidden state vector, i.e., the state at the previous time step; For complementary gating coefficients, and Complementary control state update strength; This is an event indicator vector, triggered based on traffic surge detection, with values... (No event) or (Event detected); is the hyperbolic tangent activation function, with an output range of [-1, 1]. Enhance the event weight matrix, with dimension 1. It is used to inject additional non-linear state transitions when an event occurs.

[0079] 3. Output Mapping:

[0080] In the formula: For a moment The predicted output vector contains water quality predictions for multiple future time steps, with dimensions of [missing information]. ; The Sigmoid activation function maps the output to the [0,1] interval; Let be the linear projection weight matrix from input to output, with dimension . ; For a moment The input feature vector has a dimension of It includes the sensor data and environmental variables at the current moment; This is the vector addition operator; Let be the linear projection weight matrix from the hidden state to the output, with dimension . ; For a moment The hidden state vector, with dimension . Encoding the cutoff time The timing dependency information.

[0081] in For event indicators, Enhance the weight matrix for events. This is the event detection threshold (usually set to 2.0-3.0).

[0082] Multi-scale time modeling (e.g.) Figure 5 (As shown) Using three Mamba modules with different dilation coefficients in parallel, multi-scale periodic patterns ranging from minute to day are captured: Fine-scale branch (dilation=1): Processes the sequence with the original sampling interval to capture rapid water quality fluctuations on a minute-by-minute (15-30 minute) scale, maintaining high sensitivity to sudden pollution events. Mesoscale branch (dilation=2): Sampling at twice the interval to capture hourly (1-3 hour) tidal water use cycle patterns. Coarse-scale branch (dilation=4): Sampling at 4-fold intervals to capture the diurnal water quality evolution pattern at the daily (24-hour) level. The outputs of the three branches are adaptively fused using learnable weights:

[0083] In the formula: The hidden state vector after multi-scale fusion has a dimension of . ; This is the summation operator; For scale index variables, the set of values ​​is: These correspond to three Mamba branches with expansion coefficients of 1, 2, and 4, respectively. For set membership operators; For the first The fusion weight scalars for each scale are obtained by softmax normalization, and their values ​​range from [0,1]. ; This is the vector-scalar multiplication operator; For the first Each scale (expansion coefficient is) The Mamba branch at time The output hidden state vector has dimension . ; For exponential function operators; For the first Learnable weight parameters at each scale (unnormalized logit values); To sum the index variable, iterate through all scales.

[0084] in These are learnable scale weight parameters. The representation after multi-scale fusion. The event triggering state in step S4.2 Integration through gating mechanisms:

[0085] In the formula: The final hidden state vector after gating fusion has a dimension of . ; The gating weight vector has dimension 1. The value range is [0,1], which controls the fusion ratio of the two inputs; This is Hadamard's element-wise multiplication operator; For the event-triggered state space model at time... The output hidden state vector has dimension . It is sensitive to sudden pollution incidents; This is the vector addition operator; For complementary gating weight vectors, and Element-wise complementary; The hidden state vector is a multi-scale fused vector with dimension . Encodes multi-scale periodic pattern information; The Sigmoid activation function maps the output to the [0,1] interval; The gated projection weight matrix has dimensions of . ; For vector concatenation operators, and spliced ​​as Dimensional vector.

[0086] This design enables the model to rely on multi-scale periodic patterns during stable periods, while automatically increasing the weight of event-triggered branches when events occur.

[0087] Predicted output and adaptive sampling Multi-step prediction head: Employs independent fully connected layers to output predicted water quality values ​​for the future T∈{15min, 30min, 1h, 3h}.

[0088] Quantification of uncertainty The prediction uncertainty is estimated by ensemble learning, and the prediction interval [μ-1.96σ, μ+1.96σ] is output.

[0089] Deep Ensemble Configuration (Detailed Technical Solution) 1. Integrated Architecture Design Number of members: An independent model (verified by ablation experiments), It can achieve the best balance between the accuracy of uncertainty estimation and computational cost; when The time variance estimation bias is +18.5%, when (Inference time increased by 95%, but accuracy only improved by 2.1%) Member model structure: Each member adopts a complete PC-DGSSN architecture (dynamic graph construction → physical constraint encoding → BCGP pre-training → event-triggered SSM → multi-scale fusion) to ensure the consistency of the integration. 2. Member diversity guarantee mechanism (key innovation) To avoid "ensemble collapse" (where all members converge to the same solution), This invention employs a three-layer diversity guarantee strategy:

[0090] 3. Predictive Fusion and Uncertainty Estimation Point prediction value: (Arithmetic mean), where: Predicted values ​​for integration points; The reciprocal of the number of members; This is the summation operator; For member index variables; The total number of members is set to 5. For the first Predictions from each member model Uncertainty estimation: (Unbiased sample variance, denominator is) Instead To eliminate deviation), where: To predict variance, measure the uncertainty of the forecast; The normalization coefficients for the unbiased variance estimate; This is the summation operator; For member index variables; The total number of members is set to 5. For the first The predicted values ​​of the individual member models; It is the subtraction operator; Predicted values ​​for integration points; Square operator Prediction Interval: Based on the normal distribution assumption, the 95% confidence interval is... In the formula: For range representation; Predicted values ​​for integration points; It is the subtraction operator; The two-tailed quantiles corresponding to a 95% confidence level of the standard normal distribution; To predict the standard deviation, i.e. The square root of; For addition operators Optional confidence level: Supports user-defined confidence levels The corresponding multiplier factors are respectively , ,

[0091] 4. Robust outlier handling mechanism Outlier detection rule: When a member predicts... satisfy Time (of which) (The original standard deviation including outliers) is used to determine if the prediction is an outlier. Recalculation strategy: After removing outliers, use the remaining values. Each member is recalculated. and

[0092] Minimum number of members constraint: If after removal If so, all members will be retained without being removed (to avoid information loss due to over-filtering). 5. Uncertainty calibration and verification experiment (three-city dataset)

[0093] Analysis Conclusion The coverage rate of Jinan City (local data) is 94.8%, which is highly consistent with the theoretical 95%, proving that the uncertainty quantification calibration is effective. Coverage decreased slightly during cross-domain testing (93.7% in Beijing), which is an inherent challenge of out-of-distribution generalization, but it remains within an acceptable range (>90%). The peak error (prediction error > 2 × MAE) coverage rate reached 87%-91%, demonstrating the model's ability to identify regions of high uncertainty. 6. Computational Cost Analysis Single inference time: 38ms (single model) → 167ms (5-member ensemble, including data parallel acceleration) Memory usage: 1.2GB (single model) → 5.8GB (5 members, optimized to 4.1GB after sharing the static graph structure) Meets the requirements for real-time early warning (<200ms / prediction). Adaptive sampling frequency control Dynamically adjust the sensor sampling frequency based on the prediction variance:

[0094] In the formula: The next sensor sampling interval, in minutes; The minimum allowable sampling interval, in minutes, is 5 minutes in this invention; It is the addition operator; The maximum allowable sampling interval, in minutes, is 30 minutes in this invention; It is the subtraction operator; This is the scalar multiplication operator; The Sigmoid activation function maps the output to the [0,1] interval; The sampling sensitivity coefficient controls the degree to which the prediction variance affects the sampling interval; in this invention, it is set to 0.5. For a moment The prediction variance is calculated from the predictions of the five member models of Deep Ensemble, in units of (mg / L)².

[0095] Algorithm Flowchart graph TD subgraph data input S1 [Sensor Data] --> S2 [Weather / Calendar Data] S2 --> S3 [Data Preprocessing] end Subgraph Step S1: Dynamic Graph Construction S3 --> S4 [Node Type Encoding] S4 --> S5 [Calculate the time-varying adjacency matrix A(t)] S5 --> S6 [Graph Signal Construction] end Subgraph Step S2: Physical Constraint Encoding S6 --> S7 [Calculate the loss in the convection-diffusion equation] S7 --> S8 [Calculate mass conservation loss] S8 --> S9 [Physical Parameter Update] end Subgraph Step S3: BCGP Pre-training S6 --> S10 [MAE generates branch] S6 --> S11 [Comparative Learning Branch] S10 --> S12 [Calculate the total BCGP loss] S11 -->S12 end Subgraph Step S4: Mamba Sequence Modeling S12 -->S13 [Event Detection] S13 -->S14 {Did the event trigger?} S14 --> | Yes | S15 [Event Enhancement Status Update] S14 --> |No| S16 [Standard Mamba State Update] S15 --> S17 [Multi-scale fusion] S16 -->S17 end Subgraph step S5: Predict the output S17 --> S18 [Multi-step prediction head] S18 --> S19 [Uncertainty Quantification] S19 --> S20 [Adaptive Sampling Control] end S9 -->S10 S9 -->S14 style S10 fill:#f96,stroke:#333 style S11 fill:#96f,stroke:#333 style S15 fill:#f9f,stroke:#333 style S19 fill:#9f9,stroke:#333 Time complexity analysis

[0096] Comparison of reasoning efficiency

[0097] Note: The test environment was NVIDIA A100 40GB, batch size = 32, input sequence length = 32, and number of nodes = 50.

[0098] Space complexity optimization strategies 1. Sparse Adjacency Matrix Storage: The adjacency matrix is ​​stored in CSR format, reducing space complexity from... Down to

[0099] 2. Low-rank decomposition: The 16 types of interaction matrices are decomposed using low-rank decomposition, reducing the number of parameters by approximately 8 times. 3. Mamba State Reuse: Reuse historical states during inference, eliminating the need to store the complete sequence of hidden states. Model Deployment and Real-Time Inference Solutions Deployment architecture (such as) Figure 7 (As shown)

[0100] Real-time inference pipeline Model Quantization and Acceleration

[0101] Edge-Cloud Collaborative Inference

[0102] Model update strategy 1. Incremental learning: Incremental training is performed every 72 hours using the latest data. Physical parameters are frozen, and only the neural network weights are updated. 2. Full Retraining: Perform full retraining every 6 months, including joint optimization of physical parameters. 3. A / B Testing: The new model is first validated on 10% of the nodes, and then fully deployed after performance improvements are achieved. Model update strategy effectiveness verification To verify the effectiveness of the model update strategy, a long-term tracking experiment was designed. The experiment period was from January 2023 to December 2023 (12 months in total). The dataset was divided into three parts according to time: training set (January 2022 to June 2023), validation set (July 2023 to September 2023), and test set (October 2023 to December 2023).

[0103]

[0104] Experimental results show that the incremental + full update strategy can achieve the best prediction accuracy (MAE reduced by 25.6%) with limited computational overhead (+12%), and A / B testing can ensure that the new model is stable before full deployment, reducing update risk.

[0105] To make the objectives, technical solutions, and advantages of this invention clearer, the invention will be further described in detail below with reference to the accompanying drawings and specific embodiments. This chapter is organized in the logical order of "general experimental configuration → application scenario embodiments → large-scale experimental verification". First, it introduces the general settings such as system deployment architecture, data preprocessing, model training strategies, and hyperparameter configuration. Then, it demonstrates the actual effects of this invention through four typical application scenarios. Finally, it conducts systematic experimental verification on the USN-MSHD large-scale dataset.

[0106] General Experiment Configuration System Deployment Architecture This invention adopts an edge-cloud collaborative deployment architecture (such as...) Figure 7 (As shown). Edge devices are deployed at each pump station and key inspection wells, responsible for data preprocessing and real-time updates of the adjacency matrix; cloud devices are deployed in the data center, responsible for model inference and early warning notifications. The edge devices and the cloud communicate via the MQTT protocol to ensure the real-time performance and reliability of data transmission.

[0107] Standardized Data Preprocessing Process 1. Missing value imputation: Single missing values ​​are imputed using linear interpolation (<3 consecutive steps); for long-term missing values ​​(≥3 steps), the weighted average of the previous and next 3 steps is used for imputation. 2. Outlier Handling: Outlier measurements are identified based on the 3σ criterion, marked as missing values, and then refilled. 3. Feature Standardization: Z-score standardization is used, with parameters calculated based on the training set. In the formula: These are the standardized eigenvalues, which are dimensionless. These are the original eigenvalues; It is the subtraction operator; The mean of the features in the training set; This is the division operator; Standard deviation of training set features Model training strategy and hyperparameter configuration Training strategy 1. Pre-training phase: Perform BCGP pre-training on unlabeled historical data to learn noise-robust and structurally consistent node representations. 2. Fine-tuning phase: Supervised fine-tuning is performed on the labeled data, while retaining the physical constraint loss term. 3. Inference Phase: Exponential Moving Average (EMA) is used to replicate model weights, improving inference stability. Hyperparameter configuration

[0108] Evaluation Indicator Definition The following metrics are used to comprehensively evaluate model performance:

[0109] In the formula: The unit is the mean absolute error, and the unit is the same as the target value. The coefficient is the reciprocal of the sample size, which is used to calculate the average. The total number of test samples; This is the summation operator; This is a sample index variable, with a value range of [value range missing]. ; This is the absolute value operator; For the first Predicted values ​​for each sample; It is the subtraction operator; For the first The true value of each sample.

[0110]

[0111] In the formula: The root mean square error is expressed in the same unit as the target prediction. This is the square root operator; It is the reciprocal of the sample size; The total number of test samples; This is the summation operator; For sample index variables; For the first Predicted values ​​for each sample; It is the subtraction operator; For the first The true value of each sample; This is the square operator.

[0112]

[0113] In the formula: Mean absolute percentage error, in % This is the percentage conversion factor; It is the reciprocal of the sample size; The total number of test samples; This is the summation operator; For sample index variables; This is the absolute value operator; For the first Predicted values ​​for each sample; It is the subtraction operator; For the first The true value of each sample; This is the division operator; the denominator is... For the first The true value of each sample is used as the normalization benchmark.

[0114] In the formula: It is a continuous sort probability score that measures the overall consistency between the predicted distribution and the observed values, and its unit is the same as the square of the predicted target. This is the operator for definite integrals from negative infinity to positive infinity; For the first The cumulative distribution function predicted for each sample, i.e. ; It is the subtraction operator; This is an indicator function that takes the value 1 when the condition inside the square brackets is true, and 0 otherwise. For comparison operators, evaluate variables Is it greater than or equal to the actual value? ; For square operators; Let be the integral variable, representing the value of the predicted value; It is an integral infinitesimal element; For the first Predicted values ​​for each sample; For the first The true value of each sample; This represents the total number of test samples.

[0115] Application Scenario Examples The following four embodiments, based on measured data from the raw sewage pipe network in Jinan City, Shandong Province, verify the effectiveness of the present invention in different application scenarios. Jinan City utilizes three main raw sewage pipe networks to provide heat sources for six heat pump energy stations. The pipe network includes 47 key nodes (8 pump stations, 19 inspection wells, 12 confluence points, and 8 water intakes). The data time range is from January 2022 to December 2023, with a sampling interval of 15 minutes.

[0116] Example 1: Prediction of Regional-Level Raw Wastewater Pipeline Networks in Northern Cities Scenario: Jinan City, Shandong Province, uses three main raw sewage pipe networks to provide heat sources for six heat pump energy stations. The pipe network contains 47 key nodes (the same pipe network system as in Section 6.1 of the Beneficial Effects section).

[0117] Dataset Overview

[0118] Predictive effect 1. High turbidity warning (>500 NTU) advance time: 45 minutes (traditional methods only 15 minutes) 2. Arrival time prediction error: <3 minutes 3. Peak concentration prediction error: <8% Adaptive sampling control effect verification To verify the effectiveness of the adaptive sampling control module, a comparative experiment was designed to compare fixed sampling and adaptive sampling. Experimental settings: baseline sampling interval 15 minutes (Δt=15min), adaptive sampling range [5min, 30min], hyperparameter α=0.5.

[0119]

[0120] Experimental results show that the adaptive sampling strategy can effectively reduce data transmission volume (saving 32%) and improve the operating efficiency of the energy station (+2.1%) while maintaining prediction accuracy, thus verifying the practical application value of the module.

[0121] Example 2: Prediction of water quality impact during sporting events (e.g., Figure 3 (As shown) Scenario: During a large-scale football match held at the Jinan Olympic Sports Center, the influx and dispersal of 60,000 people within a short period caused "tidal" fluctuations in water quality in the surrounding pipeline network. The Jinan Olympic Sports Center provides heat sources for five heat pump energy stations, and the pipeline network includes 32 key nodes.

[0122] Data source: Sensor data from 12 large-scale football matches held in the urban area where the Jinan Olympic Sports Center is located between March and December 2023 (each match was held approximately 2-4 weeks apart).

[0123] Technology adaptation 1. External covariate encoding: The event schedule (start of match at 19:00, halftime break from 19:45 to 20:00, and end of match at 21:00) is encoded as a highly periodic feature, and hourly sine / cosine encoding (24-hour period) is added. 2. Rapid Evolution of Dynamic Graphs: The adjacency matrix update cycle is shortened to 5 minutes (compared to the baseline of 15 minutes) to capture rapidly changing hydraulic topology. 3. Multi-scale forecasting: Focus on 15-minute and 30-minute forecast windows. 4. Event-Specific Physical Constraint Adjustment: To address the impact of high pedestrian traffic on sedimentation (SS), the deposition rate coefficient has been increased. Upper limit constraint Predictive effect

[0124] Comparison with benchmark methods

[0125] Typical prediction results during the event One hour before the show: The model accurately predicted the gradual increase in audience attendance. Halftime break: Accurately capture the "cliff-like" drop in traffic from peak to 60% within 15 minutes. Closing period: Successfully predicted a 30-minute "long-tail" emission pattern. Heating water quality assurance: The five heat pump power stations did not shut down due to water quality issues during the 12 matches. Example 3: Source Tracing of Sudden Pollution from Industrial Leaks Scenario: An upstream factory illegally discharged pollutants, causing high concentrations of pollutants to spread downstream along the pipeline network. This case validates the PC-DGSSN's capabilities in tracing the source of sudden pollution events.

[0126] Data source: Five confirmed industrial leak incidents recorded in Jinan City between 2022 and 2023. Each incident has a clear record of the leak location, the type of leaked substance, and the amount of leaked substance.

[0127] Traceability technology solution 1. Backflow Source Analysis: Utilizing the reverse topological propagation characteristics of dynamic graphs, combined with changes in pollutant concentration gradients, to pinpoint the source of anomalies. 2. Physical Constraint Assistance: By inversely calculating the deposition-resuspension source term, the initial concentration and time of the leaked substance are estimated. 3. Enhanced event triggering: Leak event triggering event indicator Enhance the state-space model's ability to capture anomalous patterns. Source traceability capability verification

[0128] Typical Case 1. In April 2023, a chemical plant illegally discharged wastewater containing high concentrations of benzene compounds. 2. The model accurately predicted that pollutants would reach the water intake 8 kilometers downstream after 12 hours. 3. The predicted peak concentration was 480 mg / L, while the measured value was 512 mg / L (error 6.3%). 4. The management department promptly activated the emergency plan based on the early warning, thus preventing pollution of the water intake. Example 4: Extreme Weather Event (Rainstorms and Floods) Scenario Scenario: A northern city experiences continuous heavy rainfall during the rainy season, leading to a surge in pipeline load and an increased risk of raw sewage overflow. This scenario validates the predictive robustness of PC-DGSSN under extreme conditions.

[0129] Data source: Sensor data from Jinan City during the flood season in July-August 2022 (July 15 - August 20, a total of 37 days), including 4 rainstorm events (rainfall > 50 mm / hour).

[0130] Technology adaptation 1. External Covariate Deep Integration: Take rainfall data as an external forced input and map it to each module of the model through a dedicated encoding layer 2. Dynamic Graph Fast Reconstruction: During heavy rainstorms, the pipe network topology changes violently. Shorten the update period of the adjacency matrix to 5 minutes and increase the smoothing coefficient To balance stability and response speed 3. Adaptive Adjustment of Physical Constraints: During heavy rainstorms, the SS concentration surges, and raise the upper limit constraint of to 1200 mg / L 4. Multi-scale Event Detection: Introduce hourly cumulative rainfall as an auxiliary event detection signal and jointly judge with the flow mutation signal Prediction effect

[0131] Comparison with the baseline method (heavy rain scenario)

[0132] Typical prediction effect in heavy rain scenario 1. One hour before heavy rain: Accurately predict the progressive upward trend of the pipe network load 2. 15 minutes after the start of heavy rain: Capture the "shock wave" effect caused by the sudden increase in flow 3. During the continuous heavy rain: Continuously track the water quality changes under high load conditions 4. After the heavy rain ends: Predict the recession curve and the time for pollutant concentration to fall back Heating System Safety Guarantee: During 4 heavy rain events, the turbidity of the water intake of 6 heat pump energy stations was always controlled within the allowable range (<50 NTU), and the heating system operated normally. <00​​​​​​​​​​​

[0135] 5.2 Node Type Distribution

[0136] 5.3 Data Partitioning Strategy Time-aware data partitioning is used to prevent future information leaks.

[0137] 5.4 Anomaly Event Labeling The following events in the test set were manually labeled by water sector experts: Sudden pollution incidents: Abnormal concentration peaks caused by industrial leaks, illegal emissions, etc. (47 incidents in total) Heavy rain events: Sudden changes in flow rate caused by rainfall >50mm / 24h (23 events in total) Equipment failure events: sensor failure, pump station shutdown, etc. (31 events in total) Example 6: Experimental Setup and Comparison Methods 6.1 Experimental Environment Configuration

[0138] 6.2 Comparison Benchmark Method The present invention selects the following representative methods as comparative benchmarks:

[0139] Example 7: Comparative Experimental Results 7.1 Overall Performance Comparison The overall performance of each method on the USN-MSHD test set is shown in the table below. Note: In this embodiment, the MAE / RMSE values ​​are dimensionless values ​​after Z-score standardization (each water quality parameter is standardized to a mean of 0 and a standard deviation of 1) to facilitate comprehensive comparison of parameters with different dimensions (SS, COD, NH3-N, etc.). For the comparison results of the original physical scale (mg / L), please refer to Section 6.5 of the Beneficial Effects section (Jinan City dataset, MAE=5.8±0.2 mg / L).

[0140]

[0141] 7.2 Performance Improvement Analysis Compared to the existing state-of-the-art (SOTA) method (PatchTST), the performance improvement of this invention on the USN-MSHD normalized dataset is as follows (see Section 6.5 for the original physical scale results, MAE reduced by 52.8%):

[0142] 7.3 Performance of different prediction step sizes

[0143] Example 8: Ablation Experiment To verify the effectiveness of each technical module of this invention, the following ablation experiment was designed: 8.1 Analysis of the Contribution of Each Module

[0144] Conclusion Analysis 1. The physical constraint in the solid-liquid two-phase system contributes the most (MAE decreases by 16.6%), demonstrating the importance of prior physical knowledge. 2. The heterogeneous node classification perception was significantly improved (MAE decreased by 23.0%), validating the necessity of node type differentiation. 3. Dual-branch hybrid pre-training further improves performance (MAE reduced by 10.6%), demonstrating the advantages of generation-contrast synergy. 4. Event-triggered SSM and multi-scale fusion provide marginal improvements and enhance long-range dependency modeling capabilities. 8.2 Sensitivity Analysis of Physical Constraint Weights

[0145] Conclusion: There is an optimal balance point for the physical constraint weights. If they are too small, the physical consistency will be insufficient, and if they are too large, the expressive power of the model may be limited.

[0146] 8.3 Impact of Dynamic Graph Update Frequency

[0147] Conclusion: A 15-minute update cycle strikes the best balance between accuracy and efficiency.

[0148] 8.4 Comparison of Pre-training Strategies

[0149] Conclusion: The BCGP hybrid pre-training framework of this invention not only improves performance but also reduces the requirement for labeled data by 70%, significantly alleviating the problem of label scarcity.

[0150] Example 9: Case Study on the Detection of Sudden Pollution Incidents 9.1 Case Background Between 2:00 and 4:00 AM on November 15, 2024, an abnormal increase in COD concentration, suspected to be caused by an industrial leak, was detected.

[0151] 9.2 Comparison of Prediction Results

[0152] 9.3 Case Analysis The advantages of this invention in this emergency are as follows: 1. The physical constraint mechanism enables the model to understand that "it takes time for pollutants to be transported with the water flow," thus avoiding the instantaneous jump prediction bias of purely data-driven methods. 2. The dynamic topology captures the flow direction changes after the leak occurred (backflow was observed in some inspection wells). 3. Event-triggered SSM automatically enhances its sensitivity to sudden traffic changes, enabling earlier warnings. Example 10: Verification of Cross-Pipeline Generalization Capability 10.1 Transfer Experiment Setup The model trained on the city A dataset of this invention is directly transferred to the city B and city C datasets for zero-shot / few-shot testing:

[0153] 10.2 Migration Performance Results

[0154] Conclusion: Thanks to the cross-domain generalization capability provided by physical constraints and the structural inductive bias of heterogeneous node classification, the cross-network migration performance of the present invention is significantly better than that of the pure data-driven method.

[0155] Example 11: Quality Assessment of Uncertainty Estimation 11.1 Calibration Curve Analysis The actual coverage of the 90% confidence interval output by this invention:

[0156] 11.2 Application of Uncertainty in Early Warning Adaptive threshold adjustment based on uncertainty interval: During normal periods: Use a fixed threshold (e.g., trigger an alert if COD > 150 mg / L) High uncertainty period: Lower the threshold and increase manual review Experiments show that after introducing uncertainty estimation, the false alarm rate decreased by 32.1% and the false negative rate decreased by 18.7%.

[0157] In summary, the technical solution described in this invention has the following advantages compared with the prior art: 1. Physical consistency assurance: By constraining the solid-liquid two-phase coupling convection-diffusion equation, the model prediction results satisfy the mass conservation and sedimentation-resuspension dynamics, significantly improving the prediction robustness in off-distribution scenarios (rainstorms, leaks), and increasing the high turbidity warning lead time from 15 minutes to 45 minutes.

[0158] 2. Dynamic topology capture capability: The heterogeneous node classification and perception dynamic graph construction mechanism can distinguish the essential hydraulic interaction differences between different node types. After low-rank decomposition and regularization control, the number of parameters of the 16 type interaction matrices is reduced by about 8 times, while maintaining accurate capture of the dynamic topology reconstruction of the pipeline network of "convergence-diversion-reverse flow".

[0159] 3. Improved pre-training efficiency: The BCGP dual-branch hybrid pre-training framework achieves efficient pre-training under conditions of scarce labeled data through a collaborative mechanism of noise-resistant generation and event-aware comparison. Ablation experiments show that BCGP pre-training improves MAE by 33.3% relative to the benchmark.

[0160] 4. Real-time inference efficiency: The event-triggered selective state-space model processes long sequences with O(N) linear complexity, and the single-step inference time is only 38ms, which is more than 60% lower than the Transformer architecture, meeting the real-time early warning requirements of heat pump heating systems.

[0161] 5. Overall Performance Improvement: Based on 24 months of real-world testing data in Jinan, the PC-DGSSN achieved a MAE of 5.8 ± 0.2, a 52.8% reduction compared to existing state-of-the-art methods, with the shortest inference time. Its effectiveness was verified in four typical scenarios (regular, event, data leak, and heavy rain). Note: p<0.001 indicates highly significant, p<0.01 indicates significant, and p<0.05 indicates marginally significant.

[0162] This invention also provides a native wastewater quality prediction system based on a physically constrained dynamic graph state-space network, comprising: Data acquisition module: used to collect time-series data from multi-source heterogeneous sensors in the urban raw sewage pipe network, and to preprocess the data and label the node types; Dynamic graph construction module: connected to the data acquisition module, used to receive preprocessed time-series data and node type labels, construct a heterogeneous node classification perception dynamic graph, and use 16 type-specific interaction matrices to model the hydraulic association between different node type pairs to achieve adaptive learning of time-varying adjacency matrices; Physical constraint encoding module: connected to the dynamic graph construction module, used to receive dynamic graph structure information, transform the solid-liquid two-phase coupled convection-diffusion equation into a differentiable physical constraint loss function, and constrain the physical consistency of the neural network output; Hybrid pre-training module: connected to the data acquisition module and the dynamic graph construction module respectively, used to receive preprocessed time series data and dynamic graph structure, including noise-resistant generation branch and event-aware comparison branch, perform self-supervised pre-training on unlabeled data, learn noise-robust and structurally consistent node representations, and pass the pre-training parameters to the state space modeling module; State space modeling module: connected to the dynamic graph construction module, the physical constraint encoding module and the hybrid pre-training module respectively, used to receive time-varying adjacency matrix, physical constraint loss and pre-trained node representation, and capture long-range temporal dependencies of water quality parameters through event-triggered selective state space model and multi-scale time fusion mechanism; Prediction output module: connected to the state space modeling module, used to receive the state space modeling output, generate predicted values ​​and uncertainty ranges of water quality parameters, and generate early warning signals based on the prediction results.

[0163] The system is deployed in an edge-cloud collaborative architecture: the edge is responsible for data preprocessing and adjacency matrix updates, while the cloud is responsible for model inference and early warning notifications; the edge and cloud communicate via the MQTT protocol. The system's predictive output is linked to the heat pump heating system: when water quality parameters are predicted to deteriorate or have excessive uncertainty, the system automatically adjusts the heat pump operating parameters or switches to a backup water source to ensure the stable operation of the heating system.

[0164] In another aspect, the present invention also discloses a computer-readable storage medium storing a computer program, which, when executed by a processor, causes the processor to perform the steps of the method described above.

[0165] In another aspect, the present invention also discloses a computer device, including a memory and a processor, wherein the memory stores a computer program, and when the computer program is executed by the processor, the processor performs the steps of the method described above.

[0166] In another embodiment provided in this application, a computer program product containing instructions is also provided, which, when run on a computer, causes the computer to execute any of the mobile source emission prediction methods based on time-series feature migration described in the above embodiments.

[0167] It is understood that the systems, devices, and storage media provided in the embodiments of the present invention correspond to the methods provided in the embodiments of the present invention, and the explanations, examples, and beneficial effects of the relevant content can be referred to the corresponding parts of the above methods.

[0168] In the above embodiments, implementation can be achieved entirely or partially through software, hardware, firmware, or any combination thereof. When implemented using software, it can be implemented entirely or partially as a computer program product. The computer program product includes one or more computer instructions. When the computer program instructions are loaded and executed on a computer, all or part of the processes or functions described in the embodiments of this application are generated. The computer can be a general-purpose computer, a special-purpose computer, a computer network, or other programmable device. The computer instructions can be stored in a computer-readable storage medium or transmitted from one computer-readable storage medium to another. For example, the computer instructions can be transmitted from one website, computer, server, or data center to another website, computer, server, or data center via wired (e.g., coaxial cable, fiber optic, digital subscriber line (DSL)) or wireless (e.g., infrared, wireless, microwave, etc.) means. The computer-readable storage medium can be any available medium that a computer can access or a data storage device such as a server or data center that integrates one or more available media. The available medium can be a magnetic medium (e.g., floppy disk, hard disk, magnetic tape), an optical medium (e.g., DVD), or a semiconductor medium (e.g., solid state disk (SSD)).

Claims

1. A method for predicting the water quality of primary wastewater based on a physically constrained dynamic graph state-space network, characterized in that, Perform the following steps using a computer device: S1: Collect time-series data from multi-source heterogeneous sensors in the urban raw sewage pipe network, including water quality parameters and hydraulic parameters; perform missing value filling, outlier detection, and feature standardization on the collected data; The network nodes are labeled with their types based on the sensor installation locations. Node types include: pump stations, inspection wells, confluence points, and water intakes. S2. Based on the node attributes and topology obtained in step S1, construct a heterogeneous graph of the original sewage pipe network; Then, a classification-aware time-varying adjacency matrix learning mechanism is adopted. The specific interaction weight matrix and node embedding are used with the learnable node type. Combined with Gaussian kernel distance weights, the elements of the time-varying adjacency matrix are calculated, and an adaptive update mechanism is used for time-series smoothing. S3: The solid-liquid two-phase coupled convection-diffusion equation is embedded as a physical prior constraint into the neural network training process; S4: A collaborative pre-training framework is adopted, consisting of a noise-resistant generation branch and an event-aware contrast branch. S5: Use an event-triggered selective state-space model to capture long-term temporal dependencies and sudden pollution sensitivity characteristics; when the flow rate change rate is detected to exceed the threshold, activate the event enhancement path; adopt a multi-scale time modeling strategy to extract features at different time scales with different expansion coefficients, and after adaptive fusion with learnable weights, integrate them with the event-triggered state through a gating mechanism to obtain the final hidden state vector. S6: Outputs predicted water quality parameters and uncertainty range for future periods; triggers an early warning signal when the predicted value exceeds the preset threshold or the uncertainty range is too large.

2. The method for predicting the quality of primary wastewater based on a physically constrained dynamic graph state-space network according to claim 1, characterized in that, The heterogeneous graph expression for the native sewage network in step S2 is: ； In the formula: For a moment Heterogeneous diagram of the original sewage pipe network; A set of nodes; For a moment The set of edges; For a moment The node feature matrix; A vector of node type labels; In step S2, the expression for the elements of the time-varying adjacency matrix is: ； In the formula: For a moment node With nodes The time-varying adjacency matrix elements between them take values ​​in the range [0,1]. The Sigmoid activation function maps the output to the [0,1] interval; For a moment node The first set of time-varying node embedding vectors has a dimension of ; superscript This is the matrix transpose operator; For node type pairs The corresponding learnable type interaction weight matrix has the following dimensions: ,common Each matrix corresponds to a hydraulic interaction mode between four node types: pumping station, inspection well, confluence point, and water intake. For a moment node The second set of time-varying node embedding vectors has a dimension of ; This is the scalar multiplication operator; For exponential function operators; For nodes With nodes Geographical distance or pipe segment length between them; The Gaussian kernel width parameter controls the sensitivity to distance decay. The type interaction weight matrix is ​​parameterized using low-rank decomposition. ,in , , As a bottleneck dimension, For the hidden layer dimension; The adjacency matrix employs an adaptive update mechanism with a basic update cycle of 15 minutes. An immediate update is triggered when the flow change rate |ΔQ / Q|>0.3 is detected, and a smoothing coefficient is used. Perform time-series smoothing: ； In the formula: For a moment The adjacency matrix elements after time-series smoothing have values ​​in the range [0,1]. The smoothing coefficient has a value range of (0,1), and in this invention, it is set to 0.

7. This is the scalar multiplication operator; For a moment Historical adjacency matrix elements; It is the addition operator; For the new information weighting coefficient, and Complementary and sum to 1; For a moment The newly calculated elements of the original adjacency matrix.

3. The method for predicting the quality of primary wastewater based on a physically constrained dynamic graph state-space network according to claim 1, characterized in that, The physical constraint loss function in step S3 is: In the formula: The loss function is the physical constraint loss function for solid-liquid two-phase coupling convection-diffusion. This represents the total number of time steps. This is the summation operator; This is the time step index variable; L2 norm operator; The partial derivative of solid concentration with respect to time describes the transient rate of change of concentration; For spatial gradient operators; This is the vector dot product operator; For the fluid velocity vector field; This is the concentration of solids, in mg / L. The spatial divergence of the convection term describes the transport process of solid particles carried by the fluid. It is the subtraction operator; This is the solid-state diffusion coefficient tensor, in m² / s. The spatial divergence of the diffusion term describes the molecular diffusion and turbulent diffusion processes of solid particles; For the sedimentation-resuspension source term function, in terms of solid concentration Calculated as an independent variable; In the formula: For the sedimentation-resuspension source term function, in terms of solid concentration As an independent variable; is the deposition rate coefficient, with a value range of [0.01, 0.5]. This is the scalar multiplication operator; Solid concentration; It is the subtraction operator; The resuspension rate coefficient has a value range of [0.001, 0.1]. The term represents the square of the solid concentration, reflecting the nonlinear enhancement effect of resuspension at high concentrations; The saturated deposition concentration has a value range of [500, 2000].

4. The method for predicting the quality of primary wastewater based on a physically constrained dynamic graph state-space network according to claim 1, characterized in that, Step S4 includes: S41, Noise Reduction Generation Branch: Employs a segmented masking strategy, prioritizing masking abnormal sensor regions at a 30% ratio and randomly masking normal regions at a 15% ratio, then performs a mask auto-encoding reconstruction task. S42, Event-Aware Comparison Branch: Triggers enhanced view generation through traffic mutation detection, and performs a comparison learning task based on InfoNCE; Steps S41 and S42 together constitute the pre-training total loss function, the expression of which is: ； In the formula: Adaptive dual-branch comparison for native wastewater domain - generating hybrid pre-trained total loss function; To mitigate the noise-resistant mask autoencoder reconstruction loss, the reconstruction accuracy of mask features is measured. It is the addition operator; To compare the loss weighting coefficient, a value of 1.0 was used; This is the scalar multiplication operator; For event-aware contrastive learning loss, maximize the similarity of source views based on InfoNCE; This is the event consistency loss weighting coefficient, with a value of 0.5; This is the event consistency loss function, ensuring the consistency of the representation across nodes at the time the event is triggered; The structural reconstruction loss weighting coefficient has a value of 0.

3. The loss function is reconstructed to fit the adjacency matrix structure, measuring the difference between the predicted and actual adjacency matrices. The type-sharing regularization loss weight coefficient is set to 0.2, which controls the strength of common constraint on the interaction matrix of nodes of the same type. The type-sharing regularization loss function forces all interaction matrices of nodes of the same type to converge to a shared basis vector.

5. The method for predicting the quality of primary wastewater based on a physically constrained dynamic graph state-space network according to claim 1, characterized in that, The final hidden state vector expression in step S5 is: In the formula: The final hidden state vector after gating fusion has a dimension of . ; The gating weight vector has dimension 1. The value range is [0,1], which controls the fusion ratio between the event trigger state and the multi-scale fusion state; This is Hadamard's element-wise multiplication operator; For the event-triggered state space model at time... The output hidden state vector has dimension . It is sensitive to sudden pollution incidents; This is the vector addition operator; For complementary gating weight vectors, and Each element is complementary and their sum is 1; The hidden state vector is a multi-scale fused vector with dimension . Encodes multi-scale periodic pattern information; The Sigmoid activation function maps the output to the [0,1] interval; The gated projection weight matrix has dimensions of . ; For vector concatenation operators, and spliced ​​as Dimensional vector.

6. The method for predicting the quality of primary wastewater based on a physically constrained dynamic graph state-space network according to claim 1, characterized in that, In step S2, type-sharing regularization loss is also introduced: In the formula: The type-sharing regularization loss function forces all interaction matrices of nodes of the same type to converge to the shared basis vectors; This is a summation operator for four node types (pump station, inspection well, confluence point, and intake); This is a node type index variable, with values ​​ranging from 1 to 4; For belonging to type All node pairs The summation operator; For type The set of nodes; The Frobenius norm operator measures the degree of difference between matrices; For node type pairs The corresponding learnable type interaction weight matrix has the following dimensions: ; This is the matrix subtraction operator; For type The average interaction matrix, by classifying the types The arithmetic mean of all interaction weight matrices is obtained.

7. The original wastewater quality prediction method based on physical constraint dynamic graph state space network according to claim 6, wherein the multi-scale modeling in step S5 adopts three parallel branches: fine-scale branch, mesoscale branch, and coarse-scale branch; the output of each branch is weighted and fused through learnable weights and then input into a gating mechanism, which together with the event state determines the final output.

8. A raw wastewater quality prediction system based on a physically constrained dynamic graph state-space network, comprising the raw wastewater quality prediction method based on a physically constrained dynamic graph state-space network as described in any one of claims 1-7, characterized in that, include: Data acquisition module: used to collect time-series data from multi-source heterogeneous sensors in the urban raw sewage pipe network, and to preprocess the data and label the node types; Dynamic graph construction module: connected to the data acquisition module, used to receive preprocessed time-series data and node type labels, construct a heterogeneous node classification perception dynamic graph, and use 16 type-specific interaction matrices to model the hydraulic association between different node type pairs to achieve adaptive learning of time-varying adjacency matrices; Physical constraint encoding module: connected to the dynamic graph construction module, used to receive dynamic graph structure information, transform the solid-liquid two-phase coupled convection-diffusion equation into a differentiable physical constraint loss function, and constrain the physical consistency of the neural network output; Hybrid pre-training module: connected to the data acquisition module and the dynamic graph construction module respectively, used to receive preprocessed time series data and dynamic graph structure, including noise-resistant generation branch and event-aware comparison branch, perform self-supervised pre-training on unlabeled data, learn noise-robust and structurally consistent node representations, and pass the pre-training parameters to the state space modeling module; State space modeling module: connected to the dynamic graph construction module, the physical constraint encoding module and the hybrid pre-training module respectively, used to receive time-varying adjacency matrix, physical constraint loss and pre-trained node representation, and capture long-range temporal dependencies of water quality parameters through event-triggered selective state space model and multi-scale time fusion mechanism; Prediction output module: connected to the state space modeling module, used to receive the state space modeling output, generate predicted values ​​and uncertainty ranges of water quality parameters, and generate early warning signals based on the prediction results.

9. The native wastewater quality prediction system based on a physically constrained dynamic graph state-space network according to claim 8, characterized in that, The system is deployed in an edge-cloud collaborative architecture: the edge is responsible for data preprocessing and adjacency matrix updates, while the cloud is responsible for model inference and early warning notifications; the edge and cloud communicate via the MQTT protocol.

10. The native wastewater quality prediction system based on a physically constrained dynamic graph state-space network according to claim 8, characterized in that, The system's predictive output is linked to the heat pump heating system: when water quality parameters are predicted to deteriorate or have excessive uncertainty, the system automatically adjusts the heat pump operating parameters or switches to a backup water source to ensure the stable operation of the heating system.

Images