Methods for identifying leakage in water supply networks used in urban lifeline engineering

Abstract

This invention belongs to the field of pipeline leakage detection technology, and relates to a method for identifying leaks in water supply networks used in urban lifeline projects. The method involves collecting pressure fluctuation data and synchronous flow data; calculating the difference between the pressure fluctuation data and a benchmark fluctuation template point by point to obtain a fluctuation residual sequence; extracting features from the fluctuation residual sequence to obtain spatial distribution characteristics; combining the spatial distribution characteristics and synchronous flow data to calculate the negative correlation between pressure drop and flow increase, as well as the time difference between pressure drop and flow increase; and marking suspected leakage features based on the calculation results. Finally, the method performs correlation analysis on the suspected leakage features based on pipe aging parameters and historical maintenance records to identify leaks. This invention uses a benchmark fluctuation template to strip away normal fluctuations to highlight weak abnormal signals, and integrates pipe aging parameters and historical maintenance records for multi-factor correlation confirmation. It progressively separates minute leakage features from complex background noise, effectively improving the accuracy of leak identification and reducing the false alarm rate.

Description

Technical Field

[0001] This invention relates to the field of pipeline leak detection technology, and in particular to a method for identifying leaks in water supply networks used in urban lifeline projects. Background Technology

[0002] Water supply networks are a core component of urban lifeline engineering, and reducing leakage in water distribution networks is crucial for improving water safety. Network pressure fluctuations are influenced by a combination of factors, including environmental disturbances, material aging, and structural fatigue, resulting in a blurred boundary between normal operating fluctuations and early signs of minor leaks in terms of data characteristics, making accurate identification difficult.

[0003] Acoustic monitoring technology is used to detect leaks in water supply networks by continuously monitoring pipeline noise using noise monitoring devices. However, leak signals are difficult to separate from complex background noise, and the extraction of minute leak characteristics is challenging. Other detection technologies locate leaks by calculating the time difference between abnormal signals reaching different probes, but under the combined influence of multiple factors, minute leak characteristics and normal fluctuation characteristics highly overlap in the time-frequency domain, resulting in blurred characteristic boundaries and a high risk of missed or false alarms. Summary of the Invention

[0004] Therefore, the purpose of this invention is to overcome the problem in the prior art that it is difficult to accurately separate and identify weak abnormal features caused by minor leaks from water supply network monitoring data under strong environmental noise and operational fluctuations, resulting in low leak identification accuracy. This invention provides a leak identification method for water supply networks used in urban lifeline projects. By using a benchmark fluctuation template to strip away normal fluctuations to highlight weak abnormal signals, and by integrating pipe aging parameters and historical maintenance records for multi-factor correlation confirmation, the invention progressively separates minor leak features from complex background noise, effectively improving the accuracy of leak identification and reducing the false alarm rate.

[0005] To address the aforementioned technical problems, this invention provides a method for identifying leakage in water supply networks used in urban lifeline engineering projects, comprising:

[0006] Collect pressure fluctuation data and synchronous flow data from sensor nodes in the water supply network;

[0007] The difference between the pressure fluctuation data and the benchmark fluctuation template constructed based on historical data is calculated point by point to obtain the fluctuation residual sequence;

[0008] Feature extraction is performed on the fluctuation residual sequence to obtain its spatial distribution characteristics;

[0009] Combining the spatial distribution characteristics and the synchronous flow data, the negative correlation between pressure drop and flow increase, as well as the time difference between pressure drop and flow increase, are calculated; based on the calculation results, the spatial distribution characteristics that meet the conditions are marked as suspected leakage characteristics.

[0010] The suspected leakage features are correlated based on pipe aging parameters and historical maintenance records, and leakage is identified based on the analysis results. The correlation analysis includes: determining the deviation between the energy attenuation rate of the suspected leakage feature propagating along the pipeline and the theoretical attenuation rate corresponding to the pipe aging parameters, as well as the spatial location of the suspected leakage feature and the distance between the historical maintenance points.

[0011] Preferably, constructing a benchmark fluctuation template based on historical data includes: extracting historical pressure fluctuation data from a historical database that have the same seasonal type, same weekday attribute, and same time period as the current period to form a historical sample set; performing a fast Fourier transform on each historical pressure fluctuation data in the historical sample set to obtain the frequency domain characteristics of each historical data, removing historical data with abnormal harmonic components in the frequency domain characteristics to obtain a clean historical sample set; clustering the clean historical sample set to divide it into multiple typical fluctuation pattern categories, selecting the category with the largest number of samples as the benchmark fluctuation pattern; performing point-by-point statistical averaging on the historical pressure fluctuation data in the benchmark fluctuation pattern to obtain the mean curve at each time point; and smoothing the mean curve to obtain the benchmark fluctuation template.

[0012] Preferably, the method of calculating the difference between the pressure fluctuation data and the benchmark fluctuation template constructed based on historical data point by point to obtain the fluctuation residual sequence includes: calculating the difference between the pressure fluctuation data and the benchmark fluctuation template point by point at the same time point to obtain an initial residual sequence; performing sliding window standard deviation analysis on the initial residual sequence to calculate the local standard deviation at each time point, and marking the time points where the local standard deviation exceeds a threshold as suspected abrupt change points; dividing the initial residual sequence into several continuous segments using the suspected abrupt change points as boundary points, performing polynomial fitting on each continuous segment to obtain the fitting baseline corresponding to each segment; subtracting the corresponding fitting baseline from the initial residual of each continuous segment to obtain the corrected residual sequence; and performing zero-mean processing on the corrected residual sequence to obtain the fluctuation residual sequence.

[0013] Preferably, the spatial distribution characteristics include the duration of the residual, the cumulative energy value, and the spatial propagation characteristics.

