# Hydrologic Time Series Anomaly Detection Based on Flink

Authors: Feng Ye; Zihao Liu; Qinghua Liu; Zhijian Wang
Journal: Mathematical Problems in Engineering
 (2020)

---

## Abstract

The data mining and calculation of time series in critical application is still worth studying. Currently, in the field of hydrological time series, most of the detection of outliers focus on improving the specificity. To efficiently detect outliers in massive hydrologic sensor data, an anomaly detection method for hydrological time series based on Flink is proposed. Firstly, the sliding window and the ARIMA model are used to forecast data stream. Then, the confidence interval is calculated for the prediction result, and the results outside the interval range are judged as alternative anomaly data. Finally, based on the historical batch data, theK-Means++ algorithm is used to cluster the batch data. The state transition probability is calculated, and the anomaly data are evaluated in quality. Taking the hydrological sensor data obtained from the Chu River as experimental data, experiments on the detection time and outlier detection performance are carried out, respectively. The results show that when calculating the tens of millions of data, the time costed by two slaves is less than that by one slave, and the maximum reduction is 17.43%. The sensitivity of the evaluation is increased from 72.91% to 92.98%. In terms of delay, the average delay of different slaves is roughly the same, which is maintained within 20 ms. It shows that, under big data platform, the proposed algorithm can effectively improve the computational efficiency of hydrologic time series detection for tens of millions of data and has a significant improvement in sensitivity.


---

## Body

## 1. Introduction

Hydrological data are divided into various types of hydrological time series according to their physical quantities. At present, many experts believe that hydrological time series is generally composed of determined and random components. The definite component has certain physical concept, and the random component is produced by the irregular oscillation and the stochastic influence [1]. Hydrological time series mainly shows the complex characteristics of randomness, fuzziness, nonlinearity, nonstationary, and multitime scale change [2].In practice, as the world gets more instrumented and connected, we are witnessing a flood of digital data generated from diversified hardware (e.g., sensors) or software in the format of big data. With the development of informatization, the hydrologic stations accumulate a great deal of important data which contains many outliers. For hydrological time series, it can be judged as an anomaly with a large difference from the general law [3]. Outliers often contain important information, and it is greatly significant for subsequent analysis decisions by accurately finding the hidden value behind the data. At present, for hydrological time series, traditional methods are only applicable to small datasets, not to the current big data environment. Moreover, the accuracy only reaches the level of 99% in specificity [4], and the sensitivity still has room for improvement. With the increase in the amount of data, how to calculate efficiently has become a problem that cannot be ignored. The anomaly detection and calculation of time series in critical application is still worth studying. This paper presents an anomaly detection method for hydrological time series based on Flink. Firstly, the sliding window and the ARIMA model are used to forecast data stream on the Flink platform. Then, the confidence interval is calculated for the prediction result, and the results outside the interval range are judged as temporary anomaly data. Finally, based on the historical batch data, the K-Means++ algorithm is used to cluster the batch data, the state transition probability is calculated, and the anomaly data are evaluated in quality. This method can effectively improve the computational efficiency and gives a reliable confidence degree to enhance the overall sensitivity. The outliers are detected quickly and accurately in massive hydrological time series.The following contents are organized as follows: Section2 discusses the research work related to this paper; Section 3 introduces the methodology of hydrologic time series anomaly detection in detail; in Section 4, we continue to use the real hydrological sensor data as experimental data to verify the effectiveness of the proposed method. Finally, the summary and prospect are given.


## 2. Related Work

### 2.1. Anomaly Detection

Outliers [5] are data that deviate most of the data in the dataset, which is not suspected of being a random error, but arises from a completely different mechanism. The following are some of the main methods of anomaly detection.Niu [6] proposes a short-term electricity price hybrid forecasting model based on wavelet transform and ARIMA, which can detect the mutation point. The model can indeed detect the situation of abrupt change point, but it is insufficient for the nonlinear part or the data with too long time series. Gil [7] presents an anomaly detection based on vector machine and principal element analysis. Firstly, the principal element analysis method is used to reduce the latitude, and then, the SVM is used to model and examine the anomaly data. But if there are more kinds of outliers, the detection accuracy is not good. Sun [3] proposes a hydrological time series anomaly value detection based on ARIMA-SVR, which uses ARIMA to predict the linear parts, and SVM predicts the nonlinear part. It adds the predicted result and evaluates the value not in the confidence interval as the anomaly value. This kind of algorithm has a good effect on small-scale datasets, but it is not able to deal with massive data or multivariate data. It is also difficult to determine the threshold value.The method based on distance detection is to set some distance function to calculate the distance of the data point. When the distance between one point and the other point is too large, it is regarded as the anomaly point. Vy and Anh [8] present an anomaly detection algorithm with variable length in time series. Firstly, it segments a time series and then calculates the anomaly factors of each mode. After that, the distance between them is calculated. At last, the anomaly is judged by the distance of the abnormal factors. The advantage of this method is that it is easy to use; the time complexity is relatively small, but it is not sensitive to local anomaly points.Ali et al. [9] propose the concept of local outlier factor (LOF) to compute the dataset density. The possibility that the object is an outlier is positively correlated with the LOF. However, the mix of different densities will result in the detection error. Although some related improvement schemes are proposed, the overall time complexity is higher.The clustering algorithm [10] divides the points in the time series into several clusters, and these points which do not belong to any cluster will be regarded as the anomaly, but the time series has a trend characteristic and cannot be classified as cluster analysis simply. Therefore, the clustering algorithm is too dependent on cluster quality, which leads to low accuracy and efficiency.Hypothesis testing is a method for discovering anomaly samples. The dataset obeys a known distribution or probability model. If a point in the dataset is inconsistent with its distribution, the anomaly is judged. Twitter opens source the traffic anomaly detection algorithm S-H-ESD in 2015 [11], and the algorithm is to use STL to decompose the sequence and investigate the residuals. Assuming this is in accordance with the normal distribution, the outlier can be extracted using generalized ESD. However, if the feature distribution is unknown, the priori hypothesis does not necessarily effect. Then, the error detection rate of this method is high, and it cannot adapt to the multivariate time series well.Yu [4] proposes the hydrological time series anomaly detection based on sliding window prediction, but the computational complexity is high. Yang et al. [12] use the knowledge granularity method to find the abnormal data in time series, reduce the time cost of the detection process, and improve the detection efficiency. Liu and Wang [13] propose an anomaly factor detection method based on the extremum difference, the slope, and the mean value. Good results are obtained.In recent years, the accuracy of model prediction is also the goal of many researchers, and Zeng et al. [14] propose fuzzy forecasting based on linear combinations of independent variables, subtractive clustering algorithm, and artificial bee colony algorithm to predict time series. The proposed method gets higher forecasting accuracy rates than the existing methods for forecasting the Taiwan Stock Exchange Capitalization Weighted Stock Index (TAIEX), the enrolment of the University of Alabama, and the daily percentage of CO2. In addition, Zeng et al. [15] still propose interval-valued intuitionistic fuzzy multiple attribute decision-making based on nonlinear programming methodology and TOPSIS method. The method can overcome the drawbacks of the MADM methods. It offers us a very useful way to deal with MADM problems in IVIF environments. Although these two methods have extremely high accuracy, the computational efficiency in the big data is still a problem worth studying.


### 2.2. Flink

