Water contamination events are threatening the safety of drinking water. Free chlorine is widely used as the disinfectant in drinking water, which can be used as a surrogate parameter to provide indications of potential contaminants. In this article, the periodic fluctuation of free chlorine is studied and the fluctuation pattern is extracted by the singular vector decomposition method, and an anomaly detection scheme for free chlorine is proposed and tested. Firstly, the normal periodic pattern and current pattern of free chlorine are both extracted from the historical and online data, and then the difference between the current data pattern and the normal data pattern are compared with thresholds for anomaly declaration. The single point detection and data series detection are investigated for the purpose of short-term and long-term inspection. Further, the anomaly data treatment and the detection method using sub-patterns are discussed. Performance tests show that the proposed method is sensitive to the anomaly data, and is effective to detect anomalous condition in typical contamination scenes.

INTRODUCTION

Drinking water is vulnerable to contaminant intrusion when transported in network pipes (Laird et al. 2005). A contamination warning system (CWS) is developed to safeguard the drinking water (USEPA 2010). Generally, the aim of a CWS is to provide an integrated solution to the detection and timely warnings of occurring or potential contamination incidents. The water quality CWSs have been deployed in many countries including the USA, China, Germany and Japan (Storey et al. 2011; Hou et al. 2013).

Anomaly detection is the most important function and component of a CWS, which analyzes the water quality data and detects potential contamination events. Anomaly detection methodologies have been proposed to analyze the real-time water quality data and provide indication of contaminants, including the statistical properties method (Byer & Carlson 2005), prediction and comparison method (McKenna et al. 2008; Perelman et al. 2012), clustering method (Klise & McKenna 2006) and neural network method (Hart et al. 2010a, b). The main task is to find the normal pattern in the normal data and compare it with the testing data pattern. If the two patterns are not consistent, then the current data can be declared as an anomaly.

In this paper, the free chlorine in drinking water networks is selected as a particular concern. Experiments show that free chlorine is sensitive to many contaminants, whose concentration significantly decreases when contaminants are introduced into water (USEPA 2005; Jeffrey Yang et al. 2009), so free chlorine is suitable to be the surrogate parameter for many contaminants. The free chlorine concentration fluctuates, due to the decay of chlorine and variation of water consumption (Byer & Carlson 2005; McKenna et al. 2008). The principle of free chlorine decay is studied and the widely accepted models are the first-order model and the second-order model (Hua et al. 1999). Generally, free chlorine concentration decreases in the night and rises in the daytime, hence exhibits a periodic pattern (Polycarpou et al. 2002). In the daytime, free chlorine concentration rises as fresh water is supplied continuously from the water treatment plant. At night the water demand is small, the free chlorine will stay in the pipe and decreases to the lowest level. A detection method utilizing the fluctuation characteristics of free chlorine is studied (Eliades & Polycarpou 2013), in which the upper and lower bounds of free chlorine concentration is calculated and used for anomaly detection. However, other studies of anomaly detection utilizing the free chlorine fluctuation are rarely found.

In this paper, we aim to extract the periodic fluctuation of free chlorine as normal pattern for event detection using singular value decomposition (SVD)-based method, and a new detection scheme using pattern matching is proposed. SVD is a matrix decomposition method which decomposes a matrix into two singular vector matrices and a diagonal matrix as where A is a matrix, U and V are singular vector matrices and is a diagonal matrix containing the singular values. SVD has been proposed for pattern extraction, in which data series is configured into a matrix and SVD is carried out, with the singular vector used as the periodic pattern (Kanjilal et al. 1999).

The proposed anomaly detection scheme has a training and testing phase. In the training phase, the periodic pattern is extracted as normal pattern, with a suitable threshold trained for anomaly declaration. In the testing phase, the residual between the current data pattern and the normal pattern is calculated and compared with the threshold. Two detection modes are proposed: the single point mode and data series mode, for short-term and long-term anomaly inspection, respectively.

The paper is organized as follows. Firstly, the methodology for periodic pattern extraction and real-time anomaly detection scheme are presented. Then, a case study using real measured free chlorine data is presented and discussed. Finally, the conclusion is drawn and future work is introduced.

METHODOLOGY FOR PATTERN EXTRACTION AND ANOMALY DETECTION

Periodic pattern extraction

Definition 1 (periodic)