[0014] Preferably, feature extraction is performed on the fluctuation residual sequence to obtain spatial distribution characteristics, including: segmenting the fluctuation residual sequence into a sliding time window to obtain multiple fluctuation residual subsequences; for each fluctuation residual subsequence within a time window, constructing a sliding energy window centered on each time point, calculating the square integral of the residual amplitude within each sliding energy window as the local energy value at that time point, and identifying the set of time points where the local energy value continuously exceeds a threshold as the effective fluctuation interval; determining the duration of the residual based on the start and end times of the effective fluctuation interval; calculating the square integral of the amplitude of the fluctuation residual subsequence within the effective fluctuation interval to obtain the cumulative energy value of the residual; for residual signals collected at different nodes for the same suspected leakage event, calculating the arrival time difference of the residual signals at each node, and combining the pipeline distance between nodes to calculate the propagation speed of the residual between adjacent nodes, thereby obtaining the spatial propagation characteristics of the residual.

[0015] Preferably, calculating the negative correlation between pressure drop and flow rate increase includes: for each time window, forming data pairs of instantaneous pressure values ​​and instantaneous flow rates at each time point within the overlapping period of the pressure drop interval and the flow rate increase interval in that time window, calculating the Pearson correlation coefficient of all data pairs, and obtaining the negative correlation coefficient between pressure drop and flow rate increase.

[0016] Preferably, calculating the time difference between pressure drop and flow increase includes: comparing the start time of the pressure drop interval with the start time of the flow increase interval, and calculating the time difference between the two. When the start time of the pressure drop interval is earlier than the start time of the flow increase interval, the time difference is a positive value, and otherwise it is a negative value.

[0017] Preferably, the conditions for marking the suspected leakage characteristics include: the absolute value of the negative correlation coefficient is greater than a threshold, and the time difference is within the time delay range.

[0018] Preferably, determining the energy attenuation rate of the suspected leakage feature propagating along the pipeline includes: extracting the pressure fluctuation residual amplitude between two adjacent sensor nodes and the pipeline distance between the nodes; determining the energy attenuation rate as follows: η represents the energy decay rate; A1 represents the pressure fluctuation residual amplitude of the upstream sensor node; A2 represents the pressure fluctuation residual amplitude of the downstream sensor node; L represents the pipe distance between two adjacent sensor nodes.

[0019] Preferably, the suspected leakage feature is correlated with pipe aging parameters and historical maintenance records, and leakage is identified based on the analysis results. This includes: querying the pipe aging parameters of the pipe segments where the two adjacent sensor nodes are located, determining the theoretical attenuation rate of the pipe segment based on the pipe aging parameters; calculating the deviation between the energy attenuation rate and the theoretical attenuation rate; in response to the deviation being less than a threshold, determining the first spatial location corresponding to the suspected leakage feature based on the time information of the two adjacent sensor nodes and the pipe distance between the nodes; querying historical maintenance records whose spatial location corresponding to the suspected leakage feature is within the search radius to obtain the second spatial location of the historical maintenance point; if there exists any historical maintenance point such that the distance between the first spatial location and the second spatial location is less than a threshold, then leakage is identified.

[0020] The technical solution of the present invention has the following advantages compared with the prior art:

[0021] The water supply network leakage identification method for urban lifeline projects described in this invention uses a benchmark fluctuation template to strip away normal fluctuations to highlight weak abnormal signals, integrates pipe aging parameters and historical maintenance records for multi-factor correlation confirmation, and progressively separates minute leakage features from complex background noise, effectively improving the accuracy of leakage identification and reducing the false alarm rate.

[0022] Among them, a benchmark fluctuation template is constructed based on historical data. By comparing the pressure fluctuation data and the benchmark fluctuation template, the normal fluctuation components that conform to the periodic pattern are stripped away, thereby obtaining the pressure fluctuation residual sequence containing potential abnormal information. This effectively solves the problem of normal fluctuations drowning out weak abnormal signals, and makes the characteristics of small leaks stand out from strong background noise.

[0023] By integrating pipe aging parameters and historical maintenance records, correlation analysis is performed on suspected leakage characteristics, fully considering the impact of the pipe's own condition on leakage characteristics, thereby improving the accuracy and reliability of identification. Attached Figure Description

[0024] To make the content of this invention easier to understand, the invention will be further described in detail below with reference to specific embodiments and accompanying drawings, wherein:

[0025] Figure 1 This is a flowchart of a water supply network leakage identification method for urban lifeline engineering in a preferred embodiment of the present invention;

[0026] Figure 2 This is a flowchart of constructing a benchmark fluctuation template in a preferred embodiment of the present invention;

[0027] Figure 3 This is a flowchart illustrating the calculation of the fluctuation residual sequence in a preferred embodiment of the present invention;

[0028] Figure 4 This is a flowchart of extracting spatial distribution features in a preferred embodiment of the present invention;

[0029] Figure 5 This is a flowchart for marking suspected leaks in a preferred embodiment of the present invention;

[0030] Figure 6 This is a flowchart of the correlation analysis for identifying leakage in a preferred embodiment of the present invention. Detailed Implementation

[0031] The present invention will be further described below with reference to the accompanying drawings and specific embodiments, so that those skilled in the art can better understand and implement the present invention. However, the embodiments described are not intended to limit the present invention.

[0032] The purpose of this invention is to provide a method for identifying leakage in water supply networks for urban lifeline projects, overcoming the problem that it is difficult to accurately separate and identify subtle abnormal features caused by minor leakage from water supply network monitoring data under strong environmental noise and operational fluctuations, resulting in low leakage identification accuracy.

[0033] The following combination Figure 1 The technical solutions of the present invention will be described in detail with specific embodiments. It should be noted that the embodiments are only used to explain the present invention and not to limit the present invention.

[0034] Step 1, Data Collection:

[0035] Pressure and flow sensors are deployed at key nodes in a city's water supply network. These key nodes include, but are not limited to, pump station outlets, network junctions, pipe segment ends, and locations where problems have historically occurred. Pressure sensors collect pressure fluctuation data at a set sampling frequency, while flow sensors collect synchronous flow data at a set sampling frequency. Synchronous flow data refers to flow values ​​collected synchronously by flow sensors on the same pipe segment or adjacent nodes within the same time window as the pressure fluctuation data; it characterizes changes in water flow within the pipeline when pressure fluctuations occur. All collected data is timestamped and uploaded to a central data platform in real time via a communication module to construct an operational status dataset.

