Large pipe network measurement anomaly whole chain closed loop tracking early warning method

Abstract

The present application belongs to the technical field of pipe network anomaly monitoring, and particularly relates to a large-scale pipe network metering anomaly full-chain closed-loop tracking and early warning method, which comprises the following steps: identifying an initial abnormal metering node of a pipe network and extracting a multi-dimensional time-frequency feature vector, thereby matching an influence domain in a model library and generating a candidate tracking chain with an initial weight along a topology; subsequently, a propagation model is constructed for each chain, a time-space search window is set, and parameters are corrected and chain confidence is updated by comparing real-time data with predicted state variables; through iterative verification, if the confidence of a single chain exceeds a convergence threshold, the chain is confirmed as a main propagation path and an early warning is given; if multiple chains exceed a background threshold, an abnormal source area is located by joint reverse tracing, and the consistency of the abnormal source area and related node monitoring data is verified by using forward transient simulation to complete the closed loop. The present application can effectively track an abnormal propagation path, and improve the response speed and positioning reliability of abnormal early warning of a large-scale complex pipe network.

Description

Technical Field

[0001] This invention relates to the field of pipeline network anomaly monitoring technology. More specifically, this invention relates to a closed-loop tracking and early warning method for metering anomalies in large-scale pipeline networks. Background Technology

[0002] Water supply networks are core infrastructure ensuring the normal operation of cities, and their operational status directly affects public safety and water resource utilization efficiency. Due to the intricate interconnected topology and complex and variable hydraulic conditions of the network, when abnormal events such as leaks, pipe bursts, or meter malfunctions occur, the resulting transient hydraulic disturbances can rapidly spread along the network to surrounding nodes, triggering widespread cascading alarms. This spatiotemporal propagation effect makes the monitoring data exhibit significant non-stationary characteristics. Traditional early warning methods based on absolute thresholds of single nodes or simple time series predictions struggle to separate the true abnormal state from background noise and normal user water usage fluctuations. This leads to managers facing a common dilemma when dealing with sudden anomalies: multiple concurrent alarms, difficulty in tracing the source, and high rates of false alarms and missed alarms.

[0003] To address these issues, the industry has gradually developed a technical approach that combines pipeline topology analysis with data-driven technology and hydraulic models. This approach uses monitoring time-series data from upstream and downstream nodes to trace and locate anomalies. For example, Chinese patent document CN113446521B discloses a method for locating pipe bursts based on transient flow. This method collects pipeline pressure data through a high-frequency SCADA system, filters it, and then uses a dual discrimination mechanism of system anomalies and transient flow anomalies to identify pipe burst events. Based on the momentum theorem and the operating characteristics of ball valves, it calculates the transient flow signal of the burst pipe and establishes a transient flow amplitude calculation model for monitoring points by combining the node transmission and reflection coefficients and friction losses. By setting no fewer than three candidate burst points in each pipe section for numerical calculation, the method minimizes the error between the calculated and measured signals to locate the burst pipe section and its position. This method has been validated in a small experimental pipeline network, providing a feasible approach for pipeline anomaly detection based on transient flow.

[0004] However, the aforementioned existing technologies still have the following shortcomings in practical engineering applications: Their global traversal calculation approach requires setting burst points for each pipe segment in the network and performing complete transient flow calculations. The computational load increases linearly with the number of pipe segments, but in large municipal networks containing tens of thousands of pipe segments, the total computational load still reaches an unacceptable level, failing to meet real-time early warning requirements. In the cross-node matching process of abnormal signals, the cumulative time delay evolution characteristics of hydraulic propagation and the instrument measurement errors of each monitoring node are not fully considered, making it prone to timing misalignment or noise interference leading to broken or mismatched tracking chains. Furthermore, they use fixed model parameters for transient flow calculations and lack a mechanism for online correction of key parameters such as pipe friction coefficients and node demand using real-time monitoring data. When there is a deviation between the actual operating conditions of the network and the model's preset parameters, the positioning accuracy will significantly decrease.

[0005] Furthermore, the existing technology only supports the location and processing of a single pipe burst event. When there are multiple suspected abnormal propagation paths in the pipeline network, it lacks a closed-loop confirmation mechanism for multi-path joint reverse tracing and forward verification. It cannot effectively distinguish between real abnormal sources and hydraulic interference signals, nor can it handle non-leakage anomalies such as meter damage. The false alarm rate and missed alarm rate under complex working conditions are difficult to meet engineering requirements. Summary of the Invention

[0006] To address the shortcomings of existing transient flow-based pipeline anomaly localization methods in large municipal pipeline networks, such as insufficient computational efficiency, susceptibility to errors in cross-node matching of anomaly signals, inability to adaptively adjust model parameters online, lack of multi-path closed-loop verification mechanisms, and inability to cover non-leakage metering anomalies like meter damage, this invention provides a closed-loop tracking and early warning method for the entire chain of metering anomalies in large pipeline networks.

[0007] This invention provides a closed-loop tracking and early warning method for metering anomalies in a large-scale pipeline network, comprising: S1, identifying initial abnormal metering nodes in the pipeline network based on pressure and flow time-series monitoring data of each metering node, and extracting their multidimensional time-domain and frequency-domain feature vectors; using the multidimensional time-domain and frequency-domain feature vectors to match the propagation influence domain from the pipeline network transient response model library, and generating candidate tracking chains with initial weights by searching along the pipeline network topology within the propagation influence domain; S2, constructing a propagation model for each candidate tracking chain, and setting the spatiotemporal search time of the next node to be verified based on the cumulative hydraulic propagation time delay and the monitoring uncertainty of the node to be verified. Search window; compare real-time monitoring data with predicted state quantities within the spatiotemporal search window, correct propagation model parameters using multidimensional first-order residual vectors, and update the posterior probability of each candidate tracking chain as confidence level; S3, iteratively execute the steps of setting the spatiotemporal search window to update confidence level. When the confidence level of a single candidate tracking chain exceeds the convergence threshold, it is confirmed as an abnormal main propagation path and an early warning is issued; if no single chain converges but the confidence levels of multiple candidate tracking chains are higher than the background threshold, then jointly reverse trace to locate the abnormal source area, and use forward transient simulation to verify the consistency between the abnormal source area and the monitoring data of related nodes to complete the closed loop.

