Early warning method for blocking risk of rainwater collection pipe network

Abstract

The invention provides a rainwater collection pipe network blockage risk early warning method, which comprises the steps of collecting multi-point flow and environmental parameter data, fusing water level and rainfall intensity, and realizing unified processing of multi-source data through timestamp alignment and mean filtering; using an empirical mode decomposition algorithm to decompose the flow time sequence data into trend, period and high-frequency fluctuation components, and identifying abnormal segments through sliding window energy calculation and dynamic threshold gating; according to the method, the change point detection algorithm and the Gaussian mixture model are combined, energy accumulation sudden change and trend attenuation are verified in a layered mode, and then structured early-stage blockage risk early warning is generated, the accuracy and timeliness of anomaly detection are effectively improved, the pipe network blockage supervision and maintenance cost can be reduced, and the operation reliability of an urban drainage system is enhanced.

Description

Technical Field

[0001] This invention relates to the field of intelligent monitoring and anomaly detection technology for urban stormwater pipe networks, and in particular to a method for early warning of blockage risks in stormwater collection pipe networks. Background Technology

[0002] Currently, intelligent monitoring and anomaly detection technologies for urban stormwater drainage networks are developing rapidly and are widely used in various fields such as urban flood control, pipeline operation and maintenance, and smart cities. Mainstream technical solutions generally rely on sensor networks to collect multi-source time-series data such as flow rate and water level, combining this data with machine learning, deep learning, or statistical analysis methods to achieve status monitoring and risk warning of stormwater drainage networks. With the maturity of artificial intelligence and big data technologies, more and more systems are adopting time-series modeling tools such as Long Short-Term Memory (LSTM) networks, autoencoders, and convolutional neural networks to predict flow rate and water level data, identify anomalies, and analyze trends, thereby helping to reduce the risk of urban flooding caused by pipeline blockage. Existing representative technologies mainly include trend-prediction-based anomaly monitoring, threshold-based anomaly detection, and multi-point collaborative analysis combining pipeline topology. In practical applications, real-time monitoring data is typically compared with a baseline by fitting parameters to pipeline flow rate and its changing trends, performing periodic statistical analysis, or establishing a reference model under normal operating conditions to quickly identify prominent anomaly signals. These methods are particularly suitable for detecting blockage events or sudden anomalies that have developed to a significant stage, and they have high accuracy in identifying obvious signals such as rapid changes in flow rate and water level rises. However, in the early stages of pipeline congestion, anomalies are typically only slight fluctuations or gradual flow attenuation. Initial signals are often obscured by high-frequency noise or environmental disturbances, making it difficult for traditional predictive reconstruction, trend detection, or clustering methods to accurately detect anomalies at this stage, resulting in significant delays in risk warnings. Commonly used noise filtering strategies ignore or eliminate residual fluctuation signals, thus losing the ability to detect early congestion trends. Furthermore, while some methods attempt to introduce deep models or attention mechanisms to enhance temporal modeling capabilities, the computational complexity is too high for resource-constrained edge nodes, making large-scale real-time deployment difficult. Analysis paths relying on network topology suffer from insufficient flexibility and generalization capabilities. Summary of the Invention

[0003] In order to solve the above-mentioned technical problems, the present invention provides a method for early warning of the risk of blockage in rainwater collection pipe networks.

[0004] The technical solution of this invention is implemented as follows: A method for early warning of rainwater harvesting pipe network blockage risk, comprising: S1: Collect real-time flow time series data and environmental parameters from rainwater pipe network monitoring points. The flow time series data includes continuous time series from multiple deployment locations, and the environmental parameters include water level height and rainfall intensity. Construct an initial input dataset. S2: Perform multi-scale fluctuation decomposition processing on the traffic time series data, and decompose the original sequence into trend principal component sequence, periodic component sequence and residual fluctuation component sequence based on empirical mode decomposition algorithm, wherein the residual fluctuation component represents high-frequency non-steady-state disturbances. S3: Apply sliding window energy calculation processing to the residual fluctuation component sequence, and generate the original fluctuation energy index based on the sum of squared deviations of the data points and local mean within the window. The length of the sliding window is determined according to the typical blockage response cycle of the pipeline network. S4: Based on the statistical characteristics of residual fluctuations under historical normal operating conditions, a dynamic noise threshold model is established. An adaptive threshold is generated by calculating the mean and standard deviation of historical data. The original fluctuation energy index is then gated, and only fluctuation segments exceeding the dynamic threshold are retained as valid abnormal candidate signals. S5: Perform cumulative trajectory construction on the effective fluctuation energy index after gating, use the change point detection algorithm to identify significant jump nodes in the energy accumulation level, and generate a preliminary set of abnormal timestamps. The change point detection algorithm uses the Pelt algorithm for optimal segmentation. S6: Perform short-time linear fitting based on the trend principal component sequence, calculate the slope change rate within the current window, and generate a first-layer verification pass marker when the slope change rate is less than a preset negative threshold and overlaps with the initial anomaly timestamp. S7: Establish a Gaussian mixture model based on the historical normal energy distribution, calculate the statistical significance score of the current cumulative energy, and generate a second-level validation pass mark when the significance score of three consecutive windows exceeds the 95% confidence interval. S8: When the first-level verification pass mark and the second-level verification pass mark are obtained simultaneously, an early blockage risk warning event containing a timestamp is generated, triggering the pipeline maintenance scheduling process. The warning event contains a risk level identifier and an associated monitoring point number.

[0005] The present invention provides a method for early warning of rainwater harvesting pipe network blockage risk, which has the following beneficial effects: (1) This invention introduces a multi-scale fluctuation decomposition mechanism to decouple the original flow time series data into three components: trend, period, and residual. It innovatively mines the non-random energy accumulation characteristics in the residual fluctuation component as an early anomaly criterion, effectively overcoming the technical defect of losing potential fault information caused by directly filtering the residual in conventional practices. Since the initial blockage often leads to the gradual deterioration of the pipeline flow stability, it is manifested as the spatial accumulation and energy enhancement of the residual fluctuation in a local time period, rather than simple high-frequency noise. Therefore, this scheme realizes the quantitative capture of such weak disturbances by implementing sliding window segmentation analysis on the residual sequence and calculating the normalized fluctuation energy index. Furthermore, combined with the dynamic noise threshold gating mechanism based on historical statistical modeling, only the energy changes exceeding the background noise level are retained for accumulation, which significantly improves the signal-to-noise separation capability and avoids false triggering caused by random fluctuations. Thus, without relying on complex prediction models or global topology, a high-sensitivity technical path for early anomaly detection is constructed. (2) This invention designs a dual-layer verification mechanism that integrates trend deviation verification and statistical significance assessment, forming a closed-loop discrimination logic from preliminary candidate to confirmed warning. The first layer of verification detects whether the sudden increase in energy is accompanied by a hidden decrease in the slope of flow growth by fitting the main trend sequence in a short time linear fit. It utilizes the physical law that the early stage of blockage development is often accompanied by a slow decrease in transport efficiency to achieve contextual consistency verification of abnormal events. The second layer establishes a probability model based on the residual energy distribution under historical normal operating conditions to assess the statistical deviation of the current observation. The risk is finally confirmed only when multiple consecutive windows fall outside the 95% confidence interval, effectively suppressing false alarms caused by occasional disturbances. This dual constraint strategy not only enhances the anti-interference ability of the method in complex operating environments, but also improves the interpretability and credibility of the warning results. Compared with the single threshold method or end-to-end deep learning model widely used in the prior art, it has stronger engineering applicability and stability, and is especially suitable for actual drainage pipe network scenarios lacking a large number of labeled samples. (3) This invention abandons the mainstream approach that relies on long-sequence prediction and reconstruction, attention mechanism modeling, or complex graph neural networks for state inference. Instead, it focuses on the residual energy evolution trajectory, which is ignored by traditional processing procedures, and proposes a lightweight, low-computational-overhead early risk identification framework with good real-time performance and edge deployment adaptability. The sliding window energy calculation, change point detection, and statistical verification modules used can all be implemented efficiently. The overall algorithm complexity is moderate and does not require high-performance computing platform support. It is particularly suitable for embedding in resource-constrained field monitoring terminals, ensuring detection accuracy while meeting the requirements of stormwater pipe network systems for millisecond-level response and continuous online operation. In addition, this method does not rely on complete periodic data or prior topological knowledge, has strong versatility and scalability, and can flexibly adapt to pipe network layouts of different scales and structures, providing practical technical support for the intelligent operation and maintenance of urban drainage systems. Attached Figure Description