[0036] Step 2, obtain the fluctuation residual sequence:

[0037] First, a baseline fluctuation template is constructed based on historical data. This template is a reference curve obtained from historical data of the same period, reflecting the typical fluctuation patterns of the pipeline network under conditions of no leakage anomalies. The construction process includes: extracting historical pressure fluctuation data from the historical database that shares the same seasonal type, weekday attribute, and time period as the current period; performing statistical analysis on the extracted historical data, removing data containing abnormal events; performing cluster analysis on the filtered data, selecting the category with the largest number of samples as the baseline fluctuation pattern; and performing point-by-point statistical averaging and smoothing on the data in this pattern to obtain the baseline fluctuation template.

[0038] Pressure fluctuation data is compared with a baseline fluctuation template at the same time point to calculate the point-by-point difference, thus obtaining the pressure fluctuation residual. The pressure fluctuation residual reflects the degree of deviation of real-time fluctuations from the normal waveform. The residuals at all time points are serialized to obtain a fluctuation residual sequence. This sequence contains weak anomalous signals that may be caused by minor leaks.

[0039] Step 3: Extract multidimensional features describing anomalous events from the fluctuation residual sequence, including the duration of the residual, the cumulative energy value, and the spatial propagation characteristics.

[0040] Specifically, the fluctuation residual sequence is segmented by adaptive time windows, and the following features are extracted from the fluctuation residuals within each time window: the duration of the residual, which refers to the length of time from the beginning to the end of the effective fluctuation interval; the cumulative energy value of the residual, which refers to the integral of the square of the residual amplitude within the effective fluctuation interval, reflecting the total energy of the abnormal event; and the spatial propagation characteristics of the residual, which refers to the propagation speed of the residual signal between adjacent sensor nodes, obtained by calculating the arrival time difference of the residual signals at each node and combining it with the pipeline distance between the nodes.

[0041] Step four involves analyzing the correlation between pressure fluctuations and flow rate changes to identify suspected characteristics that conform to the physical laws of leakage. Specifically, by combining spatial distribution characteristics and synchronous flow data, the negative correlation between pressure drop and flow rate increase, as well as the time difference between pressure drop and flow rate increase, are calculated.

[0042] Pressure drop refers to the difference between the baseline pressure value and the lowest pressure value within the pressure drop range, while flow rate increase refers to the difference between the highest flow rate value and the baseline flow rate value within the flow rate increase range. Calculating the negative correlation between the two is used to determine whether the pressure drop and flow rate increase occur simultaneously and in opposite directions. The time difference between pressure drop and flow rate increase refers to the difference between the start time of pressure drop and the start time of flow rate increase, used to determine the order of pressure and flow rate changes.

[0043] Based on the calculation results, spatial distribution characteristics that meet the criteria are identified as suspected leakage characteristics. The criteria are: the absolute value of the negative correlation coefficient is greater than a preset threshold, and the time difference is within a preset time delay range. These two conditions together ensure that the labeled characteristics conform to both the numerical relationship of pressure decrease and flow increase, and the temporal pattern of pressure decrease preceding flow, effectively distinguishing leakage from normal water use events.

[0044] Step 5: Based on the pipe aging parameters and historical maintenance records, perform correlation analysis on the marked suspected leakage characteristics, and identify the leakage based on the analysis results.

[0045] Specifically, association analysis includes the following two aspects of judgment:

[0046] First, determine the deviation between the energy attenuation rate of suspected leakage characteristics propagating along the pipeline and the theoretical attenuation rate corresponding to the pipe aging parameters. The energy attenuation rate refers to the degree of energy attenuation per unit distance when a pressure wave propagates along the pipeline, and is calculated by comparing the pressure amplitudes of adjacent sensor nodes; the theoretical attenuation rate is the expected attenuation value calculated by a physical model based on the pipe aging parameters; when the deviation between the measured attenuation rate and the theoretical attenuation rate is less than a threshold, it indicates that the energy attenuation characteristics of this feature are consistent with the current aging state of the pipeline.

[0047] Second, determine the spatial location of the suspected leak feature and its distance from historical maintenance points. The spatial location is calculated using a multi-sensor time-of-arrival (TOA) localization algorithm; historical maintenance points refer to locations recorded in the pipeline status database where maintenance work has been performed in the past. When the distance between the feature location and any historical maintenance point is less than a threshold, it indicates that the feature occurs in a weak area where problems have occurred in the past.

[0048] When both of the above conditions are met simultaneously, i.e. the energy attenuation rate deviation is less than the threshold and the distance from the historical maintenance point is less than the threshold, the suspected leakage feature is confirmed to correspond to a real micro-leak, and the leakage identification is completed.

[0049] The water supply network leakage identification method for urban lifeline projects described in this invention uses a benchmark fluctuation template to strip away normal fluctuations to highlight weak abnormal signals, integrates pipe aging parameters and historical maintenance records for multi-factor correlation confirmation, and progressively separates minute leakage features from complex background noise, effectively improving the accuracy of leakage identification and reducing the false alarm rate.

[0050] In the above embodiment, a baseline fluctuation template is constructed based on historical data. When constructing the baseline fluctuation template, the historical data will be mixed with data from abnormal periods such as leakage events and pipe burst accidents. There may be multiple fluctuation patterns on different dates in the same period. In the process of obtaining the fluctuation residual by comparing point by point, there may be a phase deviation between the template and the real-time data, which will lead to the generation of pseudo residuals by direct subtraction.

[0051] To address the problems of historical data contamination, diverse fluctuation patterns, inaccurate benchmark fluctuation templates due to phase deviation, distorted fluctuation residual extraction, and consequently, reduced accuracy of leakage identification, this invention provides a preferred method for identifying leakage in water supply networks. The following is a summary... Figures 2-3 Detailed explanation of the specific steps.

[0052] The first step is to construct a historical sample set:

[0053] Historical pressure fluctuation data with the same seasonal type, weekday attribute, and time period as the current period are extracted from the historical database to form a historical sample set.

[0054] Same season type refers to the same season among spring, summer, autumn, and winter; same weekday attribute refers to the same type among weekdays, weekends, or holidays; same time period refers to the same hour and minute intervals within a day, such as 9:00-10:00 AM every day.

