A sewage water quality on-line monitoring method and system

Abstract

The application relates to the technical field of water quality monitoring, and discloses a sewage water quality online monitoring method and system, wherein the sewage water quality online monitoring method comprises the following steps: performing time alignment processing and data quality evaluation on original water quality sequences collected by each monitoring point, obtaining aligned data and data quality weights; constructing a time delay space graph; jointly estimating a global water quality state, obtaining a global state estimation and uncertainty; performing abnormality detection and identifying an abnormal propagation path and a suspected source, obtaining an abnormality score and a source confidence ranking list; and formulating an adaptive sampling strategy and determining the sampling frequency of each monitoring point; the application realizes online joint estimation of a global water quality hidden state based on a time delay observation equation and weighted optimization, and performs reverse tracing by combining residual standardization abnormality scores and the time delay space graph, so that the abnormality can be not only found, but also positioned in terms of possible sources and propagation directions, and the pertinence and interpretability of disposal are improved.

Description

Technical Field

[0001] This invention relates to the field of water quality monitoring technology, and more specifically, to a method and system for online monitoring of wastewater quality. Background Technology

[0002] Wastewater treatment plants typically consist of multiple units connected in series and parallel, including influent pipe networks, booster pump stations, equalization tanks, primary sedimentation / flotation tanks, biological reaction tanks, secondary sedimentation / membrane treatment, and advanced treatment and discharge outlets. They also involve complex hydraulic paths such as sludge return, internal recirculation, bypass diversion, and combined sewer overflow. In this multi-pipe, multi-tank, multi-path mixing system, water quality exhibits significant non-uniformity and time-varying characteristics along the spatial direction. Different branch influent loads, intermittent industrial discharges, and rainfall runoff can cause significant differences in pollutant concentrations, conductivity, turbidity, and other indicators at different locations at the same time. Short-circuit flow, dead zones, local recirculation, and insufficient mixing introduce migration delays and spatial time lags, making the arrival time of upstream disturbances at each treatment unit and effluent outlet uncertain.

[0003] A Chinese patent with authorization announcement number CN118409064B discloses a water quality change monitoring system for wastewater treatment, including a multi-parameter sensor monitoring module for monitoring different parameters of water quality during the wastewater treatment process using different sensors. First, the water quality parameters to be monitored during the wastewater treatment process are determined, including pH value, dissolved oxygen, chemical oxygen demand, ammonia nitrogen, total phosphorus, turbidity, and temperature. Then, corresponding sensors are selected according to the requirements, with each sensor targeting specific water quality parameters.

[0004] However, in continuous operation scenarios involving pipelines, multiple pools, and multiple paths, the arrival time of upstream impact loads to downstream points can drift due to the wide distribution of hydraulic residence time and the changing paths of backflow and bypass with operating conditions. Furthermore, the order and magnitude of the same anomaly at different monitoring points can change due to mixing, reaction, and settling processes. When the system only relies on water quality changes or threshold / statistical characteristics at fixed points for alarms and recording, it is often difficult to reliably infer where the anomaly comes from, along what path it propagates, and which areas are at higher risk. Consequently, it is difficult to form targeted monitoring, scheduling, and response strategies, resulting in problems such as a compressed response window after anomaly detection and incomplete post-event evidence chains. Summary of the Invention

[0005] The purpose of this invention is to provide a method and system for online monitoring of wastewater quality to solve the above-mentioned technical problems.

[0006] This invention provides a method for online monitoring of wastewater quality, comprising the following steps:

[0007] The original water quality sequences collected from each monitoring point are time-aligned and data quality is assessed to obtain the aligned data and data quality weights for the monitoring points.

[0008] Based on the aligned data, a time-delay space graph is constructed, which includes a set of nodes, a set of edges, and the migration delay and decay coefficient of each directed edge.

[0009] Based on the aligned data, data quality weights, and time-delay space map, the global water quality state is jointly estimated to obtain the global state estimate and uncertainty.

[0010] Based on the global state estimation, uncertainty, and time delay space graph, anomaly detection is performed and the propagation path and suspected source of the anomaly are identified, resulting in anomaly scores and a source confidence ranking list.

[0011] Based on the uncertainty, anomaly score, and source confidence ranking list, an adaptive sampling strategy is formulated to determine the sampling frequency for each monitoring point.

[0012] Furthermore, the aligned data and data quality weights of the monitoring points are obtained as follows:

[0013] The original water quality sequence includes the measured values ​​of water quality indicators collected at each monitoring point at each time point. The water quality indicators include pH, dissolved oxygen, oxidation-reduction potential, turbidity, conductivity, ammonia nitrogen, and chemical oxygen demand.

[0014] The monitoring points include point metadata, which includes the spatial location of the monitoring point, the unit it belongs to, the sensor type, the sampling period, and communication quality characterization parameters, including packet loss rate and clock deviation.

[0015] Time alignment processing includes estimating the clock deviation of each monitoring point and obtaining the aligned measurement sequence through a resampling method;

[0016] Data quality assessment includes calculating data quality weights, specifically: multiplying the packet loss rate by a preset packet loss rate attenuation coefficient and taking the negative value as the exponent of the first exponential function; multiplying the standardized jump intensity by a preset jump attenuation coefficient and taking the negative value as the exponent of the second exponential function; and then multiplying the two exponential functions together to obtain the data quality weights; the standardized jump intensity is the standardized result of the original jump amplitude, and the original jump amplitude is the absolute value of the difference between the currently aligned data and the mean of the aligned data within the local sliding window.

[0017] Furthermore, constructing the time-delay spatial diagram includes:

[0018] Each monitoring point is used as a node to form a set of nodes, and the hydraulic connectivity and mixing relationship in the same pool are used as directed edges to form a set of edges. The direction of water flow determines the direction of the directed edges.

[0019] For each directed edge, a learnable migration delay and attenuation coefficient are introduced, whereby the migration delay represents the time required for a water quality change to propagate from an upstream node to a downstream node, and the attenuation coefficient represents the degree of attenuation of the water quality index during the propagation process.

[0020] The initial value of the migration time delay is the pool volume of the downstream unit divided by the representative flow rate of the unit, and the initial value of the attenuation coefficient is set according to the biochemical degradation characteristics of the water quality index or historical experience.

[0021] The attenuation coefficient is continuously corrected using real-time data. Specifically, within a preset attenuation coefficient update window, the aligned data of the downstream node is divided by the aligned data of the upstream node after migration delay to obtain a ratio. Within the update window, the median or weighted average of the ratios at multiple times is taken as the updated value of the attenuation coefficient.

[0022] The migration lag is corrected online using real-time data, and updated within a sliding time window by maximizing alignment correlation or minimizing prediction error.

[0023] Furthermore, the joint estimation of the overall water quality status includes:

[0024] Set a global hidden state vector to represent the true water quality state of all units or regions at a certain moment;

[0025] The relationship between the aligned data and the global hidden state is described by a time-delay observation equation, which represents the aligned data at the current moment as a function of the global hidden state vector at the historical moment.

[0026] A weighted optimization method is used for online estimation of the global hidden state. The estimated value is obtained by solving a weighted least squares optimization problem. The estimated value of the global hidden state vector is the global state estimate. The objective function is composed of the sum of the observation fitting term and the smoothing regularization term.

[0027] Calculate the uncertainty of the global hidden state estimation.

[0028] Furthermore, the specific form of the time-delay observation equation is as follows: the aligned data at the current moment is equal to the time-delay observation matrix multiplied by the extended state vector, plus the observation noise term; the extended state vector is formed by stacking the global hidden state vectors at the current moment and several historical moments in chronological order; the time-delay observation matrix is ​​determined by the topology of the time-delay spatial graph, the migration time delay parameters of each side, and the mapping relationship from the monitoring point to the region.

[0029] The observation fitting term is the weighted sum of squared observation fitting errors of all monitoring points. The observation fitting error of each monitoring point is the difference between the aligned data of the monitoring point and the corresponding predicted observation value. The weight is the data quality weight. The smoothing regularization term is the square of the norm of the Laplacian matrix of the time-delay spatial graph applied to the global hidden state vector at the current time. The resulting value is multiplied by a preset smoothing regularization coefficient.

[0030] Furthermore, anomaly detection and identification of abnormal propagation paths and suspected sources include:

[0031] Calculate the regional residual for each monitoring point, where the regional residual is the difference between the actual observed value and the predicted observed value;

[0032] Based on the regional residuals, a standardized anomaly score is constructed. The anomaly score of a monitoring point is the absolute value of the residual of the monitoring point divided by a marginal term, which is the sum of the marginal uncertainty of the monitoring point and a preset minimum positive number.

[0033] The set of nodes with high abnormal scores is identified based on the abnormal scores. The set of nodes with high abnormal scores includes all monitoring points whose abnormal scores exceed a preset abnormal score threshold.

[0034] By combining the aforementioned time-delay spatial diagram with reverse tracing, candidate source points are obtained;

[0035] For each candidate source point, calculate its explanatory power for the set of high anomaly score nodes, and sort the candidate source points according to the explanatory power to obtain a source confidence ranking list.

[0036] Further, the explanatory power is calculated as follows: for each node in the set of high-anomaly-score nodes, the product of the time matching factor, the path uncertainty discount factor, and the anomaly score of the node is calculated, and then the product results of all nodes are summed; the time matching factor is an exponential decay function, the base of the exponential decay function is the natural constant, the time term with a negative exponent is divided by a preset time matching decay scale, the time term is the current time minus the time when the node anomaly occurs, plus the cumulative migration delay from the candidate source point to the node, and the absolute value of the obtained result is taken; the path uncertainty discount factor is obtained by summing the time delay uncertainties of each edge on the path according to the variance propagation law to obtain the total path time delay uncertainty, and then constructing a discount factor that decreases as the total path time delay uncertainty increases; the path is the path with the minimum cumulative migration delay from the candidate source point to the target node in the time delay space graph.

[0037] Furthermore, developing an adaptive sampling strategy includes:

[0038] An information value is assigned to each monitoring point. The information value is calculated by multiplying the normalized risk term, the normalized uncertainty term, and the data quality weight.

[0039] The sampling frequency is determined based on the information value. The sampling frequency of the monitoring point is the product of the preset scheduling gain and the information value plus the preset lower limit of the sampling frequency, and then limited to between the preset lower limit of the sampling frequency and the preset upper limit of the sampling frequency through amplitude limiting processing.

[0040] Optimize sampling allocation under communication or energy consumption budget constraints. When the sum of the sampling frequencies of all monitoring points exceeds the preset total communication budget, sort all monitoring points from largest to smallest information value. Starting from the last monitoring point, reduce the sampling frequency to the preset lower limit of the sampling frequency until the sum of the sampling frequencies is less than the preset total communication budget.

[0041] Furthermore, the normalized risk item is calculated as follows: first, the original risk item is calculated, which is the anomaly score; then, the original risk item is divided by the sum of the maximum value of the original risk items of all monitoring points and the preset minimum positive number to obtain the normalized risk item with a value range of 0 to 1.

[0042] The normalized uncertainty term is calculated as follows: First, the original uncertainty term is calculated, which is the marginal uncertainty of the monitoring point. The original uncertainty term is divided by the sum of the maximum value of the original uncertainty terms of all monitoring points and the preset minimum positive number to obtain the normalized uncertainty term with a value range of 0 to 1.

[0043] This invention provides an online wastewater quality monitoring system for storing computer-readable instructions, which, when read, execute the aforementioned online wastewater quality monitoring method; the system includes:

[0044] The data processing module performs time alignment processing and data quality assessment on the raw water quality sequences collected from each monitoring point to obtain the aligned data and data quality weights of the monitoring points.

