## Abstract

The aim of the present research is to quantify the maximum feasible accuracy (location and size) for the transient-based blockage detection in a water supply pipeline. Owing to the randomness of transient measurements, a performance bound for the extended blockage detection exists which estimated parameters cannot exceed. The Cramér–Rao lower bound (CRLB) theorem is utilized to compute the lower bound variance of noise-induced estimation errors. It gives the minimum mean square error of any estimator according to information obtained from measurements and quantified by Fisher information. The Fisher information matrix is computed using direct differentiation of the compatibility equations obtained by the method of characteristics. The influence of relevant physical parameters including valve closure time, measurement time length and noise level on the best possible localization of blockage is investigated. The connection between the signal bandwidth, noise level and the performance limit is quantified for a typical case study. The results demonstrate trade-off between the size and the location/length of blockage estimates subject to different maneuver times, roughly offering half the wave speed times maneuver duration as the resolution limit.

## INTRODUCTION

Pipeline systems are one of the best and safest ways to supply fluids such as water, oil, etc. To ensure that pipeline systems are efficient and safe, appropriate design and condition assessment is essential. The water supply systems are always subject to different kinds of faults such as leaks, discrete and extended blockages, or corrosions (Massari *et al.* 2013). Blockages in pipes increase operational costs by wastage of energy, reduction in the pipe carrying capacity, increased potential for contamination, and development of leaks due to generation of overpressure in pipes (Lee *et al.* 2015). Similarly, corrosion of pipe walls (reduction in pipe thickness) is a potential cause for the emergence of leaks. In this regard, punctual detection and repair of these faults are crucial since this can prevent serious future problems like wasting a great amount of fresh water or failure of the whole pipeline.

There are various methods for fault detection and condition assessment in pipeline systems, e.g., physical inspection, acoustic method, steady-state hydraulic analysis, and device-based method, which are usually too expensive or time-consuming, and thus less efficient (Rizzo 2010; Duan *et al.* 2012). Recently, transient-based defect detection methods (TBDDM) such as inverse transient analysis (ITA) have received a great deal of attention because they are efficient, nonintrusive, and economic (Colombo *et al.* 2009; Ferrante *et al.* 2014). ITA, which requires the injection of a signal into the conduit and the measurement of its response at specified locations, has been proven as a promising approach to detect blockages, leaks, and assessment of pipe wall condition. The theory of wave propagation along with the knowledge of anomaly signatures on the reflected waves enable development of such methods. The TBDDM for extended deteriorations of the pipe wall has been developed in both time and frequency domains, such as studies by Tuck *et al.* (2013), Meniconi *et al.* (2013), Duan *et al.* (2012, 2013, 2017), Massari *et al.* (2014, 2015), Louati *et al.* (2018), Jing *et al.* (2018) for blockage, and research works of Stephens *et al.* (2013) and Gong *et al.* (2013, 2018) for corrosion; a combination of blockage and corrosion detection was investigated by Gong *et al.* (2014). Field tests by Stephens *et al.* (2013) and laboratory experiments by Meniconi *et al.* (2013) showed that deteriorations have significant impacts on the system responses.

The application results from these studies revealed that the developed TBDDM is considerably more accurate in hypothetical (ideal) numerical tests than in realistic experimental or even laboratory cases. This is explained by the inevitable uncertainties of the input data (e.g., fluid and pipe system characteristics or measured signals) and the adopted mathematical model (i.e., imprecise simulation of the forward transient waves). To address the uncertainties, Duan (2015b) studied the possible inaccuracies in the ITA results due to noise or miscalibration of measuring devices. In another study, Duan (2015a) investigated uncertainty propagation to fluid pressure owing to several random input data and, conversely, the leak detection uncertainty due to random measurements or model inputs. Massari *et al.* (2014, 2015) developed a stochastic successive linear estimator that statistically provides the best unbiased estimate of diameter distribution due to partial blockages and quantifies the uncertainty associated with these estimates.