A series of data points is called periodic, if there is a basic periodic pattern and x can be divided into continuous segments with length N (assuming that , ), and for any segments , , 
formula
where is the scaling factor for the basic pattern, and denotes the absolute value. p is a normalized pattern, so is used to scale the pattern for comparison with . If is arbitrarily small, then x is called strictly periodic, otherwise x is called nearly periodic.
The periodic pattern can be extracted using SVD. The SVD of a matrix is defined as 
formula
where ( and are column vectors), with and , and assuming ( is a zero matrix). The singular values of are denoted by . and correspond to the singular value , which is called the th left singular vector and the th right singular vector, respectively.
A data series x can be configured into an matrix An, 
formula
SVD of produces , and it can be observed that (Kanjilal et al. 1999),
  1. If x is strictly periodic with period length N and row length in , then is non-zero, but .

  2. If x is nearly periodic with period length N and row length in , then will be much larger than the rest of the singular values, i.e., .

In all the cases, the vector represents the ‘basic periodic pattern’ of the data series corresponding to in Definition 1. The successive elements of the vector represent the amplitude scaling factors of successive pattern segments corresponding to in Definition 1.

So, the main steps of the extraction of periodic pattern using SVD can be concluded:

  1. x is configured into matrix with varying row length n and the SVD of is performed as .

  2. The spectrum of the ratio of first two singular values versus row length n is calculated, which is called the periodicity spectrum or ‘ spectrum’.

  3. Find the repetitive peaks at (where i is a positive integer), if there is any embedded periodic component in x with period length N.

  4. Repeating the process is performed for residual until no repetitive peaks exists in the p spectrum.

  5. Finally, can be modeled as the addition of the basic periodic pattern multiplied by the scaling factor.

This method mainly utilizes the linear correlation between periods when the row length of is equal to the period length.

Real-time anomaly detection

In real-time water quality anomaly detection, the main objective is to find the difference between the current data pattern and the normal pattern. If the current data pattern is significantly different from the normal pattern, a possible anomaly can be declared. Based on the institution, we propose the real-time anomaly detection method which uses the periodic pattern, as shown in Figure 1.

Figure 1

Anomaly detection flow chart. Normal pattern is extracted and compared with current pattern, and if the residual exceed the threshold an anomaly is declared.

Figure 1

Anomaly detection flow chart. Normal pattern is extracted and compared with current pattern, and if the residual exceed the threshold an anomaly is declared.

Firstly, the periodic fluctuation of free chlorine is extracted as the normal pattern using the SVD method. The residual of the current data pattern from the normal data pattern is calculated and compared to a pre-set threshold. If the residual exceeds the threshold, the current pattern can be declared as an anomaly. We propose two detection modes: single point mode and data series mode, for single point and data series anomaly declaration, respectively.

If x has the periodic patterns , only is used as the basic normal pattern. The main steps for detection are as follows:

The basic normal pattern with period length N is obtained by the SVD extraction method.

The last N data point of x is used as the current data pattern pc within a sliding window of length N, denoted as , .

The residual of from p is calculated as 
formula
1
where r denotes the residual vector. The scaling factor is calculated to minimize the difference between the current pattern and the normal pattern, and denotes the 2-norm for vectors.

For the residual series , , (a) in single point mode, is compared with a threshold . If , is declared as an anomaly; (b) in data series mode, is compared with threshold . If , is declared as a series of anomalies.

Event simulation

Performance is always the biggest concern for an anomaly detection method. In performance test-simulated event is usually utilized (McKenna et al. 2008), as well as water quality simulation of drinking water network using a simulation tool such as EPANET MSX (Hart et al. 2010a, b). In this article, the simulated event is utilized. Event e can be represented as 
formula
where s is the event strength, l is the event duration and is the basic shape (such as square shape and Gaussian shape). Experiments show that the decrease of free chlorine is nearly square shaped (Jeffrey Yang et al. 2009), therefore the square shape is used in the test.

Detection rate (DR) and false alarm rate (FAR) are used to judge the detection performance. DR is the proportion of event points detected, and FAR is the proportion of normal data declared as an anomaly.

CASE STUDY OF PERIODIC PATTERN EXTRACTION OF FREE CHLORINE

Data acquisition

In the case study the real-time free chlorine data collected from an on-line drinking water quality monitoring and CWS in China measured every 15 minutes from 5 September to 20 September 2013 (about 1,500 time steps) were used, as shown in Figure 2(a), with the first 10 days of data as the training data set as a benchmark, and the following 5 days of data as the testing data.

Figure 2

(a) Training data set and testing data set of free chlorine collected every 15 minutes for a total of about 1,500 time steps. (b) p-spectrum obtained from training data with versus varying row length in plotted. The marked circles are the peaks in the curve, indicating that the period of the data is roughly 96 time steps. (c) Extracted periodic pattern data (solid line) and measured data/training data (dashed line) plotted in series for each period. The pattern data fit the measured data well, indicating that the extracted pattern well represents the fluctuation of free chlorine.