[0045] The time-delay decay module constructs a time-delay space graph based on the aligned data. The time-delay space graph includes a set of nodes, a set of edges, and the migration delay and decay coefficient of each directed edge.

[0046] The global estimation module performs a joint estimation of the global water quality state based on the aligned data, data quality weights, and time-delay space map, and obtains the global state estimate and uncertainty.

[0047] The source confidence module, based on the global state estimation, uncertainty and time delay space graph, performs anomaly detection and identifies the propagation path and suspected source of the anomaly, and obtains anomaly score and source confidence ranking list;

[0048] The adaptive sampling module formulates an adaptive sampling strategy based on the uncertainty, anomaly score, and source confidence ranking list, and determines the sampling frequency for each monitoring point.

[0049] The beneficial effects of this invention are as follows: By performing time alignment and data quality assessment on the original water quality sequences from multiple locations, this invention introduces packet loss rate, clock deviation, and jump detection to form data quality weights, ensuring that cross-location data still has a reliable basis for joint analysis even under conditions of unstable communication and asynchronous sampling, significantly reducing the interference of low-quality data on monitoring conclusions; by constructing a time-delay spatial graph that includes hydraulic connectivity and pool mixing relationships, and introducing learnable migration delays and attenuation coefficients for the edges, propagation parameters can be corrected online and physical rationality can be maintained under complex operating conditions such as multi-path mixing, backflow, and bypass, thereby explicitly integrating spatial structure and temporal propagation into the monitoring model.

[0050] This invention achieves online joint estimation of the global water quality latent state based on time-delay observation equations and weighted optimization, and simultaneously outputs uncertainty, providing a quantitative basis for global situational awareness. By combining residual standardized anomaly scores with time-delay spatial maps to conduct reverse source tracing, it outputs a confidence ranking of suspected sources and prediction of propagation paths, enabling not only the detection of anomalies but also the location of their possible sources and propagation directions, improving the targeting and interpretability of responses. By constructing information value through risk and uncertainty and adaptively adjusting the sampling frequency, high-value locations are prioritized under budget constraints, achieving intelligent scheduling driven by information value, improving resource utilization efficiency and early warning capabilities. Attached Figure Description

[0051] Figure 1 This is a flowchart illustrating an online wastewater quality monitoring method according to the present invention.

[0052] Figure 2 This is an example diagram illustrating the global state estimation and uncertainty of an online wastewater quality monitoring method according to the present invention;

[0053] Figure 3 This is an example diagram illustrating the anomaly score and source confidence ranking list obtained from an online wastewater quality monitoring method of the present invention;

[0054] Figure 4 This is a module example diagram of an online wastewater quality monitoring system according to the present invention. Detailed Implementation

[0055] The subject matter described herein will now be discussed with reference to exemplary embodiments. It should be understood that these embodiments are discussed only to enable those skilled in the art to better understand and implement the subject matter described herein, and changes may be made to the function and arrangement of the elements discussed without departing from the scope of this specification. Various processes or components may be omitted, substituted, or added as needed in the examples. Furthermore, features described in some examples may be combined in other examples.

[0056] A method and system for online monitoring of wastewater quality includes the following embodiments:

[0057] Example 1:

[0058] A wastewater quality online monitoring method is proposed for use in the continuous operation environment of wastewater treatment plants or industrial park wastewater treatment stations. Several online monitoring points are deployed at the inlet, key branch inlets, main pool sections, and outlet. These monitoring points can utilize multi-parameter probes or single-parameter combinations. Basic hydraulic topology information is provided, including pipe connection relationships, pool volume, and the status of return valves and bypass valves. Data is transmitted back from the monitoring points via wired or wireless means, supporting lightweight edge computing. The method is as follows: Figure 1 As shown, it includes the following steps:

[0059] Step 100: Perform time alignment processing and data quality assessment on the original water quality sequences collected from each monitoring point to obtain the aligned data and data quality weights of the monitoring points.

[0060] The raw water quality sequences collected from each monitoring point are time-aligned and their data quality is assessed to provide a reliable data foundation for subsequent joint inference. In wastewater treatment systems, due to the independent clocks used by sensors at each monitoring point and the existence of latency and jitter in the communication network, the sampling times of different monitoring points often cannot be precisely synchronized, resulting in nominally identical observation data corresponding to different physical times. This time misalignment severely affects the accuracy of subsequent water quality propagation analysis based on time-delay spatial maps, as the propagation time of water quality changes is typically on the order of several minutes to tens of minutes, which is on the same time scale as clock deviation. Without time alignment, it is impossible to accurately identify the causal relationship and migration time lag between upstream and downstream monitoring points, leading to failure in anomaly tracing. At the same time, the data quality of field sensors is significantly affected by factors such as environmental interference, biological attachment, and electrical faults. If low-quality data is used indiscriminately in global state estimation, it will contaminate the overall inference results and reduce the reliability of anomaly detection. Therefore, step 100 eliminates the systematic impact of clock deviation through time alignment and assigns a confidence weight to each observation value through data quality assessment, thereby providing time-consistent and quality-controllable data input for subsequent steps.

[0061] The raw water quality sequence includes measurements of specific water quality indicators collected at each monitoring point at each time point. These indicators include pH, dissolved oxygen, oxidation-reduction potential (ORP), turbidity, conductivity, ammonia nitrogen, and chemical oxygen demand (COD). These indicators were chosen as inputs because they comprehensively characterize the physicochemical properties and pollution load status of wastewater. pH reflects the acid-base balance of the water body and directly affects biochemical reaction rates and microbial activity; dissolved oxygen is a key limiting factor for aerobic biochemical processes, and its concentration changes directly indicate aeration efficiency and the state of organic matter degradation; ORP characterizes the redox environment of the water body and can be used to determine anaerobic or anoxic conditions; turbidity reflects suspended solids concentration and is closely related to sludge settling performance and effluent quality; conductivity, as a conservative indicator, is unaffected by biochemical reactions and can be used as a tracer signal for water quality propagation; ammonia nitrogen and COD represent nitrogen pollution load and organic pollution load, respectively, and are core control indicators for wastewater treatment. The combination of these indicators covers the main physical, chemical, and biological processes in wastewater treatment, supporting comprehensive water quality status inference and anomaly identification.

[0062] The monitoring points include metadata, which includes the spatial location, cell, sensor type, sampling period, and communication quality characterization parameters, such as packet loss rate and clock skew. This metadata is introduced to support accurate calculations for time alignment and data quality assessment. Spatial location and cell information determine the positional relationship of the monitoring point within the hydraulic topology, providing spatial constraints for subsequent construction of the time-delay spatial map. Sensor type information identifies the response delay characteristics and measurement accuracy differences of different sensors. Electrochemical sensors, such as dissolved oxygen probes, typically have response times ranging from several seconds to tens of seconds, while chemical analytical sensors, such as ammonia nitrogen analyzers, can have response times of several minutes; these differences need to be compensated for in the construction of the time-delay observation matrix. Sampling period information is used for time index calculation in the case of asynchronous sampling of multiple indicators. The packet loss rate in the communication quality characterization parameters directly reflects data integrity, while clock skew is the core estimation target for time alignment.

[0063] Time alignment processing involves estimating the clock offset of each monitoring point and obtaining aligned data through resampling. Specifically, this includes estimating the clock offset for each monitoring point. The clock offset is the time difference between the sampling time of that monitoring point and the system reference time. Clock offset estimation is performed using either cross-correlation alignment or alignment based on common external factors. These two methods were chosen because they are suitable for different data characteristics and topological scenarios. The cross-correlation alignment method calculates the cross-correlation function between different monitoring point sequences and finds the time offset that maximizes the cross-correlation value as the clock offset. This method utilizes the similarity of water quality changes between upstream and downstream monitoring points and is suitable for monitoring point pairs with strong hydraulic connectivity and obvious water quality change patterns. Its advantage is that it can directly extract the time relationship from the water quality data without an additional reference signal. However, in cases of weak water quality changes or severe multipath mixing, the cross-correlation peak may be insignificant, leading to a decrease in estimation accuracy.

[0064] Alignment methods based on common external factors utilize systematic changes such as flow fluctuations as reference signals to identify the time offset of each monitoring point relative to this reference signal. The physical basis of this method is that flow fluctuations synchronously affect the hydraulic retention time and dilution effect of all units in the plant, thus leaving time markers in the water quality series of each monitoring point. Clock deviation can be calculated by detecting the time differences in the water quality response to flow fluctuations at each monitoring point. Alignment methods based on common external factors are suitable for scenarios with significant flow fluctuations and clear responses of each monitoring point to flow changes. Their advantage lies in their weak dependence on the type of water quality indicator, but they require additional auxiliary information such as flow rate. In practical applications, a suitable alignment method can be selected based on data availability and topological characteristics, or the results of two methods can be fused to improve robustness. Assuming the measured value of the original water quality indicator at a certain monitoring point at a certain moment is known, and this monitoring point corresponds to a specific water quality indicator, after estimating the clock deviation of this monitoring point at that moment using the above method, the measured value of the original water quality indicator is shifted over time according to the clock deviation to obtain the aligned data. Due to the discrete sampling characteristics, the measured values ​​of the aligned data cannot usually be directly obtained from the original sequence and need to be obtained through resampling methods. Resampling methods include linear interpolation, spline interpolation, or zero-order preservation. Choosing a suitable resampling method is crucial for preserving the temporal characteristics of the water quality signal. Linear interpolation obtains the aligned data by linearly interpolating between the two nearest original sampling points on either side of the target time after the time offset. This method assumes that the water quality index changes linearly between adjacent sampling points, has low computational complexity, and does not introduce values ​​outside the original data range. It is suitable for routine operating conditions with high sampling frequencies and relatively gentle water quality changes. Spline interpolation obtains the aligned data by constructing a piecewise polynomial function to fit the original sequence and then evaluating it at the target time. This method can maintain the high-order continuity of the interpolation results and more accurately capture the curvature characteristics of water quality changes. It is suitable for situations with long sampling periods or where the water quality index has a significant nonlinear trend, but it has higher computational complexity and may cause overfitting when there is significant data noise. The zero-order hold method directly takes the nearest original sampled value before the target time after the time offset as the measured value of the aligned data. This method is equivalent to the step signal assumption and is suitable for scenarios where water quality indicators change abruptly, such as valve switching or batch dosing. Its advantage is that it completely retains the original measured value without introducing interpolation error, but it will cause the aligned data to appear step-like and lose temporal continuity. For water quality monitoring applications, linear interpolation is preferred because it is simple to calculate and can better maintain the continuity of the signal, while the interpolation error is within an acceptable range under typical sampling periods. For cases with long sampling periods or drastic signal changes, spline interpolation can be used to improve alignment accuracy, but care should be taken to avoid overfitting by controlling the smoothing parameter.

[0065] Data quality assessment includes calculating data quality weights to quantify the reliability of data from a monitoring point. The data quality weights range from 0 to 1 and are determined by a combination of packet loss rate and jump detection results. Packet loss rate and jump detection are used as the two core dimensions of data quality assessment because they reflect the key quality attributes of data integrity and data consistency, respectively. Packet loss rate characterizes the reliability of the communication link and the continuity of sensor operation. A high packet loss rate means that there are a large number of missing values ​​in the observation sequence. These missing values ​​weaken the contribution of the monitoring point to the global state estimation and increase the uncertainty of the time delay estimation. Jump detection is used to identify non-physical anomalies caused by sensor failure, electrical interference, or calibration drift. If these anomalies are included in the calculation without identification, they will severely distort the state estimation results and trigger false alarms. The packet loss rate and jump intensity are mapped to quality weights via an exponential function, rather than a simple linear mapping, because the exponential function achieves a non-linear penalty effect for quality degradation. That is, when the packet loss rate or jump intensity is low, the quality weight decreases slowly to maintain tolerance for minor quality issues. However, when the packet loss rate or jump intensity exceeds a certain threshold, the quality weight rapidly decays to near zero, effectively isolating severely poor-quality data. The multiplication of two exponential functions ensures that only data that simultaneously meets the criteria of low packet loss rate and low jump intensity receives high-quality weights, reflecting the rigor of data quality assessment.

