Statistical performance analysis of transient-based extended blockage detection in a water supply pipeline

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.

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).

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.

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₀. Accordingly, ⁠, and imply that section i is intact, corroded, or blocked, respectively.

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).

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 mth time step with respect to the ith 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.

1. The sensitivities of head and flow rate at the reservoir point with respect to all elements is set to zero: 

   

   (17)
2. The derivative of flow rate at node (see Figure 1) with respect to the pipe-sections area A_i is estimated; note that the number of pipe sections (blocked-pipe candidates), N_P, is considered equal to the number of computational points (MOC elements), N_E: 

   

   (18)
3. The derivative of head at node with respect to the ith pipe area parameter is estimated through: 

   

   (19)

   where and if k=i and otherwise.

1. The sensitivities of the characteristic coefficients are found using the calculated sensitivities at the previous time step: 

   

   (22)

    

   

   (23)

    

   

   (24)

    

   

   (25)
2. The sensitivities of the state variables (Q_k and h_k) 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): 

   

   (26)

    

   

   (27)

3. 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).

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₀ = 1,000 m, wave speed a₀ = 1,000 m/s, inner diameter, D₀ = 0.3 m, and Reynolds number of 1 × 10⁵. For the investigations, an extended blockage with beginning location of z_B = 400 m and length of L_B = 200 m, and area of A_B = 0.25 A₀ is considered (the cross-sectional area of pipe along the extended blockage is A= 0.75A₀).

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.

In Figure 2(a), 0.5aT_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.5aT_c). However, area estimation is subject to the indicated amount of inaccuracy due to noisy measurements. Considering 0.5aT_c = 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.

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₀ 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).

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.5aT_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_c = 0.06 T, when d < 125 m). The trend of change of lower-bound standard deviation is mostly dictated by valve closure time (T_c) which corresponds to the length of probing wave fronts (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.5aT_c 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.

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.5A_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_c are depicted in Figure 4 (additional lines for different SNR or measurement length are not provided due to paper length restrictions).

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.