Leak Detection and Location of Water Pipes Using Vibration Sensors and Modified ML Prefilter.

This paper proposes a new leak detection and location method based on vibration sensors and generalised cross-correlation techniques. Considering the estimation errors of the power spectral densities (PSDs) and the cross-spectral density (CSD), the proposed method employs a modified maximum-likelihood (ML) prefilter with a regularisation factor. We derive a theoretical variance of the time difference estimation error through summation in the discrete-frequency domain, and find the optimal regularisation factor that minimises the theoretical variance in practical water pipe channels. The proposed method is compared with conventional correlation-based techniques via numerical simulations using a water pipe channel model, and it is shown through field measurement that the proposed modified ML prefilter outperforms conventional prefilters for the generalised cross-correlation. In addition, we provide a formula to calculate the leak location using the time difference estimate when different types of pipes are connected.

1. Introduction

This paper deals with the leakage detection and location of water pipes which is one of the main concerns in the water management field. It has been reported that the amount of non-revenue water reaches 15.2–35.1% of drinking water supply in United States, 6.5–24.6% in Europe, and 4.3–27.0% in Korea [1,2,3]. These losses can be classified into unbilled public usage, apparent losses including unauthorized consumption and metering inaccuracies, and real losses through overflows at storage tanks and burst leaks in distribution pipelines caused by bad connections, pipe corrosion, and physical damages. It is especially very important to detect the burst leaks in buried water pipes in order to reduce water production costs as well as to protect public safety [1,4,5].

Most of the water pipes are buried underground, making it difficult to find the location of leaks. For this reason, water leakage has usually been detected when water flows out of the ground due to massive leaks in pipes. As a preemptive detection method, an acoustic-based approach is used that an experienced expert scans a suspected area by listening to leak sounds, yet this scheme is extremely time-consuming and the accuracy highly depends on the detection skills of the personnel [6,7]. As a more systematic approach, a supervisory control and data acquisition system (SCADA) has been considered to periodically measure the flow rate or water pressure using permanently installed sensors [8,9]. This technique has evolved into model-based leak detection that compares the measured data taken on the water network with the predicted values from flow models [10,11], and a model falsification method is introduced for efficiently detecting leaks in large-scale water distribution networks [12,13]. These methods identify the areas with significant water leakage, so additional techniques are required to find the exact position of leaks [8,9,12]. In addition, ground penetration radar can be used to locate leaks in water pipes by inspecting the ground surface and detecting voids in the soil caused by leaking water; however, the use of this approach is limited to specific areas because of high costs and high sensitivity to underground metal objects [4,14]. In [15,16], a non-invasive detection method based on the time domain reflectometry (TDR) has been investigated. This approach finds the position of leaks using the change of the reflection coefficient near leakage points. Because the measured reflection coefficient is sensitive to the temperature and material type around the pipe, robust leak detection is not easy.

The most popular approach to water leak detection is to use acoustic/vibration sensors or pressure transducers attached to the surface of a pipe [5,17,18,19,20,21,22,23,24,25,26,27,28,29]. When two sensors are used, the leakage location can be determined by estimating the time difference through correlation of the receive signals. The performance of various fast Fourier transform (FFT)-based algorithms is compared in terms of the time difference estimation accuracy in [18,19,20,21]; the time difference estimation schemes based on spectral transform are further enhanced by the short-time Fourier transform (STFT) [22] and the wavelet transform [23]; and the correlation-based leak detector has been verified via hardware implementation [24]. In addition, a leak detector is proposed to estimate the time difference using an entropy algorithm [25]. When a single sensor is used, the leakage can be detected using the fact that the change of the frequency spectrum is highly dependent on the change of leakage volume [26]; however, it is not easy to employ this method in real environments due to the sensitivity to the change of frequency response. In [27,28,29], leak detection schemes with single sensor have been investigated based on the wavelet transform techniques; however, it is difficult to adopt in practical environments because they require reference vibration signals without leaks or the detection performance is vulnerable to the system noise. Additionally, other authors are investigating the use of the filter diagonalization method [30], the radio frequency identification (RFID)-based leakage monitoring system [31], and magnetic flux leakage detection [32]. However, because of the excessive complexity, the difficulty of installation and management, and the technical immaturity, these techniques are rarely used in specific areas.