[0066] The data quality weight is calculated as follows: The packet loss rate is multiplied by a preset packet loss rate attenuation coefficient, and the negative value is used as the exponent of the first exponential function. The standardized jump intensity is multiplied by a preset jump attenuation coefficient, and the negative value is used as the exponent of the second exponential function. The two exponential functions are then multiplied together to obtain the data quality weight. The bases of the first and second exponential functions are natural constants. The packet loss rate is the ratio of the number of lost data packets to the total number of data packets at the monitoring point within a preset packet loss rate statistical window. The preset packet loss rate statistical window length ranges from 10 to 100 sampling periods, with a default value of 50 sampling periods. A longer window length indicates a stronger smoothing effect on short-term packet loss fluctuations and can be determined based on the typical fluctuation cycle of communication quality. The preset packet loss rate attenuation coefficient ranges from 5 to 20, with a default value of 10. A larger coefficient indicates a more sensitive impact of the packet loss rate on the quality weight and can be determined based on the statistical relationship between packet loss rate and data reliability in historical data. The standardized jump intensity is the jump indication at a monitoring point at a certain moment. It is calculated as follows: First, the absolute value of the difference between the current aligned data and the mean of the aligned data within the local sliding window is calculated as the original jump amplitude. Then, the original jump amplitude is standardized by dividing it by the median absolute deviation or standard deviation of the water quality indicator in the historical aligned data for that monitoring point, yielding the standardized jump intensity. The local sliding window length ranges from 5 to 30 sampling periods, with a default of 10 sampling periods. A longer window length indicates more stable tracking of long-term trends but a slower response to sudden changes; it can be determined based on the typical rate of change of the water quality indicator. Standardization ensures the comparability of jump intensities for different water quality indicators, avoiding deviations in data quality weight calculation due to differences in indicator dimensions. The preset jump attenuation coefficient ranges from 0.1 to 5, with a default of 1. A larger coefficient indicates a more severe penalty for data quality weighting by the jump; it can be determined based on the statistical distribution of normal fluctuation amplitudes and abnormal jump amplitudes. The exponential function ensures that the data quality weight is close to 1 when the packet loss rate is 0 and there are no jumps, while the data quality weight smoothly decays to close to zero when the packet loss rate or the intensity of the standardized jump increases.

[0067] Step 200: Based on the aligned data, construct a time-delay space graph, which includes a set of nodes, a set of edges, and the migration delay and attenuation coefficient of each directed edge.

[0068] Based on the aligned data output from step 100, a time-delay spatial graph reflecting the spatial structure and water quality migration patterns of the wastewater treatment system is constructed, treating the entire plant area as a dynamic graph rather than isolated monitoring points. The time-delay spatial graph is used as the core data structure for system modeling because the wastewater treatment system is essentially a distributed dynamic system with complex spatiotemporal coupling characteristics. Traditional single-point monitoring methods treat each monitoring point as an independent entity, ignoring the spatial propagation process and temporal delay effects of water quality changes, leading to an inability to accurately track pollutant migration paths and identify anomaly sources. The time-delay spatial graph explicitly expresses the hydraulic connectivity and mixing relationships between monitoring points through a graph structure. It quantifies the time required for water quality changes to propagate from upstream to downstream using the migration delay parameter of directed edges, and characterizes the concentration decay of water quality indicators during propagation due to biochemical degradation, dilution, or adsorption using the attenuation coefficient parameter. This graph representation method unifies physical topology, hydraulic characteristics, and biochemical processes into a mathematical framework, providing structured prior knowledge for subsequent global state estimation, anomaly tracing, and sampling scheduling. Compared to purely data-driven black-box models, time-delay spatial graphs have advantages such as strong interpretability, clear physical meaning of parameters, and good adaptability to changes in operating conditions. At the same time, the sparsity of its graph structure is also conducive to efficient computation of large-scale systems.

[0069] The information acquired in this step includes hydraulic topology information, as well as operational and hydraulic information. Specifically, the hydraulic topology information includes pipe connection relationships, tank volume, and the status of return valves and bypass valves. The operational and hydraulic information includes flow rate, valve status, and return ratio. The selection of these inputs is based on the hydraulic and process principles of the wastewater treatment system. Pipe connection relationships define the reachable paths of water flow and form the topological basis for constructing the directed edge set. Tank volume determines the residence time of water in each unit and is a key parameter for estimating the initial value of migration delay. The status of return valves and bypass valves reflects the actual process configuration during operation; different valve statuses alter the flow path and mixing mode, requiring dynamic adjustment of the effectiveness of directed edges in the time-delay space diagram. Flow rate information is used to calculate the actual hydraulic residence time of each unit and the flow rate proportion at the convergence of multiple paths, serving as an important basis for online correction of migration delay and multi-source contribution decomposition. The return ratio reflects the intensity of the impact of the return path on the main process; a high return ratio enhances upstream and downstream coupling and shortens the system's equivalent residence time. By integrating this hydraulic topology and operational information, the time-delay space diagram can dynamically reflect the actual operating state of the system, rather than relying solely on static design parameters.

[0070] First, a basic graph structure is constructed. Each monitoring point is used as a node, forming a set of nodes; hydraulic connectivity and co-pool mixing relationships are used as directed edges, forming a set of edges. The node and edge sets together constitute the time-delay spatial graph. Monitoring points are chosen as graph nodes rather than physical units because they are the actual locations where observation data is acquired. Using monitoring points as nodes directly establishes the correspondence between observed values ​​and the graph structure, simplifying the subsequent construction of the time-delay observation matrix. The distinction between hydraulic connectivity and co-pool mixing relationships is made because they represent different physical processes and time-delay characteristics. Hydraulic connectivity indicates a flow path between two monitoring points. This relationship is determined by the pipe connection, and the flow direction determines the direction of the directed edge. The migration time delay corresponding to this type of directed edge is usually relatively long, depending on the pipe length, flow velocity, and residence time in the intermediate pool. Water quality changes propagate along this type of directed edge with significant time delay and concentration decay. The co-processing relationship indicates that two monitoring points are located within the same reaction tank. This relationship is determined based on the spatial location and unit of the monitoring points. A bidirectional directed edge is established for each monitoring point within the same tank to reflect the bidirectional nature of the mixing. The migration time delay corresponding to this type of directed edge is usually short, mainly determined by the hydraulic mixing time within the tank. Under the assumption of ideal complete mixing, the water quality state at each point in the same tank should instantaneously tend to be consistent. However, in reality, short-circuit flow and dead zones exist in the tank, leading to incomplete mixing. Therefore, a limited migration time delay still needs to be retained. The setting of bidirectional directed edges ensures that water quality information between any two monitoring points within the same tank can be transmitted bidirectionally, consistent with the physical nature of the mixing process.

[0071] Secondly, learnable migration delays and attenuation coefficients are introduced for each directed edge. The migration delay represents the time required for a water quality change to propagate from an upstream node to a downstream node; the attenuation coefficient represents the degree of attenuation of the water quality index during propagation. The migration delay and attenuation coefficient are designed as learnable parameters rather than fixed constants because the hydraulic and biochemical characteristics of wastewater treatment systems are continuously changing due to various factors such as influent load, temperature, and sludge concentration; fixed parameters cannot adapt to this dynamic nature. By learning these parameters online from real-time observation data, the time-delay space graph can automatically track the evolution of the system state, maintaining consistency between the model and the actual process. The learnability of the migration delay enables the system to automatically identify changes in residence time caused by flow fluctuations, changes in effective tank volume, or short-circuit flow phenomena; the learnability of the attenuation coefficient enables the system to automatically adapt to changes in degradation rate caused by temperature changes, fluctuations in microbial activity, or the presence of inhibitory substances.

[0072] For a directed edge from an upstream node to a downstream node, the initial value of the migration delay can be roughly estimated by dividing the pool volume of the downstream unit by the representative flow rate of that unit. The representative flow rate can be taken as the average flow rate within the sliding window. This calculation method is based on an ideal continuous stirred reactor model, reflecting the average residence time of water within the pool section. The physical basis for using the pool volume divided by the flow rate as the initial value of the migration delay is the definition of hydraulic residence time. This method is simple and intuitive and does not require complex fluid dynamics calculations, making it suitable for online applications. Using the average flow rate within the sliding window instead of the instantaneous flow rate is to smooth out short-term fluctuations in flow rate and avoid noise interference in the initial value estimation. The ideal continuous stirred reactor model assumes that the water in the pool is instantaneously completely mixed and that the effluent quality is equal to the average water quality in the pool. Although this model simplifies the actual flow field distribution, its predicted average residence time basically matches the statistical characteristics of the actual system, and therefore can be used as a reasonable initial guess for the migration delay. Subsequently, through an online learning mechanism, the migration delay will be gradually corrected from this initial value to the true value. The initial value of the attenuation coefficient can be set according to the biochemical degradation characteristics of water quality indicators or historical experience. For conservative indicators such as conductivity, it can be set to a value close to 1, and for easily degradable indicators such as chemical oxygen demand, it can be set to a value between 0.7 and 0.9.

[0073] The initial attenuation coefficients are set based on the physicochemical properties of water quality indicators because different indicators exhibit fundamentally different behavioral mechanisms during wastewater treatment. Conductivity is primarily contributed by dissolved inorganic salts. These ions neither participate in biochemical reactions nor undergo phase transitions during conventional wastewater treatment; therefore, their concentration remains essentially constant during water quality propagation, only affected by dilution effects. An attenuation coefficient close to 1 reflects their conservative nature. Chemical oxygen demand (COD) represents the total amount of oxidizable organic matter and reducing substances. Under aerobic or anoxic conditions, it is degraded by microorganisms or chemically oxidized, causing its concentration to gradually decrease along the propagation path. An attenuation coefficient of 0.7 to 0.9 indicates that COD may decrease by 10% to 30% after one unit, covering typical degradation efficiencies under different temperature, sludge concentration, and retention time conditions. By setting initial attenuation coefficients that conform to the physicochemical properties of different indicators, the convergence speed of subsequent online learning can be accelerated, and estimation errors in the initial stage can be reduced. The attenuation coefficient is continuously corrected using real-time data. The correction method is as follows: within a preset attenuation coefficient update window, it is estimated using the statistical characteristics of the ratio of observed values ​​between upstream and downstream nodes. Specifically, the aligned data of the downstream node is divided by the aligned data of the upstream node after migration delay. Within the update window, the median or weighted average of the ratios at multiple times is taken as the updated value of the attenuation coefficient. The weights for the weighted average can be determined based on the data quality weights. The physical basis for estimating the attenuation coefficient using the observed value ratio statistical method is that, under stable operating conditions, the concentration ratio of upstream and downstream water quality indicators should be equal to the attenuation coefficient of the propagation path.