It is vital to conduct a systematic investigation of the influences of various uncertainties on the pipe system to handle applicability and accuracy of the TBDDM. The input parameters to the model as well as collected data of transient tests can be classified as random variables, which in turn, cause random results from the ITA (Duan 2015b; Sattar & El-Beltagy 2017), thus invoking the need to perform statistical models. The estimation theory provides an effective lower bound estimation error as a function of measured signal characteristics (bandwidth, power, time length) for the unknown parameters. Under regularity conditions, the Cramér–Rao lower bound (CRLB) theorem gives a lower bound for mean square error (MSE) of any estimator (Kay 1993). The CRLB matrix is derived from the inverse of the Fisher information matrix which itself is derived from the joint likelihood function of the estimated pipe areas of extended blockage candidates given the probability density function of measurements.

This paper aims at investigating the uncertainty quantification of the extended blockage parameters due to measurement noise in TBDDM. In extended blockages, a uniform deterioration is assumed for the entire pipe element (whose boundaries usually coincide with the adopted computational nodes), whereas in the discrete case (e.g., in Meniconi *et al.* 2014) the blockage property is lumped at a point so that it is modeled as a nodal component. A simple numerical test system with single extended blockage is scrutinized in this study. The computed CRLB of the extended blockage parameters is incorporated in quantification of minimum possible standard deviation of inner area and length of blockage estimates. The effect of valve maneuver time (signal bandwidth) and time length of measurements on the accuracy of localization and blockage size estimates in a noisy environment is studied. The standard deviations of blockage candidates are depicted as error bars on the estimated parameters to represent their precision from the ITA simulation.

## WATER HAMMER IN PIPES WITH EXTENDED BLOCKAGES

*z*= distance along pipe,

*Q*= instantaneous discharge,

*H*= instantaneous piezometric head,

*t*= time,

*a*= pressure wave speed,

*g*= gravitational acceleration, and

*f*= Darcy–Weisbach friction factor considered to be quasi-steady herein.

*A*and

*D*are cross-sectional area and inner diameter of pipe, respectively. They alter due to partial blockages along the pipeline, which can eventually lead to some reflections in the transient responses. The anomaly signatures on the reflected waves allow for the identification of blockages. Measurement noise causes uncertain reflections on the signal and hence leads to inaccurate detection of blockage parameters (Duan 2015a).

Equations (1) and (2) are solved by the method of characteristics (MOC), considering appropriate initial and boundary conditions (Wylie & Streeter 1993). It is assumed that the wave speed remains constant for unblocked and blocked sections (Duan *et al.* 2012, Duan 2015b; Jing *et al.* 2018).

## ESTIMATION MODEL

The estimation model for the extended blockage detection using random measurements is described herein. To solve the ITA problem, maximum likelihood estimator (MLE) (Press *et al.* 1992) is expected to provide the best unknowns if uncertainties are normally distributed.

**h**denote a measured pressure vector at

_{m}*M*time steps and

**h**define the corresponding true model output and

**z**,

_{B}**L**, and

_{B}**A**vectors denote beginning locations, length, and uniform cross-sectional areas of some extended blockages in the pipe system

_{B}*,*respectively. Due to the presence of measurement noise (or system uncertainties), the pressure head responses of a pipe system to a deterministic excitation can be modeled as a random process,

**H**, and the measured pressures,

_{m}**h**, can be viewed as a specific realization of

_{m}**H**. The following relation indicates the mathematical representation of the estimation problem in the presence of measurement noise

_{m}**n**:

**h**(

**z**,

_{B}**L**,

_{B}**A**) in Equation (3) represent model outputs which are governed by continuity and momentum equations of the fluid as discussed in the previous section. It is assumed that each element in the MOC method can be a blockage, with its length depending on the number of elements (

_{B}*N*) considered along the pipe length. Based on this assumption, the blockage detection problem reduces to determination of cross-sectional area of blockage candidates in the MOC method. The synthesized pressure heads then only depend on cross-sectional areas of pipe as estimation parameters.in which

_{E}**h**(

**A**) is theoretical forward model results and