[0008] This invention identifies initial abnormal metering nodes in the pipeline network and extracts their time and frequency domain features. It then matches the propagation influence domain from a model library to generate candidate tracking chains with initial weights. For each chain, a propagation model is constructed, and a spatiotemporal search window is set based on the hydraulic propagation time delay and node monitoring uncertainty. Within the window, real-time data is compared with model predictions, and residual vectors are used to correct model parameters and update the chain's posterior probability as confidence. The above steps are iteratively executed: if the confidence of a single chain reaches the convergence threshold, it is confirmed as the main propagation path of the anomaly and an early warning is issued; if the confidence of multiple chains exceeds the background threshold but none converge, the anomaly source area is located through joint reverse tracing, and its consistency with monitoring data is verified through forward transient analysis, forming a closed loop. This method can effectively track anomaly propagation paths, achieve source location, and reduce false alarms and missed alarms.

[0009] Preferably, the identification of initial abnormal metering nodes in the pipeline network includes: inputting the acquired pressure and flow time-series monitoring data into an isolated forest algorithm, calculating the abnormal score of each data point, and extracting data points with abnormal scores greater than a preset abnormal threshold as first candidate abnormal points; establishing a sliding window with a preset time step to intercept the pressure and flow time-series monitoring data, performing quantile regression calculation on the intercepted data to obtain the prediction interval boundary, and extracting data points outside the prediction interval boundary as second candidate abnormal points; calculating the intersection of the first candidate abnormal point and the second candidate abnormal point, and confirming the node and time position corresponding to the data point belonging to the intersection as the initial abnormal metering node and the time of abnormal occurrence.

[0010] By combining the global anomaly detection capability of isolated forest with the local trend analysis capability of sliding window quantile regression, and employing a dual discrimination mechanism to screen outomas and take their intersection, we can effectively filter out background noise and interference from normal water usage fluctuations, reduce the false alarm rate of initial anomaly identification, and accurately determine the nodes and specific times when anomalies occur.

[0011] Preferably, the extraction of its multidimensional time-domain and frequency-domain feature vectors includes: taking the time of the anomaly occurrence as the center, extracting continuous monitoring data as target signal segments by offsetting a preset time before and after the time of the anomaly occurrence; calculating the root mean square value, skewness, and kurtosis of the target signal segments and concatenating them in order to form a one-dimensional time-domain feature matrix; performing a fast Fourier transform on the target signal segments to convert them to the frequency domain, calculating the frequency domain amplitude spectrum to extract the main frequency amplitude and the total harmonic distortion rate of the harmonics of a preset order, forming a one-dimensional frequency-domain feature matrix; and concatenating and fusing the one-dimensional time-domain feature matrix with the one-dimensional frequency-domain feature matrix to generate the multidimensional time-domain and frequency-domain feature vectors.

[0012] By extracting the time-domain statistical features and frequency-domain harmonic features before and after the abnormal signal, the signal characteristics of hydraulic transient disturbances can be comprehensively characterized. The time-domain features reflect the energy distribution and morphological distortion of the signal, while the frequency-domain features reveal the frequency composition of the signal. The fusion of the two can significantly improve the matching accuracy with the pipeline transient response model library and accurately delineate the influence domain of abnormal propagation.

[0013] Preferably, the process of generating candidate tracing chains with initial weights by searching the pipeline topology within the propagation influence domain includes: setting the initial anomaly metering node as the source node of the traversal tree, and sequentially connecting adjacent nodes along the water flow to generate multiple candidate tracing chains; obtaining the pipe inner diameter, pipe length, and Manning roughness coefficient values ​​of the water transmission pipe segments associated with each candidate tracing chain; calculating the single-segment conductivity index of each water transmission pipe segment based on the pipe inner diameter, pipe length, and Manning roughness coefficient values, and summing the single-segment conductivity indices of all pipe segments within a single candidate tracing chain to obtain the corresponding conductivity index magnitude value; calculating the proportion of the conductivity index magnitude value of each candidate tracing chain in the sum of all chains, and assigning it as the initial weight.

[0014] Preferably, the step of constructing a propagation model for each candidate tracking chain includes: sequentially loading the pressure monitoring node data contained in all candidate tracking chains to establish a two-dimensional coordinate basis system for the pipeline network; using the continuity equation and momentum equation of one-dimensional unsteady flow, calculating the wave velocity propagation coefficient and damping hysteresis matrix within the corresponding water transmission pipeline network to construct a state transition space equation matrix; based on the preset grid differential node spacing and time step differential step size constraints of the state transition space equation matrix, discretizing the state transition space equation matrix into a Jacobian difference algebraic formula system to generate a one-dimensional transient propagation model for outputting predicted state variables.

[0015] Preferably, the step of setting the spatiotemporal search window for the next node to be verified based on the cumulative hydraulic propagation time delay and the monitoring uncertainty of the node to be verified includes: accumulating the pipe length values ​​contained in the candidate tracking chain to obtain the total pipe length; dividing the total pipe length by the theoretical propagation rate of water hammer elastic waves to establish the cumulative hydraulic propagation time delay; constructing a diagonal error variance standard matrix to extract the monitoring uncertainty of the node to be verified based on the accuracy level and calibration range value of the metering equipment of the next node to be verified; establishing the time interval boundary by extending a preset time forward and backward with the actual running time corresponding to the cumulative hydraulic propagation time delay as the axis time; and expanding the state space limit range based on three times the standard deviation of the monitoring uncertainty of the node to be verified; and establishing the spatiotemporal search window by combining the time interval boundary and the limit range.