[0006] Figure 1 A flowchart of a rainwater harvesting pipe network blockage risk early warning method according to the present invention; Figure 2 This is a sub-flowchart of a rainwater harvesting pipe network blockage risk early warning method according to the present invention; Figure 3 This is another sub-flowchart of the rainwater harvesting pipe network blockage risk early warning method of the present invention. Detailed Implementation

[0007] Embodiments of the present invention are described in detail below, examples of which are shown in the accompanying drawings, wherein the same or similar reference numerals denote the same or similar elements or elements having the same or similar functions throughout. The embodiments described below with reference to the accompanying drawings are exemplary and are only used to explain the present invention, and should not be construed as limiting the present invention.

[0008] The following disclosure provides many different embodiments or examples for implementing different structures of the invention. To simplify the disclosure, specific examples of components and arrangements are described below. Of course, these are merely examples and are not intended to limit the invention. Furthermore, reference numerals and / or letters may be repeated in different examples; such repetition is for simplification and clarity and does not in itself indicate a relationship between the various embodiments and / or arrangements discussed.

[0009] like Figure 1 As shown, the present invention provides a method for early warning of rainwater harvesting pipe network blockage risk, specifically including: S1: Collect real-time flow time series data and environmental parameters from rainwater pipe network monitoring points. The flow time series data includes continuous time series from multiple deployment locations, and the environmental parameters include water level height and rainfall intensity. Construct an initial input dataset. S2: Perform multi-scale fluctuation decomposition processing on the traffic time series data, and decompose the original sequence into trend principal component sequence, periodic component sequence and residual fluctuation component sequence based on empirical mode decomposition algorithm, wherein the residual fluctuation component represents high-frequency non-steady-state disturbances. S3: Apply sliding window energy calculation processing to the residual fluctuation component sequence, and generate the original fluctuation energy index based on the sum of squared deviations of the data points and local mean within the window. The length of the sliding window is determined according to the typical blockage response cycle of the pipeline network. S4: Based on the statistical characteristics of residual fluctuations under historical normal operating conditions, a dynamic noise threshold model is established. An adaptive threshold is generated by calculating the mean and standard deviation of historical data. The original fluctuation energy index is then gated, and only fluctuation segments exceeding the dynamic threshold are retained as valid abnormal candidate signals. S5: Perform cumulative trajectory construction on the effective fluctuation energy index after gating, use the change point detection algorithm to identify significant jump nodes in the energy accumulation level, and generate a preliminary set of abnormal timestamps. The change point detection algorithm uses the Pelt algorithm for optimal segmentation. S6: Perform short-time linear fitting based on the trend principal component sequence, calculate the slope change rate within the current window, and generate a first-layer verification pass marker when the slope change rate is less than a preset negative threshold and overlaps with the initial anomaly timestamp. S7: Establish a Gaussian mixture model based on the historical normal energy distribution, calculate the statistical significance score of the current cumulative energy, and generate a second-level validation pass mark when the significance score of three consecutive windows exceeds the 95% confidence interval. S8: When the first-level verification pass mark and the second-level verification pass mark are obtained simultaneously, an early blockage risk warning event containing a timestamp is generated, triggering the pipeline maintenance scheduling process. The warning event contains a risk level identifier and an associated monitoring point number.

[0010] Step S1: Collect real-time flow time-series data and environmental parameters from stormwater pipe network monitoring points. The flow time-series data includes continuous time series data from multiple deployment locations, and the environmental parameters include water level and rainfall intensity, constructing an initial input dataset. Specifically, this includes: S1.1: Periodically collect data from multiple sets of flow sensors deployed at key nodes of the rainwater pipe network, and obtain the instantaneous flow value of each monitoring point based on the Modbus protocol to form continuous time series data containing multiple deployment locations, which serves as the input for subsequent multi-scale fluctuation decomposition processing; A periodic sampling plan is set for multiple sets of flow sensors deployed at key nodes of the rainwater pipe network, and the instantaneous flow value parameters of each monitoring point are read using the Modbus protocol (sampling interval: determined according to the dynamic response characteristics of the pipe network flow). Furthermore, based on the Modbus register address mapping table, the raw data byte stream returned by the sensor is parsed, and the original 16-bit or 32-bit register data is restored to the physical flow value using byte order conversion and data type parsing algorithms, and the corresponding sampling timestamp information (UTC format accurate to milliseconds) is attached. Furthermore, time series splicing processing is applied to the instantaneous flow values ​​obtained from the analysis, combining the flow data of each key node at the same sampling time into a multi-dimensional time series structure to maintain the spatial synchronization and temporal continuity of the data. Furthermore, a packet loss detection and completion algorithm is used to verify the integrity of the multidimensional time series. For the detected missing sampling points, the missing values ​​are filled in by nearest neighbor interpolation or linear interpolation to ensure that the input of subsequent fluctuation decomposition is complete. Through the above data acquisition and processing methods, the periodically sampled instantaneous flow values ​​are transformed into continuous time series data containing multiple deployment locations, realizing real-time quantitative characterization of the pipeline network operation status and providing high-quality input for subsequent multi-scale fluctuation decomposition steps; For example, in a municipal rainwater harvesting system, electromagnetic flow sensors supporting the Modbus RTU protocol were installed at five key discharge nodes, with a sampling interval configured to 2 seconds. The sampling task was initiated by a centralized acquisition terminal (supporting RS485 bus), with polling executed in ascending order of node ID. Each time, registers 0x0004 and 0x0005 were read to obtain the high and low bits of the instantaneous flow value. The read raw data was converted from byte order to floating-point flow values ​​in liters per second and bound with a sampling timestamp (UTC format, accurate to milliseconds). The flow values ​​from the five nodes obtained from each sampling were concatenated in timestamp order to form a 5×N dimensional matrix time series, where N is the number of sampling cycles. During the acquisition process, data was found to be missing for node 3 in the 100th and 150th sampling cycles. A linear interpolation formula was used:

[0011] in, For flow rate value, For missing sampling time points, the flow rate value at the missing points was calculated by interpolating adjacent sample values. This resulted in a continuous time series dataset covering five monitoring locations with a length of 10,000 points, meeting the input requirements of the subsequent empirical mode decomposition algorithm. During the trial operation phase, it demonstrated significantly improved flow monitoring completeness and data smoothness. S1.2: Perform synchronous sampling operation on the water level sensors deployed at the pipeline junctions and key inspection wells, and obtain real-time water level height data based on the ADC analog-to-digital conversion module to supplement the spatial status information of the flow data and improve the spatial correlation of subsequent anomaly identification; For multiple sets of water level sensors installed at pipeline junctions and key inspection wells, a synchronous sampling control method is adopted (parameter: the synchronous trigger signal source is a unified real-time clock) to achieve time-consistent acquisition of water level signals from multiple measuring points; Furthermore, the analog-to-digital converter (ADC) module (parameters: range 0-5V, resolution 12bit, sampling frequency 100Hz) is used to convert the analog voltage signal output by the sensor into digital data and obtain digital data proportional to the water level height. Furthermore, a linear calibration algorithm (parameter: calibration coefficients are derived from the sensor's factory calibration curve) is used to map digital data to actual water level height and generate an accurate water level data sequence. Furthermore, an outlier removal algorithm (parameter: removal threshold based on ±3 standard deviations) is used to filter out abnormal data caused by instantaneous fluctuations or sensor malfunctions, and obtain a stable water level time series. Furthermore, by processing data buffering and timestamp embedding (parameter: buffer queue length N=256 sample points), the time reference of water level data and flow data is unified, providing a dual guarantee of spatial location and temporal accuracy for subsequent multi-source data fusion; By using a synchronous sampling and ADC conversion link, the real-time analog output of the water level sensor is converted into a time-aligned digital water level height index that can be directly used in spatial correlation analysis, thereby effectively supplementing the spatial status information of the pipeline network. For example, in a city's stormwater drainage network monitoring scenario, ultrasonic water level sensors with a range of 0-5m are deployed at three network junctions and two key inspection wells. The sensor's simulated output voltage range is 0.5V to 4.5V, corresponding to a linear water level range of 0.2m to 4.8m. The synchronous sampling trigger signal is provided by a GPS timing module, the sampling frequency is set to 100Hz, and the ADC module uses 12-bit resolution. The voltage conversion formula is as follows:

[0012] in, The water level is in meters (m). The sensor output voltage (V). and The coefficients obtained after calibration are the measured values. =1.25, = 0.25. The raw digital value obtained through ADC conversion is converted into water level height using the above formula, and sudden noise data is removed using an outlier removal algorithm. The fluctuation range of the water level data obtained in a single acquisition cycle remains within ±0.02m. The processed water level time series includes a UTC timestamp, achieving alignment with the flow time series obtained in step S1.1 within the millisecond range. The output structured water level height data is directly called by the multi-source data time synchronization mechanism in subsequent step S1.3, significantly improving the accuracy of spatial correlation anomaly identification. S1.3: Real-time rainfall intensity data is obtained based on the meteorological API interface. The rainfall intensity sequence is synchronized with the flow and water level data in time by combining the timestamp alignment mechanism to eliminate the time lag difference of multi-source data acquisition and form a multi-dimensional input feature set under a unified time benchmark. In the multi-source data acquisition chain of the stormwater pipe network monitoring system, for the environmental parameter of rainfall intensity, a real-time data acquisition operation based on the meteorological API interface is performed (parameters: API call frequency is 60 seconds / time, data granularity is minute level) to capture the set of rainfall intensity values ​​for the target monitoring area; Furthermore, by using the data parsing module returned by the request (parameters: JSON format parsing, fields including precipitation_rate and observation_time), the rainfall intensity value is effectively bound to the corresponding meteorological station observation time, and the initial time series of rainfall intensity is obtained; Furthermore, through a timestamp alignment mechanism (parameters: alignment tolerance ±1 second, reference master clock is synchronously generated by the flow sensor data acquisition module), the time base of the rainfall intensity sequence and the flow sequence are unified, and an aligned auxiliary time sequence of rainfall intensity is generated; Furthermore, by using a bidirectional interpolation algorithm (parameters: linear interpolation to fill in missing points, cubic spline interpolation to refine the values ​​at non-integer intervals), rainfall intensity data matching the sampling frequency of flow and water level sequences is supplemented, and a sequence matrix with high sampling consistency is generated. By using a multi-source data merging process, the rainfall intensity data that has been optimized by time alignment and interpolation in the previous step is vertically spliced ​​with the flow and water level data according to a unified time index to form a multi-dimensional input feature set under a unified time benchmark, thereby realizing the time-series fusion of data from different sources to support subsequent anomaly identification. For example, within the monitoring area of ​​a 42-kilometer-long urban stormwater main pipeline network, a meteorological API interface is deployed with a call parameter every 60 seconds to acquire rainfall intensity data from an automatic weather station 8 kilometers from the center point. The observed rainfall intensity ranges from 0 to 15 mm / h, with an accuracy of 0.1 mm / h. The flow sensor sampling frequency is 1 Hz, the water level sensor sampling frequency is 0.5 Hz, and the master clock accuracy is ±0.5 seconds. During timestamp alignment, a differential matching strategy is employed to align meteorological data with flow and water level data within a time difference of no more than ±1 second. Missing rainfall intensity data points are calculated using a linear interpolation formula. For interpolation results at non-integer intervals, a cubic spline formula is used for weighted correction to ensure matching with flow and water level characteristics under high-frequency sampling conditions. The output unified time-referenced multidimensional input feature set includes three components: flow rate (unit: L / s), water level (unit: cm), and rainfall intensity (unit: mm / h). This significantly improves temporal consistency and demonstrates more stable feature response capabilities in subsequent anomaly identification steps. S1.4: Perform sliding window mean filtering on the collected raw flow, water level and rainfall intensity data, and eliminate instantaneous noise interference based on the moving average algorithm to improve the smoothness and stability of the raw data sequence, which serves as the optimized input for subsequent empirical mode decomposition processing; S1.5: Perform standardization processing on the multi-source data after timestamp alignment, and use the Z-score normalization method to perform zero-mean unit variance transformation on the data of each dimension to eliminate the differences in the dimensions of different sensors, and generate a structured initial input dataset for the risk analysis module to call.

[0013] Step S2: Perform multi-scale fluctuation decomposition processing on the traffic time series data. Based on the empirical mode decomposition algorithm, the original sequence is decomposed into a trend principal component sequence, a periodic component sequence, and a residual fluctuation component sequence, wherein the residual fluctuation component represents high-frequency non-steady-state disturbances. Specifically, this includes: S2.1: The original traffic time series data is preprocessed by using a sliding window mean filtering algorithm to perform noise suppression processing on the original sequence to obtain a smoothed input signal sequence, which is used as the input of the empirical mode decomposition algorithm; S2.2: Based on the empirical mode decomposition algorithm, local extremum point identification processing is performed on the smoothed input signal sequence, and upper and lower envelopes are constructed using the cubic spline interpolation method to obtain the signal envelope features within the local extremum region; S2.3: Calculate the mean of the upper and lower envelopes to generate a local mean sequence, and subtract the original signal from the local mean sequence to obtain a candidate intrinsic mode function sequence, so as to determine whether it meets the two criteria of intrinsic mode function: the number of extreme points is equal to or differs from the number of zero crossover points by no more than one, and the local mean approaches zero. Based on the upper and lower envelopes obtained through processing S2.2, the arithmetic mean calculation method (parameter: the amplitude of the upper and lower envelopes at each time point is taken separately) is used to obtain the point-by-point mean of the two envelope curves, resulting in the corresponding local mean sequence. The formula is as follows:

[0014] in, The upper envelope amplitude, The lower envelope magnitude is the value of this formula, which is used to calculate the value of each point in the local mean sequence. Furthermore, by using a difference operation method (parameters: the original signal amplitude sequence and the local mean sequence), the original flow time series signal is subtracted point by point from the local mean sequence to obtain the candidate intrinsic mode function sequence, as shown in the formula:

[0015] in, The original signal amplitude sequence, It is a local mean sequence. A sequence of candidate intrinsic mode functions; Furthermore, the extreme value-zero crossover consistency criterion method is adopted to determine the candidate intrinsic mode function sequence, and the number of extreme points and the number of zero crossover points are obtained respectively. The criterion formula is as follows:

[0016] in, The number of intersections is zero. The number of extreme points. The criterion requires the difference between the two. No more than ; Furthermore, by using the mean-approaching calculation method (parameter: the overall mean of the candidate intrinsic mode function sequence), the criterion for whether the sequence satisfies the condition that the local mean approaches zero is verified. Through the above judgment and processing, the original signal difference result is transformed into candidate components with intrinsic mode function properties, which helps to improve the accuracy of subsequent EMD algorithm iteration screening. By jointly verifying the extreme value-zero crossover consistency criterion and the mean-approaching criterion, the result of the previous step is transformed into candidate IMF components that conform to the definition of EMD decomposition, thereby effectively preserving the non-stationary local features of the signal. For example, during the multi-scale fluctuation decomposition of flow monitoring data at a key node of a stormwater pipe network, a signal segment of length 1024 points is extracted. After cubic spline interpolation, the upper and lower envelope sequences have amplitudes of 1.35 and 0.95 at a certain time point t, respectively. Substituting these values ​​into the local mean calculation formula... The local mean at this moment is 1.15. Subtracting the amplitude of the original signal at this moment (1.28) from the local mean (1.15) yields a candidate IMF component value of 0.13. The total number of zero-crossing points in the entire segment is 215, and the number of extreme points is 214, with a difference of 1, satisfying the extreme-zero crossover consistency criterion. Simultaneously, the overall mean of this candidate component is 0.003, satisfying the mean-to-zero condition. Therefore, this signal segment confirms that it meets the IMF definition requirements of EMD and enters the iterative screening stage of S2.4 to ensure that the high-frequency non-steady-state disturbance characteristics used for initial blockage detection are retained in subsequent decomposition results. S2.4: If the candidate intrinsic mode function sequence does not meet the intrinsic mode function criterion, then perform iterative screening processing on it, repeat the local extremum identification and mean calculation process until the intrinsic mode function definition is met, thereby obtaining the first-order intrinsic mode function component as the initial extraction result of the residual fluctuation component; S2.5: Separate the first-order intrinsic mode function component from the original signal to obtain the remaining signal sequence, and repeat the empirical mode decomposition process to extract multiple intrinsic mode function components in sequence. Finally, classify the high-frequency components as residual fluctuation components and the low-frequency components as trend principal components and periodic components.

[0017] like Figure 2 As shown, step S3 involves applying a sliding window energy calculation to the residual fluctuation component sequence, generating an original fluctuation energy index based on the sum of squared deviations between the data points within the window and the local mean. The length of the sliding window is determined according to the typical blockage response cycle of the pipeline network. Specifically, this includes: S3.1: Perform sliding window partitioning on the residual fluctuation component sequence, set the window length based on the typical blockage response cycle of the pipeline network, and ensure that the window covers the time span of the initial disturbance of the blockage to obtain local residual fluctuation subsequences within multiple time windows; For the residual fluctuation component sequence input, a sliding window segmentation method is adopted (parameter: the window length is set according to the time span of the typical blockage response cycle of the pipeline network) to split the global time domain signal into multiple local time windows; Furthermore, by using a window boundary index generation algorithm (parameter: starting index offset is half the window length), the residual fluctuation sequence is shifted and covered between adjacent windows to ensure that there is partial overlap between the data segments of each window, thereby improving the probability of capturing weak initial disturbance signals. Furthermore, by using a timestamp mapping method (parameter: linear mapping relationship between the original sequence index and the sampling time), each local window is accurately bound to the corresponding physical time range, and a set of window timestamps is generated for time positioning in subsequent energy calculations. Furthermore, by employing data slicing operations (parameters: window start and length strictly match the set period), data vectors for corresponding time periods are extracted from the residual fluctuation component sequence to form multiple local residual fluctuation subsequences; By using a window length selection strategy, the results of the previous step are transformed into data fragments that cover the initial disturbance characteristics of the blockage, thereby providing an analytical baseline with a stable time scale for the energy calculation process. For example, in a stormwater pipe network monitoring scenario, the typical blockage response cycle is 15 minutes, the system sampling frequency is 1Hz, and the window length for the residual fluctuation component sequence is set to 900 sampling points. A sliding window segmentation method is used, with a window step size of 300 sampling points, forming a local residual fluctuation subsequence with a length of 900 points and an overlap of 600 points between adjacent windows. A window boundary index generation algorithm determines that the first window covers sampling indices 0 to 899, the second window covers sampling indices 300 to 1199, and so on. Using a timestamp mapping method, the sampling index interval is linearly converted to a physical time range; for example, the first window is mapped to the monitoring start time up to the 900th second, and the second window is mapped to the 300th to the 1200th second. Under this configuration, a total of 50 sets of local residual fluctuation subsequences are generated, which can effectively cover the initial blockage disturbance signal area in subsequent energy calculations and significantly improve detection sensitivity and positioning accuracy. S3.2: Perform local mean calculation on the local residual fluctuation subsequence within each time window, and generate a local benchmark value based on the arithmetic mean of the data points within the sliding window, which serves as the reference baseline for calculating the perturbation energy of that window; S3.3: Perform deviation square calculation on each local residual fluctuation subsequence and the corresponding local benchmark value, and generate the local perturbation energy density at each point in the window based on the point-by-point deviation square value to characterize the degree of deviation of the point from the baseline; For the local residual fluctuation subsequences and corresponding local benchmark values ​​within each time window, the difference square calculation method (parameters: local residual fluctuation value, local benchmark value) is used to achieve a quantitative characterization of the deviation degree point by point. Furthermore, the local perturbation energy density for each data point is calculated using the squared deviation formula, as follows:

[0018] in, This represents the residual fluctuation value for that data point. This is the local reference value for the current window. Let be the local perturbation energy density of the i-th data point; Furthermore, by using a vectorized operation method (parameters: window length L, residual fluctuation array X, reference value μ), the batch calculation of the squared deviation of each point in the entire window is realized, and the local energy density sequence is obtained; Furthermore, by using the array index matching method, the calculated local energy density sequence is bound and mapped with the corresponding timestamp to form window feature data containing time location and perturbation energy magnitude; The local perturbation energy density sequence generated by the point-by-point deviation square value converts the degree of deviation of each point from the baseline within the window into a quantifiable energy value, thereby achieving a precise characterization of the intensity of high-frequency perturbations. For example, a residual fluctuation subsequence of 60 seconds is detected at monitoring point A. The window length L is set to 60, and the local reference value μ is calculated to be 0.045 m³ / s by arithmetic mean. The residual fluctuation value x at a certain moment in this sequence is 0.062 m³ / s. Substituting into the formula... The local perturbation energy density at that point is obtained as follows: (m³ / s)². In vectorized calculations, the residual fluctuation values ​​of all 60 points in this window are subjected to batch deviation squared operation with μ to generate an energy density sequence of length 60. Combined with timestamp binding, for example, the local energy density at the 15th second is... The local energy density at the 47th second is The resulting time-energy pair sequence can be used for subsequent energy accumulation and change point detection. In this scenario, local energy density segments higher than 0.00025 (m³ / s)² exhibit a characteristic of continuous operation for more than 20 seconds, significantly improving the sensitivity of detecting weak anomalies in the early stages of blockage. S3.4: Perform summation processing on the local disturbance energy density within each window, and generate the original fluctuation energy index of the window based on the cumulative sum of the squared deviation values ​​of all points within the window, which serves as a comprehensive measure of the disturbance intensity within that time period; S3.5: Normalize the raw fluctuation energy index of all windows. Perform standardization operation based on the mean and standard deviation of historical data to eliminate the dimensional differences between different monitoring points and obtain the energy characteristic sequence under a unified scale. For the original fluctuation energy index sequences of each window output in the previous sub-step S3.4, the Z-score standardization method (parameters: historical data mean μ and standard deviation σ) is used to eliminate the dimensional differences between different monitoring points; Furthermore, through the mean calculation module (parameter: energy index sample set of the same window length under historical normal working conditions), the original fluctuation energy index of each monitoring point window is centrally processed, and the mean offset of each window is obtained. Furthermore, through the standard deviation calculation module (parameter: historical energy index sample set of the same monitoring point), the scaling factor of the centered energy index is extracted, and the scaling parameters of each window are generated; Furthermore, the uniform-scale energy value for each window is calculated using the following standardized formula:

[0019] in, This is the original fluctuation energy index for the current window. This is the average of historical data. The standard deviation of historical data; Furthermore, by using a vectorized batch processing method, the standardized calculation results of all windows are constructed into an energy feature sequence matrix and sorted under a unified time axis to obtain an energy feature sequence at a unified scale. This standardized processing method transforms the original fluctuation energy index into data features that can be directly compared across monitoring points, achieving consistency of energy features in the spatial dimension. For example, in a rainwater pipe network operation monitoring scenario, the sliding window length corresponding to a typical blockage response cycle of the pipe network is selected as 60 seconds. The average energy index μ under historical normal operating conditions is 25.3, and the standard deviation σ is 4.7. The original fluctuation energy index E of a certain window in the current monitoring cycle is 34.5. Substituting this value into the above formula, we can calculate:

[0020] The standardization result is 1.957, which is used in the energy feature sequence at a unified scale. For different monitoring points, the historical data μ and σ values ​​vary significantly. Through this standardization process, the output energy feature sequences exhibit direct comparability at the same scale. In this embodiment, when the standardized energy feature sequence is input into the subsequent dynamic noise threshold model, the sensitivity for identifying abnormal candidate signals is significantly improved, and the false alarm rate between different monitoring points is greatly reduced.

[0021] like Figure 3 As shown, step S4 involves establishing a dynamic noise threshold model based on the statistical characteristics of residual fluctuations under historical normal operating conditions. An adaptive threshold is generated by calculating the mean and standard deviation of historical data. This threshold is then used to gate the original fluctuation energy index, retaining only fluctuation segments exceeding the dynamic threshold as valid anomaly candidate signals. Specifically, this includes: S4.1: Perform sliding window statistical analysis on the historical residual fluctuation component sequence, and generate historical fluctuation energy index based on the sum of squared deviations between the data points within the window and the local mean. The length of the sliding window is consistent with that set in S3 to ensure time scale matching. S4.2: Normalize the historical fluctuation energy index and construct a dynamic noise threshold model based on the mean and standard deviation of historical data to adapt to the background noise differences of different monitoring points in different seasons or working conditions. S4.3: Based on the dynamic noise threshold model, perform threshold gating on the current original fluctuation energy index generated in S3, retain only the energy values ​​higher than the dynamic threshold, and generate a gated fluctuation energy sequence to suppress random noise interference; For the original fluctuation energy index sequence that has been normalized in S3, a dynamic noise threshold model (parameters are derived from S4.2) is used to realize an adaptive gating function based on the noise level under historical normal operating conditions. Furthermore, by comparing the current energy value with the dynamic threshold parameter, a point-by-point threshold determination method is used to screen the validity of the energy value, and an abnormal energy dataset that meets the threshold standard is obtained, where the threshold... The calculation formula is:

[0022] in, This represents the historical average energy fluctuation. The standard deviation of historical fluctuation energy. This is the noise tolerance factor; Furthermore, through the logic operation module (parameter: energy sequence) With threshold sequence Implement the conditional expression:

[0023] And generate the corresponding Boolean mask sequence. , used for gating and filtering operations; Furthermore, a pointwise multiplication operation is performed on the original energy sequence using a Boolean mask sequence to generate a gated fluctuation energy sequence. To suppress random noise components below the threshold; By using adaptive threshold determination and mask filtering, the original fluctuation energy from the previous step is transformed into a gated effective abnormal energy sequence, achieving high-fidelity preservation of weak and persistent abnormal signals and significantly reducing the false alarm rate under noise interference. For example, at key nodes of the stormwater drainage network, the historical average energy fluctuation... The value is 2.35, and the standard deviation is... The value is 0.47, which is the noise tolerance factor. The dynamic threshold T is set to 1.8, and the calculated dynamic threshold T = 2.35 + (1.8 × 0.47) = 3.196. For the energy sequence of the current window {2.91, 3.25, 4.02, 2.85}, point-by-point comparison reveals that the data points greater than the threshold are {3.25, 4.02}, generating a Boolean mask {0,1,1,0}. After mask filtering, the gated fluctuation energy sequence is {0,3.25, 4.02, 0}. This sequence effectively preserves two abnormal energy points exceeding the threshold, successfully triggering preliminary anomaly identification in subsequent change point detection. Performance testing shows that under the same noise background, this method significantly improves the detection sensitivity of weak early blockage signals, while significantly reducing false alarms caused by isolated noise peaks. S4.4: Timestamp the non-zero values ​​in the gated fluctuation energy sequence to generate a set of time-labeled abnormal candidate segments, wherein the timestamps are aligned with the time axis of the original flow time series data; S4.5: Perform minimum duration filtering on the set of abnormal candidate segments to remove isolated fluctuation events within a single window, thereby filtering out transient interference signals and retaining potential abnormal fluctuation segments with persistent characteristics.

[0024] Step S5: Accumulated trajectory construction is performed on the effective fluctuation energy index after gating. A change point detection algorithm is used to identify significant jumps in energy accumulation levels, generating a preliminary set of abnormal timestamps. The change point detection algorithm uses the Pelt algorithm for optimal segmentation. Specifically, this includes: S5.1: Perform cumulative trajectory construction on the effective fluctuation energy index after gating, calculate its cumulative sum based on the abnormal energy segments retained within the sliding window, and generate a continuous energy accumulation sequence as input data for the change point detection algorithm; For the effective fluctuation energy index sequence after gating, a sliding window accumulation algorithm (parameter: window length = typical blockage response period of the pipeline network) is used to calculate the sum of abnormal energy segments within the window; Furthermore, by using a weighted accumulation method (parameter: the weight coefficient is set after normalizing the energy amplitude of the anomalous segments), the contribution of different anomalous segments within the same window is weighted and summed to obtain the weighted cumulative energy value; Furthermore, a continuous sliding window stitching technique is adopted (parameter: window overlap rate is set to 50% to improve time resolution) to achieve the connection and smoothing of the cumulative energy values ​​of adjacent windows, so that the cumulative sequence maintains continuity on the time axis; Furthermore, by using a numerical integration algorithm (parameter: trapezoidal integration method), the energy sequence after gating is accumulated over the entire monitoring period, and a global energy trajectory curve is generated. By using the above sliding window accumulation and smoothing methods, abnormal energy fragments are transformed into continuous and stable energy accumulation sequences, thereby providing consistent and segmentable input features for subsequent change point detection algorithms. For example, at a pipeline monitoring point, the window length is set to 6 minutes, the window overlap rate is 0.5, and the weighting coefficient is calculated based on the segment energy amplitude / maximum window amplitude. The abnormal segment energy values ​​within the window are 2.4, 3.1, and 1.8, respectively. The weighted cumulative energy value is calculated as follows:

[0025] The cumulative value is 6.1, and the total energy is 54.8 after accumulating over 10 consecutive sliding windows using the trapezoidal integration method. The first half of the energy trajectory curve remains at a low level, and then rises significantly starting from the 7th sliding window, forming a jump feature that can be recognized by the Pelt algorithm. This enables highly sensitive capture of the energy pattern in the early stage of congestion, significantly improves the early warning recognition rate, and reduces misjudgments caused by transient noise. S5.2: The energy accumulation sequence is segmented into optimal segments based on the Pelt change point detection algorithm. The dynamic programming mechanism is used to search for the optimal set of breakpoints that minimizes the cost function in order to identify the time nodes when the energy level increases significantly. S5.3: Perform time alignment and denoising processing on the optimal breakpoint set, and filter out physically significant jump points based on the time interval threshold between adjacent breakpoints to remove false breakpoints caused by short-term disturbances or noise residues. S5.4: Generate a preliminary set of abnormal timestamps based on the selected significant jump points, and convert the time positions corresponding to the jump points into a unified system time format as candidate trigger signals for the risk warning process; For the set of significant rise points after S5.3 screening, a time location extraction method (parameter: the mapping relationship between the rise point index and the cumulative trajectory time axis) is used to realize the mapping from the breakpoint index to the actual observation time location; The original timestamp is converted from its storage format to the system's unified time format using a time format conversion algorithm (parameters: system unified time format definition, time zone, time precision), and time stamp data conforming to the global time base is obtained. Furthermore, through time validity verification processing (parameters: data source identifier, time axis range constraint), the validity of the converted timestamps is verified, and abnormal markers that are not within the monitoring period are removed to generate a compliant subset of timestamps; Furthermore, by using an index compression encoding method (parameters: monitoring point number, timestamp integer encoding rules), the timestamp and the spatial location of the monitoring point are synchronously associated and encoded, so that abnormal events can be quickly located when accessed by subsequent modules; Through the above-mentioned extraction of time location, unified format conversion, validity verification and encoding processing, the time location corresponding to the significant jump point is transformed into a preliminary abnormal timestamp set in a unified system time format, so as to achieve standardized output of candidate trigger signals for the risk warning process. For example, in a rainwater pipe network monitoring scenario, the set of significant jump points contains index values ​​[120, 250, 410]. The cumulative trajectory time axis sampling frequency is 1Hz, corresponding to original time positions of 120s, 250s, and 410s. The system's unified time format is defined as "YYYY-MM-DD HH:MM:SS", with a time zone of UTC+8 and time precision in seconds. 120s is converted to 2024-04-15 10:05:20, 250s to 2024-04-15 10:07:10, and 410s to 2024-04-15 10:09:50. Time validity verification confirms that all three are within the monitoring period "2024-04-15 10:00:00~2024-04-15 11:00:00", and all markers are retained. Furthermore, an integer encoding rule is adopted to combine the monitoring point number 002 with the timestamp to generate three encoded values, namely "002_20240415100520", which are output as a preliminary abnormal timestamp set. In this embodiment, the output set not only meets the global time reference requirements, but can also be directly called by the subsequent trend decay correlation verification module, significantly improving the speed of early warning event location and data consistency. S5.5: Perform context-related annotation processing on the preliminary abnormal timestamp set, associate each timestamp with the corresponding monitoring point number, energy jump amplitude and previous and next window state information to form a structured abnormal event record for subsequent verification module to call.

[0026] Step S6: Perform short-time linear fitting based on the trend principal component sequence, calculate the slope change rate within the current window, and generate a first-layer verification pass marker when the slope change rate is less than a preset negative threshold and overlaps with the initial anomaly timestamp. Specifically, this includes: S6.1: Perform sliding window truncation on the trend principal component sequence. Based on the principle of matching the sliding window length with the response delay of trend changes, truncate a trend data segment several time units before the current time point to obtain a local trend window sequence for subsequent short-time linear fitting analysis. S6.2: Apply the least squares method to perform short-time linear fitting on the local trend window sequence, and construct a regression model based on the linear relationship between the time variable and the trend value to obtain the trend slope parameter within the window, which serves as a quantitative indicator of the trend change rate. S6.3: Normalize the obtained trend slope parameters and calculate their standardized scores based on the slope distribution characteristics under historical normal operating conditions to eliminate scale bias caused by differences in pipeline structure between different monitoring points and obtain the standardized slope change rate index. S6.4: Based on a preset negative slope threshold, the standardized slope change rate index is judged and processed. When the standardized slope change rate is less than the preset negative threshold, it is determined that there is a trend decay phenomenon in the current window, and a trend decay marker signal is generated. Using the normalized trend slope change rate index sequence as the input data object, a low-threshold judgment analysis is performed to identify trend decay phenomena. Threshold comparison method is used (parameter: preset negative slope threshold). ), to achieve the rate of change of the current standardized slope The quantification is compared with the threshold, and it is determined whether the rate of change is lower than the set limit. Furthermore, through symbol discrimination logic (parameter: symbol function) This enables the classification and processing of trend change directions, ensuring that the attenuation detection link is triggered only when the rate of change shifts negatively; Furthermore, through the difference operation method (parameter: difference) This enables precise quantification of the difference between the rate of change and the threshold, and yields an offset metric to enhance the robustness of the judgment. The threshold out-of-bounds judgment algorithm (condition:) is used. This function enables the binarized determination and output of whether the current window trend change rate meets the decay condition, and generates a trend decay marker signal. Through this comparison link, the trend slope change rate calculation result of the previous step is transformed into a trend decay event indicator, realizing the first-order verification trigger condition for the trend change of blockage risk. For example, in the trend principal component sequence analysis of a certain stormwater pipe network node, the sliding window length is configured to be 12 hours, and the slope parameter of the current window is obtained by least squares fitting. for Historical data analysis sets a preset negative slope threshold. for Perform difference calculation. get ,satisfy The attenuation condition is determined, and the sign discrimination output is a negative trend shift. Therefore, the system generates a trend attenuation marker signal and records the node and its corresponding timestamp. Under different rainfall intensities (e.g., below 5 mm / h or above 20 mm / h), this method can accurately distinguish between the trend stabilization period and the attenuation period through the same judgment link, and provides effective event input for time alignment and overlap judgment in S6.5. Ultimately, it achieves accurate verification of early blockage trends and improves the system's risk response capability on a short timescale. S6.5: Perform time alignment and overlap judgment processing on the trend decay marker signal and the preliminary abnormal timestamp set. Based on the time window matching algorithm, identify whether there is time overlap between the two. If there is overlap, generate the first layer of verification pass mark as the preliminary basis for early congestion risk verification. Time series alignment processing is performed on the trend decay marker signal and the preliminary anomaly timestamp set, using a unified system time format (parameter: millisecond precision) to achieve accurate matching of signals from different sources under the same time base; Furthermore, by using a time window matching algorithm (parameters: window length, matching tolerance), the occurrence time of the trend decay marker signal and the preliminary abnormal timestamp set are overlapped in intervals, and an overlap matching matrix is ​​generated as the basis for judgment. Furthermore, the set intersection operation method (parameter: overlap judgment threshold is 50% of the window length) is used to extract overlapping time segments and obtain a subset of timestamps that meet the time overlap condition; Furthermore, an event validity filtering algorithm (parameter: minimum overlap duration threshold) is applied to the time overlap results to remove non-persistent overlap events caused by transient fluctuations and retain time overlap signals with clear physical meaning. Through logical AND operation processing, the trend decay marker signal is combined with the effective time overlap subset to transform it into the first-level verification pass marker, realizing the preliminary judgment of early blockage risk verification; For example, at a certain rainwater pipe network monitoring node, the normalized rate of change of the local window slope of the trend principal component sequence is: Below the preset negative threshold A trend decay marker signal was generated. The timestamp of this signal was 2024-05-20 14:32:15, which was 3 seconds apart from a timestamp in the initial abnormal timestamp set, 2024-05-20 14:32:18. The matching tolerance was set to 5 seconds, meeting the time overlap condition. The matching window length was 10 seconds, and the overlap duration reached 6 seconds, exceeding the persistence threshold of 5 seconds; therefore, the event was retained. Finally, this event was marked as a first-level verification pass marker, and the event node number WX-102 and the trend slope change rate were recorded in the system log. The continuous overlap time of 6 seconds enabled the effective identification of early blockage risk at the monitoring point.