Apache Flink is a framework and distributed processing engine for stateful computation of unbounded and bounded data flows. Flink is designed to run in all common clustered environments, performing calculations at memory speed and any size. The Flink centralized swarm mode deployment on yarn is undergone through the construction of four centralized swarms in Figure1.Figure 1
Flink architecture.Internally, Apache Flink represents job definitions using directed acyclic graphs (DAGs) [16]. The nodes of the graph are either sources, sinks, or operators. Source nodes read in or generate the input data, while sink nodes produce the output. The inner vertices are operators which execute arbitrary user-defined functions (UDFs) that consume input from incident nodes and provide input for adjacent nodes. The DAGs generated from the user’s job definitions are then transformed into the more concrete execution graphs, which contain the necessary information for running the job on a cluster. Data partitioning enables the data-parallel execution of the subtasks. During the transformation, the UDFs are split up into multiple parallel subtasks. Each subtask executes the same UDF. However, each process processes a different part of the input data. This setup makes clear where overheads might emerge. Firstly, the task scheduling as well as the deployment of the tasks on the respective machines introduce overhead. Moreover, the nodes must communicate with each other in order to distribute the workload. The communication between different nodes demand for serialization and transport buffering, which adds to the overhead.On top of the possibility to define a job by using a DAG, a layer of second-order functions is implemented in order to simplify the development for the user [17]. Apache Flink provides APIs for implementing batch as well as stream processing. The exposed functions for dataset and data stream transformations are like functions known from functional programming (e.g., map, reduce, and filter). Additionally, the DataSet API provides transformations known from relational databases like joins and grouping. The DataStream API provides additional operators which are useful in the streaming context. These operators include the definition of windows and window-based aggregations.


## 2.1. Anomaly Detection

Outliers [5] are data that deviate most of the data in the dataset, which is not suspected of being a random error, but arises from a completely different mechanism. The following are some of the main methods of anomaly detection.Niu [6] proposes a short-term electricity price hybrid forecasting model based on wavelet transform and ARIMA, which can detect the mutation point. The model can indeed detect the situation of abrupt change point, but it is insufficient for the nonlinear part or the data with too long time series. Gil [7] presents an anomaly detection based on vector machine and principal element analysis. Firstly, the principal element analysis method is used to reduce the latitude, and then, the SVM is used to model and examine the anomaly data. But if there are more kinds of outliers, the detection accuracy is not good. Sun [3] proposes a hydrological time series anomaly value detection based on ARIMA-SVR, which uses ARIMA to predict the linear parts, and SVM predicts the nonlinear part. It adds the predicted result and evaluates the value not in the confidence interval as the anomaly value. This kind of algorithm has a good effect on small-scale datasets, but it is not able to deal with massive data or multivariate data. It is also difficult to determine the threshold value.The method based on distance detection is to set some distance function to calculate the distance of the data point. When the distance between one point and the other point is too large, it is regarded as the anomaly point. Vy and Anh [8] present an anomaly detection algorithm with variable length in time series. Firstly, it segments a time series and then calculates the anomaly factors of each mode. After that, the distance between them is calculated. At last, the anomaly is judged by the distance of the abnormal factors. The advantage of this method is that it is easy to use; the time complexity is relatively small, but it is not sensitive to local anomaly points.Ali et al. [9] propose the concept of local outlier factor (LOF) to compute the dataset density. The possibility that the object is an outlier is positively correlated with the LOF. However, the mix of different densities will result in the detection error. Although some related improvement schemes are proposed, the overall time complexity is higher.The clustering algorithm [10] divides the points in the time series into several clusters, and these points which do not belong to any cluster will be regarded as the anomaly, but the time series has a trend characteristic and cannot be classified as cluster analysis simply. Therefore, the clustering algorithm is too dependent on cluster quality, which leads to low accuracy and efficiency.Hypothesis testing is a method for discovering anomaly samples. The dataset obeys a known distribution or probability model. If a point in the dataset is inconsistent with its distribution, the anomaly is judged. Twitter opens source the traffic anomaly detection algorithm S-H-ESD in 2015 [11], and the algorithm is to use STL to decompose the sequence and investigate the residuals. Assuming this is in accordance with the normal distribution, the outlier can be extracted using generalized ESD. However, if the feature distribution is unknown, the priori hypothesis does not necessarily effect. Then, the error detection rate of this method is high, and it cannot adapt to the multivariate time series well.Yu [4] proposes the hydrological time series anomaly detection based on sliding window prediction, but the computational complexity is high. Yang et al. [12] use the knowledge granularity method to find the abnormal data in time series, reduce the time cost of the detection process, and improve the detection efficiency. Liu and Wang [13] propose an anomaly factor detection method based on the extremum difference, the slope, and the mean value. Good results are obtained.In recent years, the accuracy of model prediction is also the goal of many researchers, and Zeng et al. [14] propose fuzzy forecasting based on linear combinations of independent variables, subtractive clustering algorithm, and artificial bee colony algorithm to predict time series. The proposed method gets higher forecasting accuracy rates than the existing methods for forecasting the Taiwan Stock Exchange Capitalization Weighted Stock Index (TAIEX), the enrolment of the University of Alabama, and the daily percentage of CO2. In addition, Zeng et al. [15] still propose interval-valued intuitionistic fuzzy multiple attribute decision-making based on nonlinear programming methodology and TOPSIS method. The method can overcome the drawbacks of the MADM methods. It offers us a very useful way to deal with MADM problems in IVIF environments. Although these two methods have extremely high accuracy, the computational efficiency in the big data is still a problem worth studying.


## 2.2. Flink

Apache Flink is a framework and distributed processing engine for stateful computation of unbounded and bounded data flows. Flink is designed to run in all common clustered environments, performing calculations at memory speed and any size. The Flink centralized swarm mode deployment on yarn is undergone through the construction of four centralized swarms in Figure1.Figure 1
Flink architecture.Internally, Apache Flink represents job definitions using directed acyclic graphs (DAGs) [16]. The nodes of the graph are either sources, sinks, or operators. Source nodes read in or generate the input data, while sink nodes produce the output. The inner vertices are operators which execute arbitrary user-defined functions (UDFs) that consume input from incident nodes and provide input for adjacent nodes. The DAGs generated from the user’s job definitions are then transformed into the more concrete execution graphs, which contain the necessary information for running the job on a cluster. Data partitioning enables the data-parallel execution of the subtasks. During the transformation, the UDFs are split up into multiple parallel subtasks. Each subtask executes the same UDF. However, each process processes a different part of the input data. This setup makes clear where overheads might emerge. Firstly, the task scheduling as well as the deployment of the tasks on the respective machines introduce overhead. Moreover, the nodes must communicate with each other in order to distribute the workload. The communication between different nodes demand for serialization and transport buffering, which adds to the overhead.On top of the possibility to define a job by using a DAG, a layer of second-order functions is implemented in order to simplify the development for the user [17]. Apache Flink provides APIs for implementing batch as well as stream processing. The exposed functions for dataset and data stream transformations are like functions known from functional programming (e.g., map, reduce, and filter). Additionally, the DataSet API provides transformations known from relational databases like joins and grouping. The DataStream API provides additional operators which are useful in the streaming context. These operators include the definition of windows and window-based aggregations.


## 3. The Proposed Methodology

### 3.1. Anomaly Detection Based on Sliding Window

#### 3.1.1. The Definition of Slide Window

To define the sliding neighbor windowLi of the hydrological time series X to be detected point Xi. To reduce the complexity of the algorithm, the first L point of the point is used as the input parameter of the predictive model. This algorithm selects the left neighbor window of the prediction node as the algorithm input, and the unilateral definition is as follows:(1)Lxi=xi−L,xi−L+1,…,xi−1.


#### 3.1.2. Anomaly Detection