[0016] By comprehensively considering the time delay characteristics of hydraulic propagation and the measurement error of monitoring equipment, search boundaries are set for the time dimension and the state quantity dimension respectively, and combined to form a spatiotemporal search window. This can effectively avoid abnormal signal matching errors caused by time misalignment and noise interference, and improve the accuracy of cross-node abnormal propagation signal matching.

[0017] Preferably, the step of using a multidimensional first-order residual vector to correct the propagation model parameters and update the posterior probability of each candidate tracking chain as the confidence level includes: subtracting the predicted state value at this moment from the real-time monitoring data of the next node to be verified, retaining the sign of the difference to generate a multidimensional first-order residual vector; substituting the multidimensional first-order residual vector into an evaluation function that conforms to a zero-mean multidimensional Gaussian distribution to extract the corresponding monotonic likelihood correlation evaluation score; multiplying the monotonic likelihood correlation evaluation score with the forward probability stored value of the current candidate tracking chain at the previous time step to obtain the probability level at the current time step; and normalizing the calculation by dividing the probability level of the current candidate tracking chain itself by the sum of the probability levels of all parallel chains, and updating the resulting value as the confidence level of the candidate tracking chain.

[0018] By calculating the residual between real-time monitoring data and model predictions, and combining it with Bayesian inference, the confidence of candidate chains is dynamically updated. The smaller the residual, the more accurate the model prediction, and the higher the confidence of the corresponding chain. This allows the confidence of real anomaly propagation paths to gradually increase, while the confidence of erroneous paths to decay rapidly, thus improving the reliability of path identification.

[0019] Preferably, the joint reverse tracing and locating of the anomaly source area includes: reading the real geographic longitude and physical latitude values ​​of all candidate tracing chain nodes whose confidence continuously exceeds the background threshold, and storing them as a discrete two-dimensional geospatial positioning point set; using centripetal clustering to calculate the planar centroid coordinates of the tangent polygon network of the discrete two-dimensional geospatial positioning point set, and projecting the planar centroid coordinates onto the nearest pipeline segment to establish the core source coordinates of the anomaly location; using the core source coordinates of the anomaly location as the geometric center, extending a preset absolute length reference distance in all directions to construct a boundary layer block of a completely closed boundary circle covering the area, which serves as the anomaly source area.

[0020] Preferably, the step of using forward transient simulation to verify the consistency between the anomaly source region and the monitoring data of related nodes to complete the closed loop includes: configuring the nodes in the anomaly source region as newly added leak holes or nodes with sudden change requirements in the propagation model; running forward transient simulation to generate simulated pressure and flow time series curves for all related nodes in the network; calculating the mean square error between the simulated pressure and flow time series curves and the actual real-time monitoring data of related nodes; when the mean square error is lower than the preset basic measurement noise tolerance, confirming that the reverse tracing location is accurate, generating an anomaly diagnosis report to complete the closed loop.

[0021] By setting up virtual anomaly nodes in the located anomaly source area to conduct forward transient simulation, the simulation results are compared and verified with actual monitoring data. If the error between the two is within an acceptable range, the positioning accuracy is confirmed, forming a complete closed loop of reverse tracing and forward verification, further reducing the false alarm rate and ensuring the reliability of the positioning results.

[0022] Preferably, the single-segment conductivity index of each water transmission pipe segment is calculated based on the pipe inner diameter value, pipe length value, and Manning roughness coefficient value, including: multiplying the fourth power of the pipe inner diameter value of the corresponding water transmission pipe segment by the value of pi as the dividend, multiplying the pipe length value of the corresponding water transmission pipe segment by the value of Manning roughness coefficient as the divisor, and dividing the dividend by the divisor to obtain the single-segment conductivity index.

[0023] The technical solution of the present invention has the following beneficial technical effects: This invention employs a dual discriminant method of isolated forest and sliding window quantile regression to identify initial abnormal metering nodes. It uses multi-dimensional time-frequency features to match a pipeline transient response model library to delineate the propagation influence domain, generating candidate tracking chains with initial weights. A spatiotemporal search window is set by combining hydraulic propagation time delay and monitoring uncertainty. Multi-dimensional first-order residuals are used to correct propagation model parameters online and dynamically update the confidence of each path. A complete closed loop is formed through iterative verification or joint reverse tracing and forward verification. This scheme effectively narrows the anomaly search range, reduces computational burden, achieves accurate cross-node matching of anomaly signals, dynamically adapts to changes in actual pipeline operating conditions, and supports handling anomaly scenarios with multiple paths coexisting. This overcomes the shortcomings of existing technologies, such as excessive computational load hindering real-time early warning, susceptibility to operating condition deviations in positioning accuracy, lack of closed-loop verification mechanisms, and insufficient anomaly coverage. It improves the response speed and positioning reliability of anomaly early warning in large and complex pipeline networks. Attached Figure Description

[0024] Figure 1 A flowchart for a closed-loop tracking and early warning method for metering anomalies in large-scale pipeline networks; Figure 2 A schematic diagram for extracting the dominant frequency and harmonic amplitude of the target signal in the frequency domain; Figure 3 A schematic diagram showing the initial weighting of the hydraulic fluid conduction index for candidate tracking chains; Figure 4 This is a comparison chart of the anomaly identification accuracy and spatial positioning error for each module. Detailed Implementation

[0025] The technical solutions of the present invention will be clearly and completely described below with reference to the accompanying drawings of the embodiments of the present invention. Obviously, the described embodiments are some embodiments of the present invention, but not all embodiments.

[0026] The specific embodiments of the present invention will now be described in detail with reference to the accompanying drawings.

[0027] This invention discloses a closed-loop tracking and early warning method for metering anomalies in large-scale pipeline networks, referring to... Figure 1 This includes the following steps: S1. Identify initial abnormal nodes and extract features to generate candidate chains.

