Abstract
The research about online monitoring and leakage automatic location of water distribution networks (WDN) has a wide range of applications that include water resource protection, monitoring, and allocation. Variational mode decomposition (VMD) and cross-correlation (CC) based leakage location is a popular and effective method in WDN. However, the value of K intrinsic mode functions (IMFs) based on VMD decomposition needs to be determined artificially, which affects the separation effect of signal frequency band characteristics directly. Hence, this work proposes an adaptive method to determine the parameter K of leakage vibration signal’s IMFs, which will be applied to automatic leakage location in WDN. Firstly, the number of saddle points in the frequency domain envelope of the sampled signal in different step sizes is calculated. The parameter K is determined according to the curvature change of the number of saddle points and the sampled signal. Finally, the selective IMFs are reconstituted into a new signal, which can determine a leak position using CC based time-delay estimation (TDE). To verify the effectiveness of the proposed algorithm, the different methods based on EMD and Fast ICA are compared. The experimental results demonstrate that the proposed parameter K value adaptive VMD (KVA-VMD) decomposition method is more suitable for leakage location in WDN.
1. Introduction
Leakage in urban WDN results in water resource losses and potential water pollution incidents to invade harmful substances [1, 2]. Hence, it is necessary to monitor the operation status of WDN in real time and locate the leakages in time to guarantee its safe and reliable operation [3–5]. So far, much research has been conducted for leak detection and location methods in WDN. Since the 1990s, leak location based on CC analysis of vibration signal has gradually been regarded as one of the most popular technologies in WDN, which has been widely used to locate leakage in the buried pipeline [6–10]. In real leakage location practice, the signal-to-noise ratio (SNR) of sampled data by the sensor is meager due to the influence of environmental noise (car whistle, road vibration, valve noise, and other factors), which directly affects the CC based TDE method’s reliability and errors of leak detection.
TDE can enhance the SNR of leak signals and reduce leak location errors [8]. Knapp and Carter proposed a generalized cross-correlation (GCC) method, which enhances the signal in the band with higher SNR and suppresses the signal outside these bands by prefiltering leak signal before performing the CC method [11]. However, it is difficult to determine the frequency bands of leakage signals in various pipeline situations (pipeline material, pressure, leakage aperture, etc.) and surrounding environment (car whistle, road vibration, valve noise, etc.), which results in an inability to design a prefilter for GCC to improve the SNR of leakage signal. An adaptive TDE algorithm based on the least mean square (LMS) is proposed in which GCC is replaced by an adaptive filter [12, 13]. The filter can be updated iteratively by using the leakage signal’s characteristics, and its correlation function can be obtained directly from the system parameters instead of any prior knowledge. However, LMS-TDE is a biased estimation, and the error of the leak location method based on LMS-TDE is proportional to the SNR. Hence, the application of LMS-TDE is limited to low SNR.
Most of the traditional analysis methods only consider a single narrow-band analysis method instead of each component’s differences in the leak signal. However, the leakage signal produced by a jet of fluid under high pressure is more complicated than an ordinary acoustic signal. Thus, traditional methods cannot achieve high accuracy in leak detection, especially in multiple leaks coupling [14–16]. The wavelet analysis method is proposed to suppress ambient noise by selecting the appropriate wavelet cardinality adaptively according to the leakage signals and the ambient noise [17–19]. However, the selection of wavelet cardinality is still a troubling question for the uncertainty of empirical parameters, which are varied with the pipeline situation and the ambient noise. Thus, the method based on empirical mode decomposition (EMD) is proposed, which provides a selection of analysis for nonstationary random signals [20, 21]. However, the EMD algorithm still possesses some inherent drawbacks such as the limitations in handling noise and sampling sensitivity, the limitations of mathematical methods limiting the theoretical development in this field, and the trouble of improving decomposition stability.
Furthermore, mode aliasing caused by high resolution is also a persistent problem. Dragomiretskiy and Zosso proposed the VMD methodology for the analysis of nonstationary signals in 2014 [22]. Compared to EMD, multicomponent signals can be nonrecursively decomposed by VMD into several components, which control decomposition convergence conditions. The decomposition process of VMD can effectively eliminate the mode aliasing phenomenon with good noise immunity. Therefore, it is applied in many fields such as fault detection, leakage signal analysis, signal processing scenarios, and signal separation of ambient noise [23–29]. However, VMD needs to set the decomposition number K of IMFs and penalty factor in advance, which directly affects the accuracy of decomposition [30]. The K value adaptive selection method for enhancing VMD performance is proposed in this work, where mode components with high SNR are separated to reduce the automatic leak location error in WDN. The experiment was conducted to compare the proposed parameter K value adaptive VMD leak location errors with CC and LMS-TDE.
2. Pipeline Leakage Signal Model
The vibration signal acquired by the sensor from the leaking pipeline may be composed of the parts shown in Figure 1. The acoustic signals (source signals) c1 and c2 come from upstream and downstream noise, respectively, c3 is the reflected noise produced by plugging branch pipe, c4 is the propagation noise of the bypass branch pipe and the noise generated by the pipe shunt, c5 is the noise produced by pipe throttling, c6 is the noise coupled to the pipe by external vibration, and and are the measurement noises of the sensor. Hence, the signal model is