*N*is the number of sections along the pipe to be considered as potential extended blockages (

_{P}*N*=

_{P}*N*). The vector of estimated areas

_{E}**A**implies that each pipe section is intact or deteriorated.

### Maximum likelihood estimation

**h**is a real

_{m}*M*dimensional random vector whose entries are assumed to follow Gaussian distribution with the mean vector and the covariance matrix . The joint likelihood function of cross-sectional area of pipe sections or the joint probability of measurements is:The MLE for the model parameters is estimated from the random variable

**H**(and its statistical properties). MLE corresponds to maximizing the likelihood function of the parameters:The asymptotic property of the likelihood function represents an unbeatable MSE of estimated cross-sectional area of pipe sections in the ITA model. The noise on measurements, or uncertainty in model input parameters are responsible for the uncertain blockage identification.

_{m}### Cramér–Rao lower bound

For any estimator, if the estimation process is repeated for many realizations of the uncertain variable **H _{m}**, then the statistical properties of the identified parameters provide the errors associated with the parameter estimations (Monte Carlo method, Hastings 1970). Nevertheless, the MLE has some general statistical asymptotic properties that enable us to evaluate the parameter estimation uncertainty without numerous estimations. The covariance matrix of a MLE can be asymptotically evaluated using the CRLB (Kay 1993). Also, since the MSE of MLE tends to CRLB for a large sample size of measurements, this bound can be used to anticipate the best possible performance of any estimator. In other words, it offers the lower bound variance for all estimators.

The simplest formulation of CRLB is for unbiased estimators (Vaseghi 2008) assumed in the current estimation model. The term ‘unbiasedness’ means that the expected value of the estimated cross-sectional area of pipe sections is equal to the true value of the parameters (for either intact or blocked sections). In the unblocked (intact) elements, the expected value of the estimated cross-sectional area is equal to the intact pipe area *A*_{0}. Accordingly, , and imply that section *i* is intact, corroded, or blocked, respectively.

*M*-dimensional identity matrix. Therefore, the log-likelihood function becomes:The Cramér–Rao bound states that the inverse of the Fisher information is a lower bound on the variance of unbiased pipe-area estimators of . A proof found by Van Trees (1968) and Frieden (2004) is provided herein.

_{,}therefore, its partial derivative with respect to must also be zero. By the product rule, this partial derivative is equal to:

*F*) is determined by the vector product of two time histories and . The vector (or ) means the derivative of the computed head vector at the measurement location with respect to the cross-sectional area of

_{i,j}*i*th (or

*j*th) blockage candidate. A numerical algorithm to compute the elements of these vectors is proposed in the next section. Observed in the last expression of Equation (14), the decisive experimental (simulation) parameters namely maneuver time of value (

*T*), time length of signal (

_{c}*T*), and the noise variance to determine signal to noise ratio (SNR) is captured. The computed head is obviously a function of model inputs including time closure and closure pattern of valve, time-length of analysis, pipeline and flow properties and blockage parameters. In other words, as seen in the last expression of Equation (14), lower bound of the covariance matrix of estimated parameters is independent of the measured heads and is only a function of the statistical properties of

_{T}**H**(i.e., computed head and variance, herein).

_{m}### Performance bound definition

This measure of performance has units of pipe cross-sectional area and represents the minimum average error of size estimation. Such a performance measure has also been implemented by Robinson & Milanfar (2006).

## COMPUTATIONAL APPROACH

In this section, the numerical scheme to compute the elements of the vector in Equation (14) and consequently the CRLB of estimated blockages is explained. The derivative at *m*th time step with respect to the *i*th model parameter is approximated by direct differentiation method (DDM) for all defect (intact) pipe areas as potential blockage parameters. The proposed DDM is based on exact differentiation of state variables with respect to the model parameters which is carried out on compatibility equations of MOC. The steps for numerical implementation are illustrated now.

First, the sensitivities at steady state corresponding to each point are computed using continuity and Darcy–Weisbach equations through the following steps.