[0027] Step S7: Establish a Gaussian mixture model based on historical normal energy distribution, calculate the statistical significance score of the current cumulative energy, and generate a second-layer validation pass marker when the significance score of three consecutive windows exceeds the 95% confidence interval. Specifically, this includes: S7.1: Perform data preprocessing on the residual fluctuation energy index sequence under historical normal operating conditions, including missing value imputation and standardization, to obtain a high-quality sample set that can be used for modeling; S7.2: Based on the preprocessed residual fluctuation energy index sample set, the Bayesian information criterion (BIC) is used to determine the optimal number of Gaussian mixture components in order to construct a statistical model with sufficient expressive power and without overfitting. S7.3: The parameters of the Gaussian mixture model are iteratively estimated using the expectation-maximization (EM) algorithm to obtain the mean, covariance matrix and mixing weights of each Gaussian component, so as to characterize the multimodal distribution characteristics of normal energy fluctuations. S7.4: For the residual fluctuation energy accumulation trajectory sequence after the current gating process, calculate the log likelihood score of the energy observation value of each window based on the Gaussian mixture model to quantify its deviation from the historical normal pattern; S7.5: Convert the log-likelihood score into a statistical significance score, and determine whether the significance scores of three consecutive windows all exceed the 95% confidence interval. If the condition is met, generate a second-level verification pass mark to realize the statistical confidence verification of the abnormal signal. Based on the log-likelihood scores of the energy observations for each window calculated in step S7.4, a distribution significance transformation method (parameters: mean, variance, and component weights of the Gaussian mixture model) is used to map the log-likelihood scores to statistical significance scores. Furthermore, the tail probability of each observation in the Gaussian mixture model is calculated using the cumulative distribution function of the standard normal distribution. The tail probability is used as a measure of the significance score, and the significance score sequence corresponding to each sliding window is obtained. Furthermore, an interval threshold comparison method (parameter: upper limit of confidence interval = 95%) is adopted to perform window-by-window judgment on the significance score sequence and generate a set of window labels that exceed the confidence interval; Furthermore, using the continuous window matching algorithm (parameter: number of continuous windows = 3), the above-mentioned set of labels is used to perform sequence filtering, retaining only the event sequences whose significance scores exceed the 95% confidence interval for three consecutive windows, and generating a second-level validation pass label; By using significance transformation and continuous exceedance judgment processing, the log-likelihood score result of the previous step is transformed into an anomaly confirmation label that meets dual conditions, thereby realizing the statistical confidence verification of the anomalous signal.

