Multi-stage parameter-constraining inverse transient analysis for pipeline condition assessment

Model calibration is the process of determining the model parameters to obtain a model representation of the system of interest that satisfies a certain criterion, i.e. a goodness of fit metric between the system measurements and the simulated model (Savic et al. 2009). inverse transient analysis (ITA) (Liggett & Chen 1994), one of the popular transient-based methods (Meniconi et al. 2013; Gong et al. 2014; Massari et al. 2014; Capponi et al. 2017; Duan 2017), is a process of pipe system parameter estimation using measured transient pressure data and a hydraulic transient model. The inverse analysis is achieved by adjusting pipeline system parameters (such as leak locations and sizes, friction factors, water demands or hydraulic transient wave speeds) so that the calibrated model's predicted response matches the measured response (typically the temporal measurements of pressure). In the case of hydraulic transient wave speeds, the calibrated parameter is closely related to the pipeline condition (Stephens et al. 2013), as the wave speed is related to the wall thickness and internal diameter of the pipe, which change as the pipe wall corrodes. As such, a well calibrated model provides valuable information about the pipeline condition, which serves as guidance for water authorities to assess infrastructure condition, and maintain assets strategically.

Despite the fact that ITA has been widely studied, the scale of the problems that have been considered (i.e. the number of system parameters) is relatively small. ITA was first proposed by Liggett & Chen (1994) for leakage and friction factor calibration. Within their numerical study, 15 parameters, including 11 friction factors and four leak sizes, were calibrated for a small 11-pipe network. Jung & Karney (2008) considered a numerical case study, where the wave speed, diameter and friction factor of 11 pipes, four possible leakages and one water demand were calibrated simultaneously (38 parameters overall), on the same network as Liggett & Chen (1994). The average error of the parameter estimates was, in general, greater than 30% and the authors claimed the error associated with estimated wave speeds tends to compromise the calibration of other parameters. Shamloo & Haghighi (2010) also conducted ITA on the same network to detect potential leaks. The absolute errors of the 17 parameter estimates (six leak sizes and 11 friction factors) were all below 5% by optimizing the transient generation and the choice of measurement site.

Concerning the laboratory verification of ITA, Vitkovsky et al. (2001, 2007) calibrated the leak size for a laboratory system consisting of a 37.2 m long copper pipe with an internal diameter of 22.1 mm. A total of 15 leak candidate locations were considered where the measured pressure trace contained multiple periods of the transient wave in the pipe system. Covas & Ramos (2010) applied ITA to a 272 m high-density polyethylene pipe with an internal diameter of 50.6 mm and a 1.3 km polyethylene pipe with an internal diameter of 70 mm. Multiple models were considered, where the number of calibrated parameters was less than 30 for each model. Soares et al. (2011) applied ITA to a system composed of looped PVC pipes with a total length of 203.20 m. The potential leak locations in the inverse model was 7. It is worth pointing out that both Covas & Ramos (2010) and Soares et al. (2011) adopted a strategy to minimize the number of leak candidates, starting with a set of leak candidates sparsely distributed throughout the system, and gradually reducing the set of candidates to those around the potential leak locations obtained in the earlier steps.

ITA in the frequency domain has also attracted attention. Kim (2008) developed an address-oriented impedance matrix method for the generic calibration of several parameters, such as the location and quantity of leakage, the friction factor, and the wave speed simultaneously. The method was numerically verified by a 12-pipe network. Kim (2008) also adopted a similar multi-stage strategy to Covas & Ramos (2010), starting the calibration with the initial search-space, which was the accumulated length of all of the pipeline elements, and then restraining the search-space to a candidate pipeline element. Zecchin et al. (2013) developed a frequency-domain maximum likelihood estimation method to identify parameters of fluid line networks, in which different types of parameters including wave speed, diameter, and length of pipe were calibrated simultaneously. The proposed method was applied on an 11-pipe network using turbulent steady friction and turbulent unsteady friction pipe models, which resulted in parameter estimation problems of dimension 33 and 55, respectively. Compared with time domain ITA, there are fewer laboratory or field verification works for frequency domain applications of ITA. Laboratory verification of the frequency domain ITA is found in the research by Kim et al. (2014), where location and discharge of a single leak, together with other parameters such as wave speed, were calibrated. Capponi et al. (2017) also successfully located and sized one leak, one blockage and one branch with high accuracy in a controlled laboratory environment.