Figure 2

(a) Training data set and testing data set of free chlorine collected every 15 minutes for a total of about 1,500 time steps. (b) p-spectrum obtained from training data with versus varying row length in plotted. The marked circles are the peaks in the curve, indicating that the period of the data is roughly 96 time steps. (c) Extracted periodic pattern data (solid line) and measured data/training data (dashed line) plotted in series for each period. The pattern data fit the measured data well, indicating that the extracted pattern well represents the fluctuation of free chlorine.

Extracted pattern

The p spectrum is obtained and peaks exist clearly at the 95, 192, 289, 381, 479 time steps as shown in Figure 2(b), which are roughly times 96. So it can be concluded that the free chlorine has a periodic pattern of 96 time steps which is exactly the length of 1 day. Although 2 days' pattern or longer also can be treated as a repeated pattern in the data series, we take the minimal length of 1 day as the basic period length, and the data pattern throughout a week is almost stable, not significantly different in the weekends and weekdays.

The comparison between the pattern extracted and the real data is shown in Figure 2(c), with real data presented as a dashed line. The first time step of the basic period is chosen at the 09:00, and the pattern has four peaks at 25, 46, 62 and 88 time step, respectively, corresponding to 11:15, 15:15, 21:45 and 06:00. In a fluctuation period free chlorine decreases significantly from about 21:45 at night to 06:00 in the morning, then rises to a peak at 11:15, then decreases to a low value at 15:15 and rises again to a peak at 21:45. It can be inferred that the water demand is high in the morning and in the late afternoon to midnight.

CASE STUDY OF PERFORMANCE TEST

Threshold training

As the water demand in the network fluctuates in a certain range and is relatively stable, the normal fluctuation amplitude of free chlorine can be obtained using training data, corresponding to the residual in detection. The threshold can be determined by an empirical method using the mean and standard deviation of the residuals at the same time each day, 
formula
2
where denotes the threshold at the th position, denotes the mean of residuals at the th position, ’ denotes standard deviation and denotes the threshold strength. is applied for practical threshold (USEPA 2010), and is applied for normal distributed data for anomaly detection (Byer & Carlson 2005). In this article, we use for more sensitive detection performance.

Performance test results

Performance is usually affected by the time of the event occurring and its duration (USEPA 2013), so we use a set of simulated events to test all the scenarios which has the starting time varying from the first time step to the 89th time step every eight steps (a total of 12 cases), and the duration varying from four to 92 time steps, every eight steps (a total of 12 cases). Training data (10 days) and testing data (5 days) are used for the test, respectively. As for event strength, according to the Chinese National Water Quality Standard the lower limit for free chlorine is 0.05 mg/L, and the minimum value of the training and testing data set is 0.1021 mg/L, so a suitable event strength is about 0.05 − 0.1021= −0.0521 while traditional threshold based detection cannot detect the anomaly. The test process is as follows:

  1. Obtain the thresholds from training data with set in Equation (2) ( and 2.5 are used).

  2. Set the event strength s, event duration l and starting time .

  3. Obtain an event as ( denotes the square shaped function), and add the event on the original data from the time step.

  4. Run the detection method and calculate the DR and FAR.

  5. Change event duration l, starting time and event strength s, and repeat steps 3 and 4.

  6. Calculate the averaged DR and FAR for different durations and starting time.

The test results are in Tables 16, in which event strength −0.025, −0.05 and −0.075 are used, respectively, and the labels ‘training’ and ‘testing’ in the tables imply using training data and testing data, respectively. All the tests are carried out in single point detection mode. When a fixed duration is used for the test, the averaged DR and FAR are calculated from selected cases of starting time, and the same for the starting time. The following can be concluded from the tables:

  1. DR is positively correlated with the event strength, and negatively correlated with the duration, which is discussed in the next subsection.

  2. FAR is slightly positively correlated with the duration and event strength. When duration is long, the anomaly data in the middle of the event are difficult to detect as the whole data series is shifted. So, the supplementary threshold-based method (Byer & Carlson 2005) can be applied to improve detection performance. In our test, the supplementary method makes the DR a little higher when duration is higher than 68 time steps (not shown).

  3. For different starting times, DR reaches the maximum value when the starting time is between the 40–60 time steps and has large differences. FAR is much more stable.

  4. For of threshold, the FAR decreases to half and the DR also decreases a little compared with the values for the case of .