Among them, c1–c5 may be mistaken for leakage signal in the correlation detection method. There are many forms of external noise c6 in the real environment, and its statistical characteristics are challenging to estimate. Also, and also changed according to the site environment. Although c1–c5 are correlated noises, the signal energy of which is smaller than that of leakage source. As for external noise c6 and measurement noises and , they are generally considered to be uncorrelated noises. Hence, the signal model is usually simplified as in the following equation:where s(t) is the unknown source signal and and are additive noises. Source signal s(t) and noise signals and are pairwise noncorrelation zero-mean stationary random signals. τ is the time delay of the source signal received by the two sensors, and α is the ratio of signal attenuation due to the different propagation paths of the source signal. If the distance l between the two sensors of the measured pipeline is known, τ can be calculated by LMS-TDE or GCC-TDE method to determine the leak's location (as shown in equation (3)). However, there are many leaks in WDN, or if one signal in the correlation noises c1–c5 is large enough, the simple model will not achieve satisfactory results. Consequently, a more effective leak location scheme is desirable for WDN.where c is the propagation speed of the signal in WDN.
3. K Value Adaptive VMD Method
3.1. VMD Principles
The VMD method based on Hilbert transformation (HT) and the variational problem-solving process can sufficiently suppress the defects of the EMD method. The algorithm’s principle is to decompose a signal f into K IMFs for minimizing the sum of the estimated bandwidths.
For each IMF uk, the relative analytic signal is obtained by HT, and its unilateral spectrum is obtained:where δ(t) is the impact function and represents the K IMFs which are components of the signal being decomposed. denotes convolution, and it also represents convolution in the following formula.
The estimated central frequency e−jwt is mixed with each model analytic signal, and the spectrum of each model is modulated to the corresponding fundamental frequency band:and, in formula (5), is the center frequency of each model component uk.
Calculate the norm L2 of the gradient’s square of the demodulated signal above and estimate each model component’s bandwidth. The corresponding expression of the constrained variational model is
According to the solution of the constrained variational problem, the quadratic penalty factor α and multiplication operator λ(t) are introduced, and α is a positive number large enough to guarantee the reconstruction accuracy of the signal in the presence of Gaussian noise. The Lagrange operator makes the constraint condition strict. The extended expression is given in the following equation:
The above variational problems are obtained by alternating direction multiplier algorithm (Alternate Direction Method of Multipliers, ADMM). The saddle point of the extended Lagrange expression is obtained by alternating updating (mode function of the first cycle of , center frequency, and multiplication operator). The process can be described as in the following equation:
In terms of the Plancherel theorem, the norm problem converted into a formula for the isometric Fourier transform. Among them, converted to the frequency domain and was replaced with , and the result converted to the integral form of negative frequency interval.
Following the same process, the update mode of the center frequency can be obtained:
In equation (10), τ is a time constant, which is usually 0. The multiplication operator is updated as follows:
In equation (11), ε is a sufficiently small positive number. The convergence conditions are
The challenge of this algorithm is how to choose an effective parmeter IMF to restore the information of the original signal. However, there is no principle for the selection of K value. The K value is too small to decompose the signal effectively. If the K value is too large, the signal will be discontinuous, which will affect the accuracy of subsequent processing. In this paper, a method is proposed for determining the K value adaptively by calculating the change of signal envelope curvature by considering that the leakage signal in WDN has the same spectrum distribution.
3.2. K Value Adaptive Processing
Equation (10) shows that the update of the center frequency is obtained by alternating update after equidistant transformation. In contrast, the center frequency is reflected by the instantaneous frequency mean of the component, and the number of saddle points is the K value. Given that the sampling frequency of the leakage signal is 5000 Hz and the data length is 8192, Figure 2 shows its sampling signal and unilateral spectrum.