As outlined above, the complexity of the inverse problem (i.e. the dimensionality of the parameter estimation space) is limited in most previous research. An exception is the work by Stephens et al. (2013), where ITA was applied to a mild steel cement-lined transmission mains in rural South Australia. In this work, the wave speeds of 390 reaches, along a 6 km section of pipe, were calibrated to infer internal pipe wall condition. The pattern of estimated wave speeds by ITA were generally consistent with ultrasonic measurements taken along the pipe, although this was not always the case. However, within this application, Stephens et al. (2013) claimed that it was almost certain the optimal solution was not found, due to a limitation on the number of objective function evaluations used by the optimizer (given the high computational cost of overall millions evaluations) and model error. Additionally, the authors also believed that, given the high dimensionality of the problem, it was very likely that non-unique solutions existed, meaning that there were very many local optimum solutions with similar objective function values. However, this issue was not explored in detail in their work.

Non-uniqueness of solutions for ITA applications has been observed in other work, and is a known problem for inverse problems more generally (Yeh 1986). Within the study by Jung & Karney (2008), three different ITA scenarios were analyzed on an 11-pipe network, where each had a different number of measurement sites. The authors concluded that fewer measurement sites led to a better calibration error and faster convergence to the calibrated parameter values, but to a less accurate calibration result. Kapelan et al. (2004) pointed out that some poorly estimated parameters were due to an inadequate quantity of observed information. That is, Kapelan et al. (2004) outlined that the parameter calibration problems could not be resolved by simply using an alternative search algorithm. To improve the ill-posed natural of the inverse problem, a framework for the effective incorporation of prior information into ITA for pipe networks was developed. However, these two approaches to enhance the identifiability are not always feasible in field tests, where access points to install pressure transducers are limited, and prior information is not always available. Instead of increasing the number of measurement sites or incorporating prior knowledge of a system, another strategy to enhance parameter identifiability of ITA is to reduce the complexity of the inverse problem by limiting the number of decisions (Kim 2008; Covas & Ramos 2010) or limiting the range of the search-space.

The interval for a parameter's search-space is typically the minimum and maximum feasible values of the parameter to be calibrated, and they are typically the same for every reach. However, in the case of wave speed calibration in ITA, typically sections of pipe in a normal state (e.g. sections in a generally good condition) have wave speed values within a small band about a nominal value, and it is only anomalous sections (e.g. sections with severe deterioration) that possess significant uncertainty and variability in the wave speed values they can take. Therefore, if the search-space can be properly adjusted for each reach to account for the reduced uncertainty for normal sections of pipe, better parameter estimates can be expected within the limited search-space.

To enhance parameter identifiability of ITA on a complex transmission main, a multi-stage parameter-constraining ITA is proposed in this paper. This paper is organized as follows: to highlight lack of identifiability achieved by the conventional ITA on a complex transmission mains system, a numerical example with a highly variable pipe wall condition is considered (represented by hundreds of reaches, each with different wave speeds). It is confirmed within this example that the conventional ITA is not able to uniquely identify the correct wave speed values, despite achieving parameter estimations that provide a very good fit of the predicted and measured pressure responses. To tackle this problem, a multi-stage parameter-constraining ITA is proposed. The proposed algorithm works by iteratively generating sets of parameter estimates and coincidentally modifying the search-space interval for the subsequent estimates. Based on a statistical analysis of the set of estimates, each reach is classified as normal (parameter values close to the nominal values) or anomalous (parameter values significantly different from the nominal values). Different search-space intervals are allocated to each reach based on the outcome of the classification. Parameter estimation within the properly narrowed search-space will lead to concentration of parameter estimates and elimination of outliers, thus identifiability of parameters by ITA is enhanced. This algorithm is followed by a numerical verification and discussion about the sensitivity to the relevant parameters. Finally, the method is verified by a field case study.