[0028] This step acquires continuous time series pressure and flow data of each metering node in the pipeline network through the data acquisition and monitoring control system. First, the missing values ​​are filled in using the piecewise cubic Hermit interpolation method to ensure the continuity and smoothness of the time series data. Then, four sub-steps are completed in sequence: anomaly detection, multi-dimensional time-frequency feature extraction, propagation influence domain matching, and candidate tracking chain generation. Finally, a set of candidate tracking chains with initial weights is output.

[0029] In the outlier detection stage, a strategy of parallel processing of the Isolation Forest algorithm and sliding window quantile regression is adopted, and the reliability of anomaly identification is improved through a double cross-validation mechanism. On the one hand, an Isolation Forest model is constructed, and the downsampling sample size of each decision tree and the total number of decision trees are set. After each data point is mapped by the model, an anomaly score between 0 and 1 is output based on its average path length in the binary tree. When the anomaly score exceeds a preset anomaly threshold, the data point is included in the first candidate anomaly set. On the other hand, a sliding window of data containing a preset time step is extracted at a fixed sampling interval. A regression class is called and combined with a quantile objective function to fit the conditional distribution functions corresponding to the lower and higher percentile quantiles, respectively, thereby obtaining the upper and lower bounds of the prediction interval. Data points falling outside the prediction interval are included in the second candidate anomaly set.

[0030] Finally, a Boolean AND operation is performed on the two types of candidate points to find their intersection. As a preferred solution, if the two types of candidate points within the same time window have a preset allowable error in their timestamps (such as an offset of 1 to 2 time steps), they can also be considered to meet the intersection condition. The metering node IDs and timestamps corresponding to the data points in the intersection are then identified as the initial abnormal metering nodes and the time when the abnormality occurred. This dual-channel parallel verification strategy leverages the sensitivity of isolated forests to global sparse anomalies and utilizes the adaptability of quantile regression to local operating condition fluctuations, effectively separating real anomaly signals from background noise and normal water usage fluctuations.

[0031] In this embodiment, the isolated forest model contains 100 decision trees, the downsampling size is set to 256, and the preset anomaly threshold is 0.75; the sampling interval is 1 second, the sliding window contains 200 time steps, and the quantile objective function is fitted to the conditional distribution functions corresponding to the 0.05 quantile and 0.95 quantile, respectively; assuming that at a certain moment the pressure of node A suddenly drops to 0.2 MPa, the anomaly score calculated by the isolated forest model is 0.82, the reasonable pressure range fitted by the sliding window quantile regression model is [0.28, 0.45] MPa, and the measured value of 0.25 MPa is lower than the lower bound. The two results intersect at the same node at the same time, so node A is locked as the initial anomaly measurement node, and the corresponding timestamp is marked as the time of anomaly occurrence. .

[0032] To comprehensively characterize the hydraulic physical properties of anomalous events, the time of occurrence of the anomaly is used as the starting point. Using a center of symmetry, continuous monitoring data is extracted at preset time intervals to obtain target signal segments. In the time domain, the root mean square value, skewness, and kurtosis of each signal segment are calculated sequentially and concatenated into a one-dimensional time-domain feature matrix. In the frequency domain, a Hamming window with the corresponding number of points is applied to the signal segment, followed by a Fast Fourier Transform to suppress spectral leakage. The dominant frequency amplitude is extracted from the obtained amplitude spectrum, and the total harmonic distortion (THD) is calculated based on the harmonic amplitudes of preset orders using the following formula: ; In the formula, This is the fundamental frequency amplitude, i.e., the dominant frequency amplitude. For the first Second harmonic amplitude; The highest harmonic order is represented by the matrix. The dominant frequency amplitude and total harmonic distortion rate are sequentially concatenated to form a one-dimensional frequency domain feature matrix, which is then fused with the time domain features along the column direction to obtain a multi-dimensional time-frequency feature vector. By fusing time-domain statistical features and frequency-domain harmonic features, the signal characteristics of hydraulic transient disturbances are comprehensively characterized, providing a basis for accurate matching of the subsequent propagation influence domain.

[0033] In this embodiment, with Data was extracted by shifting the center forward and backward by 5 minutes to obtain a target signal segment of 600 seconds (containing 600 sampling points at a sampling frequency of 1Hz). A 600-point Hamming window was added for Fast Fourier Transform, and the amplitudes of the 1st to 5th harmonics were extracted to calculate the total harmonic distortion rate. The time-domain feature matrix was 1×3 in size, and the frequency-domain feature matrix was 1×2 in size. Finally, they were spliced ​​together to form a 1×5 multidimensional feature vector. For the target signal segment at node A, the root mean square value was calculated to be 0.32MPa, the skewness was -1.5, and the kurtosis was 5.2.

[0034] The cosine similarity between the aforementioned multidimensional time-domain and frequency-domain feature vectors and the feature vectors of each historical sample in the pipeline transient response model library is calculated one by one. The pipeline spatial connectivity region corresponding to the historical sample with a similarity greater than a preset threshold is selected as the propagation influence domain of this abnormal event.

[0035] Starting with the initial abnormal metering node as the source, a breadth-first traversal is performed along the water flow direction based on the pipeline spatial data table. Adjacent nodes at a preset depth are read and connected sequentially to generate multiple candidate tracing chains. The initial weight of each chain is determined based on a normalized scheme of hydraulic fluid conductivity: first, the pipe inner diameters of each water transmission pipeline segment associated with the chain are retrieved from the enterprise asset management database. Pipe length and Manning roughness coefficient As a priori attribute parameter, the single-segment conductance index of each pipe segment is then calculated using the following formula: ; The hydraulic fluid conductivity index of each pipe segment within a chain is arithmetically summed to obtain the magnitude of the hydraulic fluid conductivity index for that chain. Then, the proportion of each chain's index value in the sum of all chain indices is used as the initial weight for each chain. By delineating the propagation influence domain through time-frequency feature matching, the anomaly search range can be significantly narrowed. Weighting candidate chains based on prior hydraulic conductivity indices avoids the iterative divergence problem caused by equal-weight initialization from the outset.