[0055] Taking the current time period as 9:00 AM on a summer Wednesday as an example, historical pressure fluctuation data for all summer Wednesdays from 9:00 AM to 10:00 AM over the past three years are extracted; this three-dimensional screening method ensures that the extracted historical data has similar water use patterns and environmental conditions to the current time period.

[0056] The second step is to eliminate frequency domain anomalies:

[0057] A fast Fourier transform is performed on each historical pressure fluctuation data in the historical sample set to obtain the frequency domain characteristics of each historical data; the frequency domain characteristics include the amplitude and phase information of each frequency component.

[0058] In the frequency domain, normal water usage fluctuations typically manifest as the fundamental frequency and its harmonics, such as the fundamental frequency and its harmonics with a 24-hour period. However, abnormal events such as leaks and pipe bursts can generate additional harmonic components at non-fundamental frequency harmonic positions. The frequencies and amplitudes of these components are significantly different from normal fluctuations.

[0059] The frequency domain characteristics of each historical data point are compared with the normal fluctuation frequency domain template to identify historical data containing abnormal harmonic components. Specifically, if a historical data point exhibits an energy peak at a non-fundamental frequency harmonic position, and this peak exceeds three times the average energy of the background noise in the vicinity of that frequency, the data is considered to contain an abnormal event and is removed. This process yields a clean historical sample set, ensuring that the subsequently constructed template is not contaminated by historical anomalies.

[0060] The third step is to select the baseline volatility pattern using cluster analysis:

[0061] Cluster analysis was performed on the clean historical sample set to classify several typical fluctuation patterns. The K-means algorithm was used for cluster analysis, and the clustering features included the mean, variance, peak time, trough time, rate of rise, and rate of fall of the waveform. The number of clusters K was empirically set to 3, corresponding to weekday peak patterns, weekday off-peak patterns, and holiday patterns, respectively.

[0062] For each category, the average intra-cluster distance and inter-cluster distance are calculated to verify the rationality of the clustering effect. After clustering, the number of samples in each category is counted, and the category with the largest number of samples is selected as the baseline fluctuation pattern. This pattern represents the most common and typical fluctuation pattern in the current period and can reflect the pressure change characteristics under normal conditions to the greatest extent.

[0063] Step 4: Mean curve calculation and smoothing:

[0064] First, a point-by-point statistical average is performed on all historical pressure fluctuation data in the benchmark fluctuation model to obtain the mean curve for each time point. Point-by-point statistical averaging refers to calculating the arithmetic mean of all sample values ​​at the same time point.

[0065] Then, the mean curve is smoothed by filtering to eliminate residual random noise and obtain the final reference fluctuation template. The smoothing filter uses a Savitzky-Golay filter with a window length of 11 sampling points and a polynomial order of 3. This filter can smooth noise while preserving the overall shape and peak characteristics of the waveform well.

[0066] After obtaining the baseline fluctuation template, the operating status dataset is compared point by point with the baseline fluctuation template, and the difference between the two is calculated as the fluctuation residual. Figure 3 As shown, the specific implementation steps are as follows:

[0067] The first step is to calculate the initial residual sequence:

[0068] The pressure fluctuation data in the operational status dataset is compared with the obtained benchmark fluctuation template at the same time point to calculate the point-by-point difference, thus obtaining the initial residual sequence. The calculation formula is: r0(t) = P(t) - T(t), where P(t) is the real-time pressure value, T(t) is the template value, and r0(t) is the initial residual. The initial residual sequence contains multiple components: real anomalies that may be caused by leakage, phase deviation between the template and real-time data, pipeline network trend drift, and random noise.

[0069] The second step is to identify suspected mutation points:

[0070] A sliding window standard deviation analysis is performed on the initial residual sequence to calculate the local standard deviation at each time point. A sliding window is a window of a certain length centered at the current time point, extending forward and backward; simultaneously, the global standard deviation of the entire initial residual sequence is calculated. Time points where the local standard deviation exceeds a certain multiple of the global standard deviation are marked as suspected abrupt change points; this multiple can be set to 1.5 times.

[0071] Suspected abrupt change points correspond to locations where the intensity of fluctuations changes significantly, such as the moment when water usage transitions from a peak to a trough, or the moment when a stable period enters an abnormal period; these abrupt change points are the dividing points between different fluctuation patterns.

[0072] The third step is piecewise polynomial fitting:

[0073] Using suspected mutation points as demarcation points, the initial residual sequence is divided into several continuous segments; the residuals within each segment have similar statistical properties.

[0074] For each continuous segment, a polynomial fit is performed to obtain the corresponding fitting baseline. A third-order polynomial can be used for polynomial fitting to minimize the sum of squared errors between the fitted curve and the residuals within the segment. The fitting baseline reflects the trend drift within that segment, such as slow pressure changes caused by temperature variations, or baseline shifts caused by changes in total water usage; such trend drifts are not anomalies caused by leakage and should be removed from the residuals.

[0075] Step 4: Trend removal and zero-mean normalization:

[0076] The corrected residual sequence is obtained by subtracting the corresponding fitted baseline from the initial residual of each continuous segment. The corrected residual sequence eliminates the trend drift within each segment, so that the residual mainly contains possible anomalous signals and random noise.

[0077] Calculate the overall mean of the corrected residual sequence. When the absolute value of the overall mean is greater than a preset threshold, subtract the overall mean from the corrected residual sequence to obtain the final fluctuating residual sequence. When the absolute value of the overall mean is not greater than the preset threshold, directly use the corrected residual sequence as the final fluctuating residual sequence.

[0078] The preset threshold is determined based on the pressure sensor's range and is set to 0.01 times the pressure sensor's range. For example, if the pressure sensor's range is 1 MPa, then the threshold is 0.01 MPa. This zero-meaning process ensures that the entire fluctuation residual sequence fluctuates around zero, providing a unified benchmark for subsequent feature extraction.

[0079] The embodiment of the present invention eliminates phase deviation between the template and real-time data, removes interference from trend drift, and suppresses the amplification of local noise during the point-by-point comparison process, thereby obtaining a pure fluctuation residual sequence that can truly reflect potential leakage information.