[0074] By statistically summarizing the ratios across multiple time points within the update window, the impact of noise from a single measurement can be smoothed out, and the true value of the attenuation coefficient can be extracted. The median or weighted average is chosen as the statistic to balance robustness and accuracy. The median is insensitive to outliers and is suitable for situations with occasional measurement errors or operational disturbances; the weighted average assigns greater influence to high-quality data based on data quality weights, making it suitable for situations with significant differences in data quality. Using data quality weights as weighting coefficients in the weighted average ensures that the contribution of low-quality data to the attenuation coefficient estimation is automatically suppressed, improving the reliability of parameter learning. To avoid unreasonable values ​​for the attenuation coefficient due to noise or abnormal operating conditions, physical constraints are imposed on the attenuation coefficient, requiring it to range between 0.5 and 1.2. When the calculated update value exceeds this range, it is truncated. The attenuation coefficient update cycle can be set to 2 to 5 times the migration lag update cycle to ensure the stability of parameter updates.

[0075] Migration lags are corrected online using real-time data. Within a sliding time window, migration lags are updated by maximizing alignment correlation or minimizing prediction error. Online correction, rather than offline calibration, is used to update migration lags because the hydraulic retention time of wastewater treatment systems is constantly changing due to flow fluctuations, variations in effective tank volume, and short-circuit flow phenomena; fixed parameters obtained through offline calibration cannot adapt to this dynamic nature. Online correction automatically tracks the changing trend of migration lags by continuously analyzing the temporal correlation patterns of upstream and downstream monitoring points in real-time observation data. The method of maximizing alignment correlation estimates migration lags by searching for the time offset that maximizes the correlation between the upstream delay signal and the downstream observation signal; this method is suitable for situations where water quality change patterns are obvious and noise levels are low.

[0076] The method of minimizing prediction error estimates migration delay by searching for the time offset that minimizes the mean square error between the upstream delay signal and the downstream observation signal. This method has low requirements for signal shape and is suitable for situations with gradual water quality changes or high noise levels. Using a sliding time window instead of full historical data for migration delay updates balances the statistical stability of parameter estimation with the ability to track time-varying characteristics; the choice of window length requires a trade-off between short-term noise fluctuations and long-term trend changes. To avoid time delay estimation bias in multi-path mixed scenarios, time delay constraints and multi-source decomposition mechanisms are introduced for each directed edge. In wastewater treatment systems, it is common for multiple upstream paths to converge at the same downstream node, such as multiple influent branches merging into the main treatment line or the mixing of return sludge with influent. In such multi-path mixed scenarios, the observation value of the downstream node is the superposition of the delay signals of each upstream path. If the downstream observation value is indiscriminately estimated with a single upstream path, the estimated migration delay will deviate from the true value due to the mixing effect. Time delay constraints limit the search range of migration delay within the physically feasible interval, preventing the optimization process from getting trapped in unreasonable local optima. The multi-source decomposition mechanism explicitly models the contribution weights of each upstream path, decomposing the mixed downstream observations into independent contributions from each upstream path. This provides each path with a demixed, effective signal for time delay estimation. The introduction of this mechanism significantly improves the accuracy and robustness of migration time delay estimation in multi-path scenarios. Specifically, for a directed edge, the physical feasible interval of the migration time delay is first determined based on the pool volume range and flow variation range of the corresponding unit. The lower limit of the physical feasible interval is the minimum pool volume of the unit divided by the historical maximum flow, and the upper limit is the maximum pool volume of the unit divided by the historical minimum flow. The migration time delay search range is limited to this physical feasible interval. Secondly, for downstream nodes with multiple upstream inbound directed edges, the multi-source contribution decomposition method is used to handle the mixing effect.

[0077] The multi-source contribution decomposition method constructs a linear superposition model, representing the downstream node's observations as a weighted sum of the delayed signals from each upstream node. The weights reflect the flow proportion or mixed contribution of each upstream path, and can be calculated from real-time flow data or identified from the observation data using sparse optimization methods. Based on this, for each upstream inbound directed edge, the migration delay value that minimizes the residual is searched within its physically feasible interval. During the search, the migration delay parameters of other directed edges are fixed, and the method is updated edge-by-edge using coordinate descent or alternating optimization. Furthermore, for directed edges affected by backflow or bypass, the effectiveness weight of the directed edge is dynamically adjusted based on the status of the backflow valve and bypass valve and the backflow ratio. When the backflow ratio is large or the bypass is open, the weight of the directed edge in the time delay learning is reduced to avoid noise interference from non-dominant paths. Through the above constraints and decomposition mechanisms, the migration delay update process can maintain stability and physical rationality under multi-path mixing and operating condition disturbances.

[0078] Step 300: Based on the aligned data, data quality weights, and time-delay spatial map, a joint estimation of the global water quality state is performed to obtain the global state estimate and uncertainty, as detailed below. Figure 2 As shown.

[0079] Based on the aligned data and data quality weights output in step 100 and the time-delay spatial graph output in step 200, a joint estimation of the global water quality status is performed. Global joint estimation is used instead of independent single-point estimation because the observations at each monitoring point in the wastewater treatment system are not independent but tightly coupled in time and space through hydraulic connectivity and mixing relationships. Single-point estimation methods directly treat the observations at each monitoring point as the true water quality status at that location, ignoring the influence of observation noise and sensor errors, and failing to utilize information from adjacent monitoring points for complementary correction. Global joint estimation introduces a global hidden state vector to represent the true water quality status of each region in the system. The observations at all monitoring points are considered as the result of this global hidden state vector after time-delay propagation and the superposition of observation noise. The optimal estimate of the global hidden state vector is inferred from the noisy observations by solving an optimization problem. The advantage of this method is that it can integrate information from multiple monitoring points for cross-validation and error smoothing. For monitoring points with high observation noise or missing data, compensation can be inferred from the observations of adjacent monitoring points using the topological relationships in the time-delay spatial graph, thereby improving the overall estimation accuracy and robustness. Meanwhile, global joint estimation can naturally handle the time delay propagation effect, establish an explicit connection between the state at historical moments and the observation at the current moment, and avoid the problem that the state transition model is difficult to accurately model in traditional filtering methods.

[0080] A global hidden state vector is established to represent the true water quality state of all units or regions at a given moment. The dimension of the global hidden state vector can be determined by dividing the data into regions, with each region corresponding to one component of the global hidden state vector. The concept of the global hidden state vector is introduced to distinguish between the true physical state of the system and sensor observations. In actual systems, sensor observations are affected by measurement noise, response delay, and installation location, and cannot completely and accurately reflect the average water quality state of the region. The global hidden state vector represents the true water quality state of each region under ideal noise-free conditions and is the target variable for joint estimation. The global hidden state vector is defined by region rather than by monitoring point because multiple monitoring points may be deployed within the same region. These monitoring points observe different spatial locations of the same water body, and their true states should be consistent or highly correlated.

[0081] By using regions as the basic units of state, the dimensionality of the global hidden state vector can be reduced, thereby decreasing the computational complexity of the optimization problem. This also aligns with the practical reality of wastewater treatment processes where tanks or units are the basic control objects. The relationship between the aligned data and the global hidden state vector is described by the time-delay observation equation. The time-delay observation equation represents the aligned data at the current moment as a function of the global hidden state vector at historical moments. Specifically, the aligned data at the current moment equals the time-delay observation matrix multiplied by the extended state vector, which includes global hidden state vectors from multiple historical moments, plus the observation noise term. The use of the time-delay observation equation instead of the instantaneous observation equation is to explicitly express the time delay effect of water quality propagation. In a wastewater treatment system, the current observation value at a monitoring point depends not only on the actual state at that location but also on the state from upstream historical moments after migration and time delay propagation. The time-delay observation equation incorporates global hidden state vectors from multiple historical moments into the observation model by introducing the extended state vector, allowing the time-delay propagation relationship to be uniformly handled within the state estimation framework. The time-delay observation matrix, serving as a mapping matrix connecting historical states and current observations, is directly determined by the time-delay space graph constructed in step 200, achieving seamless integration of topological information and the estimation model. The introduction of the observation noise term reflects the combined impact of sensor measurement errors and model simplification errors; its statistical characteristics can be determined based on sensor specifications or historical calibration data, providing a foundation for subsequent uncertainty quantification. Specifically, the extended state vector is formed by stacking the global hidden state vectors of the current moment and several historical moments in chronological order. The number of stacked historical moments is determined by rounding up the ratio of the system's maximum migration delay to the sampling period. The observation noise term represents the measurement error, and its statistical characteristics can be determined based on sensor specifications or historical calibration data.

[0082] The time-delay observation matrix is ​​a crucial mapping matrix connecting the observations of monitoring points to the global hidden state vector. Its construction process comprehensively considers the topological structure of the time-delay spatial graph, the migration delay parameters of each directed edge, the attenuation coefficient of each directed edge, and the mapping relationship from monitoring points to regions. The number of rows in the time-delay observation matrix equals the total number of monitoring points, and the number of columns equals the dimension of the extended state vector (i.e., the number of regions multiplied by the number of stacked historical moments). Each row of the time-delay observation matrix corresponds to a monitoring point, describing how the observations at that monitoring point are linearly combined from the states of each region at each historical moment in the extended state vector.

[0083] The method for obtaining the time-delay observation matrix is ​​as follows: For a given monitoring point, first determine the region where the monitoring point is located, denoted as the target region; then search for water quality propagation paths from all other regions to the target region in the time-delay space graph. The search for propagation paths needs to consider the connectivity and flow direction constraints of directed edges; for each propagation path found, calculate the cumulative migration delay and cumulative attenuation coefficient of the path. The cumulative migration delay is the sum of the migration delays of each directed edge on the path, and the cumulative attenuation coefficient is the product of the attenuation coefficients of each directed edge on the path; determine the historical time index corresponding to the path based on the cumulative migration delay. The historical time index is obtained by dividing the cumulative migration delay by the sampling period and rounding down. If the cumulative migration delay exceeds the maximum historical duration covered by the extended state vector, then... This path does not contribute. Based on the path's starting region, historical time index, and cumulative attenuation coefficient, the cumulative attenuation coefficient is filled into the corresponding column position of the row corresponding to the monitoring point in the time-delay observation matrix as a weighting coefficient. If multiple paths point to the same region at the same historical time, the cumulative attenuation coefficients of these paths are weighted and summed according to their flow proportions and then filled into the corresponding column position. For the target region itself, if the monitoring point is directly located within the region, the value 1 is filled into the current time target region column position of the row corresponding to the monitoring point in the time-delay observation matrix, indicating that the monitoring point directly observes the current state of the region. For elements in the time-delay observation matrix that are not assigned values ​​by the above process, their values ​​are zero, indicating that the corresponding historical time region state has no direct contribution to the monitoring point's observation value. Through the above construction process, the position of the non-zero element in each row of the time-delay observation matrix reflects which historical time and which region states the monitoring point's observation value is related to. The value of the non-zero element, i.e., the weighting coefficient, reflects the strength of the correlation. The weighting coefficient is obtained by multiplying the attenuation coefficients of each directed edge on the propagation path. The time-delay observation matrix explicitly embeds the migration time delay parameters and attenuation coefficient parameters learned in step 200 into the mathematical model of global state estimation, achieving model self-consistency between the time-delay space graph and the joint estimation.

[0084] For single water quality indicators, indicator labels can be omitted during analysis. For multiple indicators, the differences in sampling periods and sensor response delays among the different indicators need to be considered. For indicators with different sampling periods, when constructing the extended state vector, the historical states of each indicator are stacked according to its own sampling period, and the time index is adjusted accordingly in the time-delay observation matrix. Specifically, the cumulative migration delay is divided by the sampling period of the indicator rather than the unified sampling period of the system for calculating the historical time index. For chemical analysis sensors with significant response delays, an inherent delay parameter of the sensor is introduced into the observation equation corresponding to the sensor in the time-delay observation matrix. This parameter can be determined by the sensor's technical specifications or obtained through step response experiments. Specifically, the inherent delay parameter of the sensor is added when calculating the cumulative migration delay of the monitoring point corresponding to the sensor, thereby shifting the position of the non-zero element in the row corresponding to the monitoring point in the time-delay observation matrix to an earlier historical time. Through the above processing, the joint estimation in the case of multiple indicators can correctly handle the problems of asynchronous sampling and response delay, and can be processed by either independent estimation of each indicator or unified estimation of multiple indicators stacked together.