[0036] In this embodiment, the cosine similarity preset threshold is set to 0.85; candidate chains are generated by reading first-level adjacent nodes and second-level adjacent nodes along the water flow direction; the inner diameter of a certain DN500 pipe segment is 0.5m, the length is 100m, and the Manning roughness coefficient is 0.013. Substituting these values ​​into the relational formula, the numerator is approximately 0.1963 and the denominator is 1.3, resulting in a single-segment conductivity index of approximately 0.151; five candidate tracking chains are generated under a certain typical abnormal event.

[0037] As an alternative implementation, the depth of the breadth-first traversal can be adaptively adjusted within the range of [1,4] according to the actual topological complexity of the pipeline network, so as to balance the search range and computational efficiency.

[0038] S2. Construct the model, set up a spatiotemporal search window, correct parameters, and update the confidence level. (I) Construction of transient propagation model For each candidate tracing chain, the interface of the open-source hydraulic engine is invoked to construct the corresponding transient hydraulic propagation model. First, the geographic coordinates and elevation data of each pressure monitoring node included in the chain are loaded sequentially, and a two-dimensional coordinate basis system covering the current pipeline subset is established through affine transformation. Then, based on the continuity and momentum equations describing one-dimensional unsteady flow in fluid mechanics, combined with the elastic modulus of the pipe material... Water volume elastic modulus and pipe wall thickness The wave velocity propagation coefficient of the pressure wave inside the pipe is calculated using the water hammer wave velocity formula: ; In the formula, For water density, Let be the pipe's inner diameter. Simultaneously, the steady-state friction coefficient within the pipe is discretized at a set spatial and time step, transforming it into the drag coefficient in the water hammer characteristic line method (MOC) algebraic equation system, and assembled into a damped hysteresis sparse diagonal matrix according to the pipe network node topology. Then, the wave velocity propagation coefficient and the damped hysteresis matrix are substituted into the equation system to integrate a partial differential state transition space equation matrix with pressure head and pipe velocity as unknown state variables. The grid spatial difference node spacing and time difference step are set, and through finite difference discretization, the original partial differential equation matrix is ​​transformed into a difference algebraic formula system containing Jacobian matrix solution factors. Finally, a one-dimensional transient propagation model is generated that can evolve forward with a preset time series accuracy and output the predicted pressure and flow state variables at the nodes to be verified. The combination of spatial and time steps must satisfy the Courant equation. Friedrich Levi conditions are used to ensure the stability and convergence of numerical computations.

[0039] In this embodiment, the wntr library is used to construct a transient hydraulic propagation model. The elastic modulus of the pipe material corresponding to a certain chain is 160 GPa, the bulk elastic modulus of the water is 2.19 GPa, and the pipe wall thickness is 12 mm. Substituting these values ​​into the water hammer wave velocity formula, the wave velocity propagation coefficient is calculated to be approximately 1050 m / s. The grid spatial difference node spacing is set to 20 m and the time difference step size is 0.01 s. The generated one-dimensional transient propagation model can output predicted state variables with a time series accuracy of 0.01 s.

[0040] As an optional implementation method, the discretization process of the water hammer characteristic line method can be encapsulated and implemented using a professional simulation library in the water industry (such as TSNet). During engineering implementation, only the pipeline INP file and parameters such as water hammer wave velocity need to be input, and the interface can be called to output the predicted state variables.

[0041] (ii) Setting the spatiotemporal search window The spatiotemporal search window setting essentially couples the hydraulic wave propagation delay with sensor measurement noise, encompassing both time and state dimensions. In the time dimension, the total pipe length is obtained by arithmetically summing the lengths of each water pipe segment included in the candidate chain, then divided by the theoretical propagation rate of water hammer elastic waves to obtain the total water hammer travel time, i.e., the cumulative hydraulic propagation time delay. This is based on the time of anomaly occurrence. Using the theoretical arrival time obtained by superimposing this time delay as the axis, preset time intervals are extended forward and backward to define the boundaries of the time interval. In terms of state dimension, the nameplate accuracy class and calibration range of the pressure gauge and flow meter of the next node to be verified are retrieved. The standard error of a single measurement of pressure and flow is calculated according to the standard formula, and a diagonal error variance standard matrix is ​​constructed to characterize the monitoring uncertainty. Based on the predicted state variables output by the propagation model at the axis time, a Gaussian distribution is followed. The principle extends both upwards and downwards to form the boundary range of the state space.

[0042] Ultimately, the two-dimensional region formed by the intersection of the time interval boundary and the state space limit constitutes the spatiotemporal search window at the node to be verified.

[0043] In this embodiment, the time interval boundaries are defined by extending 0.8s before and after the axis center moment in the time dimension; assuming the total pipe length is 1150m, the wave velocity is 1050m / s, and the cumulative hydraulic propagation time delay is approximately 1s, the corresponding time interval boundaries are... The pressure gauge has an accuracy of 0.5 grade, a calibration range of 1 MPa, and a single standard error of 0.005 MPa; the flow meter has an accuracy of 0.2 grade, a calibration range of 500 m³ / h, and a single standard error of 1 m³ / h; if the model predicts a pressure of 0.3 MPa, then the pressure tolerance bandwidth is [0.285, 0.315] MPa, which together with the above time interval constitutes a complete spatiotemporal search window.

[0044] (iii) Confidence update When the hydraulic time advances to the next node to be verified, the real-time monitoring data matrix of that node is extracted and a vector-level subtraction operation is performed with the predicted state variables output by the propagation model at the same timestamp. The difference sign is retained to generate a multi-dimensional first-order residual vector.