The problem of lack of identifiability for applications involving the conventional ITA method are explored within this section with a numerical example application on a reservoir-pipeline-valve (RPV) system. Within this section, the numerical example illustrates that when applying the conventional ITA to a relatively complex (yet realistic) inverse problem, the identifiability of the wave speed parameters is relatively poor. This is evidenced by the fact that independent inverse runs (using a stochastic or evolutionary optimization approach with different random number sequences) yielded a large variability in the optimal parameter estimates, despite the high similarity in the estimates' associated calibration error. The data of the system are described in the next sub-section, followed by a discussion, where it is demonstrated that the conventional ITA suffers from the problem of lack of identifiability.

A typical RPV system (Figure 1) is used to simulate a transmission main throughout this paper. The wave speed distribution along the pipe is depicted in Figure 2. Wave speeds of 115 normal reaches are in the range of 950–1,000 m/s. A total of 13 reaches (nine sections in Figure 2), whose wave speeds are either less than 950 m/s or greater than 1,000 m/s, are used to represent different wall thickness pipe sections due to different pipe classes or wall condition. A transient wave is generated by instantaneously shutting the side discharge valve at node G when t = 0.2 s. The pressure responses at locations M1, M2 and G are obtained by numerical simulation using the method of characteristics (MOC) in which the time step ⁠. The length of each reach was determined by to satisfy Courant condition. This avoids interpolation errors in the forward model when synthetically generating the ‘measured’ pressure traces. Only a duration of 1.2 s after the transient generation is considered in the ITA to avoid reflections from the boundaries of the RPV system.

The ITA used in this paper is the efficient version proposed by Zhang et al. (2018), in which the head based method of characteristics (HBMOC) is employed as the transient simulator. Particle swarm optimization (PSO) (Poli et al. 2007) was employed as the optimization algorithm – to account for the stochastic nature of PSO, the ITA is run multiple times (with different starting random number seeds) to provide a statistical characterization of the algorithm's performance. The objective function is defined as the squared difference between the measured responses and the predicted responses by the calibrated model summed over all measurement stations. A time step of 0.01 s was used in the ITA. Considering the 1.2 s pressure trace used in ITA and the location of measurement sites, an inverse model with 128 reaches is required, to accommodate all reaches contributing to the 1.2 s pressure trace. As such, the number of decision variables is 128, and for PSO, the number of particles is 150 with 10,000 iterations (resulting in 1,500,000 evaluations overall), which were selected so as to allow for convergence of the algorithm's search. The search-space was set to be [800 m/s, 1,050 m/s], and the decision variables were treated as continuous within this range by the PSO.

Figure 3 presents a series of boxplots for the wave speeds estimated in each reach by ten independent runs. It can be seen that the interquartile range (IQR) of the wave speeds for each reach are generally large, which means that the wave speed estimates of the independent runs are spread diversely throughout the wide search-space (i.e. starting from different random number seeds, the PSO ended up with significantly different solutions). This high variability amongst the parameter estimates limits a clear interpretation of the results, and the ability to use these estimates for the detection of the anomalous sections.

To understand why the conventional ITA parameter estimates in Figure 3 yield such high variance, the envelopes of ten predicted pressure traces and the match between measured pressure trace and predicted pressure trace obtained by the best solution from the ten runs (ranked in terms of objective function) are given in Figure 4. Despite the existence of outliers in the estimates (O1–O4 in Figure 3), the best wave speed estimates solution still provides an excellent match between its predicted pressure traces and the measured pressure trace. The closely wrapped envelopes of predicted pressure traces indicate that reasonable match was achieved by every set of wave speed estimates. The fact that such an excellent model fit can be achieved for such vastly different parameter estimates means that the objective function is not sufficiently informative to allow for the discrimination between different parameter sets. These findings support the observations from Kapelan et al. (2004) concerning estimation issues arising from data information limitations for certain ITA applications.

In the previous example, it has been illustrated that, for the conventional ITA, independent parameter estimates are widely spread throughout the search-space. Despite the reasonable match between the modelled and measured pressure traces, distinct outliers exist in the estimates, suggesting that multiple local optimums exist in the relatively large search-space. This limitation of the conventional ITA to correctly identify parameter estimates for large-scale problems has motivated the iterative parameter-constraining approach of the proposed method.