In this paper, we focus on the time difference estimation method using two vibration sensors attached to water pipes. When the power spectral densities (PSDs) and the cross-spectral density (CSD) of the receive signals are known, the optimal time difference estimator is given by the generalised cross-correlation (GCC) with the maximum likelihood (ML) prefilter derived in [18]. Practically, the PSDs and the CSD include certain estimation errors, and thus the ML prefilter may perform worse than the smoothed coherence transform (SCOT) and Wiener processors [21]. In order to improve the performance of the ML prefilter with erroneous PSDs and CSD, we propose a modified ML method that exploits a regularisation factor to mitigate the estimation errors of the PSDs and CSD. We derive a theoretical variance of the time difference estimation error in the discrete-frequency domain through modification of the variance model in the continuous frequency domain [18], and then find the optimal regularisation factor that minimises the theoretical variance in practical water pipe channels. The time difference estimation accuracy of the proposed method is assessed through computer simulations under the channel similar to an actual water pipe. By implementing the vibration sensor in hardware and the proposed detection method in software, respectively, we carry out field measurement in buried water pipes deployed for household use to verify the performance of the proposed leak detector. Numerical simulations and field tests show that the proposed modified ML method performs better than existing correlation-based techniques, in terms of minimising the root mean square deviation (RMSD) of the time difference estimation error and maximising the peak-to-average ratio (PAR) of the cross-correlation function.

This paper is organized as follows. In Section 2, we derive a formula to find the leak location using a time difference between two vibration sensors when two different kinds of pipes are connected. In Section 3, conventional time difference estimation methods are introduced with the signal model, and the modified ML prefilter is proposed for GCC based on the theoretical variance of time difference estimation error. The performance of the proposed prefilter is compared with those of existing techniques through computer simulations and field measurement in Section 4 and Section 5, respectively. Finally, the conclusions are provided in Section 6.

2. Locating Leaks Using Time Difference

When water is leaked in pipes, internal fluid flows out to generate vibration, and it propagates in both directions from the leak point. As shown in Figure 1, vibration sensors are attached to both ends of the pipes including of a leak point, and thus a difference occurs in the vibration wave reception time due to the difference in distance between the leak point and the sensor. Considering the practical pipeline structure, the material and size of a pipe may be different on both sides of a joint. Define the time difference as where and mean the time for the vibration wave to reach sensors 1 and 2 from the leak point, respectively. The time difference can be estimated by correlating the receive signals of sensors 1 and 2. For the time being, it is assumed that is known, and the estimation of is considered in the following section. Suppose that the pipe 1 has a leak as shown in Figure 1a, i.e., . Then, and are given by where x is the distance between the sensor 1 and the leak point; and and denote the sound speed in pipes 1 and 2, respectively. From Equations (2) and (3), the time difference is expressed as

Figure 1

Leak location estimation using the time difference of signals received at sensors, when two different pipes are connected: (a) when pipe 1 has a leak; (b) when pipe 2 has a leak.

By rearranging Equation (4) with respect to x, we have where . Similarly, when pipe 2 has a leak as shown in Figure 1b, the time difference is expressed as and the distance between leak point and the sensor 1 is given by where . In summary, the leak location from the sensor 1 can be computed as where and .

The speed of sound in a water-filled pipe depends on the properties of a pipe wall that the vibration wave is traveling through. In [33], the velocity of a fluid-dominated wave is expressed as where is the bulk modulus of the contained medium; a, h, and E are the inner radius, the thickness, and the elastic modulus of the pipe, respectively; is the density of the surrounding medium; is the angular velocity; and is the free-field internal fluid velocity [34], given by where is a stiffness coefficient of the material. We use the database of the sound speed created for various types of pipes by theoretical calculations and experimental corrections provided by [35], and Table 1 presents a part of the database. Practically, the sound speed varies considerably according to the geometry and material type of the pipe. If a pipe is composed of mixed media, it is difficult to calculate the sound speed. Thus, we only consider a single material pipe.

Table 1

Sound velocity according to the pipe material and diameter [35].

Material	Diameter (mm)	Velocity (m/s)
	40	565
