Time-Series-Based Leakage Detection Using Multiple Pressure Sensors in Water Distribution Systems

Leak detection is nowadays an important task for water utilities as leakages in water distribution systems (WDS) increase economic costs significantly and create water resource shortages. Monitoring data such as pressure and flow rate of WDS fluctuate with time. Diagnosis based on time series monitoring data is thought to be more convincing than one-time point data. In this paper, a threshold selection method for the correlation coefficient based on time series data is proposed based on leak scenario falsification, to explore the advantages of data interpretation based on time series for leak detection. The approach utilizes temporal varying correlation between data from multiple pressure sensors, updates the threshold values over time, and scans multiple times for a scanning time window. The effect of scanning time window length on threshold selection is also tested. The performance of the proposed method is tested on a real, full-scale water distribution network using synthetic data, considering the uncertainty of demand and leak flow rates, sensor noise, and so forth. The case study shows that the scanning time window length of 3–6 achieves better performance; the potential of the method for leak detection performance improvement is confirmed, though affected by many factors such as modeling and measurement uncertainties.


Introduction
The urban water supply network system can be seen as a mobile data carrier which transmits the information of pipe flow rate, nodal pressure, and water quality under the condition of satisfying energy balance, pressure balance, and water quality balance. Being the city's primary infrastructure, the system provides security for urban development and residential life, which requires a guarantee of the system's safe, efficient, and economical operation. As all water distribution systems (WDS) have a certain degree of leakage, leakage management of pipe networks becomes one of the most concerning problems for water supply companies.
Over the past few decades, several researchers have conducted extensive research on leak detection in water distribution systems; a number of methods have been proposed to identify pipe bursts/leaks [1]. Leak detection techniques can be divided into two categories: external and internal [2]. External methods include radar detection, hydrophone, and acoustic logging, and are time-consuming and can only be searched locally. Internal methods are divided into four main categories: (1) transient-based methods, (2) calibration-based approaches, (3) data-driven methods, and (4) model-based techniques.
Transient-based methods determine the location of leaks through time domain or frequency domain analysis of pressure signal observations [3]. Time domain methods were introduced by Jönsson and Larson [4] and have since been derived in a number of techniques [5][6][7][8]. It consists of time domain The monitoring data of the WDS is a typical time series, and there is a certain spatial correlation between the monitoring devices distributed at different locations of the WDS. After the leakage occurs, multiple monitoring sensors will respond at the same time, producing synchronous pressure changes and showing spatial correlation; besides, the changes will often last for a period of time, and there is autocorrelation in time.
This paper proposes a threshold selection method for leak detection in water distribution systems. It is characterized by (1) fully utilizing temporal varying correlation between the data from multiple pressure sensors; (2) updating threshold values over time; and (3) scanning multiple times for a scanning time window (STW). The main objective is to explore the advantages of data mining based on time series data for leak detection, and to reduce the rate of missing detection and the rate of false alarms for the leakage events effectively. The state of the WDS is changing all the time, and the real-time update of the threshold is beneficial to reduce the influence of the model and measurement uncertainty. Section 2 describes the principle of this methodology in detail. In Section 3, a case study is presented, which uses the water supply network of J City of China to carry out the research on the performance of the leak detection method. Finally, a discussion of the relevant results is given in Section 4. Figure 1 shows the framework of leak detection and localization methodology based on monitoring data and hydraulic model. In this framework, the supervisory control and data acquisition (SCADA) system supplies raw pressure and flow rate data for the prediction and detection process. In the prediction process, it is assumed that such a tool is available to forecast the spatial and temporal demand distribution at the specified time window. Water demand forecasting is based on a combination of pattern recognition and time series models by using the historical water demand database. Then, a simulation is run to obtain the leakage residual database by adding the nodal leakage to normal water demand during the time window. Comparing the measurements from SCADA to the leakage residual database, the leak detection and identification algorithm is used to determine whether the leak occurs. If the detection alarm is not triggered, the network status is classified as a non-fault state, and the prediction hydraulic model of the next specified time window is updated using the current measurement information. Conversely, if the detection alarm is triggered, the next step is to locate the leakage position. The prediction process can be found in many studies [32][33][34]. This paper focuses on the detection process and has the following assumptions: there is only one leak that appears at one time and leaks occurs at a node in the water distribution system model; pressure sensors in the water network are placed in the selected nodes and are working well, and the measured pressures are synthetically generated from the simulated pressure value by adding random noise; the error of the prediction model is considered by adding the noise of the nodal demand to the real nodal demand; sudden special events that may produce significant relevant demand variations are not considered. To achieve the time series threshold estimation process, a few concepts should be recalled here:

Methodology
• Detection Probability (DP): the proportion of leak events detected in the total number of natural random leak events, and the rate of non-detection (false negative rate) is described as '1-DP'.

•
Rate of False Alarm (RF): the proportion of false alarms in the total number of natural random events that occur without a leak.
The authors define: • Scanning Time Window (STW): it is shown in Figure 2. It is a time container covered by continuous time steps. Its length is the number of the covered time step, namely, the scanning time window length (STWL).

•
Scanning Process: for STWL of k, one complete scan process performs a total of k scans. The data is numbered in order of reception time. The k scans use the data of time step at t = 1, t = {1, 2}..., t = {1, 2..., k}. The scanning process is shown in Figure 2. The RF of the kth scan in a STW, RF k , is assigned to calculate the threshold C k . A scanning process for an STW requires a corresponding threshold set {C k }, and the number of elements in the collection is equal to STWL. • Cumulative Rate of False alarm (CRF): for an STW containing k time steps (if STWL = k), CRF is the proportion of the total number of false alarms in the total number of natural random events without leakage when a scanning process is completed in that STW. the proportion of the total number of false alarms in the total number of natural random events without leakage when a scanning process is completed in that STW.

Leak Scenario Falsification
The core idea of the error-domain model falsification method is to falsify model instances (parameter sets). Then, the maximal plausible errors are determined by combining modeling and measurement uncertainties. There is model uncertainty (umodel) in the numerical model, as well as measurement uncertainty in the actual measured value. As shown in Equation (1), the predicted value (Pi) plus modeling error (umodel) corresponds to the real quantity (R), which is equal to the measured value (y) plus measurement error (umeasurement) [35][36][37].
The residual difference between the model prediction and the measurement is equal to the difference between the measurement error and the model error by reorganizing the terms of Equation (1). This relation is expressed in Equation (2).  assigned to calculate the threshold Ck. A scanning process for an STW requires a corresponding threshold set {Ck}, and the number of elements in the collection is equal to STWL.  Cumulative Rate of False alarm (CRF): for an STW containing k time steps (if STWL = k), CRF is the proportion of the total number of false alarms in the total number of natural random events without leakage when a scanning process is completed in that STW.

Leak Scenario Falsification
The core idea of the error-domain model falsification method is to falsify model instances (parameter sets). Then, the maximal plausible errors are determined by combining modeling and measurement uncertainties. There is model uncertainty (umodel) in the numerical model, as well as measurement uncertainty in the actual measured value. As shown in Equation (1), the predicted value (Pi) plus modeling error (umodel) corresponds to the real quantity (R), which is equal to the measured value (y) plus measurement error (umeasurement) [35][36][37].
The residual difference between the model prediction and the measurement is equal to the difference between the measurement error and the model error by reorganizing the terms of Equation (1). This relation is expressed in Equation (2).