The core is to establish the ARIMA model [16] and predict the value of the observation point by sliding window input and then get a series of predicted values. ARIMA is a time series prediction model with autoregressive (AR) and moving average (MA). It is very flexible. The model becomes MAq when AR equals 0, while when MA = 0, ARIMA becomes ARp and specifically is described as(2)yt=v+φ1yt−1+⋯+φpyt−p+εt,n≤m.In (2), εt is the error. If εt is autocorrelation, then MA (q) can be represented as(3)εt=at+Θ1at−1+⋯+Θ1at−p.In (3), ΘJj=1,2,…,q is the parameters estimated; at is the white noise.According to the above,ARIMAp,q is(4)yt=v+φ1yt−1+⋯+φpyt−p+εt+v+φ1yt−1+⋯+φpyt−p+at+Θ1at−1+⋯+Θpat−q.If the ordern is relatively large, then ARn can be approximately equivalent to ARIMAp,q; then,(5)yt=v+∑i=1nφiyt−i+ant.In equation (5), ant is the error term of order n. The estimate of ant can be achieved as follows:(6)a^nt=yt+∑i=1nφ^iyt−i−v.In equation (6), φ^i can be achieved by least square estimation. Making use of a^nt, ARIMAp,q can be established:(7)yt=v,φ1,…,φq:Θ1,…,Θq1yt−1…a^nt−1…a^nt−q+at.ARIMA parametersn, p, and q are determined by AIC criterion:(8)AIC=ln∑ap,q+2SPnum.In (8), S is the number of samples, and ∑ap,q is the determinant of Σa, which is the covariance matrix of a^t.First, unit root test is conducted for the time series . If it is a nonstationary sequence, it needs to be transformed into a stationary sequence by difference. Based on the AIC criterion, we need to determine the autoregressive orderp and the moving average order q and find the p and q combination with the minimum AIC value. The ARIMA model is for nonstationary time series. It is applicable to hydrological time series. This paper takes the confidence interval of 95%. The confidence interval calculated by the ARIMA model is compared with the original sequence, and the value not in the confidence interval is judged as outliers.


### 3.2. Anomaly Value Verification

After the exception values are identified by sliding window and ARIMA model, we also need to identify the confidence of the anomaly to determine whether the point is indeed an anomaly and reduce the error and the amount of work.


#### 3.2.1.K-Mean++ Model

K-Mean is one of the clustering algorithms. The principle of the K-Means++ algorithm is given the value of K, K represents the number of categories to divide the data into and then according to the similarity between the data to divide the data into K classes. The method of measuring the similarity of data is usually measured by the distance between data points, such as European distance, Hamming distance, and Manhattan distance.It is a common practice to use the Euclidean distance to measure the similarity between the data. For example, for two points on the two-dimensional planeAx1,y1 and Bx2,y2, the Euclidean distance between the two is(9)x1−x22+y1−y22.Generally, a cluster is described by the centre of all points in the cluster, which is also called centroid. The way of computing the centre of mass is to calculate the mean of all data points in a cluster. The advantage ofK-Means++ lies in the choice of the centre of mass. The K-Means++ model is optimized for the centroid uncertainty in the K-Means algorithm and chose the centroid K in the following strategy: hypothesis has selected the initial clustering centre n (0 < n < K). When the clustering centre is in the selection of n+1, the farther the distance from the current n clustering centre point will be a higher probability of being selected as the clustering centre n+1. However, the selection of the first cluster centre (n = 1) also adopts the random method.Specifically,K-Means takes the entire time series x1,x2,… as input and sequence T=T1,T2,… as an output; the points on the time series are converted to various cluster points, and T1,T2,… indicates which centre this time point is in, which classify the time series in different states.After providing theK-means++ model with the previous outliers and the value of their previous moment as input, the distance between each sample and the centre of the cluster is computed, and then the cluster centre is assigned to each sample.


#### 3.2.2. State Transition Probability Matrix

Markov processes in which time and state are discrete are called Markov chains, and it is assumed that(10)Xn=Xn,n=1,2,….It can be viewed as the result of successive observations of discrete Markov processes on the time set:(11)T1=0,1,2,….The state space of the chain can be described that(12)I1=a1,a2,a3,….In the state of chain, Markov property is usually expressed by conditional distribution law. That is, satisfy(13)∀n∈Z,r∈Z,0≤t1<t2<⋯<tr<m;tj,m,n+m∈T1.(14)PXm+n=ajXt1=ai1,Xt2=ai2,…,Xtr=air,Xm=ai=PXm+n=ajXm=ai,ai∈I.We sign equation (14) as follows:(15)Pijm,m+n,and describe the conditional probability as(16)Pm,m+n=PXm+n=ajXm=ai,ai∈I.It means the transition probability of Markov chain fromm+n to aj under the condition that m is in the state ai at time m. Since the chain starts from any state ai at time m to another state m+n, it must move to one of the states in a1,a2,…, so(17)∑j=1∞Pijm,m+n=1,i=1,2,….The matrix composed of transition probabilities is called the transition probability matrix of Markov chain:(18)Pm,m+n=Pijm,m+n,and the above formula (18) knows that the sum of the elements in each row of this matrix is equal to 1. When the transition probability is only related to Pijm,m+n, i, j, and time interval n, it is denoted as Pij or(19)Pm,m+n=Pijn.The transition probability is said to be stable. It is also homogeneous or time-homogeneous. In the case of homogeneous Markov chain, the transition probability is(20)Pijn=PXm+n=ajXm=ai.Equation (20) is called the n-step transition probability of Markov chains, and when n=1, it is the one-step transition probability, which is particularly important. The matrix consisting of the one-step transition probability is called the one-step transition probability matrix:(21)P=p1,1p1,1…p1,1…p1,1p1,1…p1,1…MMOMOp1,1p1,1…p1,1…MMOMO,pi,j=prji.The following is the realization method of one-step transition probability matrix:Step 1. Calculate the probability of each state:(22)PXm=ai.Step 2. Calculate the probability of stateai at moment m and state aj at moment m+1:(23)PXm+n=aj,Xm=ai.Step 3. Compute transition probability:(24)Pij=Pij1=PXm+1=ajXm=ai=PXm+n=aj,Xm=aiPXm=ai.Means model training is completed, the state sequenceT1,T2,… can be obtained, and a state transition matrix can be obtained by calculation. To facilitate the calculation, it is necessary to convert the matrix to a data frame. There are three columns: the first column represents the state i, the second column represents the state j, and the third column represents the probability that the state i transfers to state j. Suppose a time series xi,xi+1 is transformed by K-Means++ model and it becomes the state sequence Ti,Tj, xi+1 appears after xi, which in other words is equivalent to a transfer from state Ti to Tj, and the transfer probability is(25)pij=the number ofTitoTjthe number ofTito other state.


#### 3.2.3. Outliers Check

In this paper, the probability of a transition is compared with the probability of the state of the previous moment being transferred to the next most probable state as an evaluation criterion. Ifxi is to be inspected anomaly, the value of the previous moment is xi−1, and the state of Ti and Ti−1, Ti−1 is the most likely to be transferred to Tm; thus, we define an anomaly value probability:(26)pi=1−the possibility thatTi−1toTithe possibility thatTi−1toTm.From the above, it is known that the probability of stateTi to Tm is a constant, the smaller the probability of state Ti−1 transferring to Ti, the greater the pi, the greater the probability of xi being abnormal value.In the reality of hydrological data anomaly judgement standard, the data more than 2 cm are judged as outliers. Using the clustering algorithm may appear abnormal values and the moment before its values in the same condition and lead to false negatives. So after the judgement, outliers with a probability of 0 but a difference of more than 2 cm should also be judged as outliers.