[0080] In the above embodiments, feature extraction is performed on the fluctuation residual sequence to obtain spatial distribution characteristics. However, in the specific analysis process, the non-fixed fluctuation period makes feature extraction difficult, the start and end times of fluctuation are vaguely defined, and energy calculation is affected by noise. To solve this problem, this embodiment of the invention provides a preferred method for identifying leakage in water supply networks, which is summarized below. Figure 4 Detailed explanation of the specific steps.

[0081] Step 1: Sliding time window segmentation:

[0082] The fluctuation residual sequence obtained in the preceding steps is segmented by a sliding time window to obtain multiple fluctuation residual subsequences. The sliding time window segmentation adopts an adaptive window length strategy, dynamically adjusting the time window length according to the fluctuation frequency. Specifically, a Fast Fourier Transform is performed on the fluctuation residual subsequence within the current time window to obtain its spectral distribution. The frequency component with the largest amplitude is identified as the dominant frequency of the time window, and the time window length is dynamically adjusted according to the magnitude of the dominant frequency: when the dominant frequency is higher than the upper frequency threshold, it indicates that the fluctuation frequency is high and changes rapidly, and the time window needs to be shortened to avoid mixing in multiple fluctuations. In this case, the time window length is set to a shorter value; when the dominant frequency is lower than the lower frequency threshold, it indicates that the fluctuation frequency is low and changes slowly, and the time window needs to be extended to ensure that the complete fluctuation cycle is included. In this case, the time window length is set to a longer value; when the dominant frequency is between the upper and lower frequency thresholds, the standard time window length is used.

[0083] In one embodiment, the upper frequency threshold is set to 0.1Hz, corresponding to a fluctuation period of 10 seconds; the lower frequency threshold is set to 0.01Hz, corresponding to a fluctuation period of 100 seconds. Based on this, the time window length is dynamically adjusted as follows: when the main frequency > 0.1Hz, the time window length is set to 10 seconds; when the main frequency < 0.01Hz, the time window length is set to 120 seconds; and in other cases, the time window length is set to 60 seconds. This adaptive strategy ensures that each time window contains at least one complete fluctuation period.

[0084] The second step is to identify the effective fluctuation range:

[0085] For each time window, a sliding energy window is constructed with each time point as the center, and the square integral of the residual amplitude within each sliding energy window is calculated as the local energy value at that time point.

[0086] A sliding energy window refers to a window centered at the current time point t, extending forward and backward by a certain length. The window length is set to 2 seconds, from t-1 seconds to t+1 seconds. The local energy value is calculated using the formula: ∫[r(τ)²]dτ, with the integration interval [t-1s, t+1s], where r(τ) is the fluctuation residual value. The local energy value reflects the fluctuation intensity in the neighborhood at that time point. The squaring operation amplifies the contribution of large amplitude values ​​and suppresses small amplitude noise.

[0087] The energy threshold is calculated as a multiple of the historical residual energy mean. Specifically, residual data from the 24-hour period without events prior to the time window are taken, and the statistical average of their local energy values ​​is calculated. The energy threshold is set to three times the statistical average. The set of time points where the local energy value continuously exceeds the energy threshold is identified as the effective fluctuation interval, and the start and end times of this effective fluctuation interval are recorded.

[0088] Here, "continuous exceedance" refers to the existence of a continuous time interval where the local energy value at each time point within that interval is greater than the energy threshold, and the length of that interval is greater than the minimum duration threshold (set to 5 seconds). This process effectively distinguishes isolated noise points from persistent anomalies.

[0089] The third step is to determine the duration:

[0090] The duration of the residual is determined based on the start and end times of the identified valid fluctuation range; the duration reflects the complete time from the onset of the abnormal event to its return to normal.

[0091] For minor leaks, the duration typically fluctuates continuously from tens of seconds to several minutes, while the duration of transient noise is usually less than 1 second. By analyzing the duration characteristics, transient noise interference can be effectively filtered out.

[0092] Step 4: Calculate the cumulative energy value:

[0093] The amplitude square integral is calculated on the original fluctuation residual subsequence within the identified effective fluctuation interval to obtain the energy accumulation value corresponding to the effective fluctuation interval.

[0094] The cumulative energy value differs from the local energy value: the local energy value is an integral of a sliding window centered on each time point, used to identify the effective fluctuation range; the cumulative energy value is an integral of the original residual within the entire effective fluctuation range; the larger the cumulative energy value, the stronger the energy of the abnormal event, and the more severe the leakage may be.

[0095] Step 5: Obtain spatial propagation characteristics:

[0096] For the residual signals collected from different sensor nodes for the same suspected leakage event, the arrival time difference of the residual signals at each node is calculated. Combined with the pipe distance between the nodes, the propagation speed of the residual between adjacent nodes is calculated to obtain the spatial propagation characteristics of the residual.

[0097] The same suspected leakage event refers to multiple node signals confirmed to belong to the same event after spatial correlation and merging. For these node signals, the start time of the effective fluctuation interval in the residual signal of each node is extracted as the time when the event arrives at that node. Let the arrival times of the event at node i and node j be ti and tj, respectively, then the arrival time difference Δtij is: Δtij = |ti - tj|.

[0098] Query the pipeline distance Lij between node i and node j from the pipeline geographic information system. The pipeline distance refers to the curved distance along the pipeline direction, not the straight distance. Calculate the propagation speed vij of the residual between adjacent nodes as: vij=Lij / Δtij.

[0099] Propagation speed is a key parameter characterizing the spatial propagation characteristics of an event. According to the principles of fluid mechanics, the propagation speed of pressure waves in water supply pipes is usually in the range of 1000-1200 m / s. If the calculated propagation speed deviates significantly from this range, it may be due to interference or node association errors caused by other reasons. If the propagation speed is within this range, it further verifies the spatial consistency of the event.

[0100] The present invention provides a solution that, under conditions of dynamic change in fluctuation frequency, accurately extracts the duration, energy accumulation value, and spatial propagation characteristics that truly reflect the physical properties of minute leaks from the fluctuation residual sequence captured by multiple sensors.