[0085] A weighted optimization method is employed for online estimation of the global hidden state. The estimated value is obtained by solving a weighted least squares optimization problem, using the extended state vector as the optimization variable. The objective function consists of the sum of an observation fitting term and a smoothing regularization term. Weighted least squares is chosen as the optimization framework because it is optimal under the Gaussian noise assumption and has high computational efficiency, making it suitable for online applications. By minimizing the weighted sum of squares of the observation fitting error, weighted least squares ensures that the estimated global hidden state vector can interpret the actual observation data to the greatest extent. The introduction of weights allows for adjusting the influence of different monitoring points on the estimation results based on data quality weights, achieving data quality-aware state estimation. The smoothing regularization term is introduced into the objective function to constrain the estimation results using prior knowledge of the spatial continuity of the wastewater treatment system. Physically, the water quality state of adjacent areas typically does not exhibit drastic changes due to hydraulic mixing and diffusion. The smoothing regularization term penalizes solutions with excessively large state differences between adjacent areas, guiding the optimization process towards spatial smoothness. This regularization strategy is particularly important when observation points are sparse or data quality is poor, as it can prevent non-physical oscillations or singular values ​​from appearing in the estimation results. The combination of the observation fitting term and the smoothing regularization term achieves a balance between data-driven and physical constraints, ensuring both the fidelity of the estimation results to the observation data and maintaining the physical rationality of the results. The observation fitting term is the weighted sum of squared observation fitting errors for all monitoring points. The observation fitting error for each monitoring point is the difference between the aligned data and the corresponding predicted observation. The predicted observation is obtained by multiplying the time-delay observation matrix and the extended state vector, with the weights being the data quality weights obtained in step 100. The smoothing regularization term is the square of the norm of the Laplacian matrix of the time-delay spatial graph applied to the global hidden state vector at the current moment. The resulting value is multiplied by a preset smoothing regularization coefficient, which is used to constrain the smoothness of the state in adjacent regions but allows abrupt changes in abnormal situations. The preset smoothing regularization coefficient ranges from 0.01 to 10, with a default value of 0.1. The larger the coefficient, the stronger the constraint on spatial smoothness. The coefficient value that minimizes the prediction error can be determined by selecting the coefficient value from historical data using cross-validation.

[0086] The method for constructing the Laplace matrix of the time-delay spatial graph comprehensively considers the type, directionality, and coupling strength of directed edges. Specifically, for each directed edge in the time-delay spatial graph, different coupling weights are assigned according to the type of the directed edge. For directed edges with mixed relationships within the same pool, the coupling weight is set to a larger value, ranging from 0.5 to 1, with a default value of 1, reflecting that the status of each monitoring point within the same pool should be highly consistent. For directed edges with hydraulic connectivity, the coupling weight is set according to the flow ratio or propagation attenuation coefficient of the directed edge, ranging from 0.1 to 0.5, with a default value of 0.3, reflecting that the coupling strength on the transmission path is weaker than that of the mixed relationship. For directed edges affected by backflow or bypass, the coupling weight is dynamically adjusted according to the status of the backflow valve and bypass valve and the backflow ratio. When the backflow ratio increases, the coupling weight of the reverse directed edge is increased accordingly. The diagonal elements of the Laplace matrix of the time-delay space graph are equal to the sum of the coupling weights of all adjacent directed edges of that node. For the off-diagonal elements, if two nodes are connected by a directed edge, the negative value of that directed edge's coupling weight is taken; otherwise, it is 0. For asymmetric connections in the time-delay space graph, symmetry processing is used, or only the coupling weights of strong connection directions are considered. Through the above weighting and classification processing, the smoothing regularization term can more accurately reflect the actual coupling relationships between different units in the wastewater treatment system, avoiding excessive constraints on stratified or short-circuit flow regions that should not be smoothed. The optimal solution to the above optimization problem is the estimated value of the extended state vector. Extracting the state component at the current moment from it yields the estimated value of the global hidden state, and the estimated value of the global hidden state vector is the global state estimate.

[0087] The uncertainty of the global hidden state estimate is calculated synchronously to guide subsequent sampling scheduling. Uncertainty quantifies the reliability of the global hidden state vector estimate, reflecting the statistical distribution characteristics of state estimation errors caused by observation noise, missing data, and model errors. Uncertainty quantification is introduced in global state estimation because the estimate itself only provides a point estimate of the state and cannot reflect its reliability. In practical applications, the reliability of state estimates varies significantly across different regions due to differences in monitoring point density, data quality, and topological connectivity. Uncertainty quantification allows the system to identify which regions have more reliable state estimates and which regions have higher estimation uncertainties due to insufficient observation information or poor data quality. This uncertainty information is crucial for subsequent anomaly detection and sampling scheduling. In anomaly detection, uncertainty is used to standardize anomaly scores, avoiding misclassifying normal estimation errors in high-uncertainty regions as anomalies. In sampling scheduling, uncertainty is used to identify monitoring points with high information value, prioritizing sampling resources for high-uncertainty regions to reduce estimation errors. The synchronous calculation of uncertainty, rather than post-evaluation, ensures consistency between uncertainty information and state estimates, avoiding information mismatch caused by time delays. Uncertainty is represented by an uncertainty matrix, which is a square matrix whose dimension is equal to the dimension of the global hidden state vector at the current time, i.e., the number of regions. The diagonal elements of the uncertainty matrix represent the variance of the state estimates for each region, and the square root of the diagonal elements represents the marginal uncertainty of the corresponding region. The larger the value of the marginal uncertainty, the higher the uncertainty of the state estimate for that region. The off-diagonal elements of the uncertainty matrix represent the covariance between the state estimate errors of different regions, reflecting the correlation of estimation errors caused by shared observation information or topological coupling.

[0088] Methods for obtaining the uncertainty matrix include analytical methods based on covariance propagation and approximate methods based on residual statistics. The analytical method based on covariance propagation is derived from the mathematical properties of the weighted least squares optimization problem. Specifically, the weighted least squares optimization problem in step 300 is considered as a maximum a posteriori estimation problem under the Gaussian noise assumption. According to Bayesian inference theory, the optimal solution to the optimization problem corresponds to the mean of the posterior distribution, and the covariance matrix of the posterior distribution is the uncertainty matrix. The uncertainty matrix can be obtained by inverting the second derivative matrix of the objective function of the optimization problem with respect to the extended state vector, i.e., the Hessian matrix. The Hessian matrix is ​​calculated by multiplying the transpose of the time-delay observation matrix by a diagonal weight matrix formed by data quality weights, then multiplying by the time-delay observation matrix, and then adding the regularization matrix corresponding to the smoothing regularization term. The regularization matrix is ​​composed of the product of the Laplace matrix of the time-delay space graph and the preset smoothing regularization coefficients. After inverting the Hessian matrix, the covariance matrix of the extended state vector is obtained, and the submatrix corresponding to the global hidden state vector at the current time is extracted from it to obtain the uncertainty matrix. Analytical methods based on covariance propagation can accurately reflect the combined impact of observation configuration, data quality weights, and topology on uncertainty, but the computational cost of matrix inversion is high in large-scale systems.

[0089] The residual statistics-based approximation method estimates uncertainty by analyzing the residual distribution after solving the optimization problem. The specific process is as follows: Calculate the observation fitting error for all monitoring points, which is the difference between the aligned data and the estimated value of the time-delay observation matrix multiplied by the extended state vector. Estimate the variance of the observation noise based on the observation fitting error and data quality weights. The estimated variance of the observation noise can be taken as the weighted average of the products of the squares of the observation fitting errors for all monitoring points and the inverses of the data quality weights, with the weights being the data quality weights. Construct a simplified covariance propagation model using the estimated observation noise variance. Multiply the transpose of the time-delay observation matrix by the diagonal matrix formed by the observation noise variance and the inverses of the data quality weights, then multiply by the time-delay observation matrix again, and add a regularization matrix. Invert the simplified matrix and extract the submatrix corresponding to the current time step as an approximation of the uncertainty matrix. The residual statistics-based approximation method has high computational efficiency and can adaptively reflect the actual observation noise level, but its approximation accuracy may decrease when the observation points are sparse or the data quality varies greatly.

[0090] In practical applications, an appropriate uncertainty calculation method can be selected based on the system size and computational resources. For small systems with a limited number of monitoring points, an analytical method based on covariance propagation is preferred to obtain an accurate uncertainty matrix. For large systems with a large number of monitoring points, an approximate method based on residual statistics or a sparse matrix inversion algorithm for the Hessian matrix can be used to reduce computational complexity. After obtaining the uncertainty matrix, a rationality check is required. This check includes verifying the nonnegativity of diagonal elements and the positive definiteness of the matrix. If the check fails, the regularization parameters of the optimization problem need to be adjusted, or a matrix inversion method with better numerical stability needs to be adopted.

[0091] Step 400: Based on the global state estimate, uncertainty, and time delay space graph, anomaly detection is performed, and the propagation paths and suspected sources of anomalies are identified, resulting in anomaly scores and a source confidence ranking list, as detailed below. Figure 3 As shown.

[0092] Based on the global state estimate and uncertainty output in step 300, and combined with the time-delay spatial diagram from step 200, anomaly detection is performed to identify the propagation path and suspected source of the anomaly. Anomaly detection and source tracing analysis are integrated into the same step because they are logically closely related and share the same data foundation. Traditional anomaly detection methods only determine whether the observed value at a monitoring point exceeds the normal range, but cannot answer the origin location and propagation path of the anomaly, leading to maintenance personnel spending a significant amount of time on manual investigation. This invention, by combining the time-delay spatial diagram for reverse source tracing, can not only identify which monitoring points are abnormal, but also trace the most likely source location of the anomaly and predict its propagation direction, providing decision support for quickly locating pollution sources and formulating emergency response measures. The use of standardized anomaly scores based on regional residuals instead of simple threshold judgment is to consider the impact of uncertainty on anomaly judgment. In areas with high uncertainty, even a certain deviation between the observed and predicted values ​​may be a normal estimation error and should not be judged as an anomaly; while in areas with low uncertainty, a small deviation may indicate a real anomaly event. By standardizing the regional residuals by dividing them by the marginal uncertainty, the anomaly score can automatically adapt to the differences in estimation accuracy in different regions, thus achieving a unified anomaly judgment standard.

[0093] The regional residual is calculated for each monitoring point, representing the difference between the actual and predicted observations. Specifically, the regional residual equals the aligned data minus the product of the time-delay observation matrix and the global state estimate. The regional residual is the fundamental signal for anomaly detection; its physical meaning is the deviation between the actual observation and the model prediction, given the global state estimate and the time-delay spatial map. Under normal operating conditions, if the global state estimate is accurate and the time-delay spatial map parameters are correct, the regional residual should primarily be contributed by observation noise, and its statistical characteristics follow a zero-mean random distribution. When an anomaly occurs at a monitoring point, the actual water quality state at that location deviates from the normal propagation pattern, causing a significant deviation between the observed value and the predicted value based on the normal model, resulting in a significant increase in the amplitude of the regional residual. By monitoring the statistical characteristics of the regional residual, abnormal events deviating from the normal pattern can be detected in a timely manner. Using the regional residual instead of directly using the observed values ​​for anomaly detection is chosen because the regional residual has already deducted the normal time-delay propagation effect and spatial coupling influence, thus reflecting the local anomaly signal more purely and avoiding misjudging normal upstream changes as downstream anomalies.