Leak Scenario Falsification
The core idea of the error-domain model falsification method is to falsify model instances (parameter sets). Then, the maximal plausible errors are determined by combining modeling and measurement uncertainties. There is model uncertainty (u model ) in the numerical model, as well as measurement uncertainty in the actual measured value. As shown in Equation (1), the predicted value (P i ) plus modeling error (u model ) corresponds to the real quantity (R), which is equal to the measured value (y) plus measurement error (u measurement ) [35][36][37]. The residual difference between the model prediction and the measurement is equal to the difference between the measurement error and the model error by reorganizing the terms of Equation (1). This relation is expressed in Equation (2).
Using the residual of these two uncertainties (representing the expected residual between the predicted and measured values), the corresponding threshold boundaries (by taking the 95% interval of the probability density function) are calculated to evaluate the recognition performance of error-domain model falsification [27].
A set of possible leakage scenarios is traversed based on the set of the given parameter values. Each scenario represents one leakage state of the water distribution system (the leakage location and the leakage flow rate are used in this article). It is thought that the scenarios cover the extensive range of possible leakage states of systems, from 20 m 3 /h to 50 m 3 /h with increment of 1 m 3 /h for every node and from 50 m 3 /h to 350 m 3 /h with increment of 5 m 3 /h for every node. As shown in Figure 3a, using the EPANET pressure-driven model to simulate the leak scenario [38], the simulation is a 24-hourly Extended Period Simulation (EPS) based on forecasted water demands. Random noise is added to each node demand value to synthetically generate the node demand value of the predictive model,d n = d n + ∆d, ∆d is the uncertainty of nodal demands. For each leakage scenario, the estimated pressures are calculated where the sensors are placed in the water distribution network. Then, the pressure residual database of the leakage scenarios is obtained by comparing the estimated pressure with leakage to that without leakage. The total number of scenarios in the database equals the number of nodes multiplied by the number of leakage intensities considered.
As shown in Figure 3b, the pressure measurement (simulating raw pressure data supplied by the SCADA system) is synthesized by adding the noise to the simulated pressure value for the leakage scenario or no-leakage scenario. It is compared with the estimated pressure without leakage to obtain the measurement pressure residual. For large quantity leakage scenarios, their pressure residuals form a residual database of the leakage scenario set.
The predicted pressure residual vector r (r =p −p 0 ) of the leakage scenario can be obtained from the leakage residuals database. It characterizes the pressure deviation of all monitoring nodes. The size of the residual vector r depends on the number of pressure sensors in the network n s .p andp 0 are the predicted pressure vector of the nodes, where the pressure sensors are placed under the leakage and no-leakage conditions, respectively.
The measurement pressure residual vector can be obtained by comparing the SCADA data (in this paper, the synthetic pressure is substituted for the real SCADA data) to the estimated pressure without leakage, shown in the residual vector form: where the size of the residual vector r depends on the number of pressure sensors in the network n s . p is the pressure vector measured in the nodes where the pressure sensors are placed. Then, the measurement residual is compared to the leak residual database to calculate their correlation coefficient: where C r, r is the correlation coefficient; r is one of the pressure residual vectors in the leak residual database; r is the measurement residual at time step t associated with a potential leakage; cov(r, r) = E (r − r) r − r is the correlation function between two variables r and r, where r = E(r) and r = E( r).

Leak Detection
The larger the correlation coefficient is, the greater the similarity between the current scenario and the leak scenario and the greater the likelihood that the leakage is occurring. The maximum

Leak Detection
The larger the correlation coefficient is, the greater the similarity between the current scenario and the leak scenario and the greater the likelihood that the leakage is occurring. The maximum correlation value C r, rmax can be used to characterize the most similar leak scenario in the leakage residual database. The maximum correlation is usually used to diagnose the leakage events. When it is less than the threshold value, there is no leak at this time step and, conversely, an alarm is triggered. Figure 2 shows the schematic diagram of real-time leak detection when the STWL equals 3. From the viewpoint of judging whether a leakage happened at one certain time point, we not only use the current data at that time step to diagnose the leakage (first scan: t = k), but also wait to get the data at the next time step, combine it with the data at the former time step to diagnose the leakage (second scan: t = k, k + 1), and finally combine all of the data of the continuous three time steps (third scan: t = k, k + 1, k + 2). A complete scanning process at each STW performs three scans. After that, the STW moves forward in the direction of time to judge the WDS at the next time step. Leak detection follows a given rule: when the maximum correlation coefficient of all data at the kth scan are greater than the kth threshold (C k ), it is marked as abnormal condition and an alarm is triggered. If there is no alarm, after the time window moves forward, the qualified data is updated to the historical database to predict the hydraulic model of the next time window. From the viewpoint of the data stream, the data at the current time step is used to diagnose the leakage of WDS at the current time step, t = k, and the former two time steps, t = k − 1, and k − 2.