The motivation for the proposed method is outlined in the following. Transmission mains are typically installed with a similar design and are deteriorated in a similar way, unless the cumulative effect of soil or groundwater corrosion and external loaded forces increase degradation in some sections. Besides these severely deteriorated sections, undocumented class changes with different wave speeds and pipe sections replaced by plastic pipe section with lower wave speeds are also anomalous sections in transmission mains. Compared with the majority sections which are in similar conditions, these anomalous sections are localized to a limited number of short discrete sections along the pipeline. The wave speed parameter values for the normal sections are confined to within a small interval about a nominal value, as there is much less uncertainty about the values that the normal sections can take. In contrast, the anomalous sections contain wave speed values that lie within a much broader interval, as there is significant uncertainty about the values that the anomalous sections can take. Therefore, it is expected that the effectiveness and accuracy of the ITA parameter estimation process would be increased if adaptive search-space parameter intervals were used to allow for a narrowing of the interval for sections that were recognized as normal, and maintaining a broader interval for sections that were identified as anomalous. Based on this idea, a multi-stage parameter-constraining ITA (ITA_MP) method is proposed. In this approach, the search-space is iteratively limited, guiding the solutions closer to the true parameter values.

The flowchart of the ITA_MP algorithm is given in Figure 5. The algorithm comprises multiple iterative stages, and each stage involves the sub-stage of new solution generation, reach classification, search-space updating for the next stage and the termination criteria check. Details of the four sub-stages are presented in the following sub-sections.

The solution matrix provides a data set of wave speed estimates for each reach. A fundamental assumption of the proposed process is that the previous applications of ITA have been effective in their optimization strategy, and that the majority of estimates are, at least, within the neighborhood of the global optimum wave speed solution a*. Given this assumption, the classification of whether a reach is considered normal or anomalous is based on analyzing the data set for each reach to determine if the estimates are significantly different from the nominal properties of the entire set of wave speed estimates. That is, given the set of wave speed estimates for reach k, ⁠, the classification process aims to determine the likelihood that the statistical properties of this set are consistent with the statistical properties of the wave speed estimates for a normal section of pipe as characterized by the data set (the data set comprised of all elements of excluding column k, represented as the subscript ⁠). If the statistical properties of the kth reach are consistent with the properties of a normal section of pipe, then the reach is classified as normal; if not, then it is classified as anomalous.

Given the fact that the properties of the estimates of a normal section of pipe are unknown a priori, the data set comprised of the wave speed estimates for the rest of the pipeline are used to characterize a normal section. This approach to the characterization of the estimates for a normal section is justified on the basis that a fraction of reaches will be in an anomalous state (sections with severely deteriorated conditions, undocumented or replaced pipe sections with different material), and the majority reaches will be in the normal state (generally good condition or similarly deteriorated condition). This means that within this larger data set (all elements of ⁠) the influence of estimates for the anomalous sections on the statistics of this set will be small compared to the influence of the normal sections. This is particularly the case when order statistics are used to characterize the set, as these are insensitive to outliers in a data set.

An appropriate statistical test to undertake the comparison of the kth reach data set with is the Mann–Whitney U test (Hollander et al. 2013). This test provides a non-parametric approach to compare two data sets to determine the likelihood that a random sample of one data set will be less than (or greater than) a random sample of the other data set. The null hypothesis of this test is that both sets are distributed about the same median value, meaning it is equally likely that a sample from one data set will be less than (or greater than) a sample from the other. A rejection of the null hypothesis means it is more likely that both sets are not distributed about the same median value. In the classification context of this paper, a reach is classified as normal if the null hypothesis holds, and is classified as anomalous if the null hypothesis is rejected. The acceptance or rejection of the null hypothesis is dependent on the comparison of the p-value for kth reach, ⁠, returned by Mann–Whitney U, and the specified significance level ⁠. The null hypothesis is accepted if ⁠, and is rejected if ⁠.

Different search-space intervals are allocated for different reaches in the next stage, depending on the set or that the reach belongs to. The detection of an anomalous section is the focus of the parameter estimation; as a result, a reach classified as anomalous will retain the original wide wave speed search interval so that the estimation strategy is still able to search within the full search-space.

is the collection of multiple estimates for all reaches which are classified as normal, so the statistics of can be used to represent the parametric range that the wave speed values of the normal reaches are within. Considering the likely existence of outliers in the estimates in ⁠, p and q percentile of are used to determine the new search-space bounds. Compared with the original search-space ⁠, the new search-space is a narrower interval. The choice of percentile rank p and q will be discussed in the section on numerical verification.