[0094] Based on the construction of standardized anomaly scores using regional residuals, the anomaly score for a specific monitoring point at a given time is calculated as follows: the absolute value of the regional residual is divided by a marginal term, which is the sum of the marginal uncertainty of the monitoring point and a preset minimum positive number. Here, the regional residual is the component of the residual vector corresponding to the monitoring point; the marginal uncertainty of the monitoring point is obtained by projecting the uncertainty matrix output in step 300 onto the observation space and then taking the square root of the diagonal element corresponding to the monitoring point; the preset minimum positive number ranges from 10⁻⁶ to 10⁻³, with a default value of 10⁻⁴. This value is used to prevent the denominator from being zero and should be much smaller than the typical marginal uncertainty value of a monitoring point. This anomaly score is essentially a normalized form of the regional residuals. The larger the deviation between the observed and predicted values ​​or the smaller the marginal uncertainty of the monitoring point, the higher the anomaly score.

[0095] The system identifies a set of nodes with high anomaly scores, which includes all monitoring points whose anomaly scores exceed a preset anomaly score threshold. The preset anomaly score threshold ranges from 2 to 5, with a default value of 3. This threshold corresponds to the significance level under a standard normal distribution, and a threshold of 3 corresponds to a false alarm rate of approximately 0.3%. This threshold can be adjusted to balance the false alarm and false negative rates in practical applications.

[0096] Then, reverse tracing is performed using the time-delay spatial map. The goal of reverse tracing is to identify the most likely source location from all nodes in the time-delay spatial map that triggered the currently observed anomaly pattern, thus obtaining candidate source points. Reverse tracing is a key value-added function of anomaly detection; its value lies not only in answering where the anomaly occurred, but more importantly, in answering where the anomaly originated. In wastewater treatment systems, due to the propagation characteristics of water quality changes, anomalies detected at downstream monitoring points are often the result of upstream anomaly sources propagating over a period of time. If emergency responses are based solely on the downstream anomaly location, the optimal opportunity to control the source may be missed. Reverse tracing analyzes the spatiotemporal distribution pattern of anomalies in the time-delay spatial map and uses migration delay parameters for reverse reasoning to identify the most likely origin location of the anomaly. This source tracing analysis provides maintenance personnel with the causal chain of anomaly events, supporting rapid location of pollution sources, tracing of responsible units, and development of targeted emergency measures. Source tracing using the time-delay spatial map, rather than simple topological backtracking, is chosen because the time-delay spatial map contains dynamic parameters such as migration delay and attenuation coefficients, enabling quantitative time-matching analysis and significantly improving the accuracy of source tracing. Candidate source points refer to all nodes in the time-delay space graph that could potentially be the origin of anomalies. Reverse tracing includes score tracing and time tracing. Score tracing is defined as the node's anomaly score being significantly higher than the anomaly scores of its downstream nodes connected by directed edges in the time-delay space graph. Specifically, the difference between the node's anomaly score and the downstream node's anomaly score exceeds a preset anomaly score difference threshold. The preset threshold ranges from 0.5 to 2, with a default value of 1. A larger threshold indicates a stricter requirement for the significance of the source node's anomaly, and it can be determined based on the statistical distribution of source and downstream anomaly scores in historical anomaly events. Time tracing is defined as the node's anomaly occurrence time being earlier than the anomaly occurrence time of its downstream nodes. The determination of the temporal order needs to consider the migration time delay from the node to the downstream node. If the node's anomaly occurrence time plus the migration time delay is earlier than or close to the downstream node's anomaly occurrence time, then the temporal order is considered valid. The criterion for closeness is that the absolute value of the time difference is less than a preset time tolerance. The preset time tolerance ranges from 1 to 3 times the migration time delay uncertainty, with a default value of 2 times. Nodes that satisfy either the score-based or time-based source tracing criteria will be included in the candidate source point set. The construction of the candidate source point set comprehensively considers topological location, temporal sequence, and historical risk information to ensure that the source tracing analysis can cover all reasonable anomaly origin hypotheses.

[0097] For each candidate source point, its explanatory power for the set of high-anomaly-score nodes is calculated. The calculation of explanatory power requires clarifying the path selection method from the candidate source point to each high-anomaly-score node and the propagation mechanism of migration time delay uncertainty. The path selection method employs either the time-delay shortest path algorithm or a multi-path comprehensive evaluation method. The time-delay shortest path algorithm searches for the path with the minimum cumulative migration time delay from the candidate source point to the target node in the time-delay space graph. During the path search, valve status and flow direction constraints must be considered, retaining only hydraulically reachable paths under the current operating conditions. The multi-path comprehensive evaluation method searches for the top few time-delay shortest paths from the candidate source point to the target node. The number of paths ranges from 1 to 5, with a default value of 3. The explanatory power contribution of each path is calculated separately and then weighted and summed. The weights are related to the path's flow rate proportion or reliability. For each selected propagation path, the explanatory power is calculated as follows: for each node in the set of high-anomaly-score nodes, the product of the time matching factor, the path uncertainty discount factor, and the node's anomaly score is calculated, and then the product results of all nodes are summed.

[0098] The time matching factor is an exponential decay function, with the base of the exponential decay function being the natural constant. The time term with a negative exponent is divided by the preset time matching decay scale. The time term is the current time minus the time when the node anomaly occurs, plus the cumulative migration delay from the candidate source point to the node, and the absolute value of the result is taken. The cumulative migration delay is the sum of the migration delays of each directed edge along the selected path, provided by the time delay space graph maintained in step 200. The preset time matching decay scale ranges from 5 minutes to 30 minutes, with a default value of 10 minutes. The smaller the time matching decay scale, the stricter the time matching requirement, which can be determined based on the uncertainty range of the typical hydraulic residence time of the system. The path uncertainty discount factor reflects the impact of migration delay estimation uncertainty on the reliability of source tracing. It is calculated by summing the migration delay uncertainties of each directed edge on the path according to the variance propagation law to obtain the total migration delay uncertainty of the path. Then, a discount factor is constructed that decreases as the total migration delay uncertainty increases. The discount factor can be obtained by dividing the total migration delay uncertainty by a preset uncertainty reference value, adding 1, and then taking the reciprocal. The preset uncertainty reference value ranges from 10% to 30% of the typical migration delay of the system, with a default value of 20%. The smaller the uncertainty reference value, the more severe the penalty for migration delay uncertainty. The migration delay uncertainty of each directed edge can be obtained from the residual statistics of the online migration delay correction process in step 200 or the confidence interval width of the migration delay search. By introducing the path uncertainty discount factor, source tracing can comprehensively consider the time matching degree and the reliability of migration delay estimation, avoiding misjudgments caused by inaccurate migration delay parameters. The physical meaning of this calculation is: if the time when an anomaly occurs at a certain node matches the time required for the anomaly to propagate from the candidate source point to that node, and the migration delay estimation of the propagation path is relatively reliable, then the candidate source point has a high explanatory power for the anomaly; after the explanatory power of all high anomaly score nodes is accumulated, the candidate source point with the highest explanatory power is most likely to be the real source of the anomaly.

[0099] Candidate source points are ranked according to their explanatory power to obtain a source confidence ranking list, with the first one being the most likely anomaly source. Simultaneously, based on the high-weight directed edges in the time-delay spatial graph and the time-delay matching, the predicted anomaly propagation direction and path are output.

[0100] Step 500: Based on the uncertainty, anomaly score and source confidence ranking list, formulate an adaptive sampling strategy and determine the sampling frequency of each monitoring point.

[0101] Based on the uncertainty output from step 300, the anomaly score output from step 400, and the source ranking, an adaptive sampling strategy is formulated to achieve intelligent scheduling based on information value rather than simply the rate of change. The adoption of an adaptive sampling strategy instead of a fixed sampling frequency aims to maximize the information acquisition efficiency of the monitoring system within the constraints of limited communication bandwidth and energy consumption budget. Traditional fixed sampling strategies use the same or preset sampling frequency for all monitoring points, failing to consider the differences in information value at different monitoring points at different times, leading to unreasonable resource allocation. When the system is operating smoothly, the water quality status at some monitoring points changes slowly and is highly predictable, resulting in limited redundant information generated by high-frequency sampling; however, when anomalies occur, critical monitoring points require higher sampling frequencies to capture the dynamic evolution of the anomaly. The adaptive sampling strategy dynamically adjusts the sampling frequency allocation by evaluating the information value of each monitoring point in real time. This allows high-risk, high-uncertainty monitoring points to receive more sampling resources, while the sampling frequency of low-risk, low-uncertainty monitoring points is appropriately reduced, thereby reducing the overall communication load and energy consumption of the system while ensuring monitoring efficiency. Scheduling based on information value rather than simply the rate of change is important because a high rate of water quality change does not necessarily mean that high-frequency sampling is needed. For example, some rapid fluctuations may be normal process disturbances rather than abnormal events, while some slow changes may indicate a trend of system deterioration that requires close monitoring. Information value comprehensively considers anomaly risks, estimation uncertainties, and data quality, and can more accurately reflect the marginal benefits of sampling.

[0102] Information value is assigned to each monitoring point. The information value of a monitoring point comprehensively considers risk, uncertainty, and data quality, and is calculated by multiplying the normalized risk term, the normalized uncertainty term, and the data quality weight. The design goal of information value is to quantify the information gain that can be obtained by sampling a monitoring point at the current moment. Monitoring points with high information value mean that sampling them can significantly reduce the overall uncertainty of the system or promptly capture abnormal events; therefore, sampling resources should be allocated preferentially. Information value comprehensively considers three dimensions, each reflecting a different aspect of sampling value. The risk term reflects the probability of anomalies occurring at the monitoring point now or in the future; high-risk monitoring points require close monitoring to detect anomalies in a timely manner. The uncertainty term reflects the reliability of the current state estimate of the monitoring point; high uncertainty means insufficient existing observational information, and increasing sampling can effectively reduce estimation errors. The data quality weight reflects the reliability of the sensors at the monitoring point; sampling high-quality data can provide more reliable information gain, while low-quality data is unlikely to improve estimation results even with high-frequency sampling. The multiplication of these three factors ensures that only monitoring sites that simultaneously meet the criteria of high risk, high uncertainty, and high data quality can obtain the highest information value, reflecting the comprehensive optimization principle of sampling resource allocation.

[0103] To ensure comparability of factors with different dimensions and ranges, each factor needs to be normalized before multiplication. The normalized risk term is calculated as follows: First, the original risk term is calculated, which can be either the anomaly score or the probability of exceeding the standard. The anomaly score is output from step 400; if the anomaly score is negative, it is set to zero. The probability of exceeding the standard can be obtained by assuming that the global state estimate follows a normal distribution and calculating the probability that the state value exceeds the standard limit. Then, the original risk term is divided by the sum of the maximum value of the original risk terms of all monitoring points and the preset minimum positive number to obtain the normalized risk term with a value range of 0 to 1. The preset minimum positive number is 0.01 to prevent the denominator from being zero. The normalized uncertainty term is calculated as follows: First, the original uncertainty term is calculated. This original uncertainty term can be the marginal uncertainty of the monitoring point or the uncertainty after projecting the uncertainty matrix onto the monitoring point. The marginal uncertainty of the monitoring point is output in step 300. Then, the original uncertainty term is divided by the sum of the maximum value of the original uncertainty terms for all monitoring points and a preset minimum positive number, resulting in a normalized uncertainty term with a value range of 0 to 1. The preset minimum positive number is 0.01 to prevent the denominator from being zero. The data quality weight is the data quality weight obtained in step 100, and its value range is already between 0 and 1. Through the above normalization process, all three factors of information value are mapped to a dimensionless range of 0 to 1, ensuring the comparability of information value from different monitoring points and the rationality of sampling frequency allocation. The design philosophy of information value is that monitoring points with high risk, high uncertainty, and good data quality have higher sampling value.