[0101] In the above embodiments, spatial distribution features that meet the conditions are marked as suspected leakage features. However, during the marking process, the relationship between pressure fluctuations and flow rate changes is not consistent throughout the entire time window; normal and abnormal periods are mixed together; the temporal relationship between pressure and flow rate is difficult to determine; and leakage and normal event features overlap and are difficult to distinguish. To solve this problem, this invention provides a preferred method for identifying water supply network leakage, which is described below in conjunction with… Figure 5 Detailed explanation of the specific steps.

[0102] The first step is to identify the pressure drop range and the flow rate increase range:

[0103] For each time window, the start and end times of the pressure drop interval within that time window are first obtained from the extracted spatial distribution features. The pressure drop interval refers to the period during which the pressure value continuously decreases from the baseline level to the lowest point, with its start time denoted as tp1 and its end time denoted as tp2.

[0104] Meanwhile, extract the subsequence of flow data corresponding to this time window from the collected synchronous flow data, and identify the start time and end time of the rising interval of the flow data in this subsequence. The rising interval of flow is the period during which the flow value continuously rises from the reference level to the highest point, and its start time is denoted as tq1 and the end time is denoted as tq2.

[0105] The identification of the pressure drop interval and the rising interval of flow limits the effective calculation range for subsequent coupling analysis and excludes the interference of irrelevant time periods.

[0106] Step 2: Determine the overlapping period and data pairs:

[0107] Determine the overlapping period of the pressure drop interval and the rising interval of flow. The overlapping period is the part jointly covered by the two intervals on the time axis, with its start time ts = max(tp1, tq1) and end time te = min(tp2, tq2). When ts < te, there is a valid overlapping period; otherwise, it indicates that the pressure drop and the flow increase do not occur simultaneously, and this time window does not meet the conditions for further analysis.

[0108] Within the overlapping period, form data pairs (P(m), Q(m)) with the instantaneous pressure value and the instantaneous flow value at each time point, where m takes all the sampling times within the overlapping period. These data pairs constitute the original sample set for analyzing the relationship between pressure and flow.

[0109] Step 3: Calculate the Pearson correlation coefficient:

[0110] Calculate the Pearson correlation coefficient for all the data pairs constructed in Step 2 to obtain the negative correlation coefficient between the pressure drop and the flow increase; the value range of the Pearson correlation coefficient r is [-1, 1]. When r is negative, it indicates that the pressure and the flow are negatively correlated, that is, the pressure drop is accompanied by an increase in flow; when r is positive, it indicates that the pressure and the flow are positively correlated, that is, the pressure and the flow change in the same direction; when r is close to 0, it indicates that there is no obvious linear relationship between the pressure and the flow. For minor leaks, it is expected that r is negative and has a large absolute value.

[0111] Step 4: Calculate the time difference:

[0112] Compare the starting time tp1 of the pressure drop interval with the starting time tq1 of the flow increase interval, and calculate the time difference Δtpq: Δtpq = tp1 - tq1. The physical meaning of the time difference Δtpq is the amount of time by which the pressure drop precedes the flow increase. When the starting time of the pressure drop interval is earlier than the starting time of the flow increase interval, Δtpq is positive, indicating that the pressure change precedes the flow change; when the starting time of the pressure drop interval is later than the starting time of the flow increase interval, Δtpq is negative, indicating that the flow change precedes the pressure change; when both occur simultaneously, Δtpq is zero.

[0113] According to fluid mechanics principles, pressure waves caused by leakage propagate at the speed of sound in water, approximately 1000-1200 m / s, while flow rate changes depend on the water velocity and have a relatively delayed response. Therefore, a real leakage event should satisfy Δt > 0, meaning that the pressure drop occurs first.

[0114] Step 5: Mark suspected leakage characteristics:

[0115] Based on the calculated negative correlation coefficient r and the time difference Δtpq calculated in step 5-4, the spatial distribution characteristics that meet the conditions are marked to obtain suspected leakage characteristics.

[0116] The marking conditions include two criteria that must be met simultaneously:

[0117] First, the absolute value of the negative correlation coefficient is greater than a first threshold; in one embodiment, the first threshold can be set to 0.7, indicating a strong negative correlation between pressure and flow rate. When |r|>0.7, it indicates that the trends of pressure decrease and flow rate increase are highly consistent during the overlapping period.

[0118] Second, the time difference lies within the time delay range; the time delay range is determined based on the pressure wave propagation speed and typical sensor spacing, for example, it can be set to [0.1 seconds, 5 seconds]. When 0.1 seconds ≤ Δtpq ≤ 5 seconds, it indicates that the time by which the pressure drop precedes the flow rate increase conforms to the propagation law of pressure waves caused by leakage.

[0119] When both of the above conditions are met, the spatial distribution characteristics within the time window are marked as suspected leakage characteristics, and the characteristics and their corresponding time window locations are recorded.

[0120] This invention provides a preferred method for identifying water supply network leaks. By limiting the calculation scope to the overlapping period between pressure drop and flow increase intervals, the Pearson correlation coefficient is calculated to ensure that the correlation coefficient reflects the true relationship during abnormal periods. The time difference between the start time of pressure drop and the start time of flow increase is obtained and assigned positive or negative meanings, establishing a temporal criterion that conforms to the physical laws of leakage. Through the dual judgment conditions of a negative correlation coefficient absolute value greater than a first threshold and a time difference within the time delay range, a comprehensive verification is performed from both numerical and temporal dimensions, effectively distinguishing leakage events from normal water use events.

[0121] Reference Figure 6 As shown in the figure, this invention also provides a preferred method for identifying water supply network leaks, which performs correlation analysis on suspected leak characteristics based on pipe aging parameters and historical maintenance records, and identifies leaks based on the analysis results. The method includes the following steps:

[0122] The first step is to calculate the pressure amplitude extraction and energy decay rate:

[0123] From the suspected leakage features obtained in the first step, extract the pressure amplitude information of the feature at at least two adjacent sensor nodes; the pressure amplitude refers to the maximum amplitude of the suspected leakage feature within the effective fluctuation range, reflecting the intensity of the pressure fluctuation at that node.