Table 1

The averaged DR for different durations, with

 Event strength/duration41220283644526068768492
Training −0.025 0.403 0.380 0.360 0.342 0.324 0.305 0.287 0.270 0.254 0.240 0.226 0.213 
−0.05 0.743 0.720 0.701 0.680 0.652 0.622 0.594 0.568 0.541 0.514 0.487 0.460 
−0.075 0.935 0.924 0.910 0.887 0.860 0.831 0.799 0.768 0.733 0.700 0.666 0.632 
Testing −0.025 0.438 0.418 0.406 0.400 0.399 0.398 0.398 0.396 0.392 0.387 0.381 0.374 
−0.05 0.688 0.675 0.667 0.658 0.650 0.640 0.625 0.610 0.594 0.577 0.561 0.545 
−0.075 0.911 0.896 0.876 0.859 0.841 0.820 0.799 0.777 0.753 0.728 0.703 0.680 
 Event strength/duration41220283644526068768492
Training −0.025 0.403 0.380 0.360 0.342 0.324 0.305 0.287 0.270 0.254 0.240 0.226 0.213 
−0.05 0.743 0.720 0.701 0.680 0.652 0.622 0.594 0.568 0.541 0.514 0.487 0.460 
−0.075 0.935 0.924 0.910 0.887 0.860 0.831 0.799 0.768 0.733 0.700 0.666 0.632 
Testing −0.025 0.438 0.418 0.406 0.400 0.399 0.398 0.398 0.396 0.392 0.387 0.381 0.374 
−0.05 0.688 0.675 0.667 0.658 0.650 0.640 0.625 0.610 0.594 0.577 0.561 0.545 
−0.075 0.911 0.896 0.876 0.859 0.841 0.820 0.799 0.777 0.753 0.728 0.703 0.680 
Table 2

The averaged DR for different starting time, with

 Event strength/starting time1917253341495765738189
Training −0.025 0.246 0.247 0.255 0.239 0.262 0.304 0.311 0.345 0.409 0.394 0.335 0.258 
−0.05 0.434 0.466 0.535 0.537 0.558 0.610 0.672 0.772 0.773 0.717 0.648 0.559 
−0.075 0.663 0.746 0.787 0.796 0.829 0.874 0.889 0.864 0.845 0.814 0.783 0.755 
Testing −0.025 0.131 0.175 0.259 0.212 0.296 0.452 0.564 0.626 0.630 0.583 0.507 0.352 
−0.05 0.298 0.380 0.453 0.470 0.554 0.806 0.831 0.830 0.803 0.753 0.692 0.619 
−0.075 0.592 0.660 0.784 0.797 0.866 0.924 0.905 0.880 0.850 0.828 0.797 0.760 
 Event strength/starting time1917253341495765738189
Training −0.025 0.246 0.247 0.255 0.239 0.262 0.304 0.311 0.345 0.409 0.394 0.335 0.258 
−0.05 0.434 0.466 0.535 0.537 0.558 0.610 0.672 0.772 0.773 0.717 0.648 0.559 
−0.075 0.663 0.746 0.787 0.796 0.829 0.874 0.889 0.864 0.845 0.814 0.783 0.755 
Testing −0.025 0.131 0.175 0.259 0.212 0.296 0.452 0.564 0.626 0.630 0.583 0.507 0.352 
−0.05 0.298 0.380 0.453 0.470 0.554 0.806 0.831 0.830 0.803 0.753 0.692 0.619 
−0.075 0.592 0.660 0.784 0.797 0.866 0.924 0.905 0.880 0.850 0.828 0.797 0.760 
Table 3

The averaged FAR for different durations, with

 Event/duration41220283644526068768492
Training −0.025 0.061 0.063 0.065 0.066 0.067 0.069 0.070 0.072 0.073 0.074 0.075 0.075 
−0.05 0.063 0.066 0.069 0.074 0.079 0.084 0.087 0.090 0.092 0.093 0.095 0.096 
−0.075 0.063 0.069 0.077 0.085 0.092 0.098 0.101 0.104 0.107 0.109 0.110 0.110 
Testing −0.025 0.120 0.118 0.117 0.118 0.120 0.121 0.122 0.124 0.125 0.126 0.126 0.125 
−0.05 0.120 0.119 0.123 0.127 0.131 0.135 0.141 0.146 0.149 0.153 0.156 0.156 
−0.075 0.118 0.123 0.130 0.137 0.148 0.157 0.170 0.181 0.186 0.188 0.188 0.187 
 Event/duration41220283644526068768492