Polyvinyl Chloride (PVC)	80	540
	150	530
	150	1220
Cast-Iron	250	1160
	350	1120
	25	1375
	40	1350
Steel	60	1330
	90	1286
	150	1200
	250	1150

3. Time Difference Estimation

3.1. Conventional Time Difference Estimation Methods

In Figure 1, define and as the sampled receive signals of sensors 1 and 2 at time , respectively, where is the sampling period. When there is a leak between sensors 1 and 2, and are given by where is the vibration wave generated from the leak point; d is the time difference in samples between received signals at sensors 1 and 2 that satisfies ; and ⊗ means the discrete-time convolution operator. For , is the discrete-time channel impulse response between the leak point and the sensor i, and is a stationary random process representing the background noise at the sensor i. It is assumed that and are two stationary random signals uncorrelated with , and also assumed that and are uncorrelated. The time difference can be estimated using the cross-correlation of and defined as where is the time delay of , denotes the expectation operator, and means the complex conjugate.

The cross-correlation function can be computed in the time domain by directly correlating and because has the peak value at . This approach, however, is very sensitive to the non-leak vibration noises, and requires high computational load when the sample size is large. For these reasons, the time difference estimation method in the frequency domain is widely used based on the relationship between the cross-correlation function and the CSD [18,19,20,21,22,23,24]. Let us define as the N-point discrete Fourier transform (DFT) of a sequence . Then, the CSD for (13) is written as where and for ; is the PSD of ; and . By computing and taking its inverse DFT (IDFT), the cross-correlation function can be simply computed when the sample size is large.

To further enhance the frequency domain estimation method by sharpening the peak in the cross-correlation function, prefiltering is performed before the IDFT. In this case, the cross-correlation function is computed using the GCC method [18], defined as where means the N-point IDFT of a sequence and is a frequency domain prefilter. Figure 2 describes the procedure to compute the cross-correlation function in the GCC method, and Table 2 presents the prefilters for the GCC proposed in [18,19,20], where and are the PSDs of and for , respectively; and is the coherence function between and defined as

Figure 2

Block diagram to compute using the generalised cross-correlation method.

Table 2

Conventional prefilters for the generalised cross-correlation method [18,19,20].

Prefilter Type	Weight Function
Roth	1S_x_1(k)
Smoothed coherence transform (SCOT)	1S_x_1(k)S_x_2(k)
Phase transform (PHAT)	1|S_x_1x_2(k)|
Eckart	|S_x_1x_2(k)|S_w_1(k)S_w_2(k)
Wiener	|γ(k)|^2
Maximum-likelihood (ML)	|γ(k)|^21−|γ(k)|^21|S_x_1x_2(k)|

In the Wiener and ML prefilters, the background noises and are assumed to be zero-mean Gaussian random variables with variances and , respectively, to make the signal models (11) and (12) mathematically tractable. As pointed out in [12,13], is not zero due to the sensor bias, and the modeling error included in may not be expressed as a Gaussian distribution. The effect of the bias in can be removed by sensor calibration or subtracting the mean value from . In addition, for a large N, tends to be Gaussian by the central limit theorem, even if is not a Gaussian random variable. Hence, the assumption that and are Gaussian is not a strong requirement ([36], Ch.3A).

The Roth prefilter proposed by Peter R. Roth normalises the CSD by , while the smoothed coherence transform (SCOT) prefilter normalises the CSD by . The Roth and SCOT prefilters can be used only in the high signal-to-noise ratio (SNR) region because the peak of the cross-correlation function is spread out after prefiltering. To avoid the peak spreading, the phase transform (PHAT) normalises the CSD by so that the prefiltered CSD has the unit amplitude for all k. However, the peak estimation performance is degraded due to the noise enhancement in the frequency band with low SNR. The Eckart prefilter maximises the ratio of the mean correlator output power to the noise variance of the correlator output, and thus it suppresses the frequency band with low SNR. This method is difficult to use in a practical system because it is not easy to estimate the noise PSDs and from the receive signals with leaks [18,20,21]. The ML prefilter weights the CSD according to the SNR in each frequency band to minimise the variance of the time difference estimate, yet it may critically overweight the CSD in the frequency band where is close to one [18]. The Wiener prefilter multiplies the CSD by the weight , so that it mitigates the overemphasising problem at certain frequencies [19].