The algorithm is terminated from iterating in Stage 2 when the M multiple solutions generated within the updated search-space at Stage n all have a larger objective function than the previous best objective function from iteration Stage n–1. The algorithm aims to minimize the objective function, so if it fails to decrease the objective function, and thus fails to improve the match between the calibrated models and the measured data, there is sufficient reason to believe that a failure to find a better solution is attributed to an improper updated search-space, so the algorithm is terminated.

A numerical experiment was conducted to verify the proposed ITA_MP. The numerical model is the same as in Figure 1. The number of reaches is N = 128; and the number of independent ITA runs for different starting random number seeds is M = 10, where PSO was used as the optimization algorithm.

Nine cases were considered given different values of m, p and q (see Table 1). Case 3 (m = 15, p = 5 and q = 95) is discussed in detail in the section ‘Results and discussion’. The other cases are discussed in the section ‘Choice of m, p and q values’.

For Case 3, the objective functions of all independent ITA runs are given in Figure 6. From this figure it is seen that the ITA_MP terminated after Stage 3, making the total number of independent ITA runs equal to 30. The ten ITA runs in Stage 1 ended up with a diverse set of solutions due to the random nature of the optimization algorithm and the relatively large search-space, which is consistent with the diverse wave speed estimates in Figure 3. A significant reduction in the objective function value can be seen in Stage 2, after the first search-space update. Solutions in Stage 2 show more consistency in their objective function values within narrowed search space.

Figure 7 presents the boxplot of the wave speed estimates in Stage 2. Compared with the wave speed estimates in Stage 1 in Figure 3 (the result of Stage 1 is the results of the conventional ITA, since the conventional ITA is terminated after multiple independent runs is Stage 1), the IQR of wave speed of each reach become significantly less. For the proposed approach, the wave speed estimates become increasingly concentrated around the true wave speeds of the model, making the interpretation of the results and detection of anomalous sections easier.

From Figure 7, it also can be seen that the outliers existing in the Stage 1 estimates (for example, O1–O8 in Figure 3) were not present in the estimates of Stage 2. Elimination of the outliers was due to the section classification, and consequent search-space updating, in ITA_MP. Interestingly, these outliers from the single run in Stage 1 were identified as normal sections in Stage 2, because the estimates of the other independent runs yielded statistics consistent with the normal value. As a result, in the next stage, the search-space was properly narrowed and a better solution was found.

The enhanced identifiability of wave speed estimates also resulted in an improvement in the match of the measured and predicted pressure traces, for which the traces at the generator location are given in Figure 8. All independent wave speed estimates in Stage 2 provided an almost indistinguishable match between measured and predicted pressure traces, indicated by shrunk envelopes of predicted pressure traces in Stage 2, compared with those in Stage 1 in Figure 8.

The sensitivity of the ITA_MP results to the percentile rank parameters m, p and q is discussed in this section. Three parameters m, p and q are percentile ranks in Equations (2)–(4) to determine the significance level in the Mann–Whitney U test and search-space bounds for next stage ITA runs in Stage 2.

As explained in the section ‘numerical experiment’, among 128 reaches of the inverse model, 13 reaches (or a total of nine sections) are anomalous, meaning that the proportion of anomalous reaches is 10.2%. It has been explained that the percentile rank m represents the proportion of reaches which will be identified as anomalous reaches within ITA_MP. Case 1 (m = 5, p = 5, q = 95) is used as an example to illustrate the impact of a smaller m on the results. In Case 1, the ITA_MP algorithm terminated after Stage 2. Figure 9 represents the single best wave speeds of Stage 2.

It can be seen in Figure 9 that the anomalous sections D1, D4, D5, D6 and D9 were successfully identified. However, anomalous sections D2, D3, D7, which are closer to the normal wave speeds (compared with those of D1, D4, D5, D6 and D9), were identified as normal sections because a smaller m was used and only 5% of reaches were allowed to be identified as anomalous reaches. The wave speeds of these three sections reach their maximum or minimum values (the boundaries of the search-space for normal sections).