Estimation of Threshold Values
The threshold estimation can be calculated by the given RF and the cumulative probability distribution function, if only the data of one time step is used. If more data at continuous time steps is used, the corresponding threshold will be declined by the given RF t and the cumulative probability distribution function for those continuous time steps, with the CRF value of multiple time steps the same as RF value of one time step.
A large number of measurement results of simulated no-leak scenarios are generated using the Section 2.1 method; the cumulative probability distribution function of maximum correlation coefficient (C r, rmax ) can be obtained by traversing the correlation analysis between the measurement results and the leak residual database. As data is updated over time, the cumulative probability distribution function curve of the maximum correlation coefficient (C r, rmax ) at different time steps can be obtained. A curve of a certain time step is shown in Figure 4. From the cumulative probability distribution, we can get the threshold of C r, rmax by given RF t .
For single-time data, the corresponding threshold (C t ) can be obtained by assigning RF t to the cumulative probability distribution function curve. Figure 4 shows the threshold estimation C t = 0.83 corresponding to RF t = 10% for single-time data. For continuous time-serial data, Figure 5 gives the schematic diagram of estimation for a set of thresholds {C k } required for a complete scanning process, taking the STWL = 3 and RF t = 10% (t = 1, 2, 3) as an example. The complete equations and explanation involved in Figure 5 are provided in Appendix A.
The threshold estimation method of the first scan (Figure 5a) in the STW is consistent with the estimation method of the single-time threshold. The second scan (Figure 5b) in the STW uses the threshold C 2 to decide the false alarm events, and the third scan (Figure 5c) uses the threshold (C 3 ). The threshold (C 1 , C 2 , and C 3 ) is given to satisfy the RF i of different scans in an STW. The sum of the RF i for the different scans should equal CRF. The average allocation is used in this article (i.e., RF 1 = RF 2 = RF 3 ).
can be obtained by traversing the correlation analysis between the measurement results and the leak residual database. As data is updated over time, the cumulative probability distribution function curve of the maximum correlation coefficient ( , r rmax C ) at different time steps can be obtained. A curve of a certain time step is shown in Figure 4. From the cumulative probability distribution, we can get the threshold of , r rmax C by given RFt.  For an STW with STWL = k, according to the assigned rate of false alarm (RF 1 , RF 2 ...RF k ) and the cumulative probability distribution function curve of C r, rmax , the corresponding threshold set {C k } can be obtained. As the time step advances, each time new time data is received, a scan is performed, and the maximum correlation coefficient at the corresponding time is calculated to obtain a cumulative probability distribution function curve. The rule is that when the maximum correlation coefficients of all data at the kth scan are greater than the kth scan threshold value (C k ) simultaneously, the alarm is triggered (marked as abnormal). The cumulative rate of false alarm of k scans (t = 1, t = {1, 2}..., t = {1, 2..., k}) should be equal to the given CRF. (CRF = RF 1 + RF 2 + ... + RF k ). Equations (5)- (8) summarizes the formulas involved in kth scan: where CI k-k is the set of candidate scenario, the first k of the subscript represents the kth scan in an STW, and the second one represents the kth time data; N ({CI k-1 }∩{CI k-2 }...{CI k-k }) represents the number of intersection elements for candidate scenarios for the data of k different time steps; N k indicates the number of elements marked as abnormal at the kth scan; RF k is the rate of false alarm of the kth scan; N k is the total number of samples for the kth scan (the samples identified as abnormal for the kth scan will not appear in the total sample of the (k + 1)th scan, as shown in Equation (6)); N refers to the total number of simulated scenarios, that is, the total number of samples scanned for the first time.
For single-time data, the corresponding threshold (Ct) can be obtained by assigning RFt to the cumulative probability distribution function curve. Figure 4 shows the threshold estimation Ct = 0.83 corresponding to RFt = 10% for single-time data. For continuous time-serial data, Figure 5 gives the schematic diagram of estimation for a set of thresholds {Ck} required for a complete scanning process, taking the STWLion involved in Figure 5 are provided in Appendix A.

Case Study
The leak detection and location methodology are applied to a real water network model with synthetic data. The network is located in a city in Zhejiang Province, China. It consists of 509 pipes, 491 nodes, and 3 water sources, as shown in Figure 6. A total of 20 pressure sensors (red nodes in Figure 6) are arranged in the network. A model of this network is created using the software EPANET. Sensors 2019, 19, x FOR PEER REVIEW 10 of 20 Figure 6. Schematic representation of the water distribution network studied.

Uncertainties
Uncertainty consists of modeling error and measurement error. Modeling errors are caused by model simplification and errors related to model parameters. Model parameters are divided into two groups. The first group contains the main parameters characterizing the leak situation. For leak detection, as mentioned above, these are the leak location and leak strength. Other parameters (referred to as secondary parameters) do not clearly characterize the scenario, mainly including pipe diameter, pipe length, pipe roughness, node elevation, node requirements, and so forth. [27]. The sensitivity studies of these uncertainty values have been conducted to estimate the relative importance of each parameter. The results show that the relative importance of the uncertainty of node demand exceeds 99%.
Random noise is added to each node demand value to synthetically generate the node demand value of the predictive model. Random error follows the standard normal distribution N (0, σ) where σ is the standard deviation of each value, and is set to σ = pμ/3.27 [39], where p is the perturbation ratio, and μ is the true nodal demand value. Random noise N (0, 0.2 m) is added to the pressure measurement to characterize the measurement errors. Two different data sets are generated with different precision. In data set 1, the disturbance ratio is p = 5% for water demand and standard deviation is σ = 0.2 m for all pressure measurements. In data set 2, disturbance ratio is p = 10% for water demand and standard deviation is σ = 0.2 m for all pressure measurements. Thus, set 1 is expected to be more accurate than set 2.