3.2. Proposed Time Difference Estimation Scheme Using the Modified ML Method

For practical implementation of the prefilter, it is required to estimate the PSDs and as well as the CSD . Additionally, the Eckart prefilter necessitates the estimation of the noise PSDs and . Suppose that we have long sequences and measured at sensors 1 and 2. Then, the PSDs and CSD of and can be estimated by the non-overlapped Welch’s method [37] as follows: where ; M is the number of blocks; and . Using the power spectrum estimates in Equation (18), the ML prefilter is expressed as where is the estimated coherence function, given by

Note that . When a GCC technique is used to estimate the time difference, the error variance has been derived using the integration of a weighted CSD function in the frequency domain [18]. In the discrete-frequency domain, the error variance for an arbitrary prefilter is obtained through the summation of a weighted discrete CSD function as follows:

To derive the error variance of the ML prefilter using the estimated PSDs and CSD, we substitute Equation (19) into Equation (21). In addition, by approximating that , we get where , and is the estimate of given by

It can be easily shown that is minimised when for all k. Practically, is highly dependent on the estimation error of . Especially, when is close to one, is very sensitive to the variation of , thereby resulting in a rapid increase of .

In an attempt to mitigate the influence of coherence estimation error, we propose a new modified ML prefilter that employs a regularisation factor when computing , and thus it is more robust to the coherence estimation error than the conventional ML prefilter. The proposed modified ML prefilter is expressed as where is defined as

Here, is the regularisation factor. Since the denominator of Equation (25) is equal to or greater than , the estimation error included in gives less impact on the change of compared to the conventional ML prefilter. The proposed modified ML prefilter (24) is equal to the ML prefilter when , while it is equal to the Wiener prefilter when . By replacing in Equation (22) with in Equation (25), the error variance for the modified ML method is denoted as

As the regularization factor increases, becomes more robust to the estimation error of ; however, the bias of increases. In addition, the optimal minimising is varied according to SNR values of receive signals and statistical characteristics of . Therefore, we carefully design considering the variance of time difference error in numerical simulations and field measurements. The specific design method for is described in Section 4 and Section 5.

4. Simulation Results

4.1. Channel Model and Optimisation of Regularisation Factor

To perform computer simulations, we construct a channel model for vibration signals in buried water pipes based on the field measurement data. Details of field measurements are explained in Section 5. We define a channel model for propagation of vibration signals in a water pipe as below: where m is the reference distance; is the distance between the leak point and the sensor i; m is the path loss exponent; and K are the coefficient and order of a comb filter; is the coefficient of a first-order infinite impulse response (IIR) filter; and is an error term describing the channel variation. It is assumed that, for various , is an independent and identically distributed (i.i.d.) Gaussian random variable with zero mean and variance and that and are independent. In Equation (27), is the inverse comb filter to represent a sharp increase in the frequency response at harmonics of the natural frequency, and the first-order IIR filter describes the lowpass characteristics of vibration signals.

Figure 3 compares the coherence function of the signals generated by the channel model (27) and that of the measured signals in the field test, when and . For the channel model, we used m, m, , , , , and . In the field measurement, the coherence is obtained by averaging three trials of measurement data for Location 3 and sensors 1 and 3 (L3S13). The two coherence values are very similar except that the normalised frequency is around , and this shows the validity of the proposed channel model.

Figure 3

Coherence when , and : (a) signals generated by the channel model (27); (b) measured signals in the field test for Location 3 and sensors 1 and 3 (L3S13).

Using the channel model (27), we find the optimal regularisation factor for the proposed modified ML prefilter. In the simulation, the coherence is generated using the same parameters as in Figure 3, and the coherence estimate is denoted as where is the estimation error of that is assumed to be a zero-mean Gaussian random variable with variance . Figure 4 presents the normalised variance of time difference estimation error, which is defined as where is computed by Equations (25) and (26). The normalised error variance monotonically decreases as the regularisation factor increases to , whereas it grows with increment of when . Thus, the normalised error variance is minimised when , and the optimal is slightly different depending on . As decreases, the minimum of normalised error variance decreases, i.e., the performance gain of the modified ML prefilter increases compared to the conventional ML prefilter with .