### 3.3. The Anomaly Detection Method

The hydrological time series anomaly detection algorithm based on Flink is combined with two processes: the anomaly detection and result verification. First, the ARIMA model of time seriesx1,x2,…,xn is established by the idea of predictive detection by sliding window, and the predicted confidence interval is obtained, which is compared with the original data to identify the outliers. After the anomaly value is detected, the original data are clustered by K-Means++ algorithm, and the state transition matrix is computed after clustering. At last, the abnormal value is evaluated by the state transition and finally determined.The specific steps of the algorithm are provided in Algorithm1.Algorithm 1: Hydrological time series anomaly detection algorithm.
Input: hydrological time seriesX, reliability P, sliding window size L, historical batch data H.Output: the outliers in hydrological time series.Step 1: clean the sequenceH, including descending dimension, deleting duplicate value, sifting, and sortingStep 2: using the value ofL as the initial starting position of the sliding window of X, the value of the xL+1 is predicted, and as the window slides, the predicted value gradually forms a new time series xL+1,xL+2,…Step 3: the 95% confidence interval of the new time series is calculated and compared withX, the time point which is not in the confidence interval is obtained and get the exception point set e1,e2,…Step 4: taking historical dataH as input and training and establishing K-Means++ model, obtain the discrete state sequence T1,T2,…Step 5: compute the state transition matrix of the discrete state sequenceStep 6: theK-Mean++ model that is obtained in Step 4 is used for the exception point set of Step 3 and the value of its previous moment to obtain the state of anomaly and its previous momentStep 7: estimate the value of the exception and its previous moment in Step 6 by state transition data frame and then output the confidence scoreStep 8: repeat the above steps until no new data are enteredThe flow chart is shown below.Figure2 shows the overall flow of the algorithm. Initial abnormal detection is carried out on the input data stream through the sliding window in conjunction with ARIMA model.Figure 2
The whole process of anomaly detection algorithm.Figure3 shows the anomaly verification mechanism of specific process, and batch the historical data collected before. The first is to delete duplicates, dimension reduction, and sorting operation on batch data. After that, the data are clustered by the K-Means++ model. Then, original time series are converted into a Markov chain by using the clustering results. Finally, one step transfer matrix and the transition probability can be calculated. Through the principle of maximum and minimum, the hydrological sequence whose anomaly probability is higher than the threshold value is judged as the true outlier.Figure 3
Abnormal check mechanism flow chart.


## 3.1. Anomaly Detection Based on Sliding Window

### 3.1.1. The Definition of Slide Window

To define the sliding neighbor windowLi of the hydrological time series X to be detected point Xi. To reduce the complexity of the algorithm, the first L point of the point is used as the input parameter of the predictive model. This algorithm selects the left neighbor window of the prediction node as the algorithm input, and the unilateral definition is as follows:(1)Lxi=xi−L,xi−L+1,…,xi−1.


### 3.1.2. Anomaly Detection

The core is to establish the ARIMA model [16] and predict the value of the observation point by sliding window input and then get a series of predicted values. ARIMA is a time series prediction model with autoregressive (AR) and moving average (MA). It is very flexible. The model becomes MAq when AR equals 0, while when MA = 0, ARIMA becomes ARp and specifically is described as(2)yt=v+φ1yt−1+⋯+φpyt−p+εt,n≤m.In (2), εt is the error. If εt is autocorrelation, then MA (q) can be represented as(3)εt=at+Θ1at−1+⋯+Θ1at−p.In (3), ΘJj=1,2,…,q is the parameters estimated; at is the white noise.According to the above,ARIMAp,q is(4)yt=v+φ1yt−1+⋯+φpyt−p+εt+v+φ1yt−1+⋯+φpyt−p+at+Θ1at−1+⋯+Θpat−q.If the ordern is relatively large, then ARn can be approximately equivalent to ARIMAp,q; then,(5)yt=v+∑i=1nφiyt−i+ant.In equation (5), ant is the error term of order n. The estimate of ant can be achieved as follows:(6)a^nt=yt+∑i=1nφ^iyt−i−v.In equation (6), φ^i can be achieved by least square estimation. Making use of a^nt, ARIMAp,q can be established:(7)yt=v,φ1,…,φq:Θ1,…,Θq1yt−1…a^nt−1…a^nt−q+at.ARIMA parametersn, p, and q are determined by AIC criterion:(8)AIC=ln∑ap,q+2SPnum.In (8), S is the number of samples, and ∑ap,q is the determinant of Σa, which is the covariance matrix of a^t.First, unit root test is conducted for the time series . If it is a nonstationary sequence, it needs to be transformed into a stationary sequence by difference. Based on the AIC criterion, we need to determine the autoregressive orderp and the moving average order q and find the p and q combination with the minimum AIC value. The ARIMA model is for nonstationary time series. It is applicable to hydrological time series. This paper takes the confidence interval of 95%. The confidence interval calculated by the ARIMA model is compared with the original sequence, and the value not in the confidence interval is judged as outliers.


## 3.1.1. The Definition of Slide Window

To define the sliding neighbor windowLi of the hydrological time series X to be detected point Xi. To reduce the complexity of the algorithm, the first L point of the point is used as the input parameter of the predictive model. This algorithm selects the left neighbor window of the prediction node as the algorithm input, and the unilateral definition is as follows:(1)Lxi=xi−L,xi−L+1,…,xi−1.


## 3.1.2. Anomaly Detection

The core is to establish the ARIMA model [16] and predict the value of the observation point by sliding window input and then get a series of predicted values. ARIMA is a time series prediction model with autoregressive (AR) and moving average (MA). It is very flexible. The model becomes MAq when AR equals 0, while when MA = 0, ARIMA becomes ARp and specifically is described as(2)yt=v+φ1yt−1+⋯+φpyt−p+εt,n≤m.In (2), εt is the error. If εt is autocorrelation, then MA (q) can be represented as(3)εt=at+Θ1at−1+⋯+Θ1at−p.In (3), ΘJj=1,2,…,q is the parameters estimated; at is the white noise.According to the above,ARIMAp,q is(4)yt=v+φ1yt−1+⋯+φpyt−p+εt+v+φ1yt−1+⋯+φpyt−p+at+Θ1at−1+⋯+Θpat−q.If the ordern is relatively large, then ARn can be approximately equivalent to ARIMAp,q; then,(5)yt=v+∑i=1nφiyt−i+ant.In equation (5), ant is the error term of order n. The estimate of ant can be achieved as follows:(6)a^nt=yt+∑i=1nφ^iyt−i−v.In equation (6), φ^i can be achieved by least square estimation. Making use of a^nt, ARIMAp,q can be established:(7)yt=v,φ1,…,φq:Θ1,…,Θq1yt−1…a^nt−1…a^nt−q+at.ARIMA parametersn, p, and q are determined by AIC criterion:(8)AIC=ln∑ap,q+2SPnum.In (8), S is the number of samples, and ∑ap,q is the determinant of Σa, which is the covariance matrix of a^t.First, unit root test is conducted for the time series . If it is a nonstationary sequence, it needs to be transformed into a stationary sequence by difference. Based on the AIC criterion, we need to determine the autoregressive orderp and the moving average order q and find the p and q combination with the minimum AIC value. The ARIMA model is for nonstationary time series. It is applicable to hydrological time series. This paper takes the confidence interval of 95%. The confidence interval calculated by the ARIMA model is compared with the original sequence, and the value not in the confidence interval is judged as outliers.


## 3.2. Anomaly Value Verification

After the exception values are identified by sliding window and ARIMA model, we also need to identify the confidence of the anomaly to determine whether the point is indeed an anomaly and reduce the error and the amount of work.