Training −0.025 0.061 0.063 0.065 0.066 0.067 0.069 0.070 0.072 0.073 0.074 0.075 0.075 
−0.05 0.063 0.066 0.069 0.074 0.079 0.084 0.087 0.090 0.092 0.093 0.095 0.096 
−0.075 0.063 0.069 0.077 0.085 0.092 0.098 0.101 0.104 0.107 0.109 0.110 0.110 
Testing −0.025 0.120 0.118 0.117 0.118 0.120 0.121 0.122 0.124 0.125 0.126 0.126 0.125 
−0.05 0.120 0.119 0.123 0.127 0.131 0.135 0.141 0.146 0.149 0.153 0.156 0.156 
−0.075 0.118 0.123 0.130 0.137 0.148 0.157 0.170 0.181 0.186 0.188 0.188 0.187 
Table 4

The averaged FAR for different starting time, with

 Event strength/starting time1917253341495765738189
Training −0.025 0.070 0.070 0.070 0.069 0.069 0.069 0.068 0.069 0.069 0.069 0.069 0.070 
−0.05 0.084 0.083 0.083 0.082 0.082 0.082 0.081 0.081 0.082 0.082 0.083 0.084 
−0.075 0.096 0.095 0.095 0.094 0.094 0.094 0.093 0.093 0.092 0.092 0.093 0.094 
Testing −0.025 0.120 0.122 0.123 0.123 0.123 0.126 0.123 0.119 0.118 0.120 0.124 0.124 
−0.05 0.138 0.139 0.138 0.138 0.138 0.142 0.138 0.135 0.134 0.136 0.140 0.140 
−0.075 0.167 0.165 0.162 0.160 0.158 0.161 0.158 0.154 0.153 0.156 0.161 0.160 
 Event strength/starting time1917253341495765738189
Training −0.025 0.070 0.070 0.070 0.069 0.069 0.069 0.068 0.069 0.069 0.069 0.069 0.070 
−0.05 0.084 0.083 0.083 0.082 0.082 0.082 0.081 0.081 0.082 0.082 0.083 0.084 
−0.075 0.096 0.095 0.095 0.094 0.094 0.094 0.093 0.093 0.092 0.092 0.093 0.094 
Testing −0.025 0.120 0.122 0.123 0.123 0.123 0.126 0.123 0.119 0.118 0.120 0.124 0.124 
−0.05 0.138 0.139 0.138 0.138 0.138 0.142 0.138 0.135 0.134 0.136 0.140 0.140 
−0.075 0.167 0.165 0.162 0.160 0.158 0.161 0.158 0.154 0.153 0.156 0.161 0.160 
Table 5

The averaged DR for different durations, with

 Event strength/duration41220283644526068768492
Training −0.025 0.211 0.204 0.191 0.178 0.162 0.148 0.135 0.123 0.112 0.102 0.094 0.087 
−0.05 0.616 0.587 0.554 0.522 0.490 0.460 0.431 0.403 0.375 0.347 0.321 0.298 
−0.075 0.859 0.832 0.805 0.774 0.742 0.707 0.672 0.636 0.600 0.564 0.529 0.495 
Testing −0.025 0.271 0.267 0.255 0.248 0.243 0.239 0.233 0.228 0.221 0.215 0.209 0.203 
−0.05 0.568 0.559 0.547 0.534 0.519 0.502 0.485 0.471 0.455 0.441 0.426 0.412 
−0.075 0.823 0.792 0.764 0.740 0.716 0.695 0.675 0.653 0.627 0.602 0.578 0.556 
 Event strength/duration41220283644526068768492
Training −0.025 0.211 0.204 0.191 0.178 0.162 0.148 0.135 0.123 0.112 0.102 0.094 0.087 
−0.05 0.616 0.587 0.554 0.522 0.490 0.460 0.431 0.403 0.375 0.347 0.321 0.298 
−0.075 0.859 0.832 0.805 0.774 0.742 0.707 0.672 0.636 0.600 0.564 0.529 0.495 
Testing −0.025 0.271 0.267 0.255 0.248 0.243 0.239 0.233 0.228 0.221 0.215 0.209 0.203 
−0.05 0.568 0.559 0.547 0.534 0.519 0.502 0.485 0.471 0.455 0.441 0.426 0.412 
−0.075 0.823 0.792 0.764 0.740 0.716 0.695 0.675 0.653 0.627 0.602 0.578 0.556 
Table 6

The averaged FAR for different durations, with

 Event strength/duration41220283644526068768492