Figure 4

Normalised variance of time difference error as a function of , when , and .

As mentioned in [12,13], may have a non-Gaussian distribution with a bias in a practical system. However, to the authors’ best knowledge, the distribution of in practical water pipes are not available in the literature, and thus is assumed to be zero-mean Gaussian in the simulation. Since is a function of , the normalised error variance curve in Figure 4 is varied according to the distribution of and the optimal also can be changed. For more practical design, we compute the normalised error variance curve and the corresponding optimal using measured at real test sites in Section 5.2.

4.2. Performance Comparison

In this subsection, we compare the performance of the proposed modified ML method with those of existing prefilters. As mentioned in [18,21], the PHAT prefilter performs better than the Roth and SCOT prefilters with the peak spreading problem. In addition, as explained in Section 3, the Eckart prefilter requires the noise PSDs and , which are difficult to estimate in a practical system. For these reasons, we consider the PHAT, Wiener, and ML prefilters for comparison with the proposed method. To compare the performance of various prefilters, we utilize the RMSD of the time difference estimate given by

Note that when is unbiased. In addition, we employ the PAR to describe the relative magnitude of the maximum peak compared to the average of the cross-correlation function. Specifically, when the maximum peak appears at , the PAR is defined as where the cross-correlation function is given by Equation (15).

Figure 5 and Figure 6 show the RMSD and the PAR for various prefilters, respectively, when the channel model (27) is used, , and . was set to for the proposed modified ML prefilter; was determined by the SNR at sensor 1; d was randomly generated in the range between and 200; and other parameters were the same as in Figure 4. All curves were obtained by averaging the simulation results over channel realizations. The RMSD decreases with increment of SNR, and the proposed method outperforms the conventional PHAT, ML, and Wiener prefilters in the whole SNR region. The RMSD value of the PHAT prefilter approaches those of the ML and Wiener prefilters as the SNR increases, because the PHAT prefilter has good performance in high SNR region [18,21]. In Figure 6, the proposed modified ML method presents the highest PAR among all prefilters for entire SNR region, i.e., the proposed prefilter can detect the leakage point more reliably than other techniques. The PAR of the ML prefilter is saturated when SNR dB because the estimation error of influences both the numerator and denominator of Equation (23).

Figure 5

Root mean square deviation (RMSD) of time difference error for various prefilters, when , , and .

Figure 6

Peak-to-average ratio (PAR) of the time domain cross-correlation function for various prefilters, when , , and .

5. Field Measurement

5.1. System Setup for Field Testing

Figure 7 shows the leak detection system for field measurements. We developed a web-based leakage detection system including our own sensor composed of a vibration transducer, a microprocessor, and a communication module. Once the control node sends the start logging command to the sensors via the network node, the sensors record the vibration signals for 30 s and transmit the sampled data to the server through the wireless communication link and the commercial internet protocol (IP) network. Prior to data logging, the network node broadcasts a sync message, and the local clock of the sensor is synchronized using the message to avoid the time difference estimation error caused by the time offset between sensors. The measured data stored in the server can be extracted for leakage detection at the control node. At the sensor, the sampling frequency was 4 kHz ( ms) and each sample is represented as 12 bits. The wireless link between the sensor and the network node supports 40 kbps of data rate using binary phase shift-keying (BPSK) modulation.

Figure 7

Leak detection system for field measurements.

Field measurements were carried out in Seogwipo City, Korea, in order to verify the leak detection performance of the proposed method. To determine specific test areas, we utilized the leak rate information in [3] and the pipeline database provided by the local government. We have selected five locations to perform field measurements as shown in Figure 8, and specific pipe layouts in these locations are given by Figure 9. The leaks in water pipes occur mostly at the joint to connect pipes and the branch to distribute water to small pipes, and these leak points can be identified by checking manholes and the entrance of household pipes. Thus, we have found the actual leak point described as red spots in Figure 9 through the inspection of possible leak points as well as the leak detection based on the time difference estimation. The sound velocity in Figure 9 was determined by querying the database provided by [35] using the information about the material and diameter of the water supply pipe. Since two adjacent sensors are selected to log vibration signals for each field trial, there are nine kinds of test cases as shown in Table 3. For notational convenience, we use the test case names in Table 3 in Section 4 and Section 5.