### 3.2.1.K-Mean++ Model

K-Mean is one of the clustering algorithms. The principle of the K-Means++ algorithm is given the value of K, K represents the number of categories to divide the data into and then according to the similarity between the data to divide the data into K classes. The method of measuring the similarity of data is usually measured by the distance between data points, such as European distance, Hamming distance, and Manhattan distance.It is a common practice to use the Euclidean distance to measure the similarity between the data. For example, for two points on the two-dimensional planeAx1,y1 and Bx2,y2, the Euclidean distance between the two is(9)x1−x22+y1−y22.Generally, a cluster is described by the centre of all points in the cluster, which is also called centroid. The way of computing the centre of mass is to calculate the mean of all data points in a cluster. The advantage ofK-Means++ lies in the choice of the centre of mass. The K-Means++ model is optimized for the centroid uncertainty in the K-Means algorithm and chose the centroid K in the following strategy: hypothesis has selected the initial clustering centre n (0 < n < K). When the clustering centre is in the selection of n+1, the farther the distance from the current n clustering centre point will be a higher probability of being selected as the clustering centre n+1. However, the selection of the first cluster centre (n = 1) also adopts the random method.Specifically,K-Means takes the entire time series x1,x2,… as input and sequence T=T1,T2,… as an output; the points on the time series are converted to various cluster points, and T1,T2,… indicates which centre this time point is in, which classify the time series in different states.After providing theK-means++ model with the previous outliers and the value of their previous moment as input, the distance between each sample and the centre of the cluster is computed, and then the cluster centre is assigned to each sample.


### 3.2.2. State Transition Probability Matrix

Markov processes in which time and state are discrete are called Markov chains, and it is assumed that(10)Xn=Xn,n=1,2,….It can be viewed as the result of successive observations of discrete Markov processes on the time set:(11)T1=0,1,2,….The state space of the chain can be described that(12)I1=a1,a2,a3,….In the state of chain, Markov property is usually expressed by conditional distribution law. That is, satisfy(13)∀n∈Z,r∈Z,0≤t1<t2<⋯<tr<m;tj,m,n+m∈T1.(14)PXm+n=ajXt1=ai1,Xt2=ai2,…,Xtr=air,Xm=ai=PXm+n=ajXm=ai,ai∈I.We sign equation (14) as follows:(15)Pijm,m+n,and describe the conditional probability as(16)Pm,m+n=PXm+n=ajXm=ai,ai∈I.It means the transition probability of Markov chain fromm+n to aj under the condition that m is in the state ai at time m. Since the chain starts from any state ai at time m to another state m+n, it must move to one of the states in a1,a2,…, so(17)∑j=1∞Pijm,m+n=1,i=1,2,….The matrix composed of transition probabilities is called the transition probability matrix of Markov chain:(18)Pm,m+n=Pijm,m+n,and the above formula (18) knows that the sum of the elements in each row of this matrix is equal to 1. When the transition probability is only related to Pijm,m+n, i, j, and time interval n, it is denoted as Pij or(19)Pm,m+n=Pijn.The transition probability is said to be stable. It is also homogeneous or time-homogeneous. In the case of homogeneous Markov chain, the transition probability is(20)Pijn=PXm+n=ajXm=ai.Equation (20) is called the n-step transition probability of Markov chains, and when n=1, it is the one-step transition probability, which is particularly important. The matrix consisting of the one-step transition probability is called the one-step transition probability matrix:(21)P=p1,1p1,1…p1,1…p1,1p1,1…p1,1…MMOMOp1,1p1,1…p1,1…MMOMO,pi,j=prji.The following is the realization method of one-step transition probability matrix:Step 1. Calculate the probability of each state:(22)PXm=ai.Step 2. Calculate the probability of stateai at moment m and state aj at moment m+1:(23)PXm+n=aj,Xm=ai.Step 3. Compute transition probability:(24)Pij=Pij1=PXm+1=ajXm=ai=PXm+n=aj,Xm=aiPXm=ai.Means model training is completed, the state sequenceT1,T2,… can be obtained, and a state transition matrix can be obtained by calculation. To facilitate the calculation, it is necessary to convert the matrix to a data frame. There are three columns: the first column represents the state i, the second column represents the state j, and the third column represents the probability that the state i transfers to state j. Suppose a time series xi,xi+1 is transformed by K-Means++ model and it becomes the state sequence Ti,Tj, xi+1 appears after xi, which in other words is equivalent to a transfer from state Ti to Tj, and the transfer probability is(25)pij=the number ofTitoTjthe number ofTito other state.


### 3.2.3. Outliers Check

In this paper, the probability of a transition is compared with the probability of the state of the previous moment being transferred to the next most probable state as an evaluation criterion. Ifxi is to be inspected anomaly, the value of the previous moment is xi−1, and the state of Ti and Ti−1, Ti−1 is the most likely to be transferred to Tm; thus, we define an anomaly value probability:(26)pi=1−the possibility thatTi−1toTithe possibility thatTi−1toTm.From the above, it is known that the probability of stateTi to Tm is a constant, the smaller the probability of state Ti−1 transferring to Ti, the greater the pi, the greater the probability of xi being abnormal value.In the reality of hydrological data anomaly judgement standard, the data more than 2 cm are judged as outliers. Using the clustering algorithm may appear abnormal values and the moment before its values in the same condition and lead to false negatives. So after the judgement, outliers with a probability of 0 but a difference of more than 2 cm should also be judged as outliers.


## 3.2.1.K-Mean++ Model

K-Mean is one of the clustering algorithms. The principle of the K-Means++ algorithm is given the value of K, K represents the number of categories to divide the data into and then according to the similarity between the data to divide the data into K classes. The method of measuring the similarity of data is usually measured by the distance between data points, such as European distance, Hamming distance, and Manhattan distance.It is a common practice to use the Euclidean distance to measure the similarity between the data. For example, for two points on the two-dimensional planeAx1,y1 and Bx2,y2, the Euclidean distance between the two is(9)x1−x22+y1−y22.Generally, a cluster is described by the centre of all points in the cluster, which is also called centroid. The way of computing the centre of mass is to calculate the mean of all data points in a cluster. The advantage ofK-Means++ lies in the choice of the centre of mass. The K-Means++ model is optimized for the centroid uncertainty in the K-Means algorithm and chose the centroid K in the following strategy: hypothesis has selected the initial clustering centre n (0 < n < K). When the clustering centre is in the selection of n+1, the farther the distance from the current n clustering centre point will be a higher probability of being selected as the clustering centre n+1. However, the selection of the first cluster centre (n = 1) also adopts the random method.Specifically,K-Means takes the entire time series x1,x2,… as input and sequence T=T1,T2,… as an output; the points on the time series are converted to various cluster points, and T1,T2,… indicates which centre this time point is in, which classify the time series in different states.After providing theK-means++ model with the previous outliers and the value of their previous moment as input, the distance between each sample and the centre of the cluster is computed, and then the cluster centre is assigned to each sample.


## 3.2.2. State Transition Probability Matrix