[0104] The sampling frequency is determined based on the information value. The sampling frequency of the monitoring point at the next moment is calculated as follows: the product of the preset scheduling gain and the information value, plus the preset lower limit of the sampling frequency, and then limited to between the preset lower limit and the preset upper limit of the sampling frequency through amplitude limiting. The amplitude limiting process restricts the calculation result to between the lower and upper limits; if the result is less than the lower limit, the lower limit is used; if the result is greater than the upper limit, the upper limit is used. The preset lower limit of the sampling frequency ranges from 1 to 6 times per hour, with a default value of 2 times per hour. This lower limit ensures the basic monitoring needs of the system and can be determined based on the slowest timescale of water quality changes. The preset upper limit of the sampling frequency ranges from 12 to 60 times per hour, with a default value of 30 times per hour. This upper limit is constrained by the sensor response time and communication capability and can be determined according to the equipment specifications. The preset scheduling gain ranges from 0.1 to 10, with a default value of 1. A larger gain indicates a more sensitive information value to the adjustment of the sampling frequency, and the gain value that maximizes monitoring effectiveness can be determined through simulation experiments.

[0105] Optimize sampling allocation under communication or energy consumption budget constraints. The system's preset total communication budget is set to 1.5 to 3 times the sum of the lower limits of sampling frequencies at all monitoring points, with a default value of 2 times. This budget value is determined comprehensively based on the system's communication bandwidth, energy consumption limitations, and data processing capabilities, requiring that the sum of sampling frequencies at all monitoring points not exceed this budget. Introducing budget constraints reflects the limited resources of the actual system. In practical applications, the bandwidth of the communication network, the energy consumption of the sensors, and the computing power of the data processing system all have upper limits, making it impossible to support all monitoring points sampling at the highest frequency simultaneously. Budget constraints explicitly incorporate resource limitations into the sampling optimization problem, ensuring that the generated sampling strategy is executable in the actual system. Setting the preset total communication budget to 1.5 to 3 times the sum of the lower limits of sampling frequencies means that, while ensuring the basic monitoring needs of all monitoring points, the system has 50% to 200% of additional resources available for adaptive scheduling. This setting avoids both excessive resource scarcity leading to an inability to respond to abnormal events and excessive resource abundance rendering adaptive scheduling meaningless. When the sum of sampling frequencies calculated based on information value exceeds the budget, adjustments need to be made to meet the constraints. The adjustment strategy employs a priority allocation method based on information value ranking to ensure that limited resources are allocated preferentially to monitoring points with the highest information value, maximizing resource utilization efficiency. When the total sampling frequency calculated using the above method exceeds the budget, the sampling frequency of each monitoring point needs to be reduced. The adjustment strategy involves ranking all monitoring points from highest to lowest information value, prioritizing the sampling frequency calculated based on information value for the top-ranked monitoring points, and gradually reducing the sampling frequency of the lower-ranked monitoring points until the sum of the sampling frequencies of all monitoring points meets the budget constraint. Specifically, the reduction method starts with the monitoring point ranked last in information value, successively reducing its sampling frequency to a preset lower limit. If the budget constraint is still not met after reducing to the lower limit, the reduction continues for the second-to-last monitoring point, and so on, until the total sampling frequency does not exceed the budget. Through this adjustment strategy, it is ensured that, under budget constraints, monitoring points with high risk, high uncertainty, and good data quality receive sufficient sampling resources, while the sampling frequencies of monitoring points with low risk, low uncertainty, or poor data quality are appropriately compressed.

[0106] If the system is equipped with temporary or mobile probes, inspection instructions or temporary deployment instructions are issued to the top-ranked, preset number of upstream candidate areas based on the source confidence ranking in step 400. The preset number of candidate inspection areas ranges from 1 to 5, with a default value of 3. This number is determined based on the number of available mobile probes and the inspection response time.

[0107] The output of step 500 is the sampling frequency or triggering rule for each monitoring point in the next cycle, as well as temporary deployment suggestions, for the sensor node sampling controller and the mobile probe scheduler to execute.

[0108] Step 600: Based on subsequent observation data, the model established in steps 200 to 400 is verified and updated to achieve continuous optimization and adaptive adjustment of the system.

[0109] Model verification and updates based on subsequent observation data are crucial for addressing the time-varying characteristics and operational changes of wastewater treatment systems. Wastewater treatment systems are affected by various factors, including fluctuations in influent water quality, temperature changes, equipment aging, and adjustments to operating strategies, resulting in non-constant hydraulic and biochemical characteristics. If model parameters remain fixed, deviations between the model and the actual system will gradually accumulate over time, leading to decreased accuracy in global state estimation and failure in anomaly detection. Continuous model verification and updates enable the system to automatically track the evolution of the actual process, maintaining the model's long-term effectiveness. Online incremental updates, rather than periodic offline retraining, are used to achieve seamless adaptive adjustments, avoid service interruptions during retraining, and respond promptly to sudden changes in operating conditions. The verification and updates cover parameters such as migration delays and decay coefficients in the time-delay space graph, statistical characteristics of data quality weights, and parameter templates for typical operating conditions, encompassing key factors affecting system performance.

[0110] Consistency verification between predictions and observations is performed. The observed values ​​collected at subsequent time points are compared with the prediction propagation results based on the time-delay spatial graph. Statistical analysis of the prediction bias is used to determine whether model parameters need adjustment. The purpose of consistency verification is to promptly identify model mismatch issues and prevent erroneous model parameters from having a long-term impact on estimation and detection results. The prediction bias is calculated by comparing the observed values ​​collected at subsequent time points with the prediction propagation results based on the time-delay spatial graph. The prediction bias is the difference between the actual observed value of the downstream node and the predicted value after propagation based on the historical observed values ​​of the upstream node through migration delay. If the prediction bias of a directed edge continuously exceeds a preset prediction bias threshold for several consecutive time points, the migration delay of that directed edge is adjusted. The migration delay verification error is calculated as follows: within a preset verification time window, the optimal migration delay value that minimizes the mean square error between the delayed sequence of the upstream node and the observed sequence of the downstream node is searched. The difference between the optimal migration delay value and the current migration delay is taken as the migration delay verification error. The preset verification time window length ranges from 5 to 20 sampling periods, with a default value of 10 sampling periods. A longer verification time window indicates better statistical stability of the migration delay verification but a slower response to changes in migration delay. It can be determined based on the flow fluctuation cycle and system dynamic characteristics. The migration delay adjustment method is to add the current migration delay to the increment obtained by multiplying the preset online update step size by the migration delay verification error.

[0111] The preset prediction deviation threshold ranges from 1 to 3 times the standard deviation of typical observation noise, with a default value of 2 times. This threshold is used to determine whether the prediction deviation is significant and can be determined based on the statistical characteristics of observation noise in historical data. The preset online update step size ranges from 0.01 to 0.5, with a default value of 0.1. A larger online update step size indicates a faster response to migration delay adjustment but poorer stability. The online update step size value that balances convergence speed and stability can be determined through simulation experiments. To improve the robustness of migration delay updates, a migration delay change rate constraint is introduced. If the absolute value of the calculated migration delay increment exceeds the preset proportional threshold of the current migration delay, the migration delay increment is truncated. The preset proportional threshold ranges from 10% to 30%, with a default value of 20%. This preset proportional threshold constraint prevents drastic changes in migration delay parameters due to abnormal operating conditions or measurement noise. The update mechanism employs a strategy that combines incremental adjustment with constraints on the rate of change of migration delay. This ensures that the migration delay parameters can track changes in operating conditions while avoiding drastic fluctuations in migration delay parameters due to noise from a single observation, thus guaranteeing the long-term stability and convergence of the model.

[0112] The system identifies and addresses long-term low-quality monitoring points. If the data quality weight of a monitoring point remains below a preset quality weight threshold for a consecutive preset number of observation periods, the system generates a maintenance prompt for that monitoring point, reminding maintenance personnel to check the sensor status. Identifying and addressing long-term low-quality monitoring points is a crucial mechanism for ensuring the long-term reliable operation of the system. During long-term operation, sensors experience gradual data quality degradation due to factors such as biofouling, electrode aging, and calibration drift. If not detected and addressed promptly, this degradation will continuously contaminate the global state estimation results and reduce the reliability of anomaly detection. By monitoring the long-term trend of data quality weights, the system can automatically identify sensors requiring maintenance, triggering preventative maintenance procedures and avoiding monitoring blind spots caused by complete sensor failure. Setting a consecutive observation period as the judgment criterion, rather than a single instance of low quality, distinguishes between occasional communication failures or environmental interference and continuous sensor degradation, avoiding unnecessary maintenance operations triggered by short-term fluctuations. After identifying long-term low-quality monitoring points, the system takes two measures. On the one hand, maintenance prompts are generated to notify maintenance personnel to conduct physical inspections and repairs, fundamentally solving the data quality problem; on the other hand, the weight of the monitoring point is further reduced during the global state estimation process, and the impact of low-quality data on the estimation results is mitigated through algorithm-level compensation measures, ensuring that the system can still maintain basic monitoring functions before the sensor repair is completed.

[0113] The preset observation period ranges from 3 to 10 sampling periods, with a default of 5 sampling periods. A longer period indicates a higher tolerance for occasional low-quality data, and this can be determined based on the typical duration of sensor failures. The preset quality weight threshold ranges from 0.2 to 0.5, with a default of 0.3. A higher threshold indicates stricter requirements for data quality, and this can be determined based on the statistical relationship between quality weights and actual sensor failures in historical data. Simultaneously, the weight of this monitoring point is further reduced during the global state estimation process to avoid low-quality data contaminating the global inference results. Furthermore, a sensor health status assessment mechanism is introduced to distinguish between technical anomalies and process anomalies. This mechanism comprehensively judges the sensor's health status by monitoring multiple technical indicators, including but not limited to signal drift trend, response time variation, baseline offset, calibration residuals, and signal noise levels. Signal drift trend is obtained by calculating the long-term trend slope of the sensor under stable operating conditions. If the absolute value of the slope exceeds a preset drift threshold, it is judged as a drift anomaly. The preset drift threshold is determined according to the sensor type and maintenance cycle. Response time change is obtained by comparing the current response time with the factory nominal response time. If the deviation exceeds a preset percentage of the nominal value, it is judged as a response anomaly. The preset percentage ranges from 20% to 50%, with a default value of 30%. Baseline offset is obtained by evaluating the residual during zero-point or standard liquid calibration. Signal noise level is obtained by calculating the standard deviation or power spectral density of the signal under stable operating conditions. When the sensor health status assessment mechanism determines that a monitoring point has a technical anomaly, the anomaly score of that monitoring point is filtered for technical anomalies in the anomaly detection step 400. That is, the high anomaly score of that monitoring point is not included in the source analysis of propagation anomalies, but is marked as a sensor fault anomaly and the maintenance process is triggered. By separating technical anomalies from process anomalies, it is possible to avoid misjudging sensor drift as contamination propagation, thereby improving the accuracy of anomaly identification and operational stability.