[0028] in, The log-likelihood score for the i-th window. The mean parameters of the Gaussian mixture model are... The standard deviation parameter of the Gaussian mixture model. The cumulative distribution function of the standard normal distribution. This corresponds to the statistical significance score; For example, in the residual fluctuation energy accumulation trajectory analysis of a monitoring point in a stormwater pipe network, the historical normal energy distribution was modeled to obtain μ=2.35 and σ=0.45. The significance scores for the current five consecutive windows were calculated using the above formula, resulting in scores of 0.91, 0.97, 0.98, 0.94, and 0.96, respectively. The upper limit of the confidence interval was set to 0.95, and the threshold for the number of consecutive windows was set to 3. Through threshold comparison and continuous matching, the second to fourth windows were identified as consecutive abnormal sequences that met the conditions, generating a second-layer verification pass marker. This marker was combined with the first-layer verification pass marker in the subsequent S8 step to form the final early blockage risk warning event. The verification results show that the sensitivity for identifying the early stage of blockage was significantly improved under this operating condition, while maintaining a low false alarm level. Step S8: When both the first-level verification pass marker and the second-level verification pass marker are obtained simultaneously, an early congestion risk warning event containing a timestamp is generated, triggering the pipeline maintenance scheduling process. The warning event includes a risk level identifier and an associated monitoring point number. Specifically, it includes: S8.1: Perform time alignment and logical AND operation on the received first-level verification pass markers and second-level verification pass markers, and identify abnormal event nodes that simultaneously meet the two-level verification conditions based on the timestamp matching mechanism to form a joint verification pass event set; S8.2: For each event node in the joint verification pass event set, construct a structured early warning data packet based on the event occurrence timestamp and the associated monitoring point number. The structured early warning data packet includes the event timestamp, monitoring point ID, verification pass flag type, and current trend slope information. S8.3: Based on the preset risk level mapping rules, the trend slope information in the structured early warning data packet is processed by level mapping to generate the corresponding risk level identifier. The risk level identifier includes three categories: low risk, medium risk and high risk, which are used to characterize the severity of the congestion risk. Based on the trend slope information in the structured early warning data package, a risk level mapping rule table (parameter: the relationship between the slope change rate range and the risk level) is used to realize the quantitative correlation and conversion between the degree of trend decay and the risk category. Furthermore, by using an interval determination algorithm (parameters: low risk, medium risk, and high risk slope thresholds), the slope values ​​are matched with preset intervals to obtain the corresponding primary risk label data. Furthermore, through weight correction processing (parameter: historical anomaly intensity weight coefficient of monitoring points), the primary risk label is corrected and adjusted, and a corrected risk label reflecting the combined impact of current working conditions and historical patterns is generated. Furthermore, through a discretization coding algorithm (parameter: risk category enumeration mapping table), the coding conversion from corrected risk labels to system risk identifiers is realized, and data fields conforming to the system's internal communication protocol are generated; By using risk level mapping rules, the results of the previous step are transformed into structured risk level identifiers, enabling accurate representation and standardized output of the severity of congestion risk. For example, the joint validation performs risk level mapping on the trend slope information of a certain monitoring point in the event set, and sets the low-risk threshold range as the slope change rate. to The medium-risk threshold range is to The high-risk threshold is less than The rate of change of the input slope is... The interval determination algorithm matched the data to a medium-risk range, resulting in a preliminary risk label of "medium risk." This was combined with the weighting coefficients of three abnormal energy surge events that occurred at this monitoring point within the past six months. The weighted adjustment process raises the risk level to high risk. The discretization coding algorithm maps the high-risk identifier to an internal coded value. The final structured data field {Risk Level Identifier: High Risk, Encoded Value: 3} is generated and used to assemble subsequent early warning event data packets. Based on this, the early warning system prioritizes dispatching emergency response teams in pipeline maintenance scheduling, significantly improving the timeliness of handling early blockage events. S8.4: Send the early warning event data packet containing the timestamp, monitoring point number and risk level identifier to the pipeline maintenance and dispatch system through the communication interface module to trigger the corresponding early warning response mechanism.