Leak Scenario and Simulated Measurement
The parameter value boundary of the leakage residual database is set as follows: leakage position traverses all nodes (Nnodes = 491), and the leak intensity ranges from 20 m 3 /h to 350 m 3 /h (according to the historical leakage event database compiled by the Water Division, most of the leakage flow rate is concentrated in 20~350 m 3 /h, including small leakage and large burst, and when the leakage flow rate is less than 20 m 3 /h, the pressure fluctuation caused is too small, which is lower than the measurement accuracy of the sensor), with increment of 1 m 3 /h for intervals of 20 m 3 /h to 50 m 3 /h and 5 m 3 /h for intervals of 50 m 3 /h to 350 m 3 /h. The total number of scenarios is equal to the number of nodes multiplied by the number of intensities considered (491 × 91).
For the sample of the simulated leak scenarios, the parameters of the simulated measurement are set as follows: the amount of leakage is divided into six flow rate intervals: 20~30 m 3 /h,

Uncertainties
Uncertainty consists of modeling error and measurement error. Modeling errors are caused by model simplification and errors related to model parameters. Model parameters are divided into two groups. The first group contains the main parameters characterizing the leak situation. For leak detection, as mentioned above, these are the leak location and leak strength. Other parameters (referred to as secondary parameters) do not clearly characterize the scenario, mainly including pipe diameter, pipe length, pipe roughness, node elevation, node requirements, and so forth. [27]. The sensitivity studies of these uncertainty values have been conducted to estimate the relative importance of each parameter. The results show that the relative importance of the uncertainty of node demand exceeds 99%.
Random noise is added to each node demand value to synthetically generate the node demand value of the predictive model. Random error follows the standard normal distribution N (0, σ) where σ is the standard deviation of each value, and is set to σ = pµ/3.27 [39], where p is the perturbation ratio, and µ is the true nodal demand value. Random noise N (0, 0.2 m) is added to the pressure measurement to characterize the measurement errors. Two different data sets are generated with different precision. In data set 1, the disturbance ratio is p = 5% for water demand and standard deviation is σ = 0.2 m for all pressure measurements. In data set 2, disturbance ratio is p = 10% for water demand and standard deviation is σ = 0.2 m for all pressure measurements. Thus, set 1 is expected to be more accurate than set 2.

Leak Scenario and Simulated Measurement
The parameter value boundary of the leakage residual database is set as follows: leakage position traverses all nodes (N nodes = 491), and the leak intensity ranges from 20 m 3 /h to 350 m 3 /h (according to the historical leakage event database compiled by the Water Division, most of the leakage flow rate is concentrated in 20~350 m 3 /h, including small leakage and large burst, and when the leakage flow rate is less than 20 m 3 /h, the pressure fluctuation caused is too small, which is lower than the measurement accuracy of the sensor), with increment of 1 m 3 /h for intervals of 20 m 3 /h to 50 m 3 /h and 5 m 3 /h for intervals of 50 m 3 /h to 350 m 3 /h. The total number of scenarios is equal to the number of nodes multiplied by the number of intensities considered (491 × 91).
For the sample of the simulated leak scenarios, the parameters of the simulated measurement are set as follows: the amount of leakage is divided into six flow rate intervals: 20~30 m 3 /h, 30~40 m 3 /h, 40~50 m 3 /h, 50~100 m 3 /h, 100~200 m 3 /h, 200~350 m 3 /h. A node is selected as the leak position randomly, and a random leakage flow rate in the corresponding interval is selected using a random number generation function in C/C++. For every leakage event, a simulation is run to obtain the pressures at the measurement points. Each sample set corresponds to a scenario where a new leak occurs at time k, and the proportion of samples in each flow rate interval of leak is: 10%, 10%, 10%, 20%, 30%, 20%, respectively.
For samples which simulate no-leak scenarios, as described in Section 3.1, two comparison data sets are generated. On the one hand, for evaluating leak detection performance, the number of samples per data set is equal to the number of samples that simulate leak scenarios (this article uses N = 10 4 ); on the other hand, for threshold estimation, the number of samples is set as 10 5 .
The generation principle of the above three types of samples has been described in Section 2.1 in detail.

