- Email: [email protected]
Identification of multiple leaks in pipeline II: Iterative beamforming and leak number estimation
Mechanical Systems and Signal Processing 119 (2019) 346–362
Contents lists available at ScienceDirect
Mechanical Systems and Signal Processing journal homepage: www.elsevier.com/locate/ymssp
Identification of multiple leaks in pipeline II: Iterative beamforming and leak number estimation Xun Wang ⇑, Mohamed S. Ghidaoui Department of Civil and Environmental Engineering, Hong Kong University of Science and Technology, Clear Water Bay, Hong Kong, China
a r t i c l e
i n f o
Article history: Received 1 November 2017 Received in revised form 12 July 2018 Accepted 16 September 2018
Keywords: Leak detection Fluid transient Iterative beamforming EM algorithm Model selection
a b s t r a c t This paper considers the problem of detecting multiple leaks in a water-filled pipe using transient waves. A previous paper by the authors proposes a maximum likelihood (ML) method for leak detection based on an approximate linear model of wave propagation in pipe and shows that it is efficient, robust with noise, and super-resolved. However, that ML method needs to solve an N-parameter optimization problem (N is the leak number) which is high dimensional, non-convex and has many local maxima when N is large. The present paper deals with the multiple-leak detection problem where the leak number is high. A ML based iterative procedure, known as the iterative beamforming (IB) method, is proposed. Given an initial value of model parameter (locations and sizes of leaks) in each iteration, the complex inverse problem with 2N unknown parameters is decomposed into N onedimensional optimization problems, which largely decreases the computational complexity. More specifically, IB iteratively estimates the contributions from various leaks to the pressure head measurement and then performs a 1-D ML method on these estimates. In order to estimate the unknown number of leaks, model selection methods, known as Akaike information criterion (AIC) and Bayesian information criterion (BIC), are used and they are proven to be efficient to estimate the actual leak number. Numerical simulations show that five leaks can be accurately identified with super-resolution and the leak number can be correctly estimated using the methods proposed in this paper. Ó 2018 Elsevier Ltd. All rights reserved.
1. Introduction Leaks in water supply pipelines is a problem of increasing interest in the recent decades. Transient-based detection methods (TBDMs) utilize the hydraulics of the transient flows to detect leaks in the pipeline [1–18]. The reason that such methods are expected to work is that a leak in a pipeline system modifies the transient wave passing the leaks and measured by a sensor. Many TBDMs have been developed by researchers and applied to water piping systems including: (i) the transient reflection based method (TRM) [3,19–22]; (ii) the transient damping based Method (TDM) [6]; and (iii) the combination of both reflection and damping properties [2,4,5,7–10,15,16,23–32]. Recently, the problem of identifying multiple leaks has attracted increasing research attentions. Inverse transient analysis (ITA) method [2,4,29] is able to estimate multiple leaks, which assumes some discrete points in the pipe are potential leaks and estimates the corresponding leak sizes of these potential leaks by matching the measured pressure with time-domain numerical transient model. However, this approach involves a high-dimensional optimization problem (the dimension equals to the number of potential leaks), which implies a high computational cost and, more importantly, a high ⇑ Corresponding author. E-mail addresses: [email protected], [email protected] (X. Wang). https://doi.org/10.1016/j.ymssp.2018.09.020 0888-3270/Ó 2018 Elsevier Ltd. All rights reserved.
X. Wang, M.S. Ghidaoui / Mechanical Systems and Signal Processing 119 (2019) 346–362
347
Nomenclature q discharge oscillation h head oscillation xL leak location sL leak size steady-state discharge and head of leak Q L0 ; HL0 xm (m ¼ 1; ; M) sensor coordinate
Dh c n a A l D g Z R f DW
l x xth
kmin M N J T
r2
log L
head difference latent leak contribution measurement noise wave speed area of pipeline pipe length pipe diameter gravitational acceleration characteristic impedance frictional resistance Darcy–Weisbach friction factor propagation function angular frequency fundamental frequency minimum wavelength sensor number leak number frequency number sample size variance of noise log-likelihood function
Acronyms AIC Akaike information criterion BIC Bayesian information criterion EM expectation–maximization FRF frequency response function IB iterative beamforming MFP matched-field processing MFP(1) MFP with one-leak model ML maximum likelihood SNR signal-to-noise ratio
computation complexity such that the optimization problem often stops at local maxima corresponding to a wrong estimate of leak. For the same reason, ITA is not robust because the high-dimensional objective function is easy to be distorted by noise and uncertainties. Furthermore, the precision of ITA leak localization is limited to the mesh spacing of potential leaks in the pipe. A recent transient-based leak detection method, known as the matched-field processing (MFP) method [30], estimates leak location and size by matching the measurement with the frequency-domain model; it is found to be fast, accurate, robust with respect to noise and uncertainties, and can localize multiple leaks. However, MFP can only identify separate leaks provided the distance between leaks is not shorter than kmin =2, where kmin is the shortest probing wavelength. Based on the linear approximation model of wave propagation in a pipe with multiple leaks [31], a previous paper by the same authors has shown the good properties of the maximum likelihood (ML) approach for the multiple leaks [31]. That ML-based localization method has proven to be accurate, consistent, robust with respect to noise, and super-resolved. The latter property means that two leaks closer than kmin =2 can be separately estimated. In the case of one leak or two leaks, the localization can be easily accomplished by an exhaustive searching that plots the likelihood function with respect to possible leak locations. However, for a higher number of leaks, it is computationally expensive to compute the likelihood function for all the possible combinations of leak locations and it is complicated to directly solve the high dimensional optimization problem since the objective function is high-dimensional (with N parameters if N leaks are assumed) and non-convex. Therefore, the main goal of this paper is to develop an efficient method for leak detection problem where there are multiple leaks whose number, locations and sizes are all unknown. By the linearized model proposed in [31], each pressure head measured by a sensor in a pipe equals to the summation of contributions from various leaks and a term independent of all the leaks. It is clear that if these contributions could be
348
X. Wang, M.S. Ghidaoui / Mechanical Systems and Signal Processing 119 (2019) 346–362
Fig. 1. Setup of pipe system.
known, the estimation of leaks would be as straightforward as in the single leak case [30,33]. Based on this idea, an iterative method, known as the iterative beamforming (IB) [34–36], which deals with the parameter estimation of superimposed signals using the Expectation–Maximization (EM) algorithm [37–39], is proposed to solve the aforementioned problem. More specifically, IB alternates iteratively between two steps. The first one consists of computing the expected contributions from various leaks to the head measurements given the current fits of the model parameters. The location and size of each leak are then updated according to the leak contributions from the first step. The estimation of leak number is also an important issue in real applications. The ML approach cannot cope with this problem since it tends to choose the highest possible leak number which corresponds to highest possible likelihood function. In order to avoid model over-fitting, a trade-off between the goodness of fit (assessed by the likelihood function) and the model complexity (the number of leaks assumed in the model) is considered which intends to find an optimal estimate of leak number. Akaike information criterion (AIC) [40] and Bayesian information criterion (BIC) [41] are both the model selection methods for this issue and have been widely and successfully used in various domains [42–46]. In this paper, AIC and BIC are applied to decide the number of leaks. The organization of this paper is as follows. Section 2 begins with model description and problem statement. Then, estimation of multiple leaks using the IB approach is introduced in Section 3. Section 4 presents results of numerical simulation which justifies the proposed IB method is able to estimate multiple leaks with a super-resolution. Leak number determination via AIC and BIC is introduced in Section 5. Finally, conclusions are drawn in Section 6. 2. Model description and problem statement 2.1. Linearized model of wave propagation in pipeline The setup of the considered pipe system is shown in Fig. 1. A single pipeline bounded by an upstream and a downstream > node with coordinates x ¼ xU ¼ 0 and x ¼ xD is considered. The locations and sizes of N leaks, denoted as xL ¼ xL1 ; ; xLN > and sL ¼ sL1 ; ; sLN , are parameters to be estimated. Each leak size is represented by the lumped leak parameter sLn ¼ C d ALn , where C d is the discharge coefficient of the leak and ALn is the flow area of the leak orifice. The steady-state disrffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi charge of a leak is related to the lumped leak parameter by Q L0n ¼ sLn 2g HL0n zLn , in which g is the gravitational acceleration, zLn denotes the elevation of the pipe at each leak, and Q L0n and HL0n are the steady-state discharge and head1 at each leak. The leaks are estimated by pressure head measurement at M sensors. The corresponding head difference due to leaks is used, whose theoretical expression can be approximated via a linearized model [31]: N X Dh G x L s L þ n ¼ Gn xLn sLn þ n:
ð1Þ
n¼1
The derivation of Eq. (1) is detailed in Appendix A. In Eq. (1),
> Dh ¼ Dh11 ; ; DhJ1 ; ; Dh1M ; ; DhJM ; M
ð2Þ
NL
M
in which Dhjm ¼ hjm hjm (j ¼ 1; ; J and m ¼ 1; ; M) denotes the head difference between the head measurement hjm in the presence of leaks and theoretical head that does not include the leak terms NL
NL
hjm ¼ h
xj ; xm ¼ Z sinh ðlxm Þq xU þ cosh ðlxm Þh xU ;
ð3Þ
at the j-th frequency xj and at the m-th measurement station xm . Physically, Dhjm stands for the scattering of the wave at the pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi leaks. In Eq. (3), Z ¼ la2 =ð ixgAÞ is the characteristic impedance, l ¼ a1 x2 þ igAxR is the propagation function, a is the wave speed, x is the angular frequency, A is the area of pipeline, and R is the frictional resistance term. Note that 1 Discharge and head are terms in fluid dynamics and hydraulics. The discharge q (in m3 s1 ) is the volumetric flow rate of fluid and it relates to the flow velocity v (in ms1) by q ¼ Av , where A is the cross-sectional area of the pipe. The pressure head is defined as h ¼ p=ð qg Þ, where p is the pressure (in Pascal), g is the gravitational acceleration and q is the density of the fluid. For example, the atmospheric pressure is approximately patm ¼ 101 kPa. In terms of head, h ¼ 1:01 105 Pa= 103 kg m3 9:8m s2 ¼ 1:03m. This states that the atmospheric pressure is equivalent to pressure produced by 1.03 m of water. Alternatively, a 50 m of head in pipe implies that if that pipe bursts, the height of the water jet would be about 50 m.
X. Wang, M.S. Ghidaoui / Mechanical Systems and Signal Processing 119 (2019) 346–362
349
R ¼ ðf DW Q 0 Þ= gDA2 for turbulent flows where f DW is the Darcy–Weisbach friction factor, Q 0 is the steady-state discharge in the pipe and D is the pipe diameter. The matrix G in Eq. (1) is a JM N-dimensional matrix whose n-th column is
> Gn xLn ¼ G x1 ; xLn ; x1 ; ; G xJ ; xLn ; x1 ; ; G x1 ; xLn ; xM ; ; G xJ ; xLn ; xM ;
ð4Þ
in which
pffiffiffi g Z sinh l xm xLn Ln rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Z sinh lxLn q xU cosh lxLn h xU : G xj ; x ; xm ¼ 2 HL0n zLn
ð5Þ
The random noise vector
> n ¼ n11 ; ; nJ1 ; ; n1M ; ; nJM
ð6Þ
follows independent and identically distributed (i.i.d.) Gaussian distribution with 0-mean and covariance matrix r2 I, where I is an identity matrix. Note that the discharge at upstream q xU in Eqs. (3) and (5) can be estimated [31] by
U coshðlÞh xU hðx0 Þ ^ q x ¼ ; Z sinhðlðx0 xU ÞÞ
ð7Þ
if the pressure at x0 (x0 > xU ) is measured, boundary condition h xU is applied, and no leak is between xU and x0 . The derivation of Eq. (7) is detailed in Appendix B. It must be emphasized that Eq. (1) involves a number of assumptions. First, the nonlinear convective terms are neglected because they are of the order of Mach number and this Mach number is very small (103) in practice [47,48]. Second, the friction, valves, pumps, etc., in the system are given by nonlinear relations. The derivation of Eq. (1) requires that these relations are linearized. Such linearization is justified by the fact that injected transient waves for defect detection need to be small so as not to compromise the structural integrity of the system. Third, it is shown in [31] that leak detection is not affected by leak to leak interactions for small leaks. This is not a perious limitation since the motive is to detect small leaks. However, the leak to leak interaction can be accounted for but doing so makes the leak localization computationally demanding. 2.2. Problem statement Based on the linearized model Eq. (1), a ML solution for estimating multiple leaks is proposed in [31], which is briefly recalled here. Since the vector of head differences Dh follows i.i.d. Gaussian distribution, its log-likelihood function, which is a function of the leak locations xL and the leak sizes sL , has the form (after removing unnecessary terms which are independent of parameters to be estimated):
2 log L xL ; sL jDh ¼ Dh G xL sL :
ð8Þ
^L and sizes ^sL of the The MLE of xL and sL can be obtained by maximizing Eq. (8), which leads to the estimates of locations x leaks:
1 H L ^L ¼ arg max DhH G xL GH xL G xL x G x Dh xL
ð9Þ
and
1 ^sL ¼ GH x ^L G x ^L Dh: ^L GH x
ð10Þ
In the numerical examples in [31], two leaks are localized by plotting the 2D function in Eq. (9) and finding its maximum. However, in the case that the leak number N is large, this approach (exhaustive search of locations of leaks) is computationally intensive and inefficient. For example, considering a pipe with length l ¼ 2000 m and localization precision (grid density) being 1 m (as the numerical examples in [31]), the function in Eq. (9) needs to be computed 2000N times if N leaks are assumed, which implies a very high computational cost when N is large. Another approach is to directly solve the optimization problem (9). However, in the case of large leak number, this optimization problem is still difficult to solve due to the complex nature of G xL : the objective function in Eq. (9) is non-convex and has many local maxima in a highdimensional space [33]. The present paper addresses this shortcoming. More specifically, this paper proposes an indirect but simpler (in the sense of the complexity of optimization problem) method, which only needs to solve 1D optimization problems.
350
X. Wang, M.S. Ghidaoui / Mechanical Systems and Signal Processing 119 (2019) 346–362
3. Estimation of multiple leaks using the iterative beamforming method This section introduces the IB method for identifying multiple leaks. For this purpose, latent leak contribution from various leaks to pressure head measurement is introduced first. 3.1. Latent leak contributions By introducing a latent variable,2 called latent leak contribution, the optimization problem (9) can be largely simplified. The latent leak contribution is defined as the contribution from a single leak to the head difference (the data for leak detection). More specifically, the latent leak contribution from the n-th leak, denoted as cn , is
cn ¼ Gn xLn sLn þ nn :
ð11Þ
Here, the noise term nn can be obtained by arbitrarily decomposing the total noise n in Eq. (1) into N components such that their summation equals in distribution to n, i.e., N X nn ¼ n:
ð12Þ
n¼1
Therefore, the head difference Dh is related to the latent leak contributions cn by
Dh ¼
N X
cn :
ð13Þ
n¼1
In this equation, ‘‘=” means equality in distribution of random variables. By assuming that nn (n ¼ 1; ; N) are i.i.d. Gaussian random vectors with 0-mean and covariance matrix Rn ¼ r2 =N I which satisfies Eq. (12), the vector of latent leak contributions c ¼ ðc1 ; ; cN Þ> also follows a Gaussian distribution and the log-likelihood function of c is (again, removing unnecessary terms independent of model parameters): N X 2 log L xL ; sL jc ¼ cn Gn xLn sLn :
ð14Þ
n¼1
Note that each term in Eq. (14) corresponds to a single leak and is independent of other leaks. Therefore, if the contribution cn of the n-th leak is available, the corresponding leak location xLn and size sLn can be estimated by solving a single leak localization problem:
Ln Ln 2 ^x ; ^s ¼ arg min cn Gn xLn sLn ; xLn ;sLn
ð15Þ
as in [30]. 3.2. Iterative beamforming for detecting multiple leaks Unfortunately, the latent leak contributions cn cannot be directly observed. However, the IB method [34,35] is able to proceed with Eq. (14) by treating the latent contributions as random variables; this approach is based on a ML method for dealing with latent variables or missing data, known as the EM algorithm [37–39]. Here, this approach is applied to the multipleleak detection problem. The proposed method is able to maximize the observed data log-likelihood (8) by computing iteratively the expectation of the log-likelihood (14) of the latent leak contribution c, given the observed data and the current fit of the parameter. More
precisely, given an initial parameter U0 ¼ xL0 ; sL0 , the algorithm iterates back and forth between two steps until convergence: Expectation step (E-step): computing the conditional expectation
Q ðUjUk1 Þ ¼ E log L xL ; sL jc jDh; Uk1
ð16Þ
of the complete data log-likelihood with respect to the unknown variables, given the observed data Dh and a current fit
Uk1 for the model parameters;
Maximization step (M-step): determining the new fit for the parameters by maximizing Q ðUjUk1 Þ with respect to U. The algorithm guarantees that the observed data likelihood increases after each iteration [37]; whether it converges to the global maximum depends on the choice of the initial parameter U0 [38]. A strategy to deal with this problem is to 2
A latent variable is an unknown (or unobserved) information that would make the estimation process straightforward, should it be available.
X. Wang, M.S. Ghidaoui / Mechanical Systems and Signal Processing 119 (2019) 346–362
351
run the algorithm with different initial parameter values that cover a large part of the parameter space, to compute the corresponding estimates, and to finally retain the result associated with the highest likelihood function. In the following, the Eand M-steps are detailed, respectively. 3.2.1. E-step: Estimation of contributions from various leaks The E-step consists in computing N N X X 2 Eðcn jDh; Uk1 Þ Gn xLn sLn 2 ; Q ðUjUk1 Þ ¼E log L xL ; sL jc jDh; Uk1 ¼ E cn Gn xLn sLn jDh; Uk1 ¼ c n¼1
n¼1
ð17Þ where c is a constant independent of the parameters. Therefore, the E-step is to compute the expected leak contributions from various leaks given the current fit of model parameters, i.e., Eðcn jDh; Uk1 Þ. >
Note that both Dh and cn are Gaussian distributed, thus the vector ðDh; cn Þ is jointly Gaussian with 0-mean and covariance matrix
R¼
CovðDhÞ
CovðDh; cn Þ
Covðcn ; DhÞ
Covðcn Þ
¼
I 1 I N
1 I N 1 I N
! :
ð18Þ
Therefore, the conditional distribution of cn given Dh is Gaussian:
N1 1 cn Dh N cn Gn xLn sLn þ Dh G x L s L ; I : N N2
ð19Þ
Eq. (19) is derived from a property of Gaussian random variable [49], which is detailed in Appendix C. Furthermore,
n Ln 1 ^cnk ¼ Eðcn jDh; Uk1 Þ ¼ Gn xLk1 Dh G xLk1 sLk1 : sk1 þ N
ð20Þ
Finally, the conditional expectation in the E-step is obtained:
Q ðUjUk1 Þ ¼ c
N X ^cnk Gn xLn sLn 2 :
ð21Þ
n¼1
3.2.2. M-step: Single leak detection problem The M-step amounts to computing new estimates of xLn and sLn by maximizing Eq. (21), which is decomposed into N twodimensional minimization problems:
Ln Ln xk ; sk ¼ arg min k^cnk Gn xLn sLn k2 ; n ¼ 1; ; N: xLn ;sLn
ð22Þ
Eq. (22) can be further simplified as follows. For any leak location estimate xLkn , the corresponding leak size is
sLkn ¼ arg mink^cnk Gn xLkn sLn k2 ¼ sLn
GHn xLkn ^cnk : GHn xLkn Gn xLkn
ð23Þ
Then, inserting Eq. (24) into Eq. (22) results in the leak location estimate:
xLkn
2 GHn xLn Gn xLn ^ ^cnk ¼ arg mincnk H xLn Gn ðxLn ÞGn ðxLn Þ L H Ln G x Gn x n ^cnk : ¼ arg maxc^Hnk nH xLn Gn ðxLn ÞGn ðxLn Þ
ð24Þ
It is clear that the leak locations and sizes computed at each iteration amount to performing a Bartlett’s beamforming estimation (equivalently, MFP, ML or matched-filter [30,33]) given each latent leak contribution cn ; each leak location is determined via the single-leak detection method (i.e., a single-parameter optimization problem is solved) by using the corresponding estimated contributions. Finally, the IB strategy for estimating the model parameters is summarized in Algorithm 1.
352
X. Wang, M.S. Ghidaoui / Mechanical Systems and Signal Processing 119 (2019) 346–362
Algorithm 1 Identification of multiple leaks using IB > > For k ¼ 0, pick starting values xL0 ¼ xL01 ; ; xL0N and sL0 ¼ sL01 ; ; sL0N for the model parameters. For k P 1: repeat estimate the latent leak contributions ^ cnk , n ¼ 1; ; N, by Eq. (20); update the leak location xLkn , n ¼ 1; ; N, by Eq. (24); update the leak size sLkn , n ¼ 1; ; N, by Eq. (23); until the relative increase of the observed data log-likelihood (8) is less than a given threshold
j:
log LxL ; sL jDh log LxL ; sL jDh k k k1 k1 < j: log L xLk1 ; sLk1 jDh
ð25Þ
4. Numerical simulations In this section, numerical simulations are introduced to illustrate the performance of the IB algorithm. First, an example shown in [33] that cannot be well-solved by the 1D searching methods in [30,33] is studied. Then, a more complicated case with five leaks is considered. 4.1. Revisiting the example in [33]: Identifying two leaks Here, a numerical example in [33] is revisited. A single pipe which is connected by two reservoirs in a horizontal plane is considered. The pressure heads of both reservoirs are HU ¼ 25 m and HD ¼ 20 m. The pipe length is l ¼ 2 km and the diameter is D ¼ 0:5 m. The wave speed is 1200 m/s. The Darcy–Weisbach friction factor f DW ¼ 0:02. Two leaks with locations xL1 ¼ 600 m and xL2 ¼ 1600 m and sizes sL1 ¼ 1 104 m2 and xL2 ¼ 1:2 104 m2 are assumed. Three sensors are located at 50 m, 1800 m and 1960 m: the first sensor is for estimating the boundary condition q xU via Eq. (7) and the pressure heads from the latter two are then used to estimate leak locations and sizes. The resonant and anti-resonant frequencies fx : ax; a ¼ 1; 2; ; 31g are used for leak detection. The signal-to-noise ratio (SNR), which is defined by
SNR ¼ 20 log10 jDhj=r ;
ð26Þ
is 15 dB, jDhj represents the average head difference over frequencies and r stands for the standard deviation of Gaussian white noise. Fig. 2(a) shows the output function of MFP [30,33] with single leak assumption, abbreviated as MFP(1). It is clear that near each actual leak, the output function of MFP(1) has a local maximum, but the other two local maxima prevent obtaining correct estimates of leak locations. Then, the leak detection methods, i.e., directly solving the optimization problem Eq. (9) and IB, based on the model of multiple leaks in Eq. (1) are used. The former is abbreviated as ML in this section. In both methods, four leaks are assumed (i.e., N ¼ 4 in the model) according to the result of MFP: the initial values of leak locations are taken from the four local maxima in Fig. 2(a) (377 m, 576 m, 1391 m and 1604 m) and the initial leak sizes are 2 104 m2 for all the four leaks (it is arbitrarily chosen). The optimization method used in both approaches is the Broyden-Fletcher-Goldfarb-Shanno (BFGS) QuasiNewton algorithm implemented by the MATLAB function ‘‘fminunc.m”. The results of both leak detection methods are shown in Fig. 2(b) and (c), wherein the locations and sizes of the four estimated leaks (circles), along with the actual leaks (represented by the solid lines), are plotted. It is clear that both methods can accurately localize the leaks and return similar results. Taking the IB result as an example, two estimates are very close to the actual leaks, with locations ^ xL1 ¼ 600:1 m and 4 4 L L 2 L 2 ^ x 2 ¼ 1600 m and sizes ^s 1 ¼ 0:995 10 m and ^s 2 ¼ 1:197 10 m . The other two estimates have very small values of ^L4 ¼ 1238 m) and thus can be ignored. Here, the leak size (^sL3 ¼ 4:3 107 m2 at ^ xL3 ¼ 423:8 m and ^sL4 ¼ 5:2 107 m2 at x IB algorithm converges at 81st iteration with threshold of relative increase of log-likelihood
j ¼ 104 .
4.2. A five-leak example In this section, a more complicated case with five leaks is considered. The leak locations are xL1 ¼ 600 m, xL2 ¼ 640 m, L3
x
¼ 800 m, xL4 ¼ 1540 m, xL5 ¼ 1600 m; the leak sizes are sL1 ¼ 1 104 m2 ; sL2 ¼ 1:2 104 m2 ; sL3 ¼ 0:9 104
m ; sL4 ¼ 1:4 104 m2 ; sL5 ¼ 1 104 m2. As previously, the resonant and anti-resonant frequencies fx : ax; a ¼ 1; 2; ; 31g are used for leak detection. Note that with this probing wave whose maximum angular frequency 2
X. Wang, M.S. Ghidaoui / Mechanical Systems and Signal Processing 119 (2019) 346–362
353
Fig. 2. Leak detection using MFP(1) (a), ML (b), and IB (c). In subfigure (a), the solid line, dash lines and crosses stand for the objective function of MFP(1), actual leak locations and sensor locations. In subfigures (b) and (c), the solid lines and circles represent the actual and estimated leaks, respectively. The leak locations are xL1 ¼ 600 m and xL2 ¼ 1600 m and the leak sizes are sL1 ¼ 1 104 m2 and sL2 ¼ 1:2 104 m2.
is 31xth , the resolution (the minimum range of two leaks that can be separately estimated) is approximately kmin =2 ¼ 129 m according to the Nyquist-Shannon sampling theorem and the observations from the simulation results in [30,33]. However, the distance between xL1 and xL2 (40 m, 0:15kmin ) and the distance between xL4 and xL5 (60 m, 0:23kmin ) are smaller than the minimum resolvable leak distance kmin =2. In the following, IB is justified to be able to separately estimate the five leaks, which implies that IB is a super-resolution scheme. Fig. 3(a) shows the localization results using MFP(1). Two local maxima, which respectively locate between xL1 and xL2 and between xL4 and xL5 , can be found. However, other local maxima disturb the localization. Then, ML and IB are used to detect the leaks. The results in Fig. 3(a) are used for the parameter initialization. The seven highest local maxima are all considered as possible locations of leaks and thus are taken as initial points. Since the leak number is unknown, a large number of leaks is assumed (the determination of leak number is discussed in Section 5). More specifically, two initial leaks are attributed to each of the first five highest local maxima; 435 5 m, 614 5 m, 1200 5 m, 1353 5 m, 1571 5 m are taken as the initial leak locations. Another two local maxima, being 742 m and 1515 m, are also assumed as two initial leak locations, thus 12 leaks are assumed in the IB algorithm. Furthermore, the initial leak sizes are assumed to be 2 104 m2 for all the 12 leaks (it is arbitrarily chosen since no prior information is available). Fig. 3(b) and (c) and Table 1 both display the actual leak
354
X. Wang, M.S. Ghidaoui / Mechanical Systems and Signal Processing 119 (2019) 346–362
Fig. 3. Leak detection using MFP(1) (a), ML (b) and IB (c). In subfigure (a), the solid line, dash lines and crosses stand for the objective function of MFP(1), actual leak locations and sensor locations. In subfigures (b) and (c), the solid lines, crosses and circles represent the actual, initial and estimated leaks, respectively. Each estimate of leak is connected to its corresponding initial parameter by a dash line. The leak locations are xL1 ¼ 600 m, xL2 ¼ 640 m, xL3 ¼ 800 m, xL4 ¼ 1540 m, xL5 ¼ 1600 m. The leak sizes are sL1 ¼ 1 104 m2 ; sL2 ¼ 1:2 104 m2 ; sL3 ¼ 0:9 104 m2 ; sL4 ¼ 1:4 104 m2 ; sL5 ¼ 1 104 m2.
locations and sizes, their initial parameters and the estimates obtained from ML and IB. Each initial parameter and the corresponding estimates are connected by a dash line in Fig. 3(b) and (c). It can be seen that ML accurately estimate all the five leaks. However, two ‘‘ghost leaks” with high estimated leak sizes appear at around 400 m and 1400 m. By contrast, the result obtained from IB is more correct that no misleading leak estimate appears. More specifically, all the leak locations are accurately estimated except xL4 ¼ 1540 m which is replaced by two estimates at 1533 m and 1560 m. The sizes of the first three leaks are accurately estimated while the fourth and fifth are less precise. In particular, the sum of the estimated sizes at 1533 m and 1560 m approximates the actual size of the fourth leak. It is also interesting to remark that although 12 leaks is assumed, half of them almost disappears in the IB results: their estimated sizes are very small and thus neglectable. The previous result is obtained with a single initialization which is acquired from the result of MFP(1). Alternatively, a multiple-initialization strategy can be used to further improve the leak detection result of IB. Algorithm 1 is run with various starting parameters and then the result with highest likelihood function Eq. (8) is retained as the final estimate. Instead of
355
X. Wang, M.S. Ghidaoui / Mechanical Systems and Signal Processing 119 (2019) 346–362 Table 1 The actual leak locations and sizes, their initial parameters and the estimates obtained from IB and ML. Actual leaks xL [m]
L
Initial parameters 4
s [10
600 640 800
1 1.2 0.9
1540
1.4
1600
1
2
m ]
xL [m] 420 430 609 619 742 1195 1205 1348 1358 1515 1570 1580
L
s [10
4
2 2 2 2 2 2 2 2 2 2 2 2
IB Estimates 2
m ]
xL [m] 254.8 401.8 603.6 637.6 801.8 1163 1242 1322 1414 1533 1560 1606
L
ML Estimates 4
s [10
2
m ]
0.001 0.001 0.971 1.117 0.876 0.025 0.045 0.031 0.031 0.844 0.771 0.794
xL [m]
sL [104 m2]
439.2 437.8 601.3 641.3 800.1 1137 1155 1375 1376 1451 1541 1601
1.298 1.301 1.064 1.139 0.900 0.013 0.052 1.140 1.115 0.040 1.425 0.980
using the initial leak locations as in the third column in Table 1, the 12 initial leak locations are all randomly generated, each of them is taken from a uniform distribution with a support whose center is each value in the third column in Table 1 and whose radius is kmin =2 ¼ 129 m. Each initial leak size is also randomly obtained from a uniform distribution with support ½0:5; 3 104 m2. The simulation (from initial parameter generation to leak estimation) is repeated 100 times and the result with highest likelihood is taken as the final estimate, which is shown in Fig. 4(a). The actual leak locations and sizes, their estimates obtained from IB (with highest likelihood among 100 random initial parameters), and the corresponding values of the initial parameters are shown in Table 2. The first, third and fifth leaks are respectively replicated by a single leak
Fig. 4. Leak detection using IB with random initialization. The solid lines, crosses and circles represent the actual, initial and estimated leaks, respectively. Each estimate of leak is connected to its corresponding initial parameter by a dash line. The leak locations are xL1 ¼ 600 m, xL2 ¼ 640 m, xL3 ¼ 800 m, xL4 ¼ 1540 m, xL5 ¼ 1600 m. The leak sizes are sL1 ¼ 1 104 m2 ; sL2 ¼ 1:2 104 m2 ; sL3 ¼ 0:9 104 m2 ; sL4 ¼ 1:4 104 m2 ; sL5 ¼ 1 104 m2. The loglikelihood function of the four cases are: (a) 590; (b) 622; (c) 630; and (d) 2524.
356
X. Wang, M.S. Ghidaoui / Mechanical Systems and Signal Processing 119 (2019) 346–362
Table 2 The actual leak locations and sizes, their estimates obtained from IB with highest likelihood among 100 results with different random initial parameters, and the corresponding values of the initial parameters. Actual leaks L
Estimates L
x [m]
s [10
4
2
m ]
L
x [m]
Initial parameters L
s [10
4
2
m ]
x [m]
sL [104 m2]
L
600
1
600.9
0.918
559.7
2.144
640
1
639.5 634.9 635.9
0.906 0.144 0.153
653.5 392.3 427.9
1.047 0.756 1.431
800
0.9
800.8
0.865
736.5
2.667
1540
1.4
1546 1530 1524
1.034 0.318 0.108
1561 1561 1224
1.659 2.461 1.484
1600
1
1601
0.935
1623
2.505
1217 1357 1551
0.025 0.048 0.046
1154 1340 1365
1.635 0.954 0.667
estimate. The other two leaks at 640 m and 1540 m are respectively represented by three estimates. The sum of the estimated leak sizes corresponding to each of the two leaks is approximately equal to the actual leak size. As previously, since the assumed leak number is much larger than the actual leak number, some estimates have very low size and thus are neglectable. Fig. 4(b–d) shows three other IB estimates among the 100 results. Since they have lower likelihood and, thus, finally are not used, but they are discussed here to illustrate the mechanism of deciding the ultimate estimate. The log-likelihood functions of the four results in Fig. 4 are (a) 1:46 104 , (b) 1:64 104 , (c) 1:68 104 , and (d) 11:87 104 . Correspondingly, the leaks are worse estimated as the likelihood function decreases. In the first three cases, each leak is clearly identified but the error of leak size estimate increases, being (a) 0.0488, (b) 0.0835, and (c) 0.0982, respectively. Here, the error is defined as
kn 5 X L X k L L n n ^s ; e ^s ¼ s n¼1 k¼1
ð27Þ
in which ^sLn ; k ¼ 1; ; kn , are the leak size estimates corresponding to the n-th leak. In the last case, as is shown in Fig. 4(d), the estimate is inaccurate. Note that this is an extreme case whose likelihood function is much lower than the maximum likelihood. Again, this estimate can be easily excluded by comparing the likelihood function. Fig. 5 plots the likelihood function Eq. (8) for IB and ML estimates with various SNR. For each SNR the simulation is repeated 100 times with different random initial parameters as previously; the highest likelihood is picked and shown in k
Fig. 5. Plot of likelihood function Eq. (8) for IB and ML estimates with various SNR. Each likelihood function is selected among 100 simulation results (corresponding to different random initializations) with highest likelihood function.
X. Wang, M.S. Ghidaoui / Mechanical Systems and Signal Processing 119 (2019) 346–362
357
Fig. 5. This figure illustrates that the IB estimate corresponds to higher likelihood function than ML implying that IB more easily avoids stopping at local maxima and performs better than ML (maximizing the likelihood function Eq. (8) is the objective of both ML and IB). Furthermore, Eq. (8) shows that the result of IB is able to better replicate the measurements using the estimated locations and sizes of leaks. Finally, the IB estimate with various leak number is assumed. The determination of leak number is discussed in Section 5; here, only the leak detection result with different assumed leak number is shown. In each simulation, the initial parameters are randomly generated as follows. It can be observed from Fig. 3(a) that the high local maxima of MFP(1) are all in the two intervals ½300; 900 and ½1100; 1700 , thus each initial value of leak location is generated from a uniform distribution with support ½300; 900 [ ½1100; 1700 . Each initial leak size is sampled from a uniform distribution with support ½0:5; 3 104 m2. Fig. 6 displays the IB estimates with respectively 2–7 assumed leaks, each of them is obtained by repeating IB with 100 random initialization and then picking the result with highest maximum likelihood. It is interesting that when the assumed leak number is smaller than the actual leak number (subfigures (a-c)), one estimate tends to represent two or three leaks close to each other. When the assumed and actual leak numbers are identical (subfigures (d)), the estimation is precise that each estimate replicates each leak. In the case that the assumed leak number is higher than the actual leak
Fig. 6. Leak detection using IB with random initialization assuming 2–7 leaks. The solid lines and circles represent the actual and estimated leaks, respectively. The actual leak locations are xL1 ¼ 600 m, xL2 ¼ 640 m, xL3 ¼ 800 m, xL4 ¼ 1540 m, xL5 ¼ 1600 m. The actual leak sizes are sL1 ¼ 1 104 m2 ; sL2 ¼ 1:2 104 m2 ; sL3 ¼ 0:9 104 m2 ; sL4 ¼ 1:4 104 m2 ; sL5 ¼ 1 104 m2.
358
X. Wang, M.S. Ghidaoui / Mechanical Systems and Signal Processing 119 (2019) 346–362
number (subfigures (e,f)), each leak is represented by more than one estimates and/or some of the leak estimates has very low estimated size which is thus ignored. 5. Model selection: Determination of leak number In the previous sections, the multiple leaks are estimated given an assumed leak number N. This section studies the determination of the leak number N. 5.1. AIC and BIC By considering the leak number as a parameter, the log-likelihood function Eq. (8) can be rewritten as
2 N X Ln Ln L L log L N; x ; s jDh ¼ Dh Gn x s : n¼1
ð28Þ
Note that the maximum likelihood principle tends to choose the highest possible leak number. Denote the ML estimates of L ^LðNÞ ¼ ^ x 1; ; ^ xLN and ^sLðNÞ ¼ ^sL1 ; ; ^sLN , then for a higher leak leak locations and sizes under the assumption of N leaks are x ^ LðNþ1Þ ¼ ^ number assumption N þ 1, the estimate x x L1 ; ; ^ xLN ; ^ xLN and ^sLðNþ1Þ ¼ ^sL1 ; ; ^sLN ; 0 remains at the same value of the log-likelihood function, which justifies that a higher leak number does not decrease the log-likelihood function Eq. (28). As a ^LðNþ1Þ and ^sLðNþ1Þ are the optimal estimates (although it is not unique) since the matter of fact, in an noise-free environment, x ^LðNþ1Þ ; ^sLðNþ1Þ jDh log-likelihood function Eq. (28) is 0 (the maximum value). However, when noises exist, because log L N þ 1; x ^ LðNþ1Þ and ^sLðNþ1Þ in terms of having a higher log-likelihood may be found due to is in general non-zero, an estimate better than x one more available free parameter. The goal of this section is to select a reasonable number of leaks which is a trade-off between the goodness of fit (assessed by the likelihood function) and the complexity of the model, i.e., to avoid over-fitting. The solution is to modify the objective function by adding a penalty for high leak number to the log-likelihood function. Here, a commonly used criterion in model selection, known as the Akaike information criterion (AIC) [40], is introduced:
Fig. 7. Log-lilelihood (left) and AIC (right) as a function of leak number N. The sample size is T ¼ 1.
X. Wang, M.S. Ghidaoui / Mechanical Systems and Signal Processing 119 (2019) 346–362
AICðNÞ ¼
J M X X
^L ; ^sL jDhjm ¼ 2NMJ log L N; x ^L ; ^sL jDh ; 2N log L N; x
359
ð29Þ
m¼1 j¼1
in which the first term is the penalty term. The optimal leak number can be decided by minimizing Eq. (29). For a large sample size T (realizations of Dh), another criterion, known as the Bayesian information criterion (BIC) [41], can be used:
BICðN Þ ¼
J M X X
^L ; ^sL jDhjm ¼ NMJ logðT Þ log L N; x ^L ; ^sL jDh : N logðT Þ log L N; x
ð30Þ
m¼1 j¼1
The leak number N estimated from minimizing Eq. (30) has been proven to be asymptotically optimal (in the sense of T ! 1) in [41] that it maximizes the posterior likelihood function (any prior distribution of model parameter is attributed). Finally, a strategy for estimating multiple leaks which combines IB and AIC/BIC, in the case of unknown locations, sizes and number of leaks, is summarized in Algorithm 2. Algorithm 2 Identification of multiple leaks using IB and AIC/BIC in the case of unknown leak number. For N ¼ 0, compute AICð0Þ via Eq. (29) (or, BICð0Þ via Eq. (30)). For N ¼ N þ 1: repeat ^ LðNÞ and ^sLðNÞ ; Run Algorithm 1 (IB) to obtain x Compute AICðN Þ via Eq. (29) (or, BICð0Þ via Eq. (30)). until AICðN Þ > AICðN 1Þ (or, BICðN Þ > BICðN 1Þ). ^LðN1Þ and ^sLðN1Þ as the estimates of leak locations and sizes. Retain N 1 as the estimate of leak number, x 5.2. Numerical examples Here, numerical examples are introduced to illustrate the model selection methods for estimating leak number. First, the numerical example of five leaks introduced in Section 4.2 is considered again. The left plot of Fig. 7 shows the corresponding log-likelihood function with various assumed leak numbers, which increases dramatically from N ¼ 2 to N ¼ 5 (the latter is the actual leak number). After that, the log-likelihood function still increases but very slightly. The right plot of Fig. 7 displays the corresponding plot for AIC whose penalty increases as the leak number increases, such that its the optimal leak number b ¼ 5 is obtained, which is equal to its actual value. estimate N In all the previous numerical examples, only one measurement is used for leak detection, i.e., T ¼ 1. Here, a numerical simulation with a sample size T ¼ 100 is considered; the locations and sizes of leaks, noise distribution and other setups
Fig. 8. Log-lilelihood (left), AIC (middle) and BIC (right) as a function of leak number N. The sample size is T ¼ 100.
360
X. Wang, M.S. Ghidaoui / Mechanical Systems and Signal Processing 119 (2019) 346–362
are same as the previous case. Note that in this case the leak detection is easier than the case T ¼ 1 as more measurements average the noise and thus the parameter estimation is more robust. For this reason, only leak number determination is discussed here. Fig. 8 shows the log-likelihood function with respect to various assumed leak number, along with the results of AIC and BIC. The log-likelihood function has a similar feature as the previous case. Both AIC and BIC can precisely estimate the leak number. 6. Conclusions This paper investigates the problem of detecting multiple leaks in a water-filled pipe with transient waves. A pressure head measurement is assumed to be the summation of leak contributions, each of them is a function of single leak and independent of other leaks. Based on this model, a maximum likelihood (ML) based iterative procedure, known as the iterative beamforming (IB) method, is proposed. IB iteratively estimates the contributions from various leaks to the pressure head measurement and then performs a 1-D inverse method on these estimates. Therefore, this approach largely decreases the computational complexity. Furthermore, the model selection methods for estimating optimal leak number are proposed: Akaike information criterion (AIC) and Bayesian information criterion (BIC) are introduced which penalize the high leak number and thus find a trade-off between the goodness of fit and the complexity of the model. The results obtained on simulated data show the interest of the proposed methods. IB solves a two-leak localization problem which is troublesome for the 1D searching methods like the matched-field processing. Another complex example with five leaks in a pipe shows that IB is efficient, accurate and is a super-resolution scheme. A multi-initialization strategy is also proposed to guarantee the global convergence of the algorithm. Finally, various leak number assumption is discussed; AIC and BIC are both proven to be efficient to decide the leak number in a complicated multiple-leak problem. Acknowledgements This work has been supported by research grants from the Research Grant Council of the Hong Kong SAR, China (Project No. T21-602/15R). The Chinese Estates are acknowledged for their support of this research. Appendix A. A brief introduction of the model in Eq. (1) The discharge and pressure head due to a fluid transient are represented by q and h. Given the discharge q xU and U at the upstream node xU , the quantities at xm (xm > xLN > > xL1 ) are obtained from the transfer matrix head h x method [50]:
q ð xm Þ
¼M
hðxm Þ
NL
xm xLN
3 ! 1 Q L0n q xU 1 Ln 4@ 2ðH0 zLn Þ AM NL xLn xLn1 5 ; h xU n¼N 0 1
1 Y
20
in which xL0 ¼ xU ,
cosh ðlxÞ
NL
M ð xÞ ¼
Z sinh ðlxÞ
1Z sinh ðlxÞ
ðA:1Þ
! ðA:2Þ
cosh ðlxÞ
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi is the field matrix, Z ¼ la2 =ð ixgAÞ is the characteristic impedance, l ¼ a1 x2 þ igAxR is the propagation function, a is the wave speed, x is the angular frequency, A is the area of pipeline, and R is the frictional resistance term. The right hand side of Eq. (A.1) can be approximated by a linear form plus a infinitesimal term (when leak sizes tend to 0) that can be neglected in the case of small leaks [31]:
q ð xm Þ hðxm Þ
¼
MNL ðxm Þ þ
N X sLn M SL Ln ; xLn ; xm n¼1
!
! q xU U þ o max sLn ; n¼1;;N h x
ðA:3Þ
in which
! vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Z sinh lxLn cosh l xm xLn cosh lxLn cosh l xm xLn u g M SL Ln ; xLn ; xm ¼ u t Z cosh lxLn sinh l xm xLn Z 2 sinh lxLn sinh l xm xLn 2 HL0n zLn
ðA:4Þ
is a matrix related to the location xLn of the n-th leak (zLn and HL0n are decided by xLn ), but is independent of the size sLn of n-th leak and the other leaks. Note that, physically, the infinitesimal term in Eq. (A.3) implies the wave interaction (multiple reflections) between leaks, while the first term in Eq. (A.3) corresponds to the wave received by xm after a single reflection at each leak. By neglecting the infinitesimal term in Eq. (A.3), the pressure head at xm (m ¼ 1; ; M) near the downstream for the frequency xj (j ¼ 1; ; J) is
X. Wang, M.S. Ghidaoui / Mechanical Systems and Signal Processing 119 (2019) 346–362 N X NL h x j ; xm ¼ h x j ; xm þ sLn G xj ; xLn ; xm ;
361
ðA:5Þ
n¼1 NL
xj ; xm and G xj ; xLn ; xm follow Eqs. (3) and (5). NL Let Dhjm h xj ; xm h xj ; xm denote the head difference between the head measurement in the presence of leaks and the theoretical head that does not include the leak terms at the frequency xj and at the measurement station xm , and denote wherein h
> Dh ¼ Dh11 ; ; DhJ1 ; ; Dh1M ; ; DhJM ;
ðA:6Þ
then we have the following equation:
Dh ¼ G xL sL ;
ðA:7Þ
in which G is a JM N-dimensional matrix (cf. Eqs. (4) and (5)). Finally, Eq. (1) is obtained by adding a noise term. Appendix B. Upstream discharge estimation in Eq. (7) Here, the discharge q xU at the upstream node xU is estimated by adding a measurement station near the upstream boundary at xU þ ( > 0) [51,52]. Assuming that there is no leak between xU and xU þ , the discharge and head at xU and xU þ have the expression:
! ! q xU þ q xU NL ¼M ðÞ U h xU h x þ ¼
coshðlÞ
Z sinhðlÞ
1 Z
sinhðlÞ
coshðlÞ
!
! q xU þ : h xU þ
ðB:1Þ
By measuring pressure head h xU þ at xU þ ( > 0) and applying a boundary condition h xU of head at xU , the discharge U q x can be estimated by [31]:
coshðlÞh xU h xU þ ^ xU ¼ q : Z sinhðlÞ Appendix C. Derivation of Eq. (19) The following property of Gaussian random variable [49] is first reminded, which is used to derive Eq. (19). Let X and Y be n-dimensional Gaussian random vectors with expectation mX and mY , and with covariance matrix RXX and RYY . Let RXY ¼ CovðX; YÞ and RYX ¼ CovðY; XÞ, the conditional probability density function is Gaussian:
1 f YjX ðjXÞ N jmY þ RYX R1 XX ðX mX Þ; RYY RYX RXX RXY :
Therefore, the conditional distribution of cn N
ðC:1Þ
lcn ; Rcn given Dh N lDh ; RDh follows Gaussian distribution:
N1 1 1 ¼ N cn Gn xL sLn þ cn Dh N cn lcn þ Rcn R1 Dh G x L s L ; I : Dh Dh lDh ; Rcn Rcn RDh Rcs N N2
ðC:2Þ
References [1] [2] [3] [4] [5] [6] [7]
[8] [9] [10]
A.F. Colombo, P. Lee, B.W. Karney, A selective literature review of transient-based leak detection methods, J. Hydro-Environ. Res. 2 (4) (2009) 212–227. J.A. Liggett, L.-C. Chen, Inverse transient analysis in pipe networks, J. Hydraul. Eng. 120 (8) (1994) 934–955. B. Brunone, Transient test-based technique for leak detection in outfall pipes, J. Water Resour. Plann. Manage. 125 (5) (1999) 302–306. J.P. Vítkovsky`, A.R. Simpson, M.F. Lambert, Leak detection and calibration using transients and genetic algorithms, J. Water Resour. Plann. Manage. 126 (4) (2000) 262–265. W. Mpesha, S.L. Gassman, M.H. Chaudhry, Leak detection in pipes by frequency response method, J. Hydraul. Eng. 127 (2) (2001) 134–147. X.-J. Wang, M.F. Lambert, A.R. Simpson, J.A. Liggett, J.P. Vítkovsky`, Leak detection in pipelines using the damping of fluid transients, J. Hydraul. Eng. 128 (7) (2002) 697–711. P.J. Lee, J.P. Vitkovsky, M.F. Lambert, A.R. Simpson, J.A. Liggett. Frequency response coding for the location of leaks in single pipeline systems, in: The Int. Conference on Pumps, Electromechanical Devices and Systems Applied to Urban Water Management, IAHR and IHRA, Valencia, Spain, April 22–25, 2003. P.J. Lee, J.P. Vítkovsky`, M.F. Lambert, A.R. Simpson, J.A. Liggett, Frequency domain analysis for detecting pipeline leaks, J. Hydraul. Eng. 131 (7) (2005) 596–604. P.J. Lee, J.P. Vítkovsky`, M.F. Lambert, A.R. Simpson, J.A. Liggett, Leak location using the pattern of the frequency response diagram in pipelines: a numerical study, J. Sound Vib. 284 (3) (2005) 1051–1073. P.J. Lee, M.F. Lambert, A.R. Simpson, J.P. Vítkovsky`, J. Liggett, Experimental verification of the frequency response method for pipeline leak detection, J. Hydraul. Res. 44 (5) (2006) 693–707.
362
X. Wang, M.S. Ghidaoui / Mechanical Systems and Signal Processing 119 (2019) 346–362
[11] H.-F. Duan, P.J. Lee, M.S. Ghidaoui, Y.-K. Tung, Essential system response information for transient-based leak detection methods, J. Hydraul. Res. 48 (5) (2010) 650–657. [12] H.-F. Duan, P.J. Lee, M.S. Ghidaoui, Y.-K. Tung, Leak detection in complex series pipelines by using the system frequency response method, J. Hydraul. Res. 49 (2) (2011) 213–221. [13] H.-F. Duan, P.J. Lee, M.S. Ghidaoui, Y.-K. Tung, System response function–based leak detection in viscoelastic pipelines, J. Hydraul. Eng. 138 (2) (2011) 143–153. [14] W. Nixon, M.S. Ghidaoui, Numerical sensitivity study of unsteady friction in simple systems with external flows, J. Hydraul. Eng. 133 (7) (2007) 736– 749. [15] M.F. Ghazali, S.B.M. Beck, J.D. Shucksmith, J.B. Boxall, W.J. Staszewski, Comparative study of instantaneous frequency based methods for leak detection in pipeline networks, Mech. Syst. Sig. Process. 29 (2012) 187–200. [16] A.C. Zecchin, L.B. White, M.F. Lambert, A.R. Simpson, Parameter identification of fluid line networks by frequency-domain maximum likelihood estimation, Mech. Syst. Sig. Process. 37 (1) (2013) 370–387. [17] S. Meniconi, B. Brunone, M. Ferrante, C. Capponi, C.A. Carrettini, C. Chiesa, D. Segalini, E.A. Lanfranchi, Anomaly pre-localization in distribution– transmission mains by pump trip: preliminary field tests in the milan pipe system, J. Hydroinf. 17 (3) (2015) 377–389. [18] S. Meniconi, B. Brunone, M. Ferrante, C. Massari, Potential of transient tests to diagnose real supply pipe systems: what can be done with a single extemporary test, J. Water Resour. Plann. Manage. 137 (2) (2010) 238–241. [19] B. Brunone, M. Ferrante, Detecting leaks in pressurised pipes by means of transients, J. Hydraul. Res. 39 (5) (2001) 539–547. [20] D. Covas, H. Ramos, N. Graham, C. Maksimovic, Application of hydraulic transients for leak detection in water supply systems, Water Sci. Technol.: Water Supply 4 (5–6) (2005) 365–374. [21] S.T.N. Nguyen, J. Gong, M.F. Lambert, A.C. Zecchin, A.R. Simpson, Least squares deconvolution for leak detection with a pseudo random binary sequence excitation, Mech. Syst. Sig. Process. 99 (2018) 846–858. [22] M. Ferrante, B. Brunone, S. Meniconi, Wavelets for the analysis of transient pressure signals for leak detection, J. Hydraul. Eng. 133 (11) (2007) 1274–1282. [23] J.C. Liou, Pipeline leak detection by impulse response extraction, J. Fluids Eng. 120 (4) (1998) 833–838. [24] M. Ferrante, B. Brunone, Pipe system diagnosis and leak detection by unsteady-state tests. 1. harmonic analysis, Adv. Water Resour. 26 (1) (2003) 95–105. [25] D. Covas, H. Ramos, A.B. De Almeida, Standing wave difference method for leak detection in pipeline systems, J. Hydraul. Eng. 131 (12) (2005) 1106–1116. [26] A.M. Sattar, M.H. Chaudhry, Leak detection in pipelines by frequency response method, J. Hydraul. Res. 46 (EI1) (2008) 138–151. [27] M. Taghvaei, S.B.M. Beck, J.B. Boxall, Leak detection in pipes using induced water hammer pulses and cepstrum analysis, Int. J. Comadem. 13 (1) (2010) 19. [28] M.L. Stephens, Transient Response Analysis for Fault Detection and Pipeline Wall Condition Assessment in Field Water Transmission and Distribution Pipelines and Networks PhD Thesis, University of Adelaide, 2008. [29] D. Covas, H. Ramos, Case studies of leak detection and location in water pipe systems by inverse transient analysis, J. Water Resour. Plann. Manage. 136 (2) (2010) 248–257. [30] X. Wang, M.S. Ghidaoui, Pipeline leak detection using the matched-field processing method, J. Hydraul. Eng. 144 (6) (2018) 04018030. [31] X. Wang, M.S. Ghidaoui, Identification of multiple leaks in pipeline: linearized model, maximum likelihood, and super-resolution localization, Mech. Syst. Sig. Process. 107 (2018) 529–548. [32] B. Zhou, A. Liu, X. Wang, Y. She, V. Lau, Compressive sensing-based multiple-leak identification for smart water supply systems, IEEE IoT J. 5 (2) (2018) 1228–1241. [33] X. Wang, D.P. Palomar, L. Zhao, M.S. Ghidaoui, R.D. Murch, Spectral-based methods for pipeline leakage detection, J. Hydraul. Eng. (2019), https://doi. org/10.1061/(ASCE)HY.1943–7900.0001572, In Press. [34] X. Wang, B. Quost, J.-D. Chazot, J. Antoni, Estimation of multiple sound sources with data and model uncertainties using the EM and evidential EM algorithms, Mech. Syst. Sig. Process. 66–67 (2016) 159–177. [35] X. Wang, B. Quost, J.-D. Chazot, J. Antoni, Iterative beamforming for identification of multiple broadband sound sources, J. Sound Vib. 365 (2016) 260–275. [36] X. Wang, S. Khazaie, L. Margheri, P. Sagaut, Shallow water sound source localization using the iterative beamforming method in an image framework, J. Sound Vib. 395 (2017) 354–370. [37] A.P. Dempster, N.M. Laird, D.B. Rubin, Maximum likelihood from incomplete data via the EM algorithm, J. R. Stat. Soc., Ser. B (Method.) 39 (1) (1977) 1–38. [38] J.C.F. Wu, On the convergence properties of the EM algorithm, Ann. Stat. 11 (1983) 95–103. [39] M. Feder, E. Weinstein, Parameter estimation of superimposed signals using the EM algorithm, IEEE Trans. Acoust., Speech Sig. Process. 36 (4) (1988) 477–489. [40] H. Akaike, A new look at the statistical identification model, IEEE Trans. Autom. Control 19 (1974) 716–723. [41] G. Schwarz, Estimating the dimension of a model, Ann. Stat. 6 (2) (1978) 461–464. [42] K.P. Burnham, D.R. Anderson, Model Selection and Multimodel Inference: A Practical Information-theoretic Approach, Springer Science & Business Media, 2003. [43] C.M. Hurvich, C.-L. Tsai. Regression and time series model selection in small samples. Biometrika, pages 297–307, 1989. [44] X. Wang, Z. Zhang, S. Li, Set-valued and interval-valued stationary time series, J. Multivariate Anal. 145 (2016) 208–223. [45] K. Aho, D. Derryberry, T. Peterson, Model selection for ecologists: the worldviews of AIC and BIC, Ecology 95 (3) (2014) 631–636. [46] J.P. Vítkovsky`, M. F Lambert, A.R. Simpson, J.A. Liggett, Experimental observation and analysis of inverse transients for pipeline leak detection, J. Water Resour. Plann. Manage. 133 (6) (2007) 519–530. [47] M.S. Ghidaoui, M. Zhao, D.A. McInnis, D.H. Axworthy, A review of water hammer theory and practice, Appl. Mech. Rev. 58 (1) (2005) 49–76. [48] M.S. Ghidaoui, On the fundamental equations of water hammer, Urban Water J. 1 (2) (2004) 71–83. [49] I.B. Rhodes, A tutorial introduction to estimation and filtering, IEEE Trans. Automatic Control AC-16 (6) (1971) 688–706. [50] M.H. Chaudhry, Applied Hydraulic Transients, Third Edition., Springer, 2014. [51] A. Kashima, P.J. Lee, R. Nokes, Numerical errors in discharge measurements using the KDP method, J. Hydraul. Res. 50 (1) (2012) 98–104. [52] A. Kashima, P.J. Lee, M.S. Ghidaoui, M. Davidson, Experimental verification of the kinetic differential pressure method for flow measurements, J. Hydraul. Res. 51 (6) (2013) 634–644.
Contents lists available at ScienceDirect
Mechanical Systems and Signal Processing journal homepage: www.elsevier.com/locate/ymssp
Identification of multiple leaks in pipeline II: Iterative beamforming and leak number estimation Xun Wang ⇑, Mohamed S. Ghidaoui Department of Civil and Environmental Engineering, Hong Kong University of Science and Technology, Clear Water Bay, Hong Kong, China
a r t i c l e
i n f o
Article history: Received 1 November 2017 Received in revised form 12 July 2018 Accepted 16 September 2018
Keywords: Leak detection Fluid transient Iterative beamforming EM algorithm Model selection
a b s t r a c t This paper considers the problem of detecting multiple leaks in a water-filled pipe using transient waves. A previous paper by the authors proposes a maximum likelihood (ML) method for leak detection based on an approximate linear model of wave propagation in pipe and shows that it is efficient, robust with noise, and super-resolved. However, that ML method needs to solve an N-parameter optimization problem (N is the leak number) which is high dimensional, non-convex and has many local maxima when N is large. The present paper deals with the multiple-leak detection problem where the leak number is high. A ML based iterative procedure, known as the iterative beamforming (IB) method, is proposed. Given an initial value of model parameter (locations and sizes of leaks) in each iteration, the complex inverse problem with 2N unknown parameters is decomposed into N onedimensional optimization problems, which largely decreases the computational complexity. More specifically, IB iteratively estimates the contributions from various leaks to the pressure head measurement and then performs a 1-D ML method on these estimates. In order to estimate the unknown number of leaks, model selection methods, known as Akaike information criterion (AIC) and Bayesian information criterion (BIC), are used and they are proven to be efficient to estimate the actual leak number. Numerical simulations show that five leaks can be accurately identified with super-resolution and the leak number can be correctly estimated using the methods proposed in this paper. Ó 2018 Elsevier Ltd. All rights reserved.
1. Introduction Leaks in water supply pipelines is a problem of increasing interest in the recent decades. Transient-based detection methods (TBDMs) utilize the hydraulics of the transient flows to detect leaks in the pipeline [1–18]. The reason that such methods are expected to work is that a leak in a pipeline system modifies the transient wave passing the leaks and measured by a sensor. Many TBDMs have been developed by researchers and applied to water piping systems including: (i) the transient reflection based method (TRM) [3,19–22]; (ii) the transient damping based Method (TDM) [6]; and (iii) the combination of both reflection and damping properties [2,4,5,7–10,15,16,23–32]. Recently, the problem of identifying multiple leaks has attracted increasing research attentions. Inverse transient analysis (ITA) method [2,4,29] is able to estimate multiple leaks, which assumes some discrete points in the pipe are potential leaks and estimates the corresponding leak sizes of these potential leaks by matching the measured pressure with time-domain numerical transient model. However, this approach involves a high-dimensional optimization problem (the dimension equals to the number of potential leaks), which implies a high computational cost and, more importantly, a high ⇑ Corresponding author. E-mail addresses: [email protected], [email protected] (X. Wang). https://doi.org/10.1016/j.ymssp.2018.09.020 0888-3270/Ó 2018 Elsevier Ltd. All rights reserved.
X. Wang, M.S. Ghidaoui / Mechanical Systems and Signal Processing 119 (2019) 346–362
347
Nomenclature q discharge oscillation h head oscillation xL leak location sL leak size steady-state discharge and head of leak Q L0 ; HL0 xm (m ¼ 1; ; M) sensor coordinate
Dh c n a A l D g Z R f DW
l x xth
kmin M N J T
r2
log L
head difference latent leak contribution measurement noise wave speed area of pipeline pipe length pipe diameter gravitational acceleration characteristic impedance frictional resistance Darcy–Weisbach friction factor propagation function angular frequency fundamental frequency minimum wavelength sensor number leak number frequency number sample size variance of noise log-likelihood function
Acronyms AIC Akaike information criterion BIC Bayesian information criterion EM expectation–maximization FRF frequency response function IB iterative beamforming MFP matched-field processing MFP(1) MFP with one-leak model ML maximum likelihood SNR signal-to-noise ratio
computation complexity such that the optimization problem often stops at local maxima corresponding to a wrong estimate of leak. For the same reason, ITA is not robust because the high-dimensional objective function is easy to be distorted by noise and uncertainties. Furthermore, the precision of ITA leak localization is limited to the mesh spacing of potential leaks in the pipe. A recent transient-based leak detection method, known as the matched-field processing (MFP) method [30], estimates leak location and size by matching the measurement with the frequency-domain model; it is found to be fast, accurate, robust with respect to noise and uncertainties, and can localize multiple leaks. However, MFP can only identify separate leaks provided the distance between leaks is not shorter than kmin =2, where kmin is the shortest probing wavelength. Based on the linear approximation model of wave propagation in a pipe with multiple leaks [31], a previous paper by the same authors has shown the good properties of the maximum likelihood (ML) approach for the multiple leaks [31]. That ML-based localization method has proven to be accurate, consistent, robust with respect to noise, and super-resolved. The latter property means that two leaks closer than kmin =2 can be separately estimated. In the case of one leak or two leaks, the localization can be easily accomplished by an exhaustive searching that plots the likelihood function with respect to possible leak locations. However, for a higher number of leaks, it is computationally expensive to compute the likelihood function for all the possible combinations of leak locations and it is complicated to directly solve the high dimensional optimization problem since the objective function is high-dimensional (with N parameters if N leaks are assumed) and non-convex. Therefore, the main goal of this paper is to develop an efficient method for leak detection problem where there are multiple leaks whose number, locations and sizes are all unknown. By the linearized model proposed in [31], each pressure head measured by a sensor in a pipe equals to the summation of contributions from various leaks and a term independent of all the leaks. It is clear that if these contributions could be
348
X. Wang, M.S. Ghidaoui / Mechanical Systems and Signal Processing 119 (2019) 346–362
Fig. 1. Setup of pipe system.
known, the estimation of leaks would be as straightforward as in the single leak case [30,33]. Based on this idea, an iterative method, known as the iterative beamforming (IB) [34–36], which deals with the parameter estimation of superimposed signals using the Expectation–Maximization (EM) algorithm [37–39], is proposed to solve the aforementioned problem. More specifically, IB alternates iteratively between two steps. The first one consists of computing the expected contributions from various leaks to the head measurements given the current fits of the model parameters. The location and size of each leak are then updated according to the leak contributions from the first step. The estimation of leak number is also an important issue in real applications. The ML approach cannot cope with this problem since it tends to choose the highest possible leak number which corresponds to highest possible likelihood function. In order to avoid model over-fitting, a trade-off between the goodness of fit (assessed by the likelihood function) and the model complexity (the number of leaks assumed in the model) is considered which intends to find an optimal estimate of leak number. Akaike information criterion (AIC) [40] and Bayesian information criterion (BIC) [41] are both the model selection methods for this issue and have been widely and successfully used in various domains [42–46]. In this paper, AIC and BIC are applied to decide the number of leaks. The organization of this paper is as follows. Section 2 begins with model description and problem statement. Then, estimation of multiple leaks using the IB approach is introduced in Section 3. Section 4 presents results of numerical simulation which justifies the proposed IB method is able to estimate multiple leaks with a super-resolution. Leak number determination via AIC and BIC is introduced in Section 5. Finally, conclusions are drawn in Section 6. 2. Model description and problem statement 2.1. Linearized model of wave propagation in pipeline The setup of the considered pipe system is shown in Fig. 1. A single pipeline bounded by an upstream and a downstream > node with coordinates x ¼ xU ¼ 0 and x ¼ xD is considered. The locations and sizes of N leaks, denoted as xL ¼ xL1 ; ; xLN > and sL ¼ sL1 ; ; sLN , are parameters to be estimated. Each leak size is represented by the lumped leak parameter sLn ¼ C d ALn , where C d is the discharge coefficient of the leak and ALn is the flow area of the leak orifice. The steady-state disrffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi charge of a leak is related to the lumped leak parameter by Q L0n ¼ sLn 2g HL0n zLn , in which g is the gravitational acceleration, zLn denotes the elevation of the pipe at each leak, and Q L0n and HL0n are the steady-state discharge and head1 at each leak. The leaks are estimated by pressure head measurement at M sensors. The corresponding head difference due to leaks is used, whose theoretical expression can be approximated via a linearized model [31]: N X Dh G x L s L þ n ¼ Gn xLn sLn þ n:
ð1Þ
n¼1
The derivation of Eq. (1) is detailed in Appendix A. In Eq. (1),
> Dh ¼ Dh11 ; ; DhJ1 ; ; Dh1M ; ; DhJM ; M
ð2Þ
NL
M
in which Dhjm ¼ hjm hjm (j ¼ 1; ; J and m ¼ 1; ; M) denotes the head difference between the head measurement hjm in the presence of leaks and theoretical head that does not include the leak terms NL
NL
hjm ¼ h
xj ; xm ¼ Z sinh ðlxm Þq xU þ cosh ðlxm Þh xU ;
ð3Þ
at the j-th frequency xj and at the m-th measurement station xm . Physically, Dhjm stands for the scattering of the wave at the pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi leaks. In Eq. (3), Z ¼ la2 =ð ixgAÞ is the characteristic impedance, l ¼ a1 x2 þ igAxR is the propagation function, a is the wave speed, x is the angular frequency, A is the area of pipeline, and R is the frictional resistance term. Note that 1 Discharge and head are terms in fluid dynamics and hydraulics. The discharge q (in m3 s1 ) is the volumetric flow rate of fluid and it relates to the flow velocity v (in ms1) by q ¼ Av , where A is the cross-sectional area of the pipe. The pressure head is defined as h ¼ p=ð qg Þ, where p is the pressure (in Pascal), g is the gravitational acceleration and q is the density of the fluid. For example, the atmospheric pressure is approximately patm ¼ 101 kPa. In terms of head, h ¼ 1:01 105 Pa= 103 kg m3 9:8m s2 ¼ 1:03m. This states that the atmospheric pressure is equivalent to pressure produced by 1.03 m of water. Alternatively, a 50 m of head in pipe implies that if that pipe bursts, the height of the water jet would be about 50 m.
X. Wang, M.S. Ghidaoui / Mechanical Systems and Signal Processing 119 (2019) 346–362
349
R ¼ ðf DW Q 0 Þ= gDA2 for turbulent flows where f DW is the Darcy–Weisbach friction factor, Q 0 is the steady-state discharge in the pipe and D is the pipe diameter. The matrix G in Eq. (1) is a JM N-dimensional matrix whose n-th column is
> Gn xLn ¼ G x1 ; xLn ; x1 ; ; G xJ ; xLn ; x1 ; ; G x1 ; xLn ; xM ; ; G xJ ; xLn ; xM ;
ð4Þ
in which
pffiffiffi g Z sinh l xm xLn Ln rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Z sinh lxLn q xU cosh lxLn h xU : G xj ; x ; xm ¼ 2 HL0n zLn
ð5Þ
The random noise vector
> n ¼ n11 ; ; nJ1 ; ; n1M ; ; nJM
ð6Þ
follows independent and identically distributed (i.i.d.) Gaussian distribution with 0-mean and covariance matrix r2 I, where I is an identity matrix. Note that the discharge at upstream q xU in Eqs. (3) and (5) can be estimated [31] by
U coshðlÞh xU hðx0 Þ ^ q x ¼ ; Z sinhðlðx0 xU ÞÞ
ð7Þ
if the pressure at x0 (x0 > xU ) is measured, boundary condition h xU is applied, and no leak is between xU and x0 . The derivation of Eq. (7) is detailed in Appendix B. It must be emphasized that Eq. (1) involves a number of assumptions. First, the nonlinear convective terms are neglected because they are of the order of Mach number and this Mach number is very small (103) in practice [47,48]. Second, the friction, valves, pumps, etc., in the system are given by nonlinear relations. The derivation of Eq. (1) requires that these relations are linearized. Such linearization is justified by the fact that injected transient waves for defect detection need to be small so as not to compromise the structural integrity of the system. Third, it is shown in [31] that leak detection is not affected by leak to leak interactions for small leaks. This is not a perious limitation since the motive is to detect small leaks. However, the leak to leak interaction can be accounted for but doing so makes the leak localization computationally demanding. 2.2. Problem statement Based on the linearized model Eq. (1), a ML solution for estimating multiple leaks is proposed in [31], which is briefly recalled here. Since the vector of head differences Dh follows i.i.d. Gaussian distribution, its log-likelihood function, which is a function of the leak locations xL and the leak sizes sL , has the form (after removing unnecessary terms which are independent of parameters to be estimated):
2 log L xL ; sL jDh ¼ Dh G xL sL :
ð8Þ
^L and sizes ^sL of the The MLE of xL and sL can be obtained by maximizing Eq. (8), which leads to the estimates of locations x leaks:
1 H L ^L ¼ arg max DhH G xL GH xL G xL x G x Dh xL
ð9Þ
and
1 ^sL ¼ GH x ^L G x ^L Dh: ^L GH x
ð10Þ
In the numerical examples in [31], two leaks are localized by plotting the 2D function in Eq. (9) and finding its maximum. However, in the case that the leak number N is large, this approach (exhaustive search of locations of leaks) is computationally intensive and inefficient. For example, considering a pipe with length l ¼ 2000 m and localization precision (grid density) being 1 m (as the numerical examples in [31]), the function in Eq. (9) needs to be computed 2000N times if N leaks are assumed, which implies a very high computational cost when N is large. Another approach is to directly solve the optimization problem (9). However, in the case of large leak number, this optimization problem is still difficult to solve due to the complex nature of G xL : the objective function in Eq. (9) is non-convex and has many local maxima in a highdimensional space [33]. The present paper addresses this shortcoming. More specifically, this paper proposes an indirect but simpler (in the sense of the complexity of optimization problem) method, which only needs to solve 1D optimization problems.
350
X. Wang, M.S. Ghidaoui / Mechanical Systems and Signal Processing 119 (2019) 346–362
3. Estimation of multiple leaks using the iterative beamforming method This section introduces the IB method for identifying multiple leaks. For this purpose, latent leak contribution from various leaks to pressure head measurement is introduced first. 3.1. Latent leak contributions By introducing a latent variable,2 called latent leak contribution, the optimization problem (9) can be largely simplified. The latent leak contribution is defined as the contribution from a single leak to the head difference (the data for leak detection). More specifically, the latent leak contribution from the n-th leak, denoted as cn , is
cn ¼ Gn xLn sLn þ nn :
ð11Þ
Here, the noise term nn can be obtained by arbitrarily decomposing the total noise n in Eq. (1) into N components such that their summation equals in distribution to n, i.e., N X nn ¼ n:
ð12Þ
n¼1
Therefore, the head difference Dh is related to the latent leak contributions cn by
Dh ¼
N X
cn :
ð13Þ
n¼1
In this equation, ‘‘=” means equality in distribution of random variables. By assuming that nn (n ¼ 1; ; N) are i.i.d. Gaussian random vectors with 0-mean and covariance matrix Rn ¼ r2 =N I which satisfies Eq. (12), the vector of latent leak contributions c ¼ ðc1 ; ; cN Þ> also follows a Gaussian distribution and the log-likelihood function of c is (again, removing unnecessary terms independent of model parameters): N X 2 log L xL ; sL jc ¼ cn Gn xLn sLn :
ð14Þ
n¼1
Note that each term in Eq. (14) corresponds to a single leak and is independent of other leaks. Therefore, if the contribution cn of the n-th leak is available, the corresponding leak location xLn and size sLn can be estimated by solving a single leak localization problem:
Ln Ln 2 ^x ; ^s ¼ arg min cn Gn xLn sLn ; xLn ;sLn
ð15Þ
as in [30]. 3.2. Iterative beamforming for detecting multiple leaks Unfortunately, the latent leak contributions cn cannot be directly observed. However, the IB method [34,35] is able to proceed with Eq. (14) by treating the latent contributions as random variables; this approach is based on a ML method for dealing with latent variables or missing data, known as the EM algorithm [37–39]. Here, this approach is applied to the multipleleak detection problem. The proposed method is able to maximize the observed data log-likelihood (8) by computing iteratively the expectation of the log-likelihood (14) of the latent leak contribution c, given the observed data and the current fit of the parameter. More
precisely, given an initial parameter U0 ¼ xL0 ; sL0 , the algorithm iterates back and forth between two steps until convergence: Expectation step (E-step): computing the conditional expectation
Q ðUjUk1 Þ ¼ E log L xL ; sL jc jDh; Uk1
ð16Þ
of the complete data log-likelihood with respect to the unknown variables, given the observed data Dh and a current fit
Uk1 for the model parameters;
Maximization step (M-step): determining the new fit for the parameters by maximizing Q ðUjUk1 Þ with respect to U. The algorithm guarantees that the observed data likelihood increases after each iteration [37]; whether it converges to the global maximum depends on the choice of the initial parameter U0 [38]. A strategy to deal with this problem is to 2
A latent variable is an unknown (or unobserved) information that would make the estimation process straightforward, should it be available.
X. Wang, M.S. Ghidaoui / Mechanical Systems and Signal Processing 119 (2019) 346–362
351
run the algorithm with different initial parameter values that cover a large part of the parameter space, to compute the corresponding estimates, and to finally retain the result associated with the highest likelihood function. In the following, the Eand M-steps are detailed, respectively. 3.2.1. E-step: Estimation of contributions from various leaks The E-step consists in computing N N X X 2 Eðcn jDh; Uk1 Þ Gn xLn sLn 2 ; Q ðUjUk1 Þ ¼E log L xL ; sL jc jDh; Uk1 ¼ E cn Gn xLn sLn jDh; Uk1 ¼ c n¼1
n¼1
ð17Þ where c is a constant independent of the parameters. Therefore, the E-step is to compute the expected leak contributions from various leaks given the current fit of model parameters, i.e., Eðcn jDh; Uk1 Þ. >
Note that both Dh and cn are Gaussian distributed, thus the vector ðDh; cn Þ is jointly Gaussian with 0-mean and covariance matrix
R¼
CovðDhÞ
CovðDh; cn Þ
Covðcn ; DhÞ
Covðcn Þ
¼
I 1 I N
1 I N 1 I N
! :
ð18Þ
Therefore, the conditional distribution of cn given Dh is Gaussian:
N1 1 cn Dh N cn Gn xLn sLn þ Dh G x L s L ; I : N N2
ð19Þ
Eq. (19) is derived from a property of Gaussian random variable [49], which is detailed in Appendix C. Furthermore,
n Ln 1 ^cnk ¼ Eðcn jDh; Uk1 Þ ¼ Gn xLk1 Dh G xLk1 sLk1 : sk1 þ N
ð20Þ
Finally, the conditional expectation in the E-step is obtained:
Q ðUjUk1 Þ ¼ c
N X ^cnk Gn xLn sLn 2 :
ð21Þ
n¼1
3.2.2. M-step: Single leak detection problem The M-step amounts to computing new estimates of xLn and sLn by maximizing Eq. (21), which is decomposed into N twodimensional minimization problems:
Ln Ln xk ; sk ¼ arg min k^cnk Gn xLn sLn k2 ; n ¼ 1; ; N: xLn ;sLn
ð22Þ
Eq. (22) can be further simplified as follows. For any leak location estimate xLkn , the corresponding leak size is
sLkn ¼ arg mink^cnk Gn xLkn sLn k2 ¼ sLn
GHn xLkn ^cnk : GHn xLkn Gn xLkn
ð23Þ
Then, inserting Eq. (24) into Eq. (22) results in the leak location estimate:
xLkn
2 GHn xLn Gn xLn ^ ^cnk ¼ arg mincnk H xLn Gn ðxLn ÞGn ðxLn Þ L H Ln G x Gn x n ^cnk : ¼ arg maxc^Hnk nH xLn Gn ðxLn ÞGn ðxLn Þ
ð24Þ
It is clear that the leak locations and sizes computed at each iteration amount to performing a Bartlett’s beamforming estimation (equivalently, MFP, ML or matched-filter [30,33]) given each latent leak contribution cn ; each leak location is determined via the single-leak detection method (i.e., a single-parameter optimization problem is solved) by using the corresponding estimated contributions. Finally, the IB strategy for estimating the model parameters is summarized in Algorithm 1.
352
X. Wang, M.S. Ghidaoui / Mechanical Systems and Signal Processing 119 (2019) 346–362
Algorithm 1 Identification of multiple leaks using IB > > For k ¼ 0, pick starting values xL0 ¼ xL01 ; ; xL0N and sL0 ¼ sL01 ; ; sL0N for the model parameters. For k P 1: repeat estimate the latent leak contributions ^ cnk , n ¼ 1; ; N, by Eq. (20); update the leak location xLkn , n ¼ 1; ; N, by Eq. (24); update the leak size sLkn , n ¼ 1; ; N, by Eq. (23); until the relative increase of the observed data log-likelihood (8) is less than a given threshold
j:
log LxL ; sL jDh log LxL ; sL jDh k k k1 k1 < j: log L xLk1 ; sLk1 jDh
ð25Þ
4. Numerical simulations In this section, numerical simulations are introduced to illustrate the performance of the IB algorithm. First, an example shown in [33] that cannot be well-solved by the 1D searching methods in [30,33] is studied. Then, a more complicated case with five leaks is considered. 4.1. Revisiting the example in [33]: Identifying two leaks Here, a numerical example in [33] is revisited. A single pipe which is connected by two reservoirs in a horizontal plane is considered. The pressure heads of both reservoirs are HU ¼ 25 m and HD ¼ 20 m. The pipe length is l ¼ 2 km and the diameter is D ¼ 0:5 m. The wave speed is 1200 m/s. The Darcy–Weisbach friction factor f DW ¼ 0:02. Two leaks with locations xL1 ¼ 600 m and xL2 ¼ 1600 m and sizes sL1 ¼ 1 104 m2 and xL2 ¼ 1:2 104 m2 are assumed. Three sensors are located at 50 m, 1800 m and 1960 m: the first sensor is for estimating the boundary condition q xU via Eq. (7) and the pressure heads from the latter two are then used to estimate leak locations and sizes. The resonant and anti-resonant frequencies fx : ax; a ¼ 1; 2; ; 31g are used for leak detection. The signal-to-noise ratio (SNR), which is defined by
SNR ¼ 20 log10 jDhj=r ;
ð26Þ
is 15 dB, jDhj represents the average head difference over frequencies and r stands for the standard deviation of Gaussian white noise. Fig. 2(a) shows the output function of MFP [30,33] with single leak assumption, abbreviated as MFP(1). It is clear that near each actual leak, the output function of MFP(1) has a local maximum, but the other two local maxima prevent obtaining correct estimates of leak locations. Then, the leak detection methods, i.e., directly solving the optimization problem Eq. (9) and IB, based on the model of multiple leaks in Eq. (1) are used. The former is abbreviated as ML in this section. In both methods, four leaks are assumed (i.e., N ¼ 4 in the model) according to the result of MFP: the initial values of leak locations are taken from the four local maxima in Fig. 2(a) (377 m, 576 m, 1391 m and 1604 m) and the initial leak sizes are 2 104 m2 for all the four leaks (it is arbitrarily chosen). The optimization method used in both approaches is the Broyden-Fletcher-Goldfarb-Shanno (BFGS) QuasiNewton algorithm implemented by the MATLAB function ‘‘fminunc.m”. The results of both leak detection methods are shown in Fig. 2(b) and (c), wherein the locations and sizes of the four estimated leaks (circles), along with the actual leaks (represented by the solid lines), are plotted. It is clear that both methods can accurately localize the leaks and return similar results. Taking the IB result as an example, two estimates are very close to the actual leaks, with locations ^ xL1 ¼ 600:1 m and 4 4 L L 2 L 2 ^ x 2 ¼ 1600 m and sizes ^s 1 ¼ 0:995 10 m and ^s 2 ¼ 1:197 10 m . The other two estimates have very small values of ^L4 ¼ 1238 m) and thus can be ignored. Here, the leak size (^sL3 ¼ 4:3 107 m2 at ^ xL3 ¼ 423:8 m and ^sL4 ¼ 5:2 107 m2 at x IB algorithm converges at 81st iteration with threshold of relative increase of log-likelihood
j ¼ 104 .
4.2. A five-leak example In this section, a more complicated case with five leaks is considered. The leak locations are xL1 ¼ 600 m, xL2 ¼ 640 m, L3
x
¼ 800 m, xL4 ¼ 1540 m, xL5 ¼ 1600 m; the leak sizes are sL1 ¼ 1 104 m2 ; sL2 ¼ 1:2 104 m2 ; sL3 ¼ 0:9 104
m ; sL4 ¼ 1:4 104 m2 ; sL5 ¼ 1 104 m2. As previously, the resonant and anti-resonant frequencies fx : ax; a ¼ 1; 2; ; 31g are used for leak detection. Note that with this probing wave whose maximum angular frequency 2
X. Wang, M.S. Ghidaoui / Mechanical Systems and Signal Processing 119 (2019) 346–362
353
Fig. 2. Leak detection using MFP(1) (a), ML (b), and IB (c). In subfigure (a), the solid line, dash lines and crosses stand for the objective function of MFP(1), actual leak locations and sensor locations. In subfigures (b) and (c), the solid lines and circles represent the actual and estimated leaks, respectively. The leak locations are xL1 ¼ 600 m and xL2 ¼ 1600 m and the leak sizes are sL1 ¼ 1 104 m2 and sL2 ¼ 1:2 104 m2.
is 31xth , the resolution (the minimum range of two leaks that can be separately estimated) is approximately kmin =2 ¼ 129 m according to the Nyquist-Shannon sampling theorem and the observations from the simulation results in [30,33]. However, the distance between xL1 and xL2 (40 m, 0:15kmin ) and the distance between xL4 and xL5 (60 m, 0:23kmin ) are smaller than the minimum resolvable leak distance kmin =2. In the following, IB is justified to be able to separately estimate the five leaks, which implies that IB is a super-resolution scheme. Fig. 3(a) shows the localization results using MFP(1). Two local maxima, which respectively locate between xL1 and xL2 and between xL4 and xL5 , can be found. However, other local maxima disturb the localization. Then, ML and IB are used to detect the leaks. The results in Fig. 3(a) are used for the parameter initialization. The seven highest local maxima are all considered as possible locations of leaks and thus are taken as initial points. Since the leak number is unknown, a large number of leaks is assumed (the determination of leak number is discussed in Section 5). More specifically, two initial leaks are attributed to each of the first five highest local maxima; 435 5 m, 614 5 m, 1200 5 m, 1353 5 m, 1571 5 m are taken as the initial leak locations. Another two local maxima, being 742 m and 1515 m, are also assumed as two initial leak locations, thus 12 leaks are assumed in the IB algorithm. Furthermore, the initial leak sizes are assumed to be 2 104 m2 for all the 12 leaks (it is arbitrarily chosen since no prior information is available). Fig. 3(b) and (c) and Table 1 both display the actual leak
354
X. Wang, M.S. Ghidaoui / Mechanical Systems and Signal Processing 119 (2019) 346–362
Fig. 3. Leak detection using MFP(1) (a), ML (b) and IB (c). In subfigure (a), the solid line, dash lines and crosses stand for the objective function of MFP(1), actual leak locations and sensor locations. In subfigures (b) and (c), the solid lines, crosses and circles represent the actual, initial and estimated leaks, respectively. Each estimate of leak is connected to its corresponding initial parameter by a dash line. The leak locations are xL1 ¼ 600 m, xL2 ¼ 640 m, xL3 ¼ 800 m, xL4 ¼ 1540 m, xL5 ¼ 1600 m. The leak sizes are sL1 ¼ 1 104 m2 ; sL2 ¼ 1:2 104 m2 ; sL3 ¼ 0:9 104 m2 ; sL4 ¼ 1:4 104 m2 ; sL5 ¼ 1 104 m2.
locations and sizes, their initial parameters and the estimates obtained from ML and IB. Each initial parameter and the corresponding estimates are connected by a dash line in Fig. 3(b) and (c). It can be seen that ML accurately estimate all the five leaks. However, two ‘‘ghost leaks” with high estimated leak sizes appear at around 400 m and 1400 m. By contrast, the result obtained from IB is more correct that no misleading leak estimate appears. More specifically, all the leak locations are accurately estimated except xL4 ¼ 1540 m which is replaced by two estimates at 1533 m and 1560 m. The sizes of the first three leaks are accurately estimated while the fourth and fifth are less precise. In particular, the sum of the estimated sizes at 1533 m and 1560 m approximates the actual size of the fourth leak. It is also interesting to remark that although 12 leaks is assumed, half of them almost disappears in the IB results: their estimated sizes are very small and thus neglectable. The previous result is obtained with a single initialization which is acquired from the result of MFP(1). Alternatively, a multiple-initialization strategy can be used to further improve the leak detection result of IB. Algorithm 1 is run with various starting parameters and then the result with highest likelihood function Eq. (8) is retained as the final estimate. Instead of
355
X. Wang, M.S. Ghidaoui / Mechanical Systems and Signal Processing 119 (2019) 346–362 Table 1 The actual leak locations and sizes, their initial parameters and the estimates obtained from IB and ML. Actual leaks xL [m]
L
Initial parameters 4
s [10
600 640 800
1 1.2 0.9
1540
1.4
1600
1
2
m ]
xL [m] 420 430 609 619 742 1195 1205 1348 1358 1515 1570 1580
L
s [10
4
2 2 2 2 2 2 2 2 2 2 2 2
IB Estimates 2
m ]
xL [m] 254.8 401.8 603.6 637.6 801.8 1163 1242 1322 1414 1533 1560 1606
L
ML Estimates 4
s [10
2
m ]
0.001 0.001 0.971 1.117 0.876 0.025 0.045 0.031 0.031 0.844 0.771 0.794
xL [m]
sL [104 m2]
439.2 437.8 601.3 641.3 800.1 1137 1155 1375 1376 1451 1541 1601
1.298 1.301 1.064 1.139 0.900 0.013 0.052 1.140 1.115 0.040 1.425 0.980
using the initial leak locations as in the third column in Table 1, the 12 initial leak locations are all randomly generated, each of them is taken from a uniform distribution with a support whose center is each value in the third column in Table 1 and whose radius is kmin =2 ¼ 129 m. Each initial leak size is also randomly obtained from a uniform distribution with support ½0:5; 3 104 m2. The simulation (from initial parameter generation to leak estimation) is repeated 100 times and the result with highest likelihood is taken as the final estimate, which is shown in Fig. 4(a). The actual leak locations and sizes, their estimates obtained from IB (with highest likelihood among 100 random initial parameters), and the corresponding values of the initial parameters are shown in Table 2. The first, third and fifth leaks are respectively replicated by a single leak
Fig. 4. Leak detection using IB with random initialization. The solid lines, crosses and circles represent the actual, initial and estimated leaks, respectively. Each estimate of leak is connected to its corresponding initial parameter by a dash line. The leak locations are xL1 ¼ 600 m, xL2 ¼ 640 m, xL3 ¼ 800 m, xL4 ¼ 1540 m, xL5 ¼ 1600 m. The leak sizes are sL1 ¼ 1 104 m2 ; sL2 ¼ 1:2 104 m2 ; sL3 ¼ 0:9 104 m2 ; sL4 ¼ 1:4 104 m2 ; sL5 ¼ 1 104 m2. The loglikelihood function of the four cases are: (a) 590; (b) 622; (c) 630; and (d) 2524.
356
X. Wang, M.S. Ghidaoui / Mechanical Systems and Signal Processing 119 (2019) 346–362
Table 2 The actual leak locations and sizes, their estimates obtained from IB with highest likelihood among 100 results with different random initial parameters, and the corresponding values of the initial parameters. Actual leaks L
Estimates L
x [m]
s [10
4
2
m ]
L
x [m]
Initial parameters L
s [10
4
2
m ]
x [m]
sL [104 m2]
L
600
1
600.9
0.918
559.7
2.144
640
1
639.5 634.9 635.9
0.906 0.144 0.153
653.5 392.3 427.9
1.047 0.756 1.431
800
0.9
800.8
0.865
736.5
2.667
1540
1.4
1546 1530 1524
1.034 0.318 0.108
1561 1561 1224
1.659 2.461 1.484
1600
1
1601
0.935
1623
2.505
1217 1357 1551
0.025 0.048 0.046
1154 1340 1365
1.635 0.954 0.667
estimate. The other two leaks at 640 m and 1540 m are respectively represented by three estimates. The sum of the estimated leak sizes corresponding to each of the two leaks is approximately equal to the actual leak size. As previously, since the assumed leak number is much larger than the actual leak number, some estimates have very low size and thus are neglectable. Fig. 4(b–d) shows three other IB estimates among the 100 results. Since they have lower likelihood and, thus, finally are not used, but they are discussed here to illustrate the mechanism of deciding the ultimate estimate. The log-likelihood functions of the four results in Fig. 4 are (a) 1:46 104 , (b) 1:64 104 , (c) 1:68 104 , and (d) 11:87 104 . Correspondingly, the leaks are worse estimated as the likelihood function decreases. In the first three cases, each leak is clearly identified but the error of leak size estimate increases, being (a) 0.0488, (b) 0.0835, and (c) 0.0982, respectively. Here, the error is defined as
kn 5 X L X k L L n n ^s ; e ^s ¼ s n¼1 k¼1
ð27Þ
in which ^sLn ; k ¼ 1; ; kn , are the leak size estimates corresponding to the n-th leak. In the last case, as is shown in Fig. 4(d), the estimate is inaccurate. Note that this is an extreme case whose likelihood function is much lower than the maximum likelihood. Again, this estimate can be easily excluded by comparing the likelihood function. Fig. 5 plots the likelihood function Eq. (8) for IB and ML estimates with various SNR. For each SNR the simulation is repeated 100 times with different random initial parameters as previously; the highest likelihood is picked and shown in k
Fig. 5. Plot of likelihood function Eq. (8) for IB and ML estimates with various SNR. Each likelihood function is selected among 100 simulation results (corresponding to different random initializations) with highest likelihood function.
X. Wang, M.S. Ghidaoui / Mechanical Systems and Signal Processing 119 (2019) 346–362
357
Fig. 5. This figure illustrates that the IB estimate corresponds to higher likelihood function than ML implying that IB more easily avoids stopping at local maxima and performs better than ML (maximizing the likelihood function Eq. (8) is the objective of both ML and IB). Furthermore, Eq. (8) shows that the result of IB is able to better replicate the measurements using the estimated locations and sizes of leaks. Finally, the IB estimate with various leak number is assumed. The determination of leak number is discussed in Section 5; here, only the leak detection result with different assumed leak number is shown. In each simulation, the initial parameters are randomly generated as follows. It can be observed from Fig. 3(a) that the high local maxima of MFP(1) are all in the two intervals ½300; 900 and ½1100; 1700 , thus each initial value of leak location is generated from a uniform distribution with support ½300; 900 [ ½1100; 1700 . Each initial leak size is sampled from a uniform distribution with support ½0:5; 3 104 m2. Fig. 6 displays the IB estimates with respectively 2–7 assumed leaks, each of them is obtained by repeating IB with 100 random initialization and then picking the result with highest maximum likelihood. It is interesting that when the assumed leak number is smaller than the actual leak number (subfigures (a-c)), one estimate tends to represent two or three leaks close to each other. When the assumed and actual leak numbers are identical (subfigures (d)), the estimation is precise that each estimate replicates each leak. In the case that the assumed leak number is higher than the actual leak
Fig. 6. Leak detection using IB with random initialization assuming 2–7 leaks. The solid lines and circles represent the actual and estimated leaks, respectively. The actual leak locations are xL1 ¼ 600 m, xL2 ¼ 640 m, xL3 ¼ 800 m, xL4 ¼ 1540 m, xL5 ¼ 1600 m. The actual leak sizes are sL1 ¼ 1 104 m2 ; sL2 ¼ 1:2 104 m2 ; sL3 ¼ 0:9 104 m2 ; sL4 ¼ 1:4 104 m2 ; sL5 ¼ 1 104 m2.
358
X. Wang, M.S. Ghidaoui / Mechanical Systems and Signal Processing 119 (2019) 346–362
number (subfigures (e,f)), each leak is represented by more than one estimates and/or some of the leak estimates has very low estimated size which is thus ignored. 5. Model selection: Determination of leak number In the previous sections, the multiple leaks are estimated given an assumed leak number N. This section studies the determination of the leak number N. 5.1. AIC and BIC By considering the leak number as a parameter, the log-likelihood function Eq. (8) can be rewritten as
2 N X Ln Ln L L log L N; x ; s jDh ¼ Dh Gn x s : n¼1
ð28Þ
Note that the maximum likelihood principle tends to choose the highest possible leak number. Denote the ML estimates of L ^LðNÞ ¼ ^ x 1; ; ^ xLN and ^sLðNÞ ¼ ^sL1 ; ; ^sLN , then for a higher leak leak locations and sizes under the assumption of N leaks are x ^ LðNþ1Þ ¼ ^ number assumption N þ 1, the estimate x x L1 ; ; ^ xLN ; ^ xLN and ^sLðNþ1Þ ¼ ^sL1 ; ; ^sLN ; 0 remains at the same value of the log-likelihood function, which justifies that a higher leak number does not decrease the log-likelihood function Eq. (28). As a ^LðNþ1Þ and ^sLðNþ1Þ are the optimal estimates (although it is not unique) since the matter of fact, in an noise-free environment, x ^LðNþ1Þ ; ^sLðNþ1Þ jDh log-likelihood function Eq. (28) is 0 (the maximum value). However, when noises exist, because log L N þ 1; x ^ LðNþ1Þ and ^sLðNþ1Þ in terms of having a higher log-likelihood may be found due to is in general non-zero, an estimate better than x one more available free parameter. The goal of this section is to select a reasonable number of leaks which is a trade-off between the goodness of fit (assessed by the likelihood function) and the complexity of the model, i.e., to avoid over-fitting. The solution is to modify the objective function by adding a penalty for high leak number to the log-likelihood function. Here, a commonly used criterion in model selection, known as the Akaike information criterion (AIC) [40], is introduced:
Fig. 7. Log-lilelihood (left) and AIC (right) as a function of leak number N. The sample size is T ¼ 1.
X. Wang, M.S. Ghidaoui / Mechanical Systems and Signal Processing 119 (2019) 346–362
AICðNÞ ¼
J M X X
^L ; ^sL jDhjm ¼ 2NMJ log L N; x ^L ; ^sL jDh ; 2N log L N; x
359
ð29Þ
m¼1 j¼1
in which the first term is the penalty term. The optimal leak number can be decided by minimizing Eq. (29). For a large sample size T (realizations of Dh), another criterion, known as the Bayesian information criterion (BIC) [41], can be used:
BICðN Þ ¼
J M X X
^L ; ^sL jDhjm ¼ NMJ logðT Þ log L N; x ^L ; ^sL jDh : N logðT Þ log L N; x
ð30Þ
m¼1 j¼1
The leak number N estimated from minimizing Eq. (30) has been proven to be asymptotically optimal (in the sense of T ! 1) in [41] that it maximizes the posterior likelihood function (any prior distribution of model parameter is attributed). Finally, a strategy for estimating multiple leaks which combines IB and AIC/BIC, in the case of unknown locations, sizes and number of leaks, is summarized in Algorithm 2. Algorithm 2 Identification of multiple leaks using IB and AIC/BIC in the case of unknown leak number. For N ¼ 0, compute AICð0Þ via Eq. (29) (or, BICð0Þ via Eq. (30)). For N ¼ N þ 1: repeat ^ LðNÞ and ^sLðNÞ ; Run Algorithm 1 (IB) to obtain x Compute AICðN Þ via Eq. (29) (or, BICð0Þ via Eq. (30)). until AICðN Þ > AICðN 1Þ (or, BICðN Þ > BICðN 1Þ). ^LðN1Þ and ^sLðN1Þ as the estimates of leak locations and sizes. Retain N 1 as the estimate of leak number, x 5.2. Numerical examples Here, numerical examples are introduced to illustrate the model selection methods for estimating leak number. First, the numerical example of five leaks introduced in Section 4.2 is considered again. The left plot of Fig. 7 shows the corresponding log-likelihood function with various assumed leak numbers, which increases dramatically from N ¼ 2 to N ¼ 5 (the latter is the actual leak number). After that, the log-likelihood function still increases but very slightly. The right plot of Fig. 7 displays the corresponding plot for AIC whose penalty increases as the leak number increases, such that its the optimal leak number b ¼ 5 is obtained, which is equal to its actual value. estimate N In all the previous numerical examples, only one measurement is used for leak detection, i.e., T ¼ 1. Here, a numerical simulation with a sample size T ¼ 100 is considered; the locations and sizes of leaks, noise distribution and other setups
Fig. 8. Log-lilelihood (left), AIC (middle) and BIC (right) as a function of leak number N. The sample size is T ¼ 100.
360
X. Wang, M.S. Ghidaoui / Mechanical Systems and Signal Processing 119 (2019) 346–362
are same as the previous case. Note that in this case the leak detection is easier than the case T ¼ 1 as more measurements average the noise and thus the parameter estimation is more robust. For this reason, only leak number determination is discussed here. Fig. 8 shows the log-likelihood function with respect to various assumed leak number, along with the results of AIC and BIC. The log-likelihood function has a similar feature as the previous case. Both AIC and BIC can precisely estimate the leak number. 6. Conclusions This paper investigates the problem of detecting multiple leaks in a water-filled pipe with transient waves. A pressure head measurement is assumed to be the summation of leak contributions, each of them is a function of single leak and independent of other leaks. Based on this model, a maximum likelihood (ML) based iterative procedure, known as the iterative beamforming (IB) method, is proposed. IB iteratively estimates the contributions from various leaks to the pressure head measurement and then performs a 1-D inverse method on these estimates. Therefore, this approach largely decreases the computational complexity. Furthermore, the model selection methods for estimating optimal leak number are proposed: Akaike information criterion (AIC) and Bayesian information criterion (BIC) are introduced which penalize the high leak number and thus find a trade-off between the goodness of fit and the complexity of the model. The results obtained on simulated data show the interest of the proposed methods. IB solves a two-leak localization problem which is troublesome for the 1D searching methods like the matched-field processing. Another complex example with five leaks in a pipe shows that IB is efficient, accurate and is a super-resolution scheme. A multi-initialization strategy is also proposed to guarantee the global convergence of the algorithm. Finally, various leak number assumption is discussed; AIC and BIC are both proven to be efficient to decide the leak number in a complicated multiple-leak problem. Acknowledgements This work has been supported by research grants from the Research Grant Council of the Hong Kong SAR, China (Project No. T21-602/15R). The Chinese Estates are acknowledged for their support of this research. Appendix A. A brief introduction of the model in Eq. (1) The discharge and pressure head due to a fluid transient are represented by q and h. Given the discharge q xU and U at the upstream node xU , the quantities at xm (xm > xLN > > xL1 ) are obtained from the transfer matrix head h x method [50]:
q ð xm Þ
¼M
hðxm Þ
NL
xm xLN
3 ! 1 Q L0n q xU 1 Ln 4@ 2ðH0 zLn Þ AM NL xLn xLn1 5 ; h xU n¼N 0 1
1 Y
20
in which xL0 ¼ xU ,
cosh ðlxÞ
NL
M ð xÞ ¼
Z sinh ðlxÞ
1Z sinh ðlxÞ
ðA:1Þ
! ðA:2Þ
cosh ðlxÞ
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi is the field matrix, Z ¼ la2 =ð ixgAÞ is the characteristic impedance, l ¼ a1 x2 þ igAxR is the propagation function, a is the wave speed, x is the angular frequency, A is the area of pipeline, and R is the frictional resistance term. The right hand side of Eq. (A.1) can be approximated by a linear form plus a infinitesimal term (when leak sizes tend to 0) that can be neglected in the case of small leaks [31]:
q ð xm Þ hðxm Þ
¼
MNL ðxm Þ þ
N X sLn M SL Ln ; xLn ; xm n¼1
!
! q xU U þ o max sLn ; n¼1;;N h x
ðA:3Þ
in which
! vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Z sinh lxLn cosh l xm xLn cosh lxLn cosh l xm xLn u g M SL Ln ; xLn ; xm ¼ u t Z cosh lxLn sinh l xm xLn Z 2 sinh lxLn sinh l xm xLn 2 HL0n zLn
ðA:4Þ
is a matrix related to the location xLn of the n-th leak (zLn and HL0n are decided by xLn ), but is independent of the size sLn of n-th leak and the other leaks. Note that, physically, the infinitesimal term in Eq. (A.3) implies the wave interaction (multiple reflections) between leaks, while the first term in Eq. (A.3) corresponds to the wave received by xm after a single reflection at each leak. By neglecting the infinitesimal term in Eq. (A.3), the pressure head at xm (m ¼ 1; ; M) near the downstream for the frequency xj (j ¼ 1; ; J) is
X. Wang, M.S. Ghidaoui / Mechanical Systems and Signal Processing 119 (2019) 346–362 N X NL h x j ; xm ¼ h x j ; xm þ sLn G xj ; xLn ; xm ;
361
ðA:5Þ
n¼1 NL
xj ; xm and G xj ; xLn ; xm follow Eqs. (3) and (5). NL Let Dhjm h xj ; xm h xj ; xm denote the head difference between the head measurement in the presence of leaks and the theoretical head that does not include the leak terms at the frequency xj and at the measurement station xm , and denote wherein h
> Dh ¼ Dh11 ; ; DhJ1 ; ; Dh1M ; ; DhJM ;
ðA:6Þ
then we have the following equation:
Dh ¼ G xL sL ;
ðA:7Þ
in which G is a JM N-dimensional matrix (cf. Eqs. (4) and (5)). Finally, Eq. (1) is obtained by adding a noise term. Appendix B. Upstream discharge estimation in Eq. (7) Here, the discharge q xU at the upstream node xU is estimated by adding a measurement station near the upstream boundary at xU þ ( > 0) [51,52]. Assuming that there is no leak between xU and xU þ , the discharge and head at xU and xU þ have the expression:
! ! q xU þ q xU NL ¼M ðÞ U h xU h x þ ¼
coshðlÞ
Z sinhðlÞ
1 Z
sinhðlÞ
coshðlÞ
!
! q xU þ : h xU þ
ðB:1Þ
By measuring pressure head h xU þ at xU þ ( > 0) and applying a boundary condition h xU of head at xU , the discharge U q x can be estimated by [31]:
coshðlÞh xU h xU þ ^ xU ¼ q : Z sinhðlÞ Appendix C. Derivation of Eq. (19) The following property of Gaussian random variable [49] is first reminded, which is used to derive Eq. (19). Let X and Y be n-dimensional Gaussian random vectors with expectation mX and mY , and with covariance matrix RXX and RYY . Let RXY ¼ CovðX; YÞ and RYX ¼ CovðY; XÞ, the conditional probability density function is Gaussian:
1 f YjX ðjXÞ N jmY þ RYX R1 XX ðX mX Þ; RYY RYX RXX RXY :
Therefore, the conditional distribution of cn N
ðC:1Þ
lcn ; Rcn given Dh N lDh ; RDh follows Gaussian distribution:
N1 1 1 ¼ N cn Gn xL sLn þ cn Dh N cn lcn þ Rcn R1 Dh G x L s L ; I : Dh Dh lDh ; Rcn Rcn RDh Rcs N N2
ðC:2Þ
References [1] [2] [3] [4] [5] [6] [7]
[8] [9] [10]
A.F. Colombo, P. Lee, B.W. Karney, A selective literature review of transient-based leak detection methods, J. Hydro-Environ. Res. 2 (4) (2009) 212–227. J.A. Liggett, L.-C. Chen, Inverse transient analysis in pipe networks, J. Hydraul. Eng. 120 (8) (1994) 934–955. B. Brunone, Transient test-based technique for leak detection in outfall pipes, J. Water Resour. Plann. Manage. 125 (5) (1999) 302–306. J.P. Vítkovsky`, A.R. Simpson, M.F. Lambert, Leak detection and calibration using transients and genetic algorithms, J. Water Resour. Plann. Manage. 126 (4) (2000) 262–265. W. Mpesha, S.L. Gassman, M.H. Chaudhry, Leak detection in pipes by frequency response method, J. Hydraul. Eng. 127 (2) (2001) 134–147. X.-J. Wang, M.F. Lambert, A.R. Simpson, J.A. Liggett, J.P. Vítkovsky`, Leak detection in pipelines using the damping of fluid transients, J. Hydraul. Eng. 128 (7) (2002) 697–711. P.J. Lee, J.P. Vitkovsky, M.F. Lambert, A.R. Simpson, J.A. Liggett. Frequency response coding for the location of leaks in single pipeline systems, in: The Int. Conference on Pumps, Electromechanical Devices and Systems Applied to Urban Water Management, IAHR and IHRA, Valencia, Spain, April 22–25, 2003. P.J. Lee, J.P. Vítkovsky`, M.F. Lambert, A.R. Simpson, J.A. Liggett, Frequency domain analysis for detecting pipeline leaks, J. Hydraul. Eng. 131 (7) (2005) 596–604. P.J. Lee, J.P. Vítkovsky`, M.F. Lambert, A.R. Simpson, J.A. Liggett, Leak location using the pattern of the frequency response diagram in pipelines: a numerical study, J. Sound Vib. 284 (3) (2005) 1051–1073. P.J. Lee, M.F. Lambert, A.R. Simpson, J.P. Vítkovsky`, J. Liggett, Experimental verification of the frequency response method for pipeline leak detection, J. Hydraul. Res. 44 (5) (2006) 693–707.
362
X. Wang, M.S. Ghidaoui / Mechanical Systems and Signal Processing 119 (2019) 346–362
[11] H.-F. Duan, P.J. Lee, M.S. Ghidaoui, Y.-K. Tung, Essential system response information for transient-based leak detection methods, J. Hydraul. Res. 48 (5) (2010) 650–657. [12] H.-F. Duan, P.J. Lee, M.S. Ghidaoui, Y.-K. Tung, Leak detection in complex series pipelines by using the system frequency response method, J. Hydraul. Res. 49 (2) (2011) 213–221. [13] H.-F. Duan, P.J. Lee, M.S. Ghidaoui, Y.-K. Tung, System response function–based leak detection in viscoelastic pipelines, J. Hydraul. Eng. 138 (2) (2011) 143–153. [14] W. Nixon, M.S. Ghidaoui, Numerical sensitivity study of unsteady friction in simple systems with external flows, J. Hydraul. Eng. 133 (7) (2007) 736– 749. [15] M.F. Ghazali, S.B.M. Beck, J.D. Shucksmith, J.B. Boxall, W.J. Staszewski, Comparative study of instantaneous frequency based methods for leak detection in pipeline networks, Mech. Syst. Sig. Process. 29 (2012) 187–200. [16] A.C. Zecchin, L.B. White, M.F. Lambert, A.R. Simpson, Parameter identification of fluid line networks by frequency-domain maximum likelihood estimation, Mech. Syst. Sig. Process. 37 (1) (2013) 370–387. [17] S. Meniconi, B. Brunone, M. Ferrante, C. Capponi, C.A. Carrettini, C. Chiesa, D. Segalini, E.A. Lanfranchi, Anomaly pre-localization in distribution– transmission mains by pump trip: preliminary field tests in the milan pipe system, J. Hydroinf. 17 (3) (2015) 377–389. [18] S. Meniconi, B. Brunone, M. Ferrante, C. Massari, Potential of transient tests to diagnose real supply pipe systems: what can be done with a single extemporary test, J. Water Resour. Plann. Manage. 137 (2) (2010) 238–241. [19] B. Brunone, M. Ferrante, Detecting leaks in pressurised pipes by means of transients, J. Hydraul. Res. 39 (5) (2001) 539–547. [20] D. Covas, H. Ramos, N. Graham, C. Maksimovic, Application of hydraulic transients for leak detection in water supply systems, Water Sci. Technol.: Water Supply 4 (5–6) (2005) 365–374. [21] S.T.N. Nguyen, J. Gong, M.F. Lambert, A.C. Zecchin, A.R. Simpson, Least squares deconvolution for leak detection with a pseudo random binary sequence excitation, Mech. Syst. Sig. Process. 99 (2018) 846–858. [22] M. Ferrante, B. Brunone, S. Meniconi, Wavelets for the analysis of transient pressure signals for leak detection, J. Hydraul. Eng. 133 (11) (2007) 1274–1282. [23] J.C. Liou, Pipeline leak detection by impulse response extraction, J. Fluids Eng. 120 (4) (1998) 833–838. [24] M. Ferrante, B. Brunone, Pipe system diagnosis and leak detection by unsteady-state tests. 1. harmonic analysis, Adv. Water Resour. 26 (1) (2003) 95–105. [25] D. Covas, H. Ramos, A.B. De Almeida, Standing wave difference method for leak detection in pipeline systems, J. Hydraul. Eng. 131 (12) (2005) 1106–1116. [26] A.M. Sattar, M.H. Chaudhry, Leak detection in pipelines by frequency response method, J. Hydraul. Res. 46 (EI1) (2008) 138–151. [27] M. Taghvaei, S.B.M. Beck, J.B. Boxall, Leak detection in pipes using induced water hammer pulses and cepstrum analysis, Int. J. Comadem. 13 (1) (2010) 19. [28] M.L. Stephens, Transient Response Analysis for Fault Detection and Pipeline Wall Condition Assessment in Field Water Transmission and Distribution Pipelines and Networks PhD Thesis, University of Adelaide, 2008. [29] D. Covas, H. Ramos, Case studies of leak detection and location in water pipe systems by inverse transient analysis, J. Water Resour. Plann. Manage. 136 (2) (2010) 248–257. [30] X. Wang, M.S. Ghidaoui, Pipeline leak detection using the matched-field processing method, J. Hydraul. Eng. 144 (6) (2018) 04018030. [31] X. Wang, M.S. Ghidaoui, Identification of multiple leaks in pipeline: linearized model, maximum likelihood, and super-resolution localization, Mech. Syst. Sig. Process. 107 (2018) 529–548. [32] B. Zhou, A. Liu, X. Wang, Y. She, V. Lau, Compressive sensing-based multiple-leak identification for smart water supply systems, IEEE IoT J. 5 (2) (2018) 1228–1241. [33] X. Wang, D.P. Palomar, L. Zhao, M.S. Ghidaoui, R.D. Murch, Spectral-based methods for pipeline leakage detection, J. Hydraul. Eng. (2019), https://doi. org/10.1061/(ASCE)HY.1943–7900.0001572, In Press. [34] X. Wang, B. Quost, J.-D. Chazot, J. Antoni, Estimation of multiple sound sources with data and model uncertainties using the EM and evidential EM algorithms, Mech. Syst. Sig. Process. 66–67 (2016) 159–177. [35] X. Wang, B. Quost, J.-D. Chazot, J. Antoni, Iterative beamforming for identification of multiple broadband sound sources, J. Sound Vib. 365 (2016) 260–275. [36] X. Wang, S. Khazaie, L. Margheri, P. Sagaut, Shallow water sound source localization using the iterative beamforming method in an image framework, J. Sound Vib. 395 (2017) 354–370. [37] A.P. Dempster, N.M. Laird, D.B. Rubin, Maximum likelihood from incomplete data via the EM algorithm, J. R. Stat. Soc., Ser. B (Method.) 39 (1) (1977) 1–38. [38] J.C.F. Wu, On the convergence properties of the EM algorithm, Ann. Stat. 11 (1983) 95–103. [39] M. Feder, E. Weinstein, Parameter estimation of superimposed signals using the EM algorithm, IEEE Trans. Acoust., Speech Sig. Process. 36 (4) (1988) 477–489. [40] H. Akaike, A new look at the statistical identification model, IEEE Trans. Autom. Control 19 (1974) 716–723. [41] G. Schwarz, Estimating the dimension of a model, Ann. Stat. 6 (2) (1978) 461–464. [42] K.P. Burnham, D.R. Anderson, Model Selection and Multimodel Inference: A Practical Information-theoretic Approach, Springer Science & Business Media, 2003. [43] C.M. Hurvich, C.-L. Tsai. Regression and time series model selection in small samples. Biometrika, pages 297–307, 1989. [44] X. Wang, Z. Zhang, S. Li, Set-valued and interval-valued stationary time series, J. Multivariate Anal. 145 (2016) 208–223. [45] K. Aho, D. Derryberry, T. Peterson, Model selection for ecologists: the worldviews of AIC and BIC, Ecology 95 (3) (2014) 631–636. [46] J.P. Vítkovsky`, M. F Lambert, A.R. Simpson, J.A. Liggett, Experimental observation and analysis of inverse transients for pipeline leak detection, J. Water Resour. Plann. Manage. 133 (6) (2007) 519–530. [47] M.S. Ghidaoui, M. Zhao, D.A. McInnis, D.H. Axworthy, A review of water hammer theory and practice, Appl. Mech. Rev. 58 (1) (2005) 49–76. [48] M.S. Ghidaoui, On the fundamental equations of water hammer, Urban Water J. 1 (2) (2004) 71–83. [49] I.B. Rhodes, A tutorial introduction to estimation and filtering, IEEE Trans. Automatic Control AC-16 (6) (1971) 688–706. [50] M.H. Chaudhry, Applied Hydraulic Transients, Third Edition., Springer, 2014. [51] A. Kashima, P.J. Lee, R. Nokes, Numerical errors in discharge measurements using the KDP method, J. Hydraul. Res. 50 (1) (2012) 98–104. [52] A. Kashima, P.J. Lee, M.S. Ghidaoui, M. Davidson, Experimental verification of the kinetic differential pressure method for flow measurements, J. Hydraul. Res. 51 (6) (2013) 634–644.