[0045] Subsequently, an unscented Kalman filter class is instantiated using the filter library to construct an unscented Kalman filter model. The aforementioned multidimensional first-order residual vector is used as the observation input, and the pipe friction coefficient and node demand parameters in the propagation model are corrected online based on the Kalman gain matrix. Simultaneously, the residual vector is substituted into a zero-mean multidimensional Gaussian distribution evaluation function to extract the corresponding monotonic likelihood correlation score. The closer the absolute value of the residual is to zero, the larger the likelihood score. This likelihood score is multiplied by the forward probability stored at the previous time step of the current chain to obtain the probability degree at the current time step. The sum of the probability degrees of all parallel candidate chains at the current time step is calculated, and the probability degree of the current chain is divided by this sum for normalization. The resulting real percentage is overwritten and stored in the original storage unit, which is the updated confidence degree of the chain. The initial forward probability value is taken from the initial chain weights calculated based on the hydraulic conductivity index in the previous steps, giving the iteration starting point physical meaning. This closed-loop dynamic inference mechanism allows the confidence of erroneous paths to naturally decay, while the confidence of true paths continuously increases.

[0046] In this embodiment, the `UnscentedKalmanFilter` class is instantiated using the `filterpy` library to construct an unscented Kalman filter model. The measured pressure at a certain node is 0.32 MPa and the flow rate is 150 m³ / h. The corresponding model prediction is 0.31 MPa and 152 m³ / h. The generated multidimensional first-order residual vector is... Substituting the values ​​into the likelihood evaluation function, the likelihood score is calculated to be 0.85. The forward probability of the chain at the previous time step is stored as 0.3, and multiplying them together gives the current probability score of 0.255. The probability scores of the three parallel chains at the current time are 0.255, 0.120, and 0.025, respectively, with a total of 0.4. After normalization, the latest confidence score of the current chain is 63.75%.

[0047] S3. Iteratively assess confidence level, and use single-chain early warning or multi-chain closed-loop tracing.

[0048] The system continuously reads monitoring data for the next time step using a loop, iteratively executing the model correction and confidence update process in S2. A preset convergence threshold is used as the confirmation criterion for a single anomaly propagation path: when the confidence of a candidate tracing chain exceeds the convergence threshold for the first time during the loop, the iteration terminates, and the chain is confirmed as the main anomaly propagation path. The anomaly propagation path information is recorded in the system log, and an early warning message containing the anomaly node ID, occurrence timestamp, and propagation path details is pushed.

[0049] In this embodiment, a while loop is used to continuously read monitoring data, with a preset convergence threshold of 0.9. When the confidence level of a candidate chain first exceeds 0.9, the iteration is terminated and it is confirmed as the main propagation path of the anomaly. It should be noted that, based on Bayesian inference theory, a posterior probability exceeding 0.9 corresponds to a high level of confidence, so this embodiment uses 0.9 as the convergence threshold. For critical water source networks with higher response sensitivity requirements, the threshold can be lowered to 0.85, and for scenarios with high false alarm costs, it can be raised to 0.95.

[0050] As an alternative implementation method, the warning information can be pushed through the Simple Mail Transfer Protocol (SMTP) service, or through other methods such as SMS, system pop-ups, etc.

[0051] If no single chain converges after the preset maximum number of iterations, but multiple candidate chains have confidence levels exceeding the background threshold for a preset number of consecutive time steps, the joint reverse tracing module is activated to locate the anomaly source area: First, all anomaly monitoring nodes covered by the high-susceptibility chains are read, and the real geographic longitude and latitude data of each node are extracted from the spatial GIS database to construct a discrete two-dimensional geospatial location point set; then, centripetal clustering is called to construct an envelope tangent polygonal network with all discrete location points as endpoints, and the plane centroid coordinates of the core cluster containing the most pipe point samples and the highest density are calculated; if the centroid is located in a non-pipeline area, the nearest physical pipeline segment to the centroid is further retrieved through geospatial index, and the centroid is mapped to the central axis of the pipeline segment using a projection algorithm, and the resulting landing point is the coordinate of the core source of the anomaly location.

[0052] Finally, in the GIS rendering layer, using the coordinates of the core source as the geometric center, a smooth and continuous boundary circle is drawn outwards with a preset radius to construct a completely closed boundary layer covering the area. This boundary layer is then visualized as the anomaly source area on the monitoring interface. This reverse tracing method converges the anomaly source area at a global level, solving the problem of difficulty in locating the source under multiple concurrent alarms.

[0053] To form a diagnostic closed loop, a forward transient simulation verification mechanism is introduced: In the hydraulic model, the nodes in the above-mentioned anomaly source area are configured as newly added leakage holes or nodes with sudden change demand. A forward transient simulation is run to generate simulated pressure and flow time series curves for all relevant nodes in the network. The mean square error (MSE) between the simulated curve and the actual monitoring curve is calculated. When the MSE is lower than the preset basic measurement noise tolerance, it is confirmed that the reverse tracing location result is consistent with the actual operating conditions. An anomaly diagnosis report is generated and archived, thus completing the full-chain closed loop of identification → tracking → location → verification.

[0054] In this embodiment, the maximum number of iterations is preset to 50. Backward tracing is initiated when the confidence levels of multiple candidate chains exceed the background threshold for three consecutive time steps. The DBSCAN algorithm is used for centripetal clustering, with a boundary circle extension radius of 45m. The extracted anomaly node latitude and longitude data includes multiple points such as (116.345672°E, 39.987654°N). After clustering, the planar centroid coordinates are approximately (116.346001°E, 39.988002°N). A DN300 municipal pipeline exists approximately 12m from this centroid. Projecting the centroid onto the pipeline's central axis yields the core source coordinates. A closed circular region with these coordinates as the center and a radius of 45m is the anomaly source area. Anomaly source area nodes are configured in the hydraulic model of the wntr library, and a forward transient simulation is run. The preset basic measurement noise tolerance is 0.05; accurate positioning is confirmed when the mean square error is below this value.