[0124] Let the pressure amplitude of the upstream sensor node be A1 and the pressure amplitude of the downstream sensor node be A2. Since the energy gradually attenuates as it propagates along the pipeline, A1 > A2. Query the pipeline distance between these two adjacent sensor nodes from the pipeline network geographic information system and denote it as L.

[0125] Based on the pressure amplitude information and the pipe distance between nodes, the energy attenuation rate per unit distance propagating along the pipe is calculated using the following formula:

[0126] ;

[0127] Where η is the energy decay rate, representing the energy decay rate per unit distance.

[0128] The second step is to calculate the theoretical attenuation rate and obtain the deviation:

[0129] The pipeline status database is used to query the pipe aging parameters of the pipe sections where two adjacent sensor nodes are located. The pipe aging parameters include information such as pipe type, percentage of remaining wall thickness, service life of the pipeline, and corrosion rate. Among them, the percentage of remaining wall thickness is the most critical indicator, reflecting the current health status of the pipeline.

[0130] Based on the pipe aging parameters, the theoretical attenuation rate η0 of the pipe section is calculated using a pre-constructed physical model. The physical model is established based on laboratory test data and field measurement data, creating a mapping relationship between pressure wave attenuation rates for different pipe materials and different aging degrees. For ductile iron pipes, the theoretical attenuation rate model can be expressed as:

[0131] η0 = k × (100 / δ) × (100 / D);

[0132] Where δ is the percentage of remaining wall thickness, D is the pipe diameter, and k is the pipe material coefficient, which is obtained through experimental calibration; for ductile iron pipes, the value of k is 0.015-0.025 dB·mm.

[0133] Calculate the deviation Δη between the energy decay rate η and the theoretical decay rate η0:

[0134] Δη = |η - η0| / η0 × 100%;

[0135] When the deviation Δη is less than the second preset threshold, it is determined that the energy decay characteristic of this feature matches the aging degree of the pipe. The second preset threshold is set to 15% based on engineering experience, indicating that the difference between the measured decay rate and the theoretical decay rate is within a reasonable error range.

[0136] The third step is to determine the spatial location of the features:

[0137] Based on the time information of two adjacent sensor nodes and the distance between the nodes in the pipeline, the spatial position corresponding to the feature is calculated by the time difference of arrival localization algorithm, and the first spatial position of the feature on the pipeline is obtained.

[0138] Let t1 be the time when the feature arrives at the upstream sensor node, t2 be the time when it arrives at the downstream sensor node, L be the distance between the nodes, and v be the propagation speed of the pressure wave. Then the time when the feature occurs, t0, can be calculated as: t0 = (t1 + t2 - L / v) / 2.

[0139] The distance d from the feature location to the upstream node is: d = (t1 - t0) × v;

[0140] Then the first spatial position coordinate x of the feature can be expressed as: x = x 上游 +d;x 上游 The location coordinates of the upstream sensor node.

[0141] Step 4: Calculate the distance to the historical repair point:

[0142] Query historical maintenance records from the pipeline status database that are within a preset search radius of the first spatial location x, and extract the second spatial location of the historical maintenance point corresponding to each historical maintenance record.

[0143] The preset search radius is set to 100 meters based on engineering experience; searching within this radius ensures coverage of potentially related areas while avoiding the introduction of irrelevant points that are too far away. Historical maintenance records include information such as the location, time, type, and reason of each maintenance operation, with location coordinates stored in the form of pipeline station numbers or GIS coordinates.

[0144] For each historical maintenance point found, calculate the distance along the pipeline between its second spatial location and the first spatial location x of the feature. The distance along the pipeline refers to the curved distance along the pipeline, not the straight-line distance. The shortest path length along the pipeline between two points can be obtained by querying the pipeline topology table.

[0145] Step 5, Leakage Identification and Confirmation:

[0146] Based on the obtained energy attenuation matching and location overlap determination results, a comprehensive decision is made to identify leaks.

[0147] Specifically, when there exists any historical maintenance point such that the distance d along the pipeline 管道 If the distance is less than a threshold, the spatial location of the feature is determined to coincide with a historical maintenance point. This threshold is set to 50 meters based on engineering experience, indicating that the feature location and the historical maintenance point can be considered to be in the same area if they are within an acceptable range for engineering purposes.

[0148] When both of the following conditions are met, it is confirmed that this pure anomaly characteristic corresponds to a real micro-leak:

[0149] First, the energy decay rate deviation Δη is less than 15%, meaning that the energy decay characteristics match the aging degree of the pipe.

[0150] Second, there exists at least one historical repair point that makes d 管道 Less than 50 meters, meaning the spatial location coincides with the historical maintenance point.

[0151] After confirming the leak, output the first spatial location as the coordinates of the minor leak.

[0152] The preferred water supply network leakage identification method provided in this invention establishes a physical model with the percentage of remaining wall thickness as the core parameter, transforming static pipe aging parameters into dynamic theoretical attenuation rates; accurately calculates the feature spatial location using a time-of-arrival (TOA) positioning algorithm; associates historical maintenance points by calculating the distance along the pipeline and using preset search radii and overlap thresholds; and comprehensively verifies the identification of minute leaks from three dimensions: physical characteristics, spatial location, and historical patterns, using the dual judgment conditions of energy attenuation rate deviation being less than a threshold and distance from historical maintenance points being less than a threshold.

[0153] Obviously, the above embodiments are merely illustrative examples for clear explanation and are not intended to limit the implementation. Those skilled in the art will recognize that other variations or modifications can be made based on the above description. It is neither necessary nor possible to exhaustively list all possible implementations here. However, obvious variations or modifications derived therefrom are still within the scope of protection of this invention.

Claims