Markov processes in which time and state are discrete are called Markov chains, and it is assumed that(10)Xn=Xn,n=1,2,….It can be viewed as the result of successive observations of discrete Markov processes on the time set:(11)T1=0,1,2,….The state space of the chain can be described that(12)I1=a1,a2,a3,….In the state of chain, Markov property is usually expressed by conditional distribution law. That is, satisfy(13)∀n∈Z,r∈Z,0≤t1<t2<⋯<tr<m;tj,m,n+m∈T1.(14)PXm+n=ajXt1=ai1,Xt2=ai2,…,Xtr=air,Xm=ai=PXm+n=ajXm=ai,ai∈I.We sign equation (14) as follows:(15)Pijm,m+n,and describe the conditional probability as(16)Pm,m+n=PXm+n=ajXm=ai,ai∈I.It means the transition probability of Markov chain fromm+n to aj under the condition that m is in the state ai at time m. Since the chain starts from any state ai at time m to another state m+n, it must move to one of the states in a1,a2,…, so(17)∑j=1∞Pijm,m+n=1,i=1,2,….The matrix composed of transition probabilities is called the transition probability matrix of Markov chain:(18)Pm,m+n=Pijm,m+n,and the above formula (18) knows that the sum of the elements in each row of this matrix is equal to 1. When the transition probability is only related to Pijm,m+n, i, j, and time interval n, it is denoted as Pij or(19)Pm,m+n=Pijn.The transition probability is said to be stable. It is also homogeneous or time-homogeneous. In the case of homogeneous Markov chain, the transition probability is(20)Pijn=PXm+n=ajXm=ai.Equation (20) is called the n-step transition probability of Markov chains, and when n=1, it is the one-step transition probability, which is particularly important. The matrix consisting of the one-step transition probability is called the one-step transition probability matrix:(21)P=p1,1p1,1…p1,1…p1,1p1,1…p1,1…MMOMOp1,1p1,1…p1,1…MMOMO,pi,j=prji.The following is the realization method of one-step transition probability matrix:Step 1. Calculate the probability of each state:(22)PXm=ai.Step 2. Calculate the probability of stateai at moment m and state aj at moment m+1:(23)PXm+n=aj,Xm=ai.Step 3. Compute transition probability:(24)Pij=Pij1=PXm+1=ajXm=ai=PXm+n=aj,Xm=aiPXm=ai.Means model training is completed, the state sequenceT1,T2,… can be obtained, and a state transition matrix can be obtained by calculation. To facilitate the calculation, it is necessary to convert the matrix to a data frame. There are three columns: the first column represents the state i, the second column represents the state j, and the third column represents the probability that the state i transfers to state j. Suppose a time series xi,xi+1 is transformed by K-Means++ model and it becomes the state sequence Ti,Tj, xi+1 appears after xi, which in other words is equivalent to a transfer from state Ti to Tj, and the transfer probability is(25)pij=the number ofTitoTjthe number ofTito other state.


## 3.2.3. Outliers Check

In this paper, the probability of a transition is compared with the probability of the state of the previous moment being transferred to the next most probable state as an evaluation criterion. Ifxi is to be inspected anomaly, the value of the previous moment is xi−1, and the state of Ti and Ti−1, Ti−1 is the most likely to be transferred to Tm; thus, we define an anomaly value probability:(26)pi=1−the possibility thatTi−1toTithe possibility thatTi−1toTm.From the above, it is known that the probability of stateTi to Tm is a constant, the smaller the probability of state Ti−1 transferring to Ti, the greater the pi, the greater the probability of xi being abnormal value.In the reality of hydrological data anomaly judgement standard, the data more than 2 cm are judged as outliers. Using the clustering algorithm may appear abnormal values and the moment before its values in the same condition and lead to false negatives. So after the judgement, outliers with a probability of 0 but a difference of more than 2 cm should also be judged as outliers.


## 3.3. The Anomaly Detection Method

The hydrological time series anomaly detection algorithm based on Flink is combined with two processes: the anomaly detection and result verification. First, the ARIMA model of time seriesx1,x2,…,xn is established by the idea of predictive detection by sliding window, and the predicted confidence interval is obtained, which is compared with the original data to identify the outliers. After the anomaly value is detected, the original data are clustered by K-Means++ algorithm, and the state transition matrix is computed after clustering. At last, the abnormal value is evaluated by the state transition and finally determined.The specific steps of the algorithm are provided in Algorithm1.Algorithm 1: Hydrological time series anomaly detection algorithm.
Input: hydrological time seriesX, reliability P, sliding window size L, historical batch data H.Output: the outliers in hydrological time series.Step 1: clean the sequenceH, including descending dimension, deleting duplicate value, sifting, and sortingStep 2: using the value ofL as the initial starting position of the sliding window of X, the value of the xL+1 is predicted, and as the window slides, the predicted value gradually forms a new time series xL+1,xL+2,…Step 3: the 95% confidence interval of the new time series is calculated and compared withX, the time point which is not in the confidence interval is obtained and get the exception point set e1,e2,…Step 4: taking historical dataH as input and training and establishing K-Means++ model, obtain the discrete state sequence T1,T2,…Step 5: compute the state transition matrix of the discrete state sequenceStep 6: theK-Mean++ model that is obtained in Step 4 is used for the exception point set of Step 3 and the value of its previous moment to obtain the state of anomaly and its previous momentStep 7: estimate the value of the exception and its previous moment in Step 6 by state transition data frame and then output the confidence scoreStep 8: repeat the above steps until no new data are enteredThe flow chart is shown below.Figure2 shows the overall flow of the algorithm. Initial abnormal detection is carried out on the input data stream through the sliding window in conjunction with ARIMA model.Figure 2
The whole process of anomaly detection algorithm.Figure3 shows the anomaly verification mechanism of specific process, and batch the historical data collected before. The first is to delete duplicates, dimension reduction, and sorting operation on batch data. After that, the data are clustered by the K-Means++ model. Then, original time series are converted into a Markov chain by using the clustering results. Finally, one step transfer matrix and the transition probability can be calculated. Through the principle of maximum and minimum, the hydrological sequence whose anomaly probability is higher than the threshold value is judged as the true outlier.Figure 3
Abnormal check mechanism flow chart.


## 4. Results and Discussion

### 4.1. Experimental Environment and Dataset

The runtime environment is deployed in a cluster using four PCs, and the hardware environment is as follows: CPU is Intel (R) Xeon (R) CPU E5645@2.40 GHz dual-core 24 CPU; memory is Kingston DDR3 1333 MHz 8G and 500 GB SSD flash memory. Operating system tools are Ubuntu 16.04 64-bit and Linux 3.11.0 kernel. The relevant software versions are as follows: Java 1.8 and Flink 1.5.2.The experimental data are from the monitoring data from 2015 to 2017 for more than 70 hydrological stations located on Chu river, with a data size of 18910864 rows.


### 4.2. Forecast Testing

Since it is difficult to describe the flow data directly in the form of a graph, Figure4 shows the change in water level in a whole certain period of time. Because of the time format of this experiment, the horizontal mark overlaps and the observation effect is poor. So, the number of days is used as the horizontal axis. Then, we introduce a time-continuous data stream for anomaly detection.Figure 4
The part of initial hydrologic time series.From Figure4, there are indeed data in the data that deviates significantly from its neighbors’ node, which is the anomaly we are trying to detect. The test results are given in Figure 5:Figure 5
The part of exception detection results.Figure5 shows the measured values, predictive values, and confidence intervals with a confidence level of 95% in the case of a sliding window with a length of 6 and an anomaly detected on a given dataset when the reliability is 95%. From Figure 5, most of the points are very close to the normal value, but there are some points outside the interval, so they are judged as the suspect anomaly point.


### 4.3. Batch Data Clean