Training −0.025 0.000 0.000 0.000 0.001 0.002 0.003 0.003 0.004 0.005 0.005 0.006 0.006 
−0.05 0.000 0.001 0.003 0.006 0.009 0.012 0.016 0.019 0.021 0.023 0.024 0.025 
−0.075 0.000 0.002 0.007 0.013 0.021 0.027 0.032 0.035 0.038 0.039 0.041 0.041 
Testing −0.025 0.034 0.034 0.035 0.036 0.036 0.037 0.039 0.039 0.040 0.041 0.042 0.042 
−0.05 0.035 0.036 0.037 0.041 0.046 0.051 0.055 0.057 0.060 0.064 0.066 0.066 
−0.075 0.034 0.037 0.043 0.052 0.058 0.066 0.074 0.082 0.090 0.094 0.097 0.097 
 Event strength/duration41220283644526068768492
Training −0.025 0.000 0.000 0.000 0.001 0.002 0.003 0.003 0.004 0.005 0.005 0.006 0.006 
−0.05 0.000 0.001 0.003 0.006 0.009 0.012 0.016 0.019 0.021 0.023 0.024 0.025 
−0.075 0.000 0.002 0.007 0.013 0.021 0.027 0.032 0.035 0.038 0.039 0.041 0.041 
Testing −0.025 0.034 0.034 0.035 0.036 0.036 0.037 0.039 0.039 0.040 0.041 0.042 0.042 
−0.05 0.035 0.036 0.037 0.041 0.046 0.051 0.055 0.057 0.060 0.064 0.066 0.066 
−0.075 0.034 0.037 0.043 0.052 0.058 0.066 0.074 0.082 0.090 0.094 0.097 0.097 

Further, an anomaly data replacement treatment is proposed in the next section for better performance.

Performance theoretical analysis

In this section we analyze the DR and FAR. In single point mode, assuming the normal pattern p, current pattern , the residual and the scaling factor in Equation (1), and with event added we get current data pattern , residual and the scaling factor , if the simulated event has a strength s which lasts from th to th time step, then , and it can be derived that , where 
formula
So for the anomaly part with event added, , in which 
formula
if event duration is short, then which significantly exceeds the threshold because of the term s, then the DR is high; if the event duration is nearly as long as p, then , and which is small, then DR is low. Similar analysis can reveal that FAR has a positive correlation with event duration.

DISCUSSION

Although the proposed detection scheme is tested by several event scenes, it should be noted that the proposed detection scheme has some limitations. For example, it is suitable for free chlorine with a stable periodic pattern which is apparent compared with the noise data, and the detection results may be affected by the treatment method of the time series. Here, we discuss the problems of treatment of anomaly data and the detection scheme using the sub-pattern.

Unstable free chlorine data pattern

As we have discussed, the proposed method is effective for the stable free chlorine data pattern, so when the free chlorine has some anomalous fluctuation how does proposed method perform in detection? For example, in weekends the water demand is generally greater than weekdays, which may cause the data to increase to a higher value than in weekdays, or the chlorine input is changed deliberately due to some control actions which may cause it to become lower or higher.

For the first case, the water demand variation may cause the pattern extracted to be not suitable for everyday use in a week, but we can extract the patterns for weekdays and weekends, respectively, for more accurate detection. Furthermore, patterns extracted in different situations can be stored and a suitable pattern can be chosen from the stored patterns for anomaly detection. New patterns are added when dealing with new scenarios, which is the future work needed. For the second case, the control actions, including the tank operations and so on, usually have a regular operation time in a day and the chlorine would also show a regular change at that time, which can also be extracted into the pattern with a regular change, so the method will not treat that normal operation or control actions as an anomaly. Conversely, if the control action is not regular, some other network operation-based method can be combined in the detection to reduce false alarms, such as utilizing the tank level data and conductivity data (Hart et al. 2010a, b).

Overall we focus on proposing an idea and some testing cases for utilizing the fluctuation pattern of chlorine for anomaly detection. Further work for dealing with varying fluctuation patterns can be done based on our work.

Anomaly data treatment in single point detection

Anomaly data will make the residuals of the following points larger than expected in single point mode detection. We test a replacement method in which anomaly data are replaced with a value close to the normal pattern to reduce that bad impact. Let denote the data pattern at the th time step, then the current pattern after replacement is, 
formula
where is the scaling factor obtained in the previous detection, i.e., the anomaly data are replaced with normal pattern data multiplied by the scaling factor in the last detection, and denotes the anomaly declared at the th time step and denotes no anomaly declared. In Table 7 DR and FAR differences calculated from the test with replacement and without treatment are presented as DR (treatment)-DR (without treatment) and FAR (treatment)-FAR (without treatment), which shows that the DR is significantly improved with treatment especially when the duration is long and FAR is almost the same.
Table 7