[0029] S8.5: Record and archive all generated early congestion risk warning events, and build a warning event log database based on event timestamps and risk level identifiers to provide data support for subsequent system performance evaluation and risk evolution analysis.

[0030] The technical solution of the present invention has been described above with reference to the preferred embodiments shown in the accompanying drawings. However, it will be readily understood by those skilled in the art that the scope of protection of the present invention is obviously not limited to these specific embodiments. Without departing from the principles of the present invention, those skilled in the art can make equivalent changes or substitutions to the relevant technical features, and the technical solutions after these changes or substitutions will all fall within the scope of protection of the present invention.

[0031] The above description is merely a preferred embodiment of the present invention and is not intended to limit the invention. Various modifications and variations can be made to the present invention by those skilled in the art. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and rules of the present invention should be included within the scope of protection of the present invention.

Claims

1. A method for early warning of rainwater harvesting pipe network blockage risk, characterized in that, Includes the following steps: S1: Collect time-series flow data and environmental parameters from rainwater pipe network monitoring points to construct the initial input dataset; S2: Perform multi-scale fluctuation decomposition processing on the traffic time series data, decomposing it into trend principal component sequence, periodic component sequence and residual fluctuation component sequence; S3: Apply sliding window energy calculation to the residual fluctuation component sequence to generate the original fluctuation energy index; S4: Based on the statistical characteristics of residual fluctuations under historical normal operating conditions, a dynamic noise threshold model is established, and the original fluctuation energy index is gating processed based on the dynamic noise threshold model. S5: Perform cumulative trajectory construction on the effective fluctuation energy index after gating, use the change point detection algorithm to identify significant jump nodes in the energy accumulation level, and generate preliminary abnormal timestamps; S6: Perform short-time linear fitting based on the trend principal component sequence, calculate the slope change rate within the current window, and generate a first-layer verification pass marker when the slope change rate is less than a preset negative threshold and overlaps with the initial anomaly timestamp; S7: Establish a Gaussian mixture model based on the historical normal energy distribution, calculate the statistical significance score of the current cumulative energy, and generate a second-level validation pass mark when the significance score of three consecutive windows exceeds the confidence interval; S8: When the first layer verification pass flag and the second layer verification pass flag are obtained simultaneously, an early blockage risk warning event is generated, triggering the pipeline maintenance scheduling process.