(a)

(b)
Figure 3 presents the saddle point transformation relationship of the leakage signal’s frequency domain signal after Hilbert transformation at different step sizes. The trend of the number of saddle points after transformation tends to converge with the increase of step size. However, the changing trend of saddle points number has certain volatility rather than monotonic decline. With the increase of step size, the number of saddle points gradually tends to a stable value. If the step size increases, the number of saddle points caused by envelope distortion will increase first and then decrease. Finally, with the further increase of step size, the number of saddle points tends to 1. Hence, as long as the saddle points number’s stable interval is solved, the appropriate step size and K are determined.

Figure 4 shows the relationship among the step size, the number of saddle points, and continuous interval. The duration interval of saddle points also tends to increase with the increase of step size. When the step size starts from 107, the saddle points number’s continuous interval achieves the longest, and the corresponding saddle points number is 6. Once the step size increases, the stable range of the number of saddle points decreases sharply (when the step size is 142, the continuous interval is the minimum), which indicates that the envelope shows distortion after transformation. Figure 5 shows the HT transform of the frequency-domain signal at step sizes of 107 and 142. It can be seen from the figure that when the step size is 142, the envelope after HT transformation performs noticeable distortion. Hence, the step size between 107 and 142 is the best. Because the smaller the step size is, the smaller the envelope error is, the optimal step size is 107, and the corresponding K value is determined to be 6. This process of finding the longest stable interval can be converted to finding the maximum rate of change (derivative) of the stable interval. When the stable interval’s length increases monotonically, the step with a sharp decrease becomes the maximum distortion location, and the step with the longest stable interval before this step is the optimal step. The calculation method is shown in the following equation:


After the K value is determined, the VMD model is used to decompose the signal, as shown in Figure 6.

After model decomposition, each mode is sorted in order of frequency from low to high. The leakage signal sampled by another sensor can be decomposed automatically in the same way. According to the corresponding time-domain signal reconstructed by each IMF, leakage location can be carried out by CC, GCC-TDE, LMS-TDE, and other methods.
4. Experiments of Leak Location in WDN Using GCC-TDE, LMS-TDE, and KVA-VMD
4.1. Introduction to the Experimental Environment
A section of buried water supply pipeline with leak points on the campus is selected to verify the method’s effectiveness. The pipeline diagram is shown in Figure 7(a). The pipeline is a galvanized steel pipe, with a total length of 91.12 meters, a hole diameter of 75 mm, and a buried depth of 65 cm. The pipe’s water supply pressure is about 3.0 MPa, and valves are installed at both ends. The acquisition device for leak signal is the SS-01 sensor developed by Wizepipes Company. The device has a sampling frequency of 5000 Hz, a sampling accuracy of 12 bits, and sampling data of 65536 at a time. It can sample vibration signal ranges of 0∼0.5 G in real time. The two sensors with the synchronous acquisition bias of ±30 ns were installed on the valves at both ends of the pipeline (as shown in Figure 7(b)). The collected vibration signals are uploaded to the host computer through the wireless network. The host computer is equipped with data acquisition and analysis software to receive and store the collected vibration data. At the same time, the host computer can control two wireless acceleration sensor nodes to collect data synchronously and perform the data (as shown in Figure 7(c)).