DR and FAR difference between replacement treatment and without treatment

Event strength/duration41220283644526068768492
DR difference −0.025 −0.019 −0.008 0.003 0.015 0.026 0.035 0.043 0.050 0.056 0.061 0.067 0.072 
−0.05 0.007 0.014 0.024 0.039 0.062 0.093 0.118 0.139 0.161 0.180 0.197 0.213 
−0.075 −0.009 0.004 0.015 0.035 0.061 0.088 0.117 0.146 0.178 0.211 0.245 0.278 
FAR difference −0.025 0.017 0.017 0.018 0.018 0.018 0.017 0.016 0.014 0.012 0.011 0.012 0.013 
−0.05 0.016 0.018 0.020 0.017 0.014 0.010 0.006 0.004 0.001 0.001 0.001 0.002 
−0.075 0.016 0.014 0.011 0.005 −0.001 −0.005 −0.009 −0.012 −0.015 −0.016 −0.016 −0.016 
Event strength/duration41220283644526068768492
DR difference −0.025 −0.019 −0.008 0.003 0.015 0.026 0.035 0.043 0.050 0.056 0.061 0.067 0.072 
−0.05 0.007 0.014 0.024 0.039 0.062 0.093 0.118 0.139 0.161 0.180 0.197 0.213 
−0.075 −0.009 0.004 0.015 0.035 0.061 0.088 0.117 0.146 0.178 0.211 0.245 0.278 
FAR difference −0.025 0.017 0.017 0.018 0.018 0.018 0.017 0.016 0.014 0.012 0.011 0.012 0.013 
−0.05 0.016 0.018 0.020 0.017 0.014 0.010 0.006 0.004 0.001 0.001 0.001 0.002 
−0.075 0.016 0.014 0.011 0.005 −0.001 −0.005 −0.009 −0.012 −0.015 −0.016 −0.016 −0.016 

Anomaly detection in sub-pattern mode

A series of anomaly points may indicate a contamination event. Therefore, data series mode detection can provide an indication of a potential event, but one period (1 day) is too long to locate the event precisely. Detection using part of the period pattern may locate the event more precisely. Using sub-pattern the residuals r is 
formula
Sub-pattern is part of the whole pattern taken in a moving sliding window, corresponding to the starting time of the current pattern. Corresponding thresholds for sub-pattern are also required, which can be obtained similarly from the training data.

Tests are performed using 5 days of training data with events starting at the (96 + 50)th time step with event strength and length of 48 time steps. The result is shown in Figure 3, in which the 96, 48, 24 and 12 time steps long sub-pattern is used, respectively. With a shorter sub-pattern, the detection is more sensitive, but more false alarms are reported. In Figure 3, the event range is marked by a rectangle, and the stems represent the percentage of residuals obtained in detection over the threshold, that is. So the clusters of vertical lines within the event range are accurate detection of the anomaly and those outside the event range are false alarms. Owing to small scaling factor used in the threshold, the false alarms are frequent. In Figure 3(a) and (b) no accurate anomaly points are declared, and in Figure 3(c) and (d) the beginning points of the event are detected. So a shorter sub-pattern is more sensitive to the event, and more false alarms are reported. Another noteworthy result is that the data points immediately after the event are also declared as an anomaly, as the residuals of those points are also large for sudden change in shape.

Figure 3

Detection results of percentage of residual exceeding threshold using sub-pattern detection with pattern length: (a) 96, (b) 48, (c) 24, and (d) 12. A simulated event is added at the 146th time step with event strength −0.05 represented by a rectangle. Event is detected at the beginning in (c) (d), and some data points following the event are also declared as anomaly.

Figure 3

Detection results of percentage of residual exceeding threshold using sub-pattern detection with pattern length: (a) 96, (b) 48, (c) 24, and (d) 12. A simulated event is added at the 146th time step with event strength −0.05 represented by a rectangle. Event is detected at the beginning in (c) (d), and some data points following the event are also declared as anomaly.

CONCLUSION