Estimation of Threshold Values and Leak Detection
The threshold value used was obtained from the method of Section 2.2 and the synthetic data from Section 3.2 for leak detection testing. The time series is normalized before the data at different times is processed to obtain a threshold, as shown in Equations (9) and (10).
where R nor(t) is the normalization coefficient at time t; N represents the total number of samples; C r, rmax (N) is maximum correlation coefficient of the Nth sample at time t; C r, rmax (N) nor indicates the maximum correlation coefficient of the Nth sample after normalization. The online monitoring k-hour scenario is simulated, and then the corresponding k-group thresholds are obtained for a given STWL of 1~24 by threshold analysis. Table 1 shows threshold set for each STW with entire time horizon obtained by threshold analysis when CRF = 10%, STWL = 3. The threshold set corresponding to the scanning process varies over time.

Results and Discussion
Two different predicted water demand errors (p = 5% and 10%) are considered. The cumulative rate of false alarm (CRF) is set from 5% to 30%, with an increment of 2.5%. The RF k is set to average allocation in this paper. Noticeably, it has been confirmed that there is a better allocation ratio scheme (see the Appendix B for details); it is necessary to automatically optimize the allocation ratio by using optimization methods in the future. For every CRF, the DP versus STWL curves are obtained. The curves shown in Figure 7a are the scenarios of CRF = 10%, 15%, 20%, and 30%, with leak intensity ranges from 20 m 3 /h to 350 m 3 /h. Compared to the single-time STW, we can clearly see that the proposed method using multitime STW can effectively improve the detection probability for a given cumulative rate of false alarm.
The detection probability increases with the STWL and then descends with the STWL. In the initial STWLs, when the STWL is set from 1 to 4 (the time step of one hour is used in this paper), the detection probability is significantly improved. After that, the detection probability shows a downward trend with the increasing STWL. There is an optimal range of STWL for detection probability. A wide STW means a large number of scans in an STW. The last number of scans cover continuous monitoring data over the STW. For successful leakage detection, all correlation coefficients of the data should support the same judgments. However, the data includes random demand noise and pressure measurement noise. The longer the data series is used, the harder it is to make consistent judgment. Thus, the detection probability descends with the increase of the STWL when the STWL is over a certain value.
In comparison with the clustering-based method using cosine distance to evaluate dissimilarity between data from multiple flow sensors to detect bursts [40], we can find a similar trend where RF increases with the decrease in window size when the DP remains constant, although different analysis methods are used for multitime series data. Using the optimal STWL value can enhance DP by 3%, compared to using single-time STW (STWL = 1) for leak detection. When STWL increases beyond 16, the detection probability is lower than that of single-time STW. The discussion of the influence of STWL on detection performance can be found in many other leak detection methods. Abokifa et al. [41] demonstrated that an optimal window period can be selected to achieve the best performance, while smaller and larger windows would generally yield less accurate results. Wu et al. [42] drew the conclusion that detection performance of the method changes only slightly when STWL varies from 1 to 6 days. Figure 7b shows that the predicted water demand error has an obvious effect on the proposed leakage detection method. The water demand disturbance of set 1 (p = 5%) is smaller than for set 2 (p = 10%), and the leakage detection probability of set 1 is higher than for set 2. On the other hand, the leakage detection probability shows a similar trend for different cumulative rates of false alarm In order to exploit the advantages of the proposed method using multitime STW, the DP versus STWL curves for different leakage flow rate intervals are shown in Figure 8 when the CRF is set at 10%. Six intervals are divided: 20~30 m 3 /h, 30~40 m 3 /h, 40~50 m 3 /h, 50~100 m 3 /h, 100~200 m 3 /h, 200~350 m 3 /h. The optimal solution points indicate the optimal STWL value corresponding to the highest detection probability within the corresponding flow rate interval. The curves of each flow rate interval show a similar trend for differently tested data sets. The detection performance of the the affected node and may be overlooked. In short, fluctuation caused by larger leakages is greater and less likely to be masked under uncertainty at a given level. The STWLs that have optimal detection performance for different flow rate intervals are concentrated from 3 to 6 for the tested data set. Similarly, Wu et al. [40] chose the STWL of approximately 6 h for burst detection. However, in his paper, every STW only scans once, while a total of k scans (t = 1, t = {1,2}...,t = {1,2,...,k}) are performed for every STW in this paper.  In order to exploit the advantages of the proposed method using multitime STW, the DP versus STWL curves for different leakage flow rate intervals are shown in Figure 8 when the CRF is set at 10%. Six intervals are divided: 20~30 m 3 /h, 30~40 m 3 /h, 40~50 m 3 /h, 50~100 m 3 /h, 100~200 m 3 /h, 200~350 m 3 /h. The optimal solution points indicate the optimal STWL value corresponding to the highest detection probability within the corresponding flow rate interval. The curves of each flow rate interval show a similar trend for differently tested data sets. The detection performance of the large flow rate interval is better than that of the small flow rate interval. This is primarily due to the relationship between the fluctuation caused by the leakage and the degree of uncertainty. Small leakages located in areas with high uncertainty are not detectable due to the small variations caused by them; the fluctuation caused by the leakage is smaller than the fluctuation of the uncertainty of the affected node and may be overlooked. In short, fluctuation caused by larger leakages is greater and less likely to be masked under uncertainty at a given level. The STWLs that have optimal detection performance for different flow rate intervals are concentrated from 3 to 6 for the tested data set. Similarly, Wu et al. [40] chose the STWL of approximately 6 h for burst detection. However, in his paper, every STW only scans once, while a total of k scans (t = 1, t = {1,2}...,t = {1,2,...,k}) are performed for every STW in this paper.
In comparison with the single-time model-based leak detection methodology [24], the method in this paper fully uses the temporal varying (multiple time rather than single time) correlation between the data from multiple pressure sensors, and the detection performance is improved. In addition, Casillas Ponce et al. [43] proposed a model-based approach for leak detection and location, which considers an extended time horizon analysis of pressure sensitivities. To extract data information in the time horizon, the authors look for the mean value or mode in an STW for different metrics (binarization, correlation, angle between vectors, and Euclidean distance). The advantage of the method using the correlation coefficient metric in this paper is that multitime series residual analysis can be more sensitive to leakage than using mode or mean value in an STW. Figure 9 shows the trade-off curve for the given cumulative rate of false alarm and rate of missing detection of the tested data set when selecting the optimal solution of STWL. The users can make the decisions to choose a combination between the false positive rate (rate of false alarm) and false negative rate (rate of missing detection) according to the actual situation, and every point of the curve corresponds to the solution of the optimal STWL. For the leakage detection performance at the optimal STWL, the demand error of p = 10% decreased the DP of 1% more than that of p = 5%. In comparison with the single-time model-based leak detection methodology [24], the method in this paper fully uses the temporal varying (multiple time rather than single time) correlation between the data from multiple pressure sensors, and the detection performance is improved. In addition, Casillas Ponce et al. [43] proposed a model-based approach for leak detection and location, which considers an extended time horizon analysis of pressure sensitivities. To extract data information in the time horizon, the authors look for the mean value or mode in an STW for different metrics (binarization, correlation, angle between vectors, and Euclidean distance). The advantage of the method using the correlation coefficient metric in this paper is that multitime series residual analysis can be more sensitive to leakage than using mode or mean value in an STW. Figure 9 shows the trade-off curve for the given cumulative rate of false alarm and rate of missing detection of the tested data set when selecting the optimal solution of STWL. The users can make the decisions to choose a combination between the false positive rate (rate of false alarm) and false negative rate (rate of missing detection) according to the actual situation, and every point of the curve corresponds to the solution of the optimal STWL. For the leakage detection performance at the optimal STWL, the demand error of p = 10% decreased the DP of 1% more than that of p = 5%. Figure 9. Relationship between cumulative rate of false alarm and rate of missing detection. (Each STWL was tested for 24 different start times, and the vertical coordinate value of the point in Figure  7 and Figure 8 is the average of 24 different DPs. Figure 9 utilizes a similar processing method to obtain an average of 24 1-DPs when selecting the optimal solution of STWL.)