- The derivative of flow rate at node (see Figure 1) with respect to the pipe-sections area
*A*is estimated; note that the number of pipe sections (blocked-pipe candidates),_{i}*N*, is considered equal to the number of computational points (MOC elements),_{P}*N*:_{E}

- The sensitivities of the state variables (
*Q*and_{k}*h*) at the new time step are found by first taking a derivative of simultaneous solution of the compatibility relations for flow rate , and then a derivative of a characteristics equation, e.g., Equation (21):_{k} The sensitivity of the head at the reservoir is zero since pressure head of the reservoir is a specified quantity and remains constant during transients. The negative characteristic line allows for the computation of flow-rate sensitivity at the reservoir point. Likewise, the positive compatibility equation in conjunction with the orifice relation is employed to find sensitivities at the valve point.

A similar approach was previously used by Nash & Karney (1999) to estimate sensitivities with respect to friction factors of pipe sections and later on by Keramat *et al.* (2019) for gradients with respect to leak areas in the ITA. They both reported that the DDM is accurate and computationally efficient when adopted to MOC equations. The proposed computational approach has also the great advantage of versatility, allowing for the application of models which account for advanced models of viscoelasticity (Covas *et al.* 2005; Keramat *et al.* 2013; Zanganeh *et al.* 2015), turbulence, and fluid–structure interaction (Keramat & Tijsseling 2012; Keramat *et al.* 2012), non-Newtonian fluids (Majd *et al.* 2016), and various types of boundary conditions (Ahmadi & Keramat 2010; Meniconi *et al.* 2014, 2017; Ferreira *et al.* 2018).

## RESULTS AND DISCUSSION

The main aim of this paper is to provide performance bounds of any blockage detection method when one measurement signal is used to recover blockage characteristics. Due to noise in the measured signal, error in estimation of blockage unknowns is inevitable. The performance limit is estimated in terms of mean square error of size estimation, which is quantified by Cramér–Rao (CR) inequalities as described earlier.

The performance bounds under various values of input data are estimated. This is carried out by computing CRLB when different time closure of valve (*T _{c}*) and length of potential blockages (

*d*) are tested. The use of CRLB to find minimum standard deviation of parameter estimates is also studied. Subsequently, it is illustrated that the lower bound errors of potential blockage sizes provide considerable insight into spatial resolution of blockage localization.

The hypothetical pipe system under consideration has length *L*_{0} = 1,000 m, wave speed *a*_{0} = 1,000 m/s, inner diameter, *D*_{0} = 0.3 m, and Reynolds number of 1 × 10^{5}. For the investigations, an extended blockage with beginning location of *z _{B}* = 400 m and length of

*L*= 200 m, and area of

_{B}*A*= 0.25

_{B}*A*

_{0}is considered (the cross-sectional area of pipe along the extended blockage is

*A*

*=*0.75

*A*

_{0}).

### Standard deviation in cross-sectional area detection

The CRLB of the covariance matrix of the cross-sectional pipe-area estimates is computed based on Equation (14) and the scheme described in the section ‘Computational approach’. The estimated standard deviations are often shown along with the results of the ITA as complementary information, which visualize the amount of variation of estimated parameter.

*N*=

_{P}*N*= 20. Consequently, the length of each blocked pipe is

_{E}*d*

*=*50 m which is equivalent to the numerical reach of the MOC scheme. The computed standard deviations correspond to system characteristics indicated above each figure part. The ratio of the signal to noise (SNR) energy is used to characterize the signal strength in the localization exercise. Figure 2(a) shows standard deviations in an experimental case where the blockage location can be identified. In fact, the following criteria can be defined for a successful blockage localization:where subscript

*i*stands for intact section and

*j*indicates a blocked section. Consequently, the key parameters of the experiment being maneuver time of value (

*T*), time length of signal (

_{c}*T*), and noise level (SNR) allow for favorable localization as observed in Figure 2(a). However, with a higher closure time, the accuracy of estimated parameters (pipe-section areas) reduces so that it is highly probable to perform wrong detections. This is observed in Figure 2(b) when closure time increased to 0.1

_{T}*T*,