Before the anomaly detection, the data obtained should be cleaned and the data before cleaning are as follows.We can see from Table1, there are many problems in the original data, such as duplication, sort confusion, date format not conforming to data mining requirement, and existence of unrelated series. To solve the above problems, we have cleaned the initial 18910864 rows hydrological data based on Flink and compared with the cleaning time in traditional single mode, and the results are as follows.Table 1
Data before cleaning.
Station IDTimeThe water levelSource1291052027/4/2016 19:45:0036.780Null1291052027/4/2016 19:45:0036.780Null1291056027/4/2016 19:45:0029.600Null1291058027/4/2016 19:40:0033.010Null6040320027/4/2016 19:35:005.940Nj34As shown in Figure6, when 15 stations are selected, the running speed of double node is slower than that of single node. However, with the increase in data size, the result calculated by double node has a small fluctuation, while the time of single node is significantly increased. Some Flink data after speed cleaning are shown in Table 2.Figure 6
Comparison of running time in Flink in different slaves.Table 2
Data after cleaning.
Station IDTimeThe water level129102802016-01-29 14:15:0022.3129102802016-01-29 14:20:0022.3129102802016-01-29 14:25:0022.3129102802016-01-29 14:30:0022.3129102802016-01-29 14:35:0022.3Table2 retains important data, eliminates invalid data and duplicate data, and unifies the data format. These data play a key role in the later abnormal check.


### 4.4. Anomaly Verification Model Runtime

This experiment aims at the problem of that if predicting time series based on sliding window, high delay, and long running calculation time will appear. Flink is adopted for calculation, and the time of using Flink to execute the algorithm under different computing resources is compared. The results are as follows.As can be seen from Figure7, with the data of 15 and 35 hydrological stations selected, the running speed of the double node is not ideal, but when the data size rises to the level of tens of millions, the advantages of the double node can be reflected with the advantages of time growth speed and calculation time. It can be seen that, under the larger dataset, the double-node operation speed is faster, with a decrease in 17.43% in the fastest case. The relevant experimental results for the delay will be combined with the algorithm effectiveness experiment in the last part of the experiment.Figure 7
Comparison of Flink running time in different slaves.


### 4.5. Outliers Assessment

After anomaly detection, it is necessary to evaluate the detected outliers and use the state transition matrix to calculate the probability of its true anomaly value. The display section in the form of a data frame is as follows.From Table3, t0 represents the state before the transfer and t1 represents the state after the transfer. By classifying the initial hydrological time series, selecting the state of the anomaly value and their previous moments and looking in Table 3 to get the probability of the real anomaly value, some of the results are as follows.Table 3
State transition data frame.
t0t1Probability110.9882835381100.0000000001110.0000000001120.001054482120.000000000Table4 shows that the value detected by the experiment in Forecast Testing section is really the probability of abnormal value, it can be seen that the probability of real abnormal value for some detected value is 0, so we will remove these values whose probability of real abnormal value is less than 50% from the detected values.Table 4
Anomaly evaluation.
Station IDTimeThe water levelPre_rzPro_ex1291052027/4/2016 19:45:0036.78022.3198.2701291052027/4/2016 19:45:0036.78022.320.00001291056027/4/2016 19:45:0029.60022.3299.3641291058027/4/2016 19:40:0033.01022.3299.3646040320027/4/2016 19:35:005.94022.3299.364


### 4.6. Effectiveness and Accuracy

After anomaly detection, it is necessary to evaluate the detected outliers and use the state transition matrix to calculate the probability of its true anomaly value. In order to verify the validity and accuracy of this mechanism, this paper divides the experimental results into 4 categories. The first category is TP (true positive), and the actual abnormal is judged to be abnormal; the second category is FN (false negative), and the actual abnormal is judged to be normal; the third category is FP (false positive), and the actual normal is judged to be abnormal; the last Category is TN (true negative), and the actual normal is judged to be normal. TP and TN are ideal situations where FN and FP are not desired. This paper defines sensitivity,sensitivity=TP/TP+FP, and specificity, specificity=TN/TN+FN. In this paper, the ARIMA model based on sliding window, MARS, and the algorithm in this paper are compared with the same hydrological time series, and the results are as follows (Table 5).Table 5
Comparison of the sensitivity and specificity on the ARIMA model based on sliding window, MARS, and the algorithm in this paper.
Window sizeNode numberProposed algorithmMARSARIMA based on slide windowSpecificity (%)Sensitivity (%)Specificity (%)Sensitivity (%)Specificity (%)Sensitivity (%)101 node82.1680.9887.4482.5581.5754.552 nodes82.6381.4387.6983.1282.6157.573 nodes82.4180.8487.5782.5181.4455.82301 node94.7490.8194.9987.8493.1267.282 nodes94.7790.6794.5687.7794.3468.163 nodes93.9890.9894.7387.9893.4568.45501 node99.4692.4097.6590.4199.1372.062 nodes99.3292.0298.0490.3399.2072.513 nodes99.6192.9898.1190.6399.2972.91Table5 shows the optimal results of the three models for comparison. The results show that the traditional window method is not significantly different from the method proposed in this paper in terms of specificity, maintaining high accuracy. In terms of sensitivity, as a nonlinear prediction model, MARS performs better in the case of a small sample, reaching a maximum of 83.12%. However, with the increase in the size of window, the increase is not as good as the algorithm proposed in this paper because the algorithm in this paper includes a historical verification mechanism, which uses historical data to verify real-time data. With the increase in window size, the algorithm proposed in this paper has a significant improvement compared with the traditional algorithm, with the sensitivity increasing from 72.91% to 92.98%. Comparing with the MARS model, the improvement is up to 2.75%.As can be seen from Figure8, with the increase in the number of windows, the specificity and sensitivity increased rapidly when the window size is about 30. During the window size from 10 to 30, the specificity increased by 29.47%, reaching 99.47%, the sensitivity increased by 42.08%, reaching 92.44%, and the delay increased by 93.32%, reaching 19.92 ms. However, when the window size is larger than 30, the specificity and sensitivity did not change much, but the delay rate increased to 57.21 ms and increased by 187.2%. Therefore, it is an ideal choice to set the window length to about 30, which can make the average delay less than 20 ms.Figure 8
Effects of window size on specificity, sensitivity, and latency.


## 4.1. Experimental Environment and Dataset

The runtime environment is deployed in a cluster using four PCs, and the hardware environment is as follows: CPU is Intel (R) Xeon (R) CPU E5645@2.40 GHz dual-core 24 CPU; memory is Kingston DDR3 1333 MHz 8G and 500 GB SSD flash memory. Operating system tools are Ubuntu 16.04 64-bit and Linux 3.11.0 kernel. The relevant software versions are as follows: Java 1.8 and Flink 1.5.2.The experimental data are from the monitoring data from 2015 to 2017 for more than 70 hydrological stations located on Chu river, with a data size of 18910864 rows.


## 4.2. Forecast Testing

Since it is difficult to describe the flow data directly in the form of a graph, Figure4 shows the change in water level in a whole certain period of time. Because of the time format of this experiment, the horizontal mark overlaps and the observation effect is poor. So, the number of days is used as the horizontal axis. Then, we introduce a time-continuous data stream for anomaly detection.Figure 4
The part of initial hydrologic time series.From Figure4, there are indeed data in the data that deviates significantly from its neighbors’ node, which is the anomaly we are trying to detect. The test results are given in Figure 5:Figure 5
The part of exception detection results.Figure5 shows the measured values, predictive values, and confidence intervals with a confidence level of 95% in the case of a sliding window with a length of 6 and an anomaly detected on a given dataset when the reliability is 95%. From Figure 5, most of the points are very close to the normal value, but there are some points outside the interval, so they are judged as the suspect anomaly point.


## 4.3. Batch Data Clean