1. A method for water supply network leak identification for urban lifeline engineering, characterized by, include: Collect pressure fluctuation data and synchronous flow data from sensor nodes in the water supply network; The difference between the pressure fluctuation data and the benchmark fluctuation template constructed based on historical data is calculated point by point to obtain the fluctuation residual sequence; Feature extraction is performed on the fluctuation residual sequence to obtain its spatial distribution characteristics; Combining the spatial distribution characteristics and the synchronous flow data, the negative correlation between pressure drop and flow increase, as well as the time difference between pressure drop and flow increase, are calculated; based on the calculation results, the spatial distribution characteristics that meet the conditions are marked as suspected leakage characteristics. The suspected leakage features are correlated based on pipe aging parameters and historical maintenance records, and leakage is identified based on the analysis results. The correlation analysis includes: determining the deviation between the energy attenuation rate of the suspected leakage feature propagating along the pipeline and the theoretical attenuation rate corresponding to the pipe aging parameters, as well as the spatial location of the suspected leakage feature and the distance between the historical maintenance points. Constructing a benchmark fluctuation template based on historical data includes: extracting historical pressure fluctuation data from a historical database that have the same seasonal type, weekday attribute, and time period as the current period, forming a historical sample set; performing a Fast Fourier Transform on each historical pressure fluctuation data in the historical sample set to obtain the frequency domain characteristics of each historical data, removing historical data with abnormal harmonic components in the frequency domain characteristics to obtain a clean historical sample set; clustering the clean historical sample set to divide it into multiple typical fluctuation pattern categories, selecting the category with the largest number of samples as the benchmark fluctuation pattern; performing point-by-point statistical averaging on the historical pressure fluctuation data in the benchmark fluctuation pattern to obtain the mean curve at each time point; and smoothing the mean curve to obtain the benchmark fluctuation template.

2. The water supply network leakage identification method for urban lifeline engineering according to claim 1, characterized by, The difference between the pressure fluctuation data and the benchmark fluctuation template constructed based on historical data is calculated point by point to obtain the fluctuation residual sequence, including: The pressure fluctuation data and the benchmark fluctuation template are compared point by point at the same time point to obtain the initial residual sequence; Sliding window standard deviation analysis was performed on the initial residual sequence to calculate the local standard deviation at each time point, and time points where the local standard deviation exceeded the threshold were marked as suspected mutation points. Using the suspected mutation point as the dividing point, the initial residual sequence is divided into several continuous segments, and polynomial fitting is performed on each continuous segment to obtain the fitting baseline corresponding to each segment. The corresponding fitted baseline is subtracted from the initial residual of each continuous segment to obtain the corrected residual sequence. The corrected residual sequence is then zero-mean processed to obtain the fluctuation residual sequence.

3. The water supply network leak identification method for urban lifeline engineering according to claim 1, characterized by, The spatial distribution characteristics include the duration of the residual, the cumulative energy value, and the spatial propagation characteristics.

4. The water supply network leakage identification method for urban lifeline engineering according to claim 3, characterized by, Feature extraction is performed on the fluctuation residual sequence to obtain spatial distribution features, including: The fluctuation residual sequence is segmented by a sliding time window to obtain multiple fluctuation residual subsequences; For each time window, a sliding energy window is constructed with each time point as the center. The square integral of the residual amplitude within each sliding energy window is calculated as the local energy value at that time point. The set of time points where the local energy value continuously exceeds the threshold is identified as the effective fluctuation interval. The duration of the residual is determined based on the start and end times of the effective fluctuation range; The cumulative energy value of the residual is obtained by performing amplitude square integral on the fluctuation residual subsequence within the effective fluctuation range. For the residual signals collected at different nodes for the same suspected leakage event, the arrival time difference of the residual signals at each node is calculated. Combined with the pipe distance between the nodes, the propagation speed of the residual between adjacent nodes is calculated to obtain the spatial propagation characteristics of the residual.

5. The water supply network leak identification method for urban lifeline engineering according to claim 1, characterized by, The calculation of the negative correlation between pressure drop and flow rate increase includes: for each time window, forming data pairs of instantaneous pressure and instantaneous flow rate values ​​at each time point within the overlapping period of the pressure drop interval and the flow rate increase interval, calculating the Pearson correlation coefficient of all data pairs, and obtaining the negative correlation coefficient between pressure drop and flow rate increase.

6. The water supply network leak identification method for urban lifeline engineering according to claim 5, characterized by, Calculating the time difference between pressure drop and flow increase includes comparing the start time of the pressure drop interval with the start time of the flow increase interval, and calculating the time difference between the two. When the start time of the pressure drop interval is earlier than the start time of the flow increase interval, the time difference is a positive value, and otherwise it is a negative value.

7. The water supply network leak identification method for urban lifeline engineering according to claim 6, characterized by, The conditions for marking the suspected leakage features include: the absolute value of the negative correlation coefficient is greater than a threshold, and the time difference is within the time delay range.

8. The method for identifying leakage in water supply networks for urban lifeline projects according to claim 1, characterized in that, Determining the energy attenuation rate of the suspected leakage feature propagating along the pipeline includes: Extract the residual amplitude of pressure fluctuation between two adjacent sensor nodes and the pipe distance between the nodes; The energy decay rate is determined to be: ; η represents the energy decay rate; A1 represents the pressure fluctuation residual amplitude of the upstream sensor node; A2 represents the pressure fluctuation residual amplitude of the downstream sensor node; L represents the pipe distance between two adjacent sensor nodes.

9. The method for identifying leakage in water supply networks for urban lifeline projects according to claim 8, characterized in that, Based on pipe aging parameters and historical maintenance records, a correlation analysis is performed on the suspected leakage characteristics. Leakage is then identified based on the analysis results, including: The pipe aging parameters of the pipe segments where the two adjacent sensor nodes are located are queried, and the theoretical attenuation rate of the pipe segment is determined based on the pipe aging parameters; and the deviation between the energy attenuation rate and the theoretical attenuation rate is calculated. In response to the deviation being less than a threshold, the first spatial location corresponding to the suspected leakage feature is determined based on the time information of the two adjacent sensor nodes and the pipe distance between the nodes; Query historical maintenance records whose spatial location corresponds to the suspected leakage feature is within the search radius to obtain the second spatial location of the historical maintenance point; If any of the historical maintenance points exists such that the distance between the first spatial location and the second spatial location is less than a threshold, then a leak is identified.

Images