[0055] To address the aforementioned technical solution, an experimental scenario was set up, comprising 600 kilometers of municipal water supply pipeline and 500 pressure and flow monitoring nodes. The sampling frequency was uniformly set to 1Hz, and the test set covered 300 real pipeline rupture anomaly events and 200 normal hydraulic fluctuation disturbances.

[0056] The experiment set up four comparative modules: Module 1 only uses isolated forest for anomaly identification and nearest node localization; Module 2 adds sliding window quantile regression and time-frequency feature extraction on the basis of Module 1; Module 3 adds transient propagation model and spatiotemporal search window mechanism on the basis of Module 2; Module 4 is the complete technical solution of this disclosure, which includes reverse clustering and forward closed-loop verification.

[0057] After 60 hours of continuous offline testing, Module 1 achieved an anomaly recognition accuracy of only 72.3% with a false alarm rate as high as 15.4%, and an average spatial positioning error of 420m. Module 2 improved the accuracy to 86.5%, reduced the false alarm rate to 5.2%, and reduced the error to 280m. Module 3 achieved an accuracy of 94.8%, reduced the false alarm rate to 2.1%, and reduced the error to 85m. Module 4, employing the complete technical components of this invention, performed best in the experiment, achieving an anomaly recognition accuracy of 98.7%, reducing the false alarm rate to only 0.6%, and precisely controlling the average positioning error to within 17m, fully demonstrating the progressive improvement effect of the combined and superimposed technical components of this invention.

[0058] like Figure 2 As shown in the figure, this diagram illustrates the amplitude distribution characteristics of the abnormal pressure signal after being converted to the frequency domain by a Fast Fourier Transform. The figure clearly presents the amplitudes corresponding to the dominant frequency and multiple harmonics, from which the dominant frequency amplitude and total harmonic distortion rate can be extracted to form a frequency domain feature matrix characterizing the hydraulic transient disturbance, providing a basis for matching the subsequent propagation influence domain.

[0059] like Figure 3 As shown, the initial weight allocation results of multiple candidate tracing chains generated under a typical anomaly event are presented. The initial weights of each chain are significantly different, and the weight values ​​are positively correlated with the hydraulic fluid conductivity of the corresponding chain. This weight allocation method based on hydraulic characteristics can prioritize guiding the tracing direction of the anomaly propagation path, reduce invalid searches, and lower the risk of oscillations in subsequent iterations.

[0060] like Figure 4 As shown in the figure, the anomaly identification accuracy and average spatial positioning error of the four comparison modules are displayed respectively. With the gradual addition of technical components, the anomaly identification accuracy continues to improve, and the average spatial positioning error is significantly reduced. This intuitively demonstrates the progressive improvement effect of each technical aspect of the present invention on system performance, and verifies the superiority of the complete technical solution in anomaly early warning and location in large-scale pipeline networks.

[0061] The above description represents the preferred embodiments of the present invention. It should be noted that those skilled in the art can make various improvements and modifications without departing from the principles of the present invention, and these improvements and modifications should also be considered within the scope of protection of the present invention.

Claims

1. A closed-loop tracking and early warning method for metering anomalies in large-scale pipeline networks, characterized in that: include: S1. Based on the time-series monitoring data of pressure and flow at each metering node in the pipeline network, identify the initial abnormal metering nodes in the pipeline network and extract their multi-dimensional time-domain and frequency-domain feature vectors. Using the multidimensional time-domain and frequency-domain feature vectors, a propagation influence domain is matched from the pipeline network transient response model library, and candidate tracking chains with initial weights are generated by searching along the pipeline network topology within the propagation influence domain; S2, a propagation model is constructed for each candidate tracking chain. Based on the cumulative hydraulic propagation time delay and the monitoring uncertainty of the node to be verified, a spatiotemporal search window for the next node to be verified is set. The construction of the spatiotemporal search window includes: establishing the time interval boundary by extending a preset time forward and backward with the actual running time corresponding to the cumulative hydraulic propagation time delay as the axis time, and expanding the state space boundary range based on three times the standard deviation of the monitoring uncertainty of the node to be verified. The spatiotemporal search window is established by combining the time interval boundary and the boundary range. Within an empty search window, real-time monitoring data and predicted state values ​​are compared. The predicted state value at this moment is subtracted from the real-time monitoring data of the next node to be verified. The difference sign is retained to generate a multidimensional first-order residual vector. The propagation model parameters are corrected using the multidimensional first-order residual vector, and the posterior probability of each candidate tracking chain is updated as the confidence level. S3, the steps of setting the spatiotemporal search window and updating the confidence level are executed iteratively. When the confidence level of a single candidate tracking chain exceeds the convergence threshold, it is confirmed as an abnormal main propagation path and an early warning is issued. If no single chain converges but the confidence levels of multiple candidate tracking chains are higher than the background threshold, the abnormal source area is located by joint reverse tracing, and the consistency between the abnormal source area and the monitoring data of related nodes is verified by forward transient simulation to complete the closed loop.

2. The method for closed-loop tracking and early warning of metering anomalies in large-scale pipeline networks according to claim 1, characterized in that, The identification of initial abnormal metering nodes in the pipeline network includes: inputting the acquired pressure and flow time-series monitoring data into an isolated forest algorithm, calculating the abnormal score of each data point, and extracting data points with abnormal scores greater than a preset abnormal threshold as first candidate abnormal points; setting up a sliding window with a preset time step to intercept the pressure and flow time-series monitoring data, performing quantile regression calculation on the intercepted data to obtain the prediction interval boundary, and extracting data points outside the prediction interval boundary as second candidate abnormal points; calculating the intersection of the first candidate abnormal point and the second candidate abnormal point, and confirming the node and time position corresponding to the data point belonging to the intersection as the initial abnormal metering node and the time of abnormal occurrence.