With regard to Cases 2–4 where the m% values are all greater than 10.2% (the true proportion of anomalous reaches), the boxplots of estimates in Cases 2–4 combined (that is 30 estimates for each reach) are given in Figure 10. Comparison between the boxplots in Case 3 alone (Figure 7) and Cases 2–4 combined demonstrates that as long as m% is greater than the proportion of anomalous sections, the choice of the percentile rank m has little sensitivity to the wave speed calibration. In conclusion, it is recommended that m is set to be greater than the assumed proportion of anomalous sections. In field applications, m can be chosen based on the prior information, if available, or by a trial and error strategy.

A comparison of Cases 5–9 with Case 3 is made to explore the impact of variations in the percentile rank parameters p and q. Figure 11 shows the search-space interval for normal sections of different stages for the different cases. In Case 5 (p = 0, q = 100), the search-space bound and for the normal sections are determined as the minimum and maximum of normal section matrix ⁠. These Case 5 bounds are almost the same as the original search-space, and as a result, ITA_MP terminated after Stage 2 as no further improvement was achievable. This proves that limiting the search-space is the driving force behind the improvement of the wave speed estimation in ITA_MP.

As the percentile rank p becomes greater and q becomes smaller, the updated search-space becomes narrower, and ITA_MP was observed to terminate earlier. This is because for the increasingly narrow search-space, more normal reaches have wave speed values that violate the search-space bounds, yielding a reduced calibration performance. This side effect of ITA_MP is that some normal reaches are not allowed to take their true values, if their wave speeds are out of the narrowed search-space range. However, from the perspective of the detection of anomalous sections, this side effect can be considered unimportant as the accurate calibration of the normal sections does not affect the anomalous section detection. Overall, the results from Cases 5–9 demonstrate that the estimated wave speeds, especially those of anomalous sections, are not sensitive to different values for the percentile ranks p and q.

A field case study on the Morgan–Whyalla transmission pipeline (MTP) in South Australia was conducted to verify the multi-stage parameter-constraining ITA (the reader is referred to Stephens et al. (2013) and Gong et al. (2015) for more details on the field case study). The pipeline is made of mild steel with cement mortar lining (MSCL), is in good condition with intact cement mortar lining, and has an internal diameter of 727.5 mm. The configuration of the pipeline under investigation (from chainage 14,900 to 18,900 m) is given in Figure 12. For the series of field experiments, the transient waves were generated by shutting scour valve SC24 abruptly. Three transducers were installed at air valves AVFP43 and AVFP44, and scour valve SC24 to record the transient pressure response.

The measured pressure traces at AVFR43, SC24 and AVFR44, in which the underlying long-period pressure oscillations were removed by a low-pass filter, are given in Figures 16–18 (where they are compared with predicted pressure trace of the calibrated model). A time step of 0.0135 s and a duration of 4.4 s of the measured pressure traces were used in ITA_MP to cover the 4 km pipe section of interest. This setup of the inverse problem resulted in the number of reaches in the inverse model being N = 383. The number of independent ITA runs of different random number seeds in each sub-stage was set to be M = 10. In each independent ITA run, for the adopted optimization algorithm PSO, the number of particles and the number of maximum iterations were set to be 400 and 3,000, respectively. The percentile ranks were chosen as m = 20, p = 5 and q = 95.

Objective function values of all independent ITA runs are given in Figure 13. Similar to Figure 6 (the objective function of the numerical study), a significant reduction is found when the search-space was updated for the first time. The ten objective function values became concentrated gradually up to Stage 4, when multiple ITA runs failed to find a better solution within the third updated search-space.

Figure 14(a) and 14(b) represent boxplots of the wave speed estimates in Stage 1 and Stage 3. Extensive ultrasonic wall thickness measurements were manually undertaken using an ultrasonic thickness measurement instrument at 5 m intervals along the MTP between chainages 14,900 and 18,900 m (Stephens et al. 2013). Measurements were taken at eight points around the circumference of the pipe at each location. All ultrasonic wall thickness measurements are given in Figure 14(c).

Four known thicker-walled sections of pipeline A, B, C and D, confirmed by the ultrasonic measurements, are captured by wave speed estimates in both Stage 1 and Stage 3. Other unknown anomalous sections (for example, section F) are also identfied by ultrasonic measurements and wave speed estimates in both Stage 1 and 3. However, the IQR of Stage 3 is significantly less than that of Stage 1. Less variation in wave speed estimates enable less false positives identification and more accurate abnomal section detection. Section E, which is depicted in Figure 15, is used to demonstrate the benefits the ITA_MP brings.