Before the anomaly detection, the data obtained should be cleaned and the data before cleaning are as follows.We can see from Table1, there are many problems in the original data, such as duplication, sort confusion, date format not conforming to data mining requirement, and existence of unrelated series. To solve the above problems, we have cleaned the initial 18910864 rows hydrological data based on Flink and compared with the cleaning time in traditional single mode, and the results are as follows.Table 1
Data before cleaning.
Station IDTimeThe water levelSource1291052027/4/2016 19:45:0036.780Null1291052027/4/2016 19:45:0036.780Null1291056027/4/2016 19:45:0029.600Null1291058027/4/2016 19:40:0033.010Null6040320027/4/2016 19:35:005.940Nj34As shown in Figure6, when 15 stations are selected, the running speed of double node is slower than that of single node. However, with the increase in data size, the result calculated by double node has a small fluctuation, while the time of single node is significantly increased. Some Flink data after speed cleaning are shown in Table 2.Figure 6
Comparison of running time in Flink in different slaves.Table 2
Data after cleaning.
Station IDTimeThe water level129102802016-01-29 14:15:0022.3129102802016-01-29 14:20:0022.3129102802016-01-29 14:25:0022.3129102802016-01-29 14:30:0022.3129102802016-01-29 14:35:0022.3Table2 retains important data, eliminates invalid data and duplicate data, and unifies the data format. These data play a key role in the later abnormal check.


## 4.4. Anomaly Verification Model Runtime

This experiment aims at the problem of that if predicting time series based on sliding window, high delay, and long running calculation time will appear. Flink is adopted for calculation, and the time of using Flink to execute the algorithm under different computing resources is compared. The results are as follows.As can be seen from Figure7, with the data of 15 and 35 hydrological stations selected, the running speed of the double node is not ideal, but when the data size rises to the level of tens of millions, the advantages of the double node can be reflected with the advantages of time growth speed and calculation time. It can be seen that, under the larger dataset, the double-node operation speed is faster, with a decrease in 17.43% in the fastest case. The relevant experimental results for the delay will be combined with the algorithm effectiveness experiment in the last part of the experiment.Figure 7
Comparison of Flink running time in different slaves.


## 4.5. Outliers Assessment

After anomaly detection, it is necessary to evaluate the detected outliers and use the state transition matrix to calculate the probability of its true anomaly value. The display section in the form of a data frame is as follows.From Table3, t0 represents the state before the transfer and t1 represents the state after the transfer. By classifying the initial hydrological time series, selecting the state of the anomaly value and their previous moments and looking in Table 3 to get the probability of the real anomaly value, some of the results are as follows.Table 3
State transition data frame.
t0t1Probability110.9882835381100.0000000001110.0000000001120.001054482120.000000000Table4 shows that the value detected by the experiment in Forecast Testing section is really the probability of abnormal value, it can be seen that the probability of real abnormal value for some detected value is 0, so we will remove these values whose probability of real abnormal value is less than 50% from the detected values.Table 4
Anomaly evaluation.
Station IDTimeThe water levelPre_rzPro_ex1291052027/4/2016 19:45:0036.78022.3198.2701291052027/4/2016 19:45:0036.78022.320.00001291056027/4/2016 19:45:0029.60022.3299.3641291058027/4/2016 19:40:0033.01022.3299.3646040320027/4/2016 19:35:005.94022.3299.364


## 4.6. Effectiveness and Accuracy

After anomaly detection, it is necessary to evaluate the detected outliers and use the state transition matrix to calculate the probability of its true anomaly value. In order to verify the validity and accuracy of this mechanism, this paper divides the experimental results into 4 categories. The first category is TP (true positive), and the actual abnormal is judged to be abnormal; the second category is FN (false negative), and the actual abnormal is judged to be normal; the third category is FP (false positive), and the actual normal is judged to be abnormal; the last Category is TN (true negative), and the actual normal is judged to be normal. TP and TN are ideal situations where FN and FP are not desired. This paper defines sensitivity,sensitivity=TP/TP+FP, and specificity, specificity=TN/TN+FN. In this paper, the ARIMA model based on sliding window, MARS, and the algorithm in this paper are compared with the same hydrological time series, and the results are as follows (Table 5).Table 5
Comparison of the sensitivity and specificity on the ARIMA model based on sliding window, MARS, and the algorithm in this paper.
Window sizeNode numberProposed algorithmMARSARIMA based on slide windowSpecificity (%)Sensitivity (%)Specificity (%)Sensitivity (%)Specificity (%)Sensitivity (%)101 node82.1680.9887.4482.5581.5754.552 nodes82.6381.4387.6983.1282.6157.573 nodes82.4180.8487.5782.5181.4455.82301 node94.7490.8194.9987.8493.1267.282 nodes94.7790.6794.5687.7794.3468.163 nodes93.9890.9894.7387.9893.4568.45501 node99.4692.4097.6590.4199.1372.062 nodes99.3292.0298.0490.3399.2072.513 nodes99.6192.9898.1190.6399.2972.91Table5 shows the optimal results of the three models for comparison. The results show that the traditional window method is not significantly different from the method proposed in this paper in terms of specificity, maintaining high accuracy. In terms of sensitivity, as a nonlinear prediction model, MARS performs better in the case of a small sample, reaching a maximum of 83.12%. However, with the increase in the size of window, the increase is not as good as the algorithm proposed in this paper because the algorithm in this paper includes a historical verification mechanism, which uses historical data to verify real-time data. With the increase in window size, the algorithm proposed in this paper has a significant improvement compared with the traditional algorithm, with the sensitivity increasing from 72.91% to 92.98%. Comparing with the MARS model, the improvement is up to 2.75%.As can be seen from Figure8, with the increase in the number of windows, the specificity and sensitivity increased rapidly when the window size is about 30. During the window size from 10 to 30, the specificity increased by 29.47%, reaching 99.47%, the sensitivity increased by 42.08%, reaching 92.44%, and the delay increased by 93.32%, reaching 19.92 ms. However, when the window size is larger than 30, the specificity and sensitivity did not change much, but the delay rate increased to 57.21 ms and increased by 187.2%. Therefore, it is an ideal choice to set the window length to about 30, which can make the average delay less than 20 ms.Figure 8
Effects of window size on specificity, sensitivity, and latency.


## 5. Conclusions

In the era of big data, traditional detection algorithms cannot meet the current needs. Based on the characteristics of the sliding window and the defects of traditional sliding window inspection, such as high time complexity and high error detection rate disadvantages, this paper puts forward a kind of hydrological time series anomaly detection method based on Flink. By using Flink calculation, this method reduces the computing time and combines two processes. After cleaning the data, the sliding window and the ARIMA model are used to forecast on the Flink platform. Then, the confidence interval is calculated for the predicted result and evaluated it as an anomaly value outside the interval range. Based on the detection result, theK-means algorithm is used to cluster the original data and the state transition probability is calculated.Taking the data of hydrologic sensor obtained from the Chu River as an example of experimental data, experiments on the detection time and validity of outliers are carried out, respectively. The result shows that the million data using 2 slaves cost more time than 1 slave in the calculation time, but when the tens data are calculated, 2 slaves are better than 1 slave, and the maximum is reduced by 17.43%. The sensitivity of the evaluation is increased from 72.91% to 92.98%. In terms of delay, the average delay of different slaves is roughly the same, which is maintained within 20 ms. It shows that Flink can effectively improve the calculation efficiency by adding nodes when the method is used to detect tens of millions of hydrologic data. At the same time, the sensitivity of the method is significantly improved compared with the traditional method.However, there is still room for further improvement in the detection accuracy of the algorithm proposed in this paper. The follow-up work will focus on more accurate identification of which are outliers and which are normal fluctuations caused by natural factors.



---

*Source: 3187697-2020-05-28.xml*