3. The method for closed-loop tracking and early warning of metering anomalies in large-scale pipeline networks according to claim 2, characterized in that, The extraction of its multidimensional time-domain and frequency-domain feature vectors includes: taking the time of the anomaly occurrence as the center, extracting continuous monitoring data as target signal segments by offsetting by a preset time before and after the time of the anomaly occurrence; calculating the root mean square value, skewness, and kurtosis of the target signal segments and concatenating them in order to form a one-dimensional time-domain feature matrix; performing a fast Fourier transform on the target signal segments to convert them to the frequency domain, calculating the frequency domain amplitude spectrum to extract the main frequency amplitude and the total harmonic distortion rate of the harmonics of a preset order, forming a one-dimensional frequency-domain feature matrix; and concatenating and fusing the one-dimensional time-domain feature matrix with the one-dimensional frequency-domain feature matrix to generate the multidimensional time-domain and frequency-domain feature vectors.

4. The method for closed-loop tracking and early warning of metering anomalies in large-scale pipeline networks according to claim 1, characterized in that, Within the propagation influence domain, candidate tracing chains with initial weights are generated through a topological search along the pipeline network. This includes: setting the initial anomaly metering node as the source node of the traversal tree, and sequentially connecting adjacent nodes along the water flow to generate multiple candidate tracing chains; obtaining the pipe inner diameter, pipe length, and Manning roughness coefficient values ​​of the water transmission pipe segments associated with each candidate tracing chain; calculating the single-segment conductivity index of each water transmission pipe segment based on the pipe inner diameter, pipe length, and Manning roughness coefficient values, and summing the single-segment conductivity indices of all pipe segments within a single candidate tracing chain to obtain the corresponding conductivity index magnitude value; calculating the proportion of the conductivity index magnitude value of each candidate tracing chain in the sum of all chains, and assigning it as the initial weight.

5. The method for closed-loop tracking and early warning of metering anomalies in large-scale pipeline networks according to claim 1, characterized in that, The construction of a propagation model for each candidate tracking chain includes: sequentially loading the pressure monitoring node data contained in all candidate tracking chains to establish a two-dimensional coordinate basis system for the pipeline network; using the continuity equation and momentum equation of one-dimensional unsteady flow, calculating the wave velocity propagation coefficient and damping hysteresis matrix within the corresponding water transmission pipeline network to construct a state transition space equation matrix; based on the preset grid differential node spacing and time step differential step size constraints of the state transition space equation matrix, discretizing the state transition space equation matrix into a Jacobian difference algebra formula system to generate a one-dimensional transient propagation model for outputting predicted state variables.

6. The method for closed-loop tracking and early warning of metering anomalies in large-scale pipeline networks according to claim 1, characterized in that, The method based on cumulative hydraulic propagation time delay and monitoring uncertainty of the node to be verified includes: accumulating the pipe length values ​​contained in the candidate tracking chain to obtain the total pipe length, dividing the total pipe length by the theoretical propagation rate of water hammer elastic waves to establish the cumulative hydraulic propagation time delay; and constructing a diagonal error variance standard matrix to extract the monitoring uncertainty of the node to be verified based on the accuracy level and calibration range value of the metering equipment of the next node to be verified.

7. The method for closed-loop tracking and early warning of metering anomalies in large-scale pipeline networks according to claim 1, characterized in that, The step of using a multidimensional first-order residual vector to correct the propagation model parameters and update the posterior probability of each candidate tracking chain as the confidence level includes: substituting the multidimensional first-order residual vector into an evaluation function that conforms to a zero-mean multidimensional Gaussian distribution to extract the corresponding monotonic likelihood correlation score; multiplying the monotonic likelihood correlation score by the forward probability stored value of the current candidate tracking chain at the previous time step to obtain the probability level at the current time step; and normalizing the calculation by dividing the probability level of the current candidate tracking chain by the sum of the probability levels of all parallel chains, and updating the resulting value as the confidence level of the candidate tracking chain.

8. The method for closed-loop tracking and early warning of metering anomalies in large-scale pipeline networks according to claim 1, characterized in that, The joint reverse tracing and location of the anomaly source area includes: reading the real geographic longitude and physical latitude values ​​of all candidate tracking chain nodes whose confidence continuously exceeds the background threshold, and storing them as a discrete two-dimensional geospatial location point set; using centripetal clustering to calculate the planar centroid coordinates of the tangent polygon network of the discrete two-dimensional geospatial location point set, and projecting the planar centroid coordinates onto the nearest pipeline segment to establish the core source coordinates of the anomaly location; using the core source coordinates of the anomaly location as the geometric center, extending a preset absolute length reference distance in all directions to construct a boundary layer block of a completely closed boundary circle covering the area, which serves as the anomaly source area.

9. The method for closed-loop tracking and early warning of metering anomalies in large-scale pipeline networks according to claim 8, characterized in that, The step of using forward transient simulation to verify the consistency between the anomaly source region and the monitoring data of related nodes to complete the closed loop includes: configuring the nodes in the anomaly source region as newly added leak holes or nodes with sudden change requirements in the propagation model; running forward transient simulation to generate simulated pressure and flow time series curves for all related nodes in the network; calculating the mean square error between the simulated pressure and flow time series curves and the actual real-time monitoring data of related nodes; and when the mean square error is lower than the preset basic measurement noise tolerance, confirming that the reverse tracing location is accurate, generating an anomaly diagnosis report to complete the closed loop.

10. The method for closed-loop tracking and early warning of metering anomalies in large-scale pipeline networks according to claim 4, characterized in that, The single-segment conductivity index of each water transmission pipe section is calculated based on the pipe inner diameter, pipe length, and Manning roughness coefficient. This includes: multiplying the fourth power of the pipe inner diameter of the corresponding water transmission pipe section by the value of pi as the dividend, multiplying the pipe length of the corresponding water transmission pipe section by the Manning roughness coefficient as the divisor, and dividing the dividend by the divisor to obtain the single-segment conductivity index.

Images