According to the ultrasonic measurements, a small fraction of the section from chainage 16,159.2 to 16,329.2 m (Xc in Figure 15(c)) is very likely to be in the anomalous state. It is identified by ITA results in both Stage 1 and Stage 2; they are Xa in Figure 15(a) and Xb in Figure 15(b), respectively. However, several sections in the dashed box in Figure 15(a) are likely to be interpreted as the anomalous sections, although they should not be. After ITA_MP is applied, wave speed estimates of those sections are close to their nominal values, thus they are identified as normal sections, which is consistent with ultrasonic measurements. The true anomalous section Xb is highlighted as a result. To summarize, ITA_MP provided a clear interpretation of the parameter estimates, and improved the ability to use these estimates for the detection of the anomalous sections.

However, a few true anomalous sections (for example, section G in Figure 14) were missed by ITA_MP. This is likely due to the coarse spatial resolution (roughly 12 m, the 4 km pipe section was discretized into more than 300 reaches) used in the ITA_MP. The ultrasonic measurement is taken at 5 m intervals, so the length of several anomalous sections was less than the spatial resolution used in ITA_MP. As a result, those anomalous sections were missed by ITA_MP.

The comparison between the measured pressure trace and the predicted pressure trace obtained by the best estimated wave speeds (ranked in terms of objective function value) in Stage 3 is also examined. Figures 16–18 show the comparison between the measured pressure trace and the predicted pressure trace obtained by the best estimated wave speeds from Stage 3. The high frequency pressure fluctuations at SC24 are due to high frequency pressure oscillations, and resultant cavitation, within the pipework of the generator. The inverse model is not able to describe this phenomenon, so they did not appear in the predicted pressure traces in Figure 17. The pressures at ACFP43 and ACFP44 also exhibit significant oscillations in the first 0.5 s period after the transient has arrived. The oscillations were likely due to the oscillations generated at SC24, so the pressures of the first 0.5 s have been excluded from the objective function calculations. It can be seen that a reasonable match between the measured and predicted pressure trace is achieved for the time period 0.5 s after the wave front.

Overall, given the difficulty of the field problem that ITA_MP attempted, the anomalous sections identified by Stage 3 are generally consistent with the sections identified by ultrasonic measurements. The identifiability is improved compared with the conventional ITA (Stage 1 in the ITA_MP), which is indicated by less variation in the wave speed estimates and a decreased objective function (meaning a better match between the predicted and measured pressure traces).

Parameter identifiability is a problem commonly encountered in large-scale ITA applications for determining pipeline condition. To tackle this problem, multi-stage parameter-constraining ITA (ITA_MP), a new method for pipeline condition assessment, has been developed and is presented in this paper. ITA_MP iteratively generates sets of parameter estimates and modifies the search-space interval for the generation of subsequent estimates based on identifying a reach as either in a normal or anomalous state. The parameter estimates for the normal reaches become increasingly concentrated and outliers in the parameter estimates are eliminated, thus identifiability of parameter values is enhanced.

This method has been verified by a numerical study in this paper. The results show that the wave speeds estimated by the conventional ITA (equivalent to Stage 1 of ITA_MP) possess a number of anomalous wave speed estimates, despite the methods seemingly being a reasonable match to the numerical data. After Stage 2, the estimated wave speeds are much closer to the true wave speeds, and a large proportion of anomalous estimates are eliminated.

A field case study on a section of the Morgan transmission pipeline in South Australia has been conducted to verify the capability of ITA_MP to assess pipeline condition assessment by identifying anomalous sections. Wave speed estimates by ITA_MP have less variance and higher identifiability than the conventional ITA. The results show a reasonable consistency between anomalous sections detection by ITA_MP and available ultrasonic measurements, demonstrating that ITA_MP can be used as a tool for pipeline condition assessment in the field.

The research presented in this paper has been supported by the Australia Research Council through the Discovery Project Grant DP140100994 and the Linkage Project Grant LP130100567. The first author thanks the Chinese Scholarship Council and the University of Adelaide for providing a joint postgraduate scholarship.