[0114] Periodically solidified operating condition templates. Model parameters under different typical operating conditions are identified and recorded, forming a switchable parameter set. Typical operating conditions include, but are not limited to, rainy season conditions, nighttime conditions, and high-load conditions. When a change in operating condition is detected, the corresponding parameter template is automatically invoked, reducing model mismatch problems caused by the change in operating conditions. The operating condition template mechanism is an effective strategy to address the periodic and predictable changes in the operating conditions of wastewater treatment systems. Wastewater treatment systems exhibit significant differences in operating modes at different times and seasons. For example, during the rainy season, the influent flow rate increases significantly and the water quality is diluted; at night, the influent load decreases and the hydraulic retention time increases; and during high-load conditions, the system operates at full capacity. The hydraulic and biochemical characteristics under these typical operating conditions differ significantly, and the corresponding migration lag and attenuation coefficient parameters also differ. If a single parameter set is used to address all operating conditions, a long period of model mismatch and increased estimation errors will occur when operating conditions change. The operating condition template mechanism identifies and records the characteristic parameters of each typical operating condition in historical operation, forming a parameter library. When an operating condition change signal is detected, the corresponding parameter template is automatically loaded as the initial value, significantly shortening the time for the model to adapt to the new operating condition. Operating condition identification can be based on external information such as flow rate, time period, and season, or it can be based on the statistical characteristics of water quality data for automatic clustering. Periodic updates of parameter templates ensure that the templates can reflect the long-term evolution trend of the system and maintain the timeliness of the templates.

[0115] Through the coordinated execution of steps 100 to 600 above, the present invention achieves global online situational awareness of wastewater quality.

[0116] Example 2

[0117] See Figure 4 As shown, an online wastewater quality monitoring system is provided, which stores computer-readable instructions. When the computer-readable instructions are read, the system can execute the aforementioned online wastewater quality monitoring method. The system includes:

[0118] The data processing module 101 performs time alignment processing and data quality assessment on the raw water quality sequences collected from each monitoring point to obtain the aligned data and data quality weights of the monitoring points.

[0119] The time-delay attenuation module 102 constructs a time-delay space graph based on the aligned data. The time-delay space graph includes a set of nodes, a set of edges, and the migration time delay and attenuation coefficient of each directed edge.

[0120] The global estimation module 103 performs a joint estimation of the global water quality state based on the aligned data, data quality weights, and time-delay space map, and obtains the global state estimate and uncertainty.

[0121] The source confidence module 104, based on the global state estimation, uncertainty and time delay space map, performs anomaly detection and identifies the propagation path and suspected source of the anomaly, and obtains anomaly score and source confidence ranking list;

[0122] The adaptive sampling module 105 formulates an adaptive sampling strategy based on the uncertainty, anomaly score and source confidence ranking list, and determines the sampling frequency of each monitoring point.

[0123] The embodiments of the present invention have been described above. However, the embodiments are not limited to the specific implementation methods described above. The specific implementation methods described above are merely illustrative and not restrictive. Those skilled in the art can make more equivalent embodiments under the guidance of the present embodiments, and all of them are within the protection scope of the present embodiments.

Claims

1. A method for online monitoring of wastewater quality, characterized in that, Includes the following steps: The original water quality sequences collected from each monitoring point are time-aligned and data quality is assessed to obtain the aligned data and data quality weights for the monitoring points. Based on the aligned data, a time-delay space graph is constructed, which includes a set of nodes, a set of edges, and the migration delay and decay coefficient of each directed edge. Based on the aligned data, data quality weights, and time-delay space map, the global water quality state is jointly estimated to obtain the global state estimate and uncertainty. Based on the global state estimation, uncertainty, and time delay space graph, anomaly detection is performed and the propagation path and suspected source of the anomaly are identified, resulting in anomaly scores and a source confidence ranking list. Based on the uncertainty, anomaly score, and source confidence ranking list, an adaptive sampling strategy is formulated to determine the sampling frequency for each monitoring point.

2. The method for online monitoring of wastewater quality according to claim 1, characterized in that, The aligned data and data quality weights obtained from the monitoring points include: The original water quality sequence includes the measured values ​​of water quality indicators collected at each monitoring point at each time point. The water quality indicators include pH, dissolved oxygen, oxidation-reduction potential, turbidity, conductivity, ammonia nitrogen, and chemical oxygen demand. The monitoring points include point metadata, which includes the spatial location of the monitoring point, the unit it belongs to, the sensor type, the sampling period, and communication quality characterization parameters, including packet loss rate and clock deviation. Time alignment processing includes estimating the clock deviation of each monitoring point and obtaining the aligned measurement sequence through a resampling method; Data quality assessment includes calculating data quality weights, specifically: multiplying the packet loss rate by a preset packet loss rate attenuation coefficient and taking the negative value as the exponent of the first exponential function; multiplying the standardized jump intensity by a preset jump attenuation coefficient and taking the negative value as the exponent of the second exponential function; and then multiplying the two exponential functions together to obtain the data quality weights; the standardized jump intensity is the standardized result of the original jump amplitude, and the original jump amplitude is the absolute value of the difference between the currently aligned data and the mean of the aligned data within the local sliding window.

3. The method for online monitoring of wastewater quality according to claim 1, characterized in that, Constructing the time-delay spatial diagram includes: Each monitoring point is used as a node to form a set of nodes, and the hydraulic connectivity and mixing relationship in the same pool are used as directed edges to form a set of edges. The direction of water flow determines the direction of the directed edges. For each directed edge, a learnable migration delay and attenuation coefficient are introduced, whereby the migration delay represents the time required for a water quality change to propagate from an upstream node to a downstream node, and the attenuation coefficient represents the degree of attenuation of the water quality index during the propagation process. The initial value of the migration time delay is the pool volume of the downstream unit divided by the representative flow rate of the unit, and the initial value of the attenuation coefficient is set according to the biochemical degradation characteristics of the water quality index or historical experience. The attenuation coefficient is continuously corrected using real-time data. Specifically, within a preset attenuation coefficient update window, the aligned data of the downstream node is divided by the aligned data of the upstream node after migration delay to obtain a ratio. Within the update window, the median or weighted average of the ratios at multiple times is taken as the updated value of the attenuation coefficient. The migration lag is corrected online using real-time data, and updated within a sliding time window by maximizing alignment correlation or minimizing prediction error.

4. The method for online monitoring of wastewater quality according to claim 1, characterized in that, Joint estimation of the overall water quality status includes: Set a global hidden state vector to represent the true water quality state of all units or regions at a certain moment; The relationship between the aligned data and the global hidden state is described by a time-delay observation equation, which represents the aligned data at the current moment as a function of the global hidden state vector at the historical moment. A weighted optimization method is used for online estimation of the global hidden state. The estimated value is obtained by solving a weighted least squares optimization problem. The estimated value of the global hidden state vector is the global state estimate. The objective function is composed of the sum of the observation fitting term and the smoothing regularization term. Calculate the uncertainty of the global hidden state estimation.

5. The method for online monitoring of wastewater quality according to claim 4, characterized in that, The specific form of the time-delay observation equation is as follows: the aligned data at the current moment is equal to the time-delay observation matrix multiplied by the extended state vector, plus the observation noise term; the extended state vector is formed by stacking the global hidden state vectors at the current moment and several historical moments in chronological order; the time-delay observation matrix is ​​determined by the topology of the time-delay space graph, the migration time delay parameters of each side, and the mapping relationship from the monitoring point to the region. The observation fitting term is the weighted sum of squared observation fitting errors of all monitoring points. The observation fitting error of each monitoring point is the difference between the aligned data of the monitoring point and the corresponding predicted observation value. The weight is the data quality weight. The smoothing regularization term is the square of the norm of the Laplacian matrix of the time-delay spatial graph applied to the global hidden state vector at the current time. The resulting value is multiplied by a preset smoothing regularization coefficient.

6. The method for online monitoring of wastewater quality according to claim 1, characterized in that, Anomaly detection and identification of anomaly propagation paths and suspected sources include: Calculate the regional residual for each monitoring point, where the regional residual is the difference between the actual observed value and the predicted observed value; Based on the regional residuals, a standardized anomaly score is constructed. The anomaly score of a monitoring point is the absolute value of the residual of the monitoring point divided by a marginal term, which is the sum of the marginal uncertainty of the monitoring point and a preset minimum positive number. The set of nodes with high abnormal scores is identified based on the abnormal scores. The set of nodes with high abnormal scores includes all monitoring points whose abnormal scores exceed a preset abnormal score threshold. By combining the aforementioned time-delay spatial diagram with reverse tracing, candidate source points are obtained; For each candidate source point, calculate its explanatory power for the set of high anomaly score nodes, and sort the candidate source points according to the explanatory power to obtain a source confidence ranking list.

7. The method for online monitoring of wastewater quality according to claim 6, characterized in that, The explanatory power is calculated as follows: for each node in the set of high anomaly score nodes, the product of the time matching factor, the path uncertainty discount factor and the anomaly score of the node is calculated, and then the product results of all nodes are summed. The time matching factor is an exponential decay function, with the base of the exponential decay function being the natural constant. The time term with a negative exponent is divided by a preset time matching decay scale. The time term is the current time minus the time when the node anomaly occurs, plus the cumulative migration delay from the candidate source point to the node, and the absolute value of the result is taken. The path uncertainty discount factor is obtained by summing the time delay uncertainties of each edge on the path according to the variance propagation law to obtain the total path time delay uncertainty, and then constructing a discount factor that decreases as the total path time delay uncertainty increases. The path is the path with the minimum cumulative migration delay from the candidate source point to the target node in the time delay space graph.

8. The method for online monitoring of wastewater quality according to claim 1, characterized in that, Developing an adaptive sampling strategy includes: An information value is assigned to each monitoring point. The information value is calculated by multiplying the normalized risk term, the normalized uncertainty term, and the data quality weight. The sampling frequency is determined based on the information value. The sampling frequency of the monitoring point is the product of the preset scheduling gain and the information value plus the preset lower limit of the sampling frequency, and then limited to between the preset lower limit of the sampling frequency and the preset upper limit of the sampling frequency through amplitude limiting processing. Optimize sampling allocation under communication or energy consumption budget constraints. When the sum of the sampling frequencies of all monitoring points exceeds the preset total communication budget, sort all monitoring points from largest to smallest information value. Starting from the last monitoring point, reduce the sampling frequency to the preset lower limit of the sampling frequency until the sum of the sampling frequencies is less than the preset total communication budget.

9. The method for online monitoring of wastewater quality according to claim 8, characterized in that, The normalized risk item is calculated as follows: First, the original risk item is calculated, which is the anomaly score. Then, the original risk item is divided by the sum of the maximum value of the original risk items of all monitoring points and the preset minimum positive number to obtain the normalized risk item with a value range of 0 to 1. The normalized uncertainty term is calculated as follows: First, the original uncertainty term is calculated, which is the marginal uncertainty of the monitoring point. The original uncertainty term is divided by the sum of the maximum value of the original uncertainty terms of all monitoring points and the preset minimum positive number to obtain the normalized uncertainty term with a value range of 0 to 1.

10. An online wastewater quality monitoring system, characterized in that, The system is used to store computer-readable instructions, which, when read, execute the online wastewater quality monitoring method according to any one of claims 1-9; the system includes: The data processing module performs time alignment processing and data quality assessment on the raw water quality sequences collected from each monitoring point to obtain the aligned data and data quality weights of the monitoring points. The time-delay decay module constructs a time-delay space graph based on the aligned data. The time-delay space graph includes a set of nodes, a set of edges, and the migration delay and decay coefficient of each directed edge. The global estimation module performs a joint estimation of the global water quality state based on the aligned data, data quality weights, and time-delay space map, and obtains the global state estimate and uncertainty. The source confidence module, based on the global state estimation, uncertainty and time delay space graph, performs anomaly detection and identifies the propagation path and suspected source of the anomaly, and obtains anomaly score and source confidence ranking list; The adaptive sampling module formulates an adaptive sampling strategy based on the uncertainty, anomaly score, and source confidence ranking list, and determines the sampling frequency for each monitoring point.

Images