*T*= fundamental water hammer period, which is equivalent to the reduction of the signal bandwidth. As a result, the estimated uncertainty of the unknown parameters of the inverse problem (the areas of the pipe sections), will certainly lead to misdiagnosing the blocked section. Unfortunately, it is usually hard in reality to conduct experiments with short maneuver time. In return, a higher time length of measurement as observed in Figure 2(c) or data acquisition with more precision (higher SNR) like that of Figure 2(d) can be incorporated to achieve similar successful localization to Figure 2(a).

In Figure 2(a), 0.5*aT _{c}* = 20 m is the super-resolution limit (Gong

*et al.*2013; Lee

*et al.*2015) so that successful localization within

*d*

*=*50 m accuracy is predictable as this is not a super-resolution exercise (

*d*> 0.5

*aT*). However, area estimation is subject to the indicated amount of inaccuracy due to noisy measurements. Considering 0.5

_{c}*aT*= 200 m for Figure 2(b)–2(d), the accuracy of 50 m is super-resolution, among which case (b) illustrates failure, whereas cases (c) and (d) depict successful localizations.

_{c}It is worthwhile to note that due to the fundamental property of the extended blockage detection model expressed in the section ‘Estimation model’, the location and length of blockages are lumped at computational pipe sections. As a result, the length of pipe sections for blockage candidates indicates the desired resolution or accuracy of localization. Since the estimation model is unbiased, higher areas than *A*_{0} are possible from the ITA model to represent the extended-corroded pipe wall.

The computed standard deviations can be interpreted and anticipated if an inverse analysis is successful in a system with some given blockages. An example is provided in Figure 3 where the standard deviations for two different cases of blockage location/length estimations are compared. The time closure, sample time length, and SNR are the same in this figure, but the resolutions for blockage candidates are different. In Figure 2(a) (also repeated in Figure 3(a) for the simple comparison), 20 parameters for blockage detection are used, which in turn, seek the accuracy of localization within 50 m. In Figure 3(b), the number of parameters are raised to 100 which corresponds to the localization accuracy of 10 m. As seen, the higher the number of parameters in the inverse problem, which corresponds to higher accuracy in detecting the location/length of blockages, results in larger values of standard deviations; this trend can be interpreted as follows.

Localizations under the condition that the time length of the blockage reflection in the pressure signal is smaller than the valve maneuver-induced wave front (time closure) is termed as super-resolution: or . Accordingly, the blockage reflections are affected by the valve maneuver so that blockage characteristics are not accurately identified even in noise-free conditions. In the application of ITA using unbiased estimators (which allows for the identification of both blocked and corroded sections), this especially seems to occur owing to the opposing reflections of blockage and corrosion, which in turn, enables fitting to the collected pressure signal using both the types of defects. In other words, similar band-limited pressure signals can be generated at the downstream valve from various combinations of neighboring blockage or corrosions. This eventually leads to higher size localizations, which means higher standard deviations, and is more probable for blockage candidates of relatively shorter length . For immediate closures (infinite bandwidth), CRLB is independent from the length of potential blockages and does not increase with resolution because the wave-front length is zero.

The significance of signal bandwidth or maneuver time of the downstream valve in the localization is further elaborated upon by Gong *et al.* (2013), Tuck *et al.* (2013), and Jing *et al.* (2018). Regarding the duration of the maneuver time, possible methods for generating a fast transient can be addressed: the spark wave generator described by Gong *et al.* (2018); the portable pressure wave maker (PPWM) proposed by Brunone *et al.* (2008); the side discharge valve method used by Stephens *et al.* (2013) and Meniconi *et al.* (2011); and the device incorporated by Taghvaei *et al.* (2010).

### Maximum precision in localization (super-resolution limit)