2. The method for early warning of rainwater harvesting pipe network blockage risk according to claim 1, characterized in that, Step S1 specifically includes: Periodically collect data from multiple sets of flow sensors deployed at key nodes of the rainwater pipe network, obtain the instantaneous flow value of each monitoring point based on the Modbus protocol, and form continuous time series data containing multiple deployment locations; Synchronous sampling is performed on water level sensors deployed at pipeline junctions and key inspection wells, and real-time water level data is obtained based on the ADC analog-to-digital conversion module; Real-time rainfall intensity data is obtained based on the meteorological API interface. The real-time rainfall intensity data is synchronized with the instantaneous flow rate and water level data using a timestamp alignment mechanism to form a multi-dimensional input feature set under a unified time reference. The instantaneous flow rate value, the water level height data, and the real-time rainfall intensity data are subjected to sliding window mean filtering to eliminate instantaneous noise interference based on the moving average algorithm; Standardization processing is performed on the multi-source data after timestamp alignment. The Z-score normalization method is used to perform zero-mean unit variance transformation on the data of each dimension to generate a structured initial input dataset.

3. The method for early warning of rainwater harvesting pipe network blockage risk according to claim 2, characterized in that, Step S1 further includes parsing the Modbus register address mapping and byte order conversion to restore the physical flow value and bind it to the UTC timestamp. Missing data is filled in by linear interpolation and nearest neighbor interpolation and then spliced ​​into a spatially synchronized multidimensional time series.