Conclusions
This work proposes a leak detection method combined with multitime series. Leak detection is based on the comparison between real-time online data (simulated measurements) and predictive model data. The proposed method is about threshold selection method based on time series monitoring data. The detection is based on multitime series residual analysis. The influences of the uncertainty of the prediction model and measurement noise are considered for this proposed method.
The higher the CRF, the higher the leakage detection probability. Under given parameters, detection performance for the large leakage is better than for the small leakage because the pressure change caused by the low-intensity leakage may be overlooked by the modeling and measurement  Figure 9 utilizes a similar processing method to obtain an average of 24 1-DPs when selecting the optimal solution of STWL.)

Conclusions
This work proposes a leak detection method combined with multitime series. Leak detection is based on the comparison between real-time online data (simulated measurements) and predictive model data. The proposed method is about threshold selection method based on time series monitoring data. The detection is based on multitime series residual analysis. The influences of the uncertainty of the prediction model and measurement noise are considered for this proposed method.
The higher the CRF, the higher the leakage detection probability. Under given parameters, detection performance for the large leakage is better than for the small leakage because the pressure change caused by the low-intensity leakage may be overlooked by the modeling and measurement uncertainties. For the optimal STWL, the DP with demand error of p = 10% decreased by 1% over p = 5%. The future improvement of the water demand prediction method and installing of extra sensors can reduce the effects of uncertainty. For different flow rate intervals and different cumulative rate of false alarm, the leakage detection performance and the STWL show a positive correlation trend firstly, and then the trend is a negative correlation. The STWLs for achieving optimal detection performance for different conditions are concentrated from 3 to 6 for the tested data set. Using the optimal STWL can enhance DP of 3% more than using single-time STW for leak detection.
The multitime series analysis method is used to enhance the detection probability and to reduce the rate of missing detection and the rate of false alarm effectively, compared with the single-time series method. Many factors (the STWL, RF k , CRF, ∆d, etc.) will affect the performance of the detection method. The diversity combination of the parameter settings will be the next research work. The main purpose of this paper is to prove the feasibility of the method. In addition, various leakage indicators can be obtained by different metrics, and then leakage detection methods derived can be coupled with this method, which has wide applicability. Future work will consider the minimum detectable leak depending on the resolution of sensors. In addition, flow sensors will be tested and compared with pressure sensors in order to assess which is the best option. A real-case test will be performed when real data become available in the future. Finally, the applicability of the method will be evaluated by comparison with other methods.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
The equations involved in Figure 5 in this study are summarized as follows: where N 1 is the total number of samples for the 1st scan; P C r, rmax1 ≥ C 1 is the probability that the maximum correlation coefficient of the first time is greater than the threshold value (C 1 ); N refers to the total number of simulated scenarios; CI 1-1 is the set of candidate scenario indexes representing the first scan threshold (the former 1) for the first time data (the latter 1); N 1 indicates the number of elements marked as abnormal at the first scan; RF 1 is the rate of false alarm of the first scan.
P C r, r max 1 ≥ C 2 * N 2 = N{CI 2−1 } (A6) P C r, r max 2 ≥ C 2 * N 2 = N{CI 2−2 } (A7) where N 2 is the total number of samples for the second scan; P C r, rmax1 ≥ C 2 , P C r, rmax2 ≥ C 2 are the probability that the maximum correlation coefficient of the first time and second time are greater than the threshold value, respectively, when the threshold value is C 2 ; CI 2-1 , CI 2-2 are the set of candidate scenario indexes representing the second scan threshold for the first time data and second time data, respectively; N 2 indicates the number of elements marked as abnormal at the second scan; RF 2 is the rate of false alarm of the second scan.
where N 3 is the total number of samples for the third scan; P C r, rmax1 ≥ C 3 , P C r, rmax2 ≥ C 3 , P C r, rmax3 ≥ C 3 are the probability that the maximum correlation coefficient of the first time, second time, and third time are greater than the threshold value, respectively, when the threshold value is C 3 ; CI 3-1 , CI 3-2 , CI 3-3 are the set of candidate scenario indexes representing the third scan threshold for the first time data, second time data, and third time data, respectively; N 3 indicates the number of elements marked as abnormal at the third scan; RF 3 is the rate of false alarm of the third scan. RF 1 = RF 2 = RF 3 (A16) CRF = RF 1 +RF 2 +RF 3 = 30% (A17) Figure A1. The DPs of different allocation ratio schemes for different given CRF when STWL = 3: (a) σ = 0.2m, p = 10%; (b) σ = 0.2m, p = 5%.
The category axis of the radar chart represents the different schemes of the allocation ratio (e.g., '1/3 1/3 1/3' represents RF 1 = RF 2 = RF 3 = CRF/3), and the value axis represents the corresponding DP. Partial allocation ratio schemes are shown in Figure A1 for different conditions.