Figure 8

Locations for field measurements.

Figure 9

Pipe layout used in field measurements where S1, S2, and S3 means Sensors 1, 2, and 3; and red spots denote leak points: (a) Location 1; (b) Location 2; (c) Location 3; (d) Location 4; (e) Location 5.

Table 3

Test cases in the field measurement.

Test Case	L1S13	L1S23	L2S12	L3S13	L3S23	L4S12	L4S23	L5S12	L5S23
Location	1	1	2	3	3	4	4	5	5
Sensor Pair	S1, S3	S2, S3	S1, S2	S1, S3	S2, S3	S1, S2	S2, S3	S1, S2	S2, S3

5.2. Performance Evaluation

In a practical system, the real coherence is not available. Instead, is estimated using the received signals and through multiple trials of measurement. Using measured at five test locations, we compare the normalised variance of time difference error as a function of as shown in Figure 10. The coherence function was obtained by averaging three separate measurement data for each test case. We used and , and was set to . For varying , the normalised error variance presents huge variations depending on a test place because of the change of and some interference by the external environment; however, the normalised error variance changes slightly when having the minimum value between 0.35 and 0.60 as shown in Table 4. Based on the analysis in Figure 10, was set to for the proposed prefilter.

Figure 10

Normalised variance of time difference error as a function of , when measured coherence values were used, , , and .

Table 4

Optimal regularisation factor minimising the time difference error variance in field measurements.

Test Case	L1S23	L2S12	L3S23	L4S12	L5S12
Optimal α	0.35	0.40	0.60	0.35	0.45

As mentioned in Section 5.1, the field measurement was carried out for nine test cases in five locations. The field measurement was performed three times for each test case. We used and for estimating the PSDs and , and the CSD ; and was set to 0.45 according to the analysis in Section 4.1. Based on the measurement data, Figure 11 and Figure 12 present the RMSD and the PAR of various prefilters, respectively. The proposed modified ML prefilter achieves the smallest RMSD value in all test cases except L1S13 in which the PHAT prefilter has slightly less RMSD than the proposed method. Note that, in the test case L1S13, the RMSD values of all prefilters are less than ms, i.e., the proposed method achieves proper estimation accuracy for finding the leak location. The RMSD of the proposed method is less than ms for entire test cases. In Figure 12, the proposed prefilter obtains the highest PAR in all test cases except L2S12 in which the PAR value of the PHAT prefilter is dB greater than that of the proposed scheme. The PHAT and Wiener prefilters have similar PAR values; however, the RMSD of the PHAT prefilter is higher than that of the Wiener method as shown in Figure 11. In addition, while the ML and Wiener prefilters have similar RMSD values except L5S23, the ML method presents much lower PAR values than the Wiener prefilter. Thus, the Wiener prefilter is better than the PHAT and ML methods. Overall, the proposed method outperforms the conventional prefilters because it achieves the minimum RMSD value among all prefilters while maximising the PAR. In other words, when the distance between sensors are identical, the proposed prefilter has higher leak detection probability than the conventional methods. In addition, under a fixed RMSD requirement, the distance between sensors can be longer in the proposed prefilter than those in the conventional schemes, and thus the sensor deployment cost can be reduced.

Figure 11

RMSD of various prefilters in field measurements, when , , and .

Figure 12

PAR of various prefilters in field measurements, when , , and .

6. Conclusions

This paper proposed a new time difference estimation method for leakage detection based on the modified ML prefilter with the regularisation factor. Through analysis and numerical simulations, we have found the optimal regularisation factor and compared the performance of the proposed prefilter with that of conventional prefilters. The proposed leak detection method was verified by field measurements using a practical leak detection system, and it was shown that the proposed method outperforms the conventional techniques in terms of RMSD and PAR. The proposed modified ML prefilter can be applied to correlation-based leak detection systems utilizing the STFT and the wavelet transform. Furthermore, the proposed method can contribute to developing an automatic leakage management solution that collects leakage data, alarms about the risk of leaks, and informs about the specific leak locations.