In this article, we propose an anomaly detection scheme for free chlorine in drinking water network. The proposed detection scheme is based on the periodic fluctuation pattern of free chlorine, which is extracted using the SVD method, and the current data pattern is compared with the normal pattern in single point mode or data series mode for short-term or long-term anomaly inspection. A case study is presented using on-line data, with a 1-day-long periodic pattern obtained and event detection performance tests carried out. Performance tests show that the proposed method can detect the anomaly when the free chlorine data are still in the normal range of the national standard. Single point mode is sensitive to a short event. In addition, we presented the detection using a sub-pattern, which is better than single point mode and data series mode when the sub-pattern has a similar length to the event. Further study would need to be carried out to improve the detect methodology when the periodic pattern of free chlorine is not stable and more real contamination data should be analyzed to verify the detection performance in application.

ACKNOWLEDGEMENTS

This work was funded by the National Major Projects on Control and Rectification of Water-Body Pollution of China Granted by (2008ZX07420-004) ‘Research and Application of Water Quality Security Evaluation and Early-warning Technologies', the National Natural Science Foundation of China Grant (41101508) ‘Research on Water Quality Event Detection Methods based on Time-Frequency Analysis and Multi-sensor Data Fusion’, the Fundamental Research Funds for the Central Universities of China Grant (2014FZA5008) and the Program for Zhejiang Provincial Innovative Research Team (no. 2012R10037-08) ‘Research on Water Quality Early-warning Universal Methods and Technologies’.

REFERENCES

REFERENCES
Byer
D.
Carlson
K. H.
2005
Real-time detection of intentional chemical contamination in the distribution system
.
J. Am. Water Works Assoc.
97 (
7
),
130
133
.
Eliades
D. G.
Polycarpou
M. M.
2013
Contaminant detection in urban water distribution networks using chlorine measurements
. In:
Critical Information Infrastructures Security
,
Springer
,
Berlin, Heidelberg
, pp.
203
214
.
Hart
D. B.
McKenna
S. A.
Murray
R.
Haxton
T.
2010a
Combining water quality and operational data for improved event detection, paper presented at Water Distribution Systems Analysis (WDSA) Conference
.
Hart
D. B.
Hart
W. E.
McKenna
S.
Murray
R.
Phillips
C. A.
2010b
Event Detection System Operating Characteristics into Sensor Placement Optimization, paper presented at Water Distribution Systems Analysis 2010, ASCE
.
Hou
D.
Song
X.
Zhang
G.
Zhang
H.
Loaiciga
H.
2013
An early warning and control system for urban, drinking water quality protection: China's experience, Environ.
Sci. Pollut. R
2013
(
20
),
1
13
.
Hua
F.
West
J. R.
Barker
R. A.
Forster
C. F.
1999
Modelling of chlorine decay in municipal water supplies
.
Water Res.
33
(
12
),
2735
2746
.
Klise
K. A.
McKenna
S. A.
2006
Water quality change detection: multivariate algorithms, paper presented at Defense and Security Symposium, International Society for Optics and Photonics
.
Laird
C. D.
Biegler
L. T.
van Bloemen Waanders
B. G.
Bartlett
R. A.
2005
Contamination source determination for water networks
.
J. Water Resour. Plann. Manage
131
(
2
),
125
134
.
McKenna
S. A.
Wilson
M.
Klise
K. A.
2008
Detecting changes in water quality data
.
J. Am. Water Works Assoc.
100
(
1
),
74
85
.
Perelman
L.
Arad
J.
Housh
M.
Ostfeld
A.
2012
Event detection in water distribution systems from multivariate water quality time series
.
Environ. Sci. Technol.
 
46
(
15
),
8212
8219
.
Polycarpou
M. M.
Uber
J. G.
Wang
Z.
Shang
F.
Brdys
M.
2002
Feedback control of water quality
.
Control Sys., IEEE
22
(
3
),
68
87
.
Storey
M. V.
van der Gaag
B.
Burns
B. P.
2011
Advances in on-line drinking water quality monitoring and early warning systems
.
Water Res.
45
(
2
),
741
747
.
USEPA
2005
WaterSentinel Online Water Quality Monitoring as an Indicator of Drinking Water Contamination
,
Tech. Report No. EPA/817/D-05/002
,
Water Security Division, EPA
,
Washington, DC, USA
.
USEPA
2010
Water Quality Event Detection Systems for Drinking Water Contamination Warning Systems-Development, Testing, and Application of CANARY, Tech
.
Report No. EPA/600/R-10/036
,
National Homeland Security Research Center, Office of Research and Development, EPA
,
Cincinnati, OH
.
USEPA
2013
Water Quality Event Detection System Challenge: Methodology and Findings, Tech
.
Report No. EPA/817/R-13/002
,
Office of Water, EPA
,
Washington, DC, USA
.