4.2. Experimental Results and Discussion
In the experiment, the vibration signals collected from the valves at the two ends are uploaded to the host computer through wireless signals. Choose a group of signals for analysis (as shown in Figure 8). It can be seen from the figure that each segment of vibration signal lasts about 14 seconds, among which node 1 has large-signal energy and rich spectrum distribution, with energy gathering in 200–500 Hz, 800–1000 Hz, and 1300–1500 Hz frequency bands, while the signal spectrum collected by node 2 is distributed in 200–500 Hz and 2100–2200 Hz frequency bands. According to wave law, the high-frequency component is more easily absorbed and scattered than the low-frequency component. Therefore, the high-frequency signal transmitted from the leakage point to sensor 2 has been basically submerged by the pipe network’s background noise, which indicates that the leakage point should be close to node 1. Simultaneously, the peak of 2100–2200 Hz spectrum energy of node 2 may be caused by the external environment. If the TDE location is carried out directly, the error will be relatively large.

To determine the leak location by the KVA-VMD method, the number K of IMFs that each signal to be decomposed must be calculated first. Based on formula (12), the longest saddle point’s duration before the maximum distortion can be obtained (as shown in Figure 9). As can be seen from the diagram, the optimal HT transformation step sizes for node signals 1 and 2 are 1730 and 2563, respectively, and the corresponding saddle points are 5 and 4. Thus, it is determined that the optimal numbers of IMFs for VMD decomposition of node signals 1 and 2 are 5 and 4, and the decomposition results are shown in Figures 10 and 11. The two decomposed IMFs are compared in the spectrum (shown in Figure 12). The IMF1 spectrum overlaps, and the rest of the IMF spectrum does not overlap, so IMF1 can be selected as the reconstructed signal to locate the leakage location. Similarly, GCC and LMS methods can be used for location analysis (as shown in Figure 13). It can be seen from the graph that the result of leak location processed by the KVA-VMD method has the highest SNR, while the result of leak location processed by the LMS method has the lowest SNR.





The offset base point is known, and the sampling frequency is 5000 Hz, so the offset time τ can be calculated. Based on the wave equation, the sound velocity c of the galvanized steel pipe is 1020 m/s, and the distance between the two sensors is 91.12 meters. According to formula (3), the leakage location distance from valve 1 is 11.70 m, 11.49 m, and 11.59 m. The calculation results show that the leakage location errors obtained by the KVA-VMD, GCC, and LMS methods are 0.50 m, 0.71 m, and 0.61 m.
To further verify this method's effectiveness, 15 sets of leak data were collected at the valves of both ends to compare the leak locating performances of KVA-VMD, GCC, and LMS. Figure 14 compares the location errors of 50 sets of collected data for leak location using the KVA-VMD, GCC, and LMS methods. The results show that the relative error of KVA-VMD based leak location method is small, but that of the GCC method is large because the leak vibration signal contains multimodal and dispersed components. The average relative errors of leakage locating using the GCC, LMS, and KVA-VMD methods are 24.86%, 15.84%, and 3.85%, respectively. Therefore, the proposed KVA-VMD method is more suitable for leak location in WDN.

5. Conclusions
This work considers the challenge that many parameters need to be set artificially in the WDN monitoring and maintenance process. We proposed the KVA-VMD based leak location method, and its performance is compared with GCC and LMS through the experiment, suggesting the following:(1)The leakage signal in WDN has the optimal HT transformation step size. Under the optimal transformation step size, the number of saddle points is consistent with the IMF number K of VMD optimal decomposition, which can be used for VMD automatic decomposition.(2)Leak signal in WDN contains multimode and dispersive components, which can be considered narrow-band noise, leading to large leak location error. The attenuation rate of the low-frequency signal below 800 Hz is much slower than that of the high-frequency signal. This method can reconstruct the frequency band signal for leak location and reduce leak error.(3)The average relative errors of leakage locating using the GCC, LMS, and KVA-VMD methods are 24.86%, 15.84%, and 3.85%. Therefore, the proposed KVA-VMD method is more suitable for automatic leak location in WDN without artificial preset parameters.
The algorithm also can be applied to automatic leakage location and maintenance of the gas pipeline, heating pipeline, and other pipelines with certain pressure.
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This work was supported by the Scientific Research Project of the Hubei Education Department (Grant no. B2018084) and the Natural Science Foundation of Hubei Province (Grant no. 2019CFB250).