4. The method for early warning of rainwater harvesting pipe network blockage risk according to claim 1, characterized in that, Step S2 specifically includes: The raw traffic time series data is preprocessed, and a sliding window mean filtering algorithm is used to perform noise suppression processing on the raw traffic time series data to obtain a smoothed input signal sequence. Based on the empirical mode decomposition algorithm, local extremum point identification processing is performed on the smoothed input signal sequence, and upper and lower envelopes are constructed using cubic spline interpolation to obtain the signal envelope features within the local extremum region. The mean of the upper and lower envelopes is calculated to generate a local mean sequence. The smoothed input signal sequence is subtracted from the local mean sequence to obtain a candidate intrinsic mode function sequence, so as to determine whether it meets the intrinsic mode function criterion. If the candidate intrinsic mode function sequence does not satisfy the intrinsic mode function criterion, then iterative screening is performed on it, and the local extremum identification and mean calculation process is repeated until the intrinsic mode function definition is satisfied, and the first-order intrinsic mode function components are obtained. The first-order intrinsic mode function component is separated from the original signal to obtain the remaining signal sequence. The empirical mode decomposition process is repeated to extract multiple intrinsic mode function components in sequence. Finally, the high-frequency components are classified as residual fluctuation components, and the low-frequency components are classified as trend principal components and periodic components.

5. The method for early warning of rainwater harvesting pipe network blockage risk according to claim 4, characterized in that, Step S2 further includes obtaining the upper and lower envelopes using cubic spline interpolation, calculating the mean and subtracting it from the original signal to obtain candidate intrinsic mode function components, and using the number of extreme points and zero crossover points and the mean approaching zero as criteria. The specific criteria are: the number of extreme points and the number of zero crossover points are equal or differ by no more than one, and the local mean approaches zero. If these conditions are not met, the process is iterated until a qualified intrinsic mode function is obtained.

6. The method for early warning of rainwater harvesting pipe network blockage risk according to claim 1, characterized in that, Step S3 specifically includes: A sliding window partitioning process is performed on the residual fluctuation component sequence to obtain local residual fluctuation subsequences within multiple time windows; For each time window, a local mean is calculated for the local residual fluctuation subsequence, and a local benchmark value is generated based on the arithmetic mean of the data points within the sliding window. For each local residual fluctuation subsequence, the deviation squared is calculated between the local reference value and the corresponding local reference value. Based on the point-by-point deviation squared value, the local perturbation energy density of each point in the window is generated. The local disturbance energy density within each window is summed, and the original fluctuation energy index for that window is generated based on the sum of the squared deviations of all points within the window. The original fluctuation energy indices of all windows are normalized. Standardization is performed based on the mean and standard deviation of historical data to obtain an energy characteristic sequence at a uniform scale.

7. The method for early warning of rainwater harvesting pipe network blockage risk according to claim 6, characterized in that, The window length in the sliding window division is set based on the typical blockage response cycle of the pipeline network.

8. The method for early warning of rainwater harvesting pipe network blockage risk according to claim 1, characterized in that, Step S4 specifically includes: A sliding window statistical analysis is performed on the historical residual fluctuation component sequence, and a historical fluctuation energy index is generated based on the sum of squared deviations of data points within the window from the local mean. The historical fluctuation energy index is normalized, and a dynamic noise threshold model is constructed based on the mean and standard deviation of historical data. Based on the dynamic noise threshold model, threshold gating is performed on the current original fluctuation energy index to retain only the energy values ​​higher than the dynamic threshold, thereby generating a gated fluctuation energy sequence. The non-zero values ​​in the gated fluctuation energy sequence are timestamped to generate a set of time-tagged abnormal candidate segments; The minimum duration filtering is performed on the set of abnormal candidate segments to remove isolated fluctuation events within a single window and retain potential abnormal fluctuation segments with persistent characteristics.

9. The method for early warning of rainwater harvesting pipe network blockage risk according to claim 1, characterized in that, Step S5 specifically includes: The effective fluctuation energy index after gating is used to construct an accumulated trajectory. The sum of the abnormal energy segments retained within the sliding window is calculated to generate a continuous energy accumulation sequence. The energy accumulation sequence is segmented and optimally divided based on the Pelt change point detection algorithm, and the optimal set of breakpoints that minimizes the cost function is searched using a dynamic programming mechanism. Time alignment and denoising are performed on the optimal breakpoint set, and significant jump points are selected based on the time interval threshold between adjacent breakpoints. Based on the selected significant jump points, a preliminary set of abnormal timestamps is generated. The preliminary set of abnormal timestamps is subjected to context-related annotation processing, which associates each timestamp with the corresponding monitoring point number, energy jump magnitude, and previous and subsequent window state information to form a structured abnormal event record.

10. A method for early warning of rainwater harvesting pipe network blockage risk according to claim 9, characterized in that, Step S5 further includes smoothing the energy accumulation trajectory through a sliding window, weighting, and numerical integration, then using the Pelt change point detection algorithm to identify nodes with significant energy jumps, and mapping the jump point index to a unified system timestamp.

Images