The influence of signal bandwidth on the localization performance is further investigated in this section. The performance bound defined by Equation (16) is used towards this aim. Figure 4 reveals the change of mean standard deviation of cross-sectional area estimates with length of extended blockage candidates. The localization of blockages with accuracy *d* < 0.5*aT _{c}* is prone to higher uncertainty in size detection. The significant rise of CRLB for quite short lengths of extended blockage candidates reveals that no confident size detection can be made within that order of location accuracy (e.g., see the line corresponding to

*T*= 0.06

_{c}*T*, when

*d*< 125 m). The trend of change of lower-bound standard deviation is mostly dictated by valve closure time (

*T*) which corresponds to the length of probing wave fronts (

_{c}*aT*)

_{c}*.*More specifically, the travel time of a wave in a blockage candidate should be compared with the closure time of the valve. For localizations with higher resolution than , the estimated standard deviation starts to rise. In this zone, it is still likely to detect blockage but with a larger standard deviation which corresponds to more error in area detection. According to these curves, the quantity 0.5

*aT*can be suggested as the efficient mesh size for the ITA-based blockage estimation; a higher mesh size than that is less accurate in location of blockage, while lower mesh size will be less accurate in detecting cross-sectional area.

_{c}*j-*th element with pipe area

*A*, and there is no blockage at the remaining element

_{j}*i*(

*A*=

_{i}*A*

_{0},

*i*≠

*j*), then Equation (28) can be used as a criterion for successful localization of this blocked section. Based on the numerical manifestations in Figures 2 and 3, one can assume that all elements have similar standard deviations so that , thus Equation (28) reduces to:

It means that only blockages of size higher than twice the blockage-detection error can be localized. Equivalently, to detect the blocked section, the lower bound error at all sections should be smaller than 0.5*A _{B}*. As a result, is used as the threshold for localization, which for the blockage case under consideration is indicated by the dashed line in Figure 4. Given each experiment property (e.g., valve closure, SNR, measurement length, etc.), the corresponding best possible resolution in localization can be found using the intersection of the relevant line with the dashed line; different lines for different

*T*are depicted in Figure 4 (additional lines for different SNR or measurement length are not provided due to paper length restrictions).

_{c}## SUMMARY AND CONCLUSIONS

Noise on a measured signal can blur the desired reflections from extended blockages and reduce the amount of information that can be achieved from blockages using waves. It is therefore of great interest to quantify the maximum achievable information or lower bound of error in the blockage detection and localization. According to the CRLB theorem, the covariance matrix of parameter estimates of maximum likelihood approaches CRLB, so that it makes a powerful tool to quantify the lower bound error. An efficient numerical algorithm based on MOC in a multi-blockage pipe system is proposed in this research to compute pressure head sensitivities with respect to the cross-sectional area of pipe elements. The gradients' computations are carried out using direct differentiation of characteristics' equations so that a computationally fast and accurate method for the numerical evaluation of CRLB is proposed.

Some case studies to reveal the existence of a limit in the blockage localization corresponding to the applied set of experiment characteristics including time closure, time length, and SNR along with the mesh size of the ITA are provided. The numerical results for a typical reservoir–pipe–valve system considering an extended blockage in the middle of the pipe are discussed. The main contributions of CRLB in terms of lower bound error of cross-sectional area of pipe section estimates are addressed. The results demonstrate that the computed standard deviations provide additional information regarding the ITA accuracy.

The provided figures reveal a sharp rise of inaccuracy with increased time closure, and demonstrate the trade-off between the resolution in the localization of blockage and its size estimation. The study manifests practical implications regarding optimal design of time closure, measurement time duration, and noise level in conjunction with optimal mesh size for ITA when a specified accuracy is required. To quantify the success of the blockage localization, a criterion in terms of the pipe-area estimates and the lower bound error of estimates is proposed (e.g., if the error is as high as the estimated blockage size, no useful information from the ITA can be drawn). It is observed that the error of size detection rises with reducing length of potential blockage sections. A significant correlation between the blockage candidates as the unknowns can be concluded from the results since the random pressure signal can be reconstructed by several blockages/corrosions of different size. The change of the mean error of estimates versus the length of blockage candidates enables prediction of the best possible resolution in localization for a given set of physical parameters. This further allows for setting the appropriate number of MOC elements while using ITA.