基于时间序列的异常检测

简介

搜罗了网上几乎所有的基于时间序列的异常检测方法，没有包括文献，整理记录一下。

综合引用以下文章：

• 数据挖掘导论
• 时间序列异常检测机制的研究
• KPI异常检测竞赛笔记
• 异常检测之时间序列的异常检测
• Skyline timeseries异常判定算法
• 腾讯蓝鲸数据平台之告警系统
• 异常点检测算法综述
• 时间序列简介（一）

时间序列预测

• 周期性：AE/VAE
• 稳定性：ARIMA
• 不稳定性：频谱分析/小波分析

RNN貌似对以上都适用。

异常检测

异常检测的方法有很多，根据《数据挖掘导论》，有以下几个分类：

基于模型的方法

许多异常检测技术首先建立一个数据模型。异常是那些同模型不能完美拟合的对象。例如，数据分布模型可以通过估计概率分布的参数来创建。如果一个对象不能很好地同该模型拟合，即如果它很可能不服从该分布，则它是一个异常。如果模型是簇的集合，则异常是不显著属于任何簇的对象。在使用回归模型时，异常是相对远离预测值的对象。

由于异常和正常对象可以看作两个不同的类，因此可以使用分类方法来建立这两个类的模型。当然，仅当某些对象存在类标号，使得我们可以构造训练数据集时才可以使用分类方法。此外，异常相对稀少，在选择分类方法和评估度量是需要考虑这一因素。

基于模型的方法又称为统计方法，主要通过拟合单个模型或多个模型来判断该点的概率。

• 单维情况

  • 3$\sigma$-法则

    $(μ-3σ,μ+3σ)$区间内的概率为99.74。所以可以认为，当数据分布区间超过这个区间时，即可认为是异常数据。

  • 箱型图

    IQR是第三四分位数减去第一四分位数，大于Q3+1.5*IQR之外的数和小于Q1-1.5*IQR的值被认为是异常值。

  • Grubbs测试

    Grubbs测试是一种从样本中找出outlier的方法，所谓outlier，是指样本中偏离平均值过远的数据，他们有可能是极端情况下的正常数据，也有可能是测量过程中的错误数据。使用Grubbs测试需要总体是正态分布的。

    算法流程：

    1. 样本从小到大排序
    2. 求样本的mean和std.dev
    3. 计算min/max与mean的差距，更大的那个为可疑值
    4. 求可疑值的z-score (standard score)，如果大于Grubbs临界值，那么就是outlier；

    Grubbs临界值可以查表得到，它由两个值决定：检出水平$α$（越严格越小），样本数量$n$。排除outlier，对剩余序列循环做 1-4 步骤。由于这里需要的是异常判定，只需要判断tail_avg是否outlier即可。

• 高维情况

  • Mahalanobis距离
  • $\mathcal{X}^2$统计
  • 单个/混合高斯分布

在某些情况下，很难建立模型，例如，因为数据的统计分布未知或没有训练数据可用。在这些情况下，可以使用如下所述的不需要模型的技术。

基于距离的方法/基于邻近度的方法

通常可以在对象之间定义邻近性度量，并且许多异常检测方法都基于邻近度。异常对象是那些远离大部分其他对象的对象。这一领域的许多方法都基于距离，称作基于距离的离群点检测方法。当数据能够以二维或三维散布图显示时，通过寻找与大部分其他点分离的点，可以从视觉上检测出基于距离的离群点。

但是该方法计算复杂度过高，且分布不均匀的点，容易出错。

方法：

• 基于角度的离群点检测

  角度越小，说明距离越远。

• k-最近邻

  到k-最近邻的距离。评分为数据对象与最近的k个点的距离之和。很明显，与k个最近点的距离之和越小，异常分越低；与k个最近点的距离之和越大，异常分越大。设定一个距离的阈值，异常分高于这个阈值，对应的数据对象就是异常点。

• Local Outlier Factor（LOF）

  LOF得分为数据对象的k个最近邻的平均局部密度与数据对象本身的局部密度之比。

• Connectivity Outlier Factor（COF）

  如果点p的平均连接距离大于它的k最近邻的平均连接距离，则点p是异常点。COF将异常值识别为其邻域比其近邻的邻域更稀疏的点。

• stochastic outlier selection algorithm（无监督）

  基于方差，近邻点越多，方差越小；近邻点越大，方差越大。

基于密度的方法

对象的密度估计可以相对直接地计算，特别是当对象之间存在邻近性度量时，低密度区域中的对象相对远离近邻，可能被看作异常。一种更复杂的方法考虑到数据集可能有不同密度区域这一事实，仅当一个点的局部密度显著地低于它的大部分近邻时才将其分类为离群点。

定义1：基于密度的异常

异常就是那些在低密度区域的数据对象，一个数据对象的异常分就是该对象所在区域的密度的倒数，下面是基于密度的异常分的计算公式：

$$\text{densitiy}(x,k)=(\frac{\sum_{y\in N(x,k)}\text{distance}(x,y)}{|N(x,k)|})^{-1}$$

其中$N(x,k)$指的是$x$的$k$个最近的邻居的集合，$|N(x,k)|$表示该集合的大小，$y$是$x$最近的邻居。

定义2：给定半径的邻域内的数据对象数

一个数据对象的密度等于半径为d的邻域内的数据对象数。

d的选择很重要，若d太小，则会有很多正常的数据对象被认为是异常点；若d太大，则很多异常数据对象会被误判为正常点。事实上，当密度分布不均匀的时候，上述方法得到的异常点会不正确。为了克服密度不均匀的情况，我们使用下面的平均相对密度来作为异常分。

定义3：平均相对密度

$$\text{average relative densitiy}(x,k)=(\frac{\text{density}(x,k)}{\Sigma _{y\in N(x,k)}\text{density}(y,k)/|N(x,k)| })$$

基于相对密度的异常检测算法

1. {$k$是最近邻的个数}（类似KNN中的$k$）
2. for 所有的数据对象$x$ do
   1. 计算$N(x,k)$
   2. 计算$\text{density}(x,k)$
3. end for
4. for 所有的数据对象x do
   1. $\text{outlier score}(x) = \text{average relative density}(x,k)$
5. end for
6. 取$\text{outlier score}$最大的$N$个数据对象作为异常

基于聚类的方法

一种利用聚类检测离群点的方法是丢弃远离其他簇的小簇。这种方法可以与任何聚类方法一起使用，但是需要最小簇大小和小簇与其他簇之间距离的國值。通常，该过程可以简化为丢弃小于某个最小尺寸的所有簇。这种方案对簇个数的选择高度敏感。此外，使用这一方案，很难将离群点得分附加在对象上。注意，把一组对像看作离群点，将离群点的概念从个体对象扩展到对象组，但是本质上没有任何改变。

一种更系统的方法是，首先聚类所有对象，然后评估对象属于簇的程度。对于基于原型的聚类，可以用对象到它的簇中心的距离来度量对象属于簇的程度。更一般地，对于基于目标函数的聚类方法，可以使用该目标函數来评估对象属于任意簇的程度。特殊情况下，如果删除一个对象导致该目标的显著改进，则我们可以将对象分类为离群点。如，对于K均值，脷除远离其相关簇中心的对象能够显著地改进簇的误差的平方和（SSE）。总而言之，聚类创建数据的模型，而异常扭曲模型。

基于划分的方法

孤立森林

基于划分的思想，划分成树，深度越低，说明越容易被划分，即为离群点。

• https://www.jianshu.com/p/4d817084a69a
• https://www.jianshu.com/p/5af3c66e0410
• http://www.cnblogs.com/fengfenggirl/p/iForest.html

基于线性的方法

• Principal Component Analysis（PCA）

  PCA在做特征值分解之后得到的特征向量反应了原始数据方差变化程度的不同方向，特征值为数据在对应方向上的方差大小。所以，最大特征值对应的特征向量为数据方差最大的方向，最小特征值对应的特征向量为数据方差最小的方向。原始数据在不同方向上的方差变化反应了其内在特点。如果单个数据样本跟整体数据样本表现出的特点不太一致，比如在某些方向上跟其它数据样本偏离较大，可能就表示该数据样本是一个异常点。

基于非线性的方法

• Replicator Neural Networks（RNNs）
• AutoEncoder（AE）：Anomaly Detection异常检测的几种方法
• Variational AutoEncoder（VAE）：AIOps探索：基于VAE模型的周期性KPI异常检测方法

针对非数值型的方法

• Attribute Value Frequency

  针对非数值型的数据，即类别离散数据的方法。

基于深度学习的方法

• Deep Learning for Anomaly Detection: A Survey

相关综述可以阅读上面的文献。另外，随着贝叶斯与深度网络的结合，能够衡量模型的不确定性的贝叶斯深度网络在各个领域中应用广泛，也诞生了很多专门针对于神经网络计算不确定性的方法，具体文献可以参考Uncertainty in Deep Learning。另外一篇博客介绍的Uber用于预测和检测的一个框架，就是基于不确定性的思路做的。

总结

其实，可以看到，以上方法的分类之间也是相近的，说到底，概率、密度和聚类也是由空间中的数据点之间的距离所决定的。而本文主要是想记录基于时间序列的异常检测方法，当然，以上的方法，在时间序列抽取后的特征空间中也能够使用。

基于时间序列的异常检测方法

同比，是指在相邻时段中的某一相同时间点进行比较。

环比，是指在同一时段中的相邻时间点进行比较。

短期环比（SS）

对于时间序列（是指将同一统计指标的数值按其发生的时间先后顺序排列而成的数列）来说，$T$时刻的数值对于$T-1$时刻有很强的依赖性。比如流量在8:00很多，在8:01时刻的概率是很大的，但是如果07:01时刻对于8:01时刻影响不是很大。

首先，我们可以使用最近时间窗口（$T$）内的数据遵循某种趋势的现象来做文章。比如我们将$T$设置为7，则我们取检测值（$\text{now value}$）和过去7个（记为$i$）点进行比较，如果大于阈值我们将$\text{count}$加1，如果$\text{count}$超过我们设置的$\text{count num}$，则认为该点是异常点。

$$\text{count}(\sum_{i=0}^T(\text{now value}-i)>\text{threshold})>\text{count num}$$

上面的公式涉及到$\text{threshold}$和$\text{count num}$两个参数，$\text{threshold}$如何获取我们将在下节进行介绍，而$\text{count num}$可以根据的需求进行设置，比如对异常敏感，可以设置$\text{count num}$小一些，而如果对异常不敏感，可以将$\text{count num}$设置的大一些。

动态阈值

业界关于动态阈值设置的方法有很多，这里介绍一种针对360的lvs流量异常检测的阈值设置方法。通常阈值设置方法会参考过去一段时间内的均值、最大值以及最小值，我们也同样应用此方法。取过去一段时间（比如$T$窗口）的平均值、最大值以及最小值，然后取$\text{max-avg}$和$\text{avg-min}$的最小值。之所以取最小值的原因是让筛选条件设置的宽松一些，让更多的值通过此条件，减少一些漏报的事件。

$$\text{threshold} = \text{min(max-avg, avg-min)}$$

窗口特征

统计特征

3-sigma

一个很直接的异常判定思路是，拿最新3个datapoint的平均值（tail_avg方法）和整个序列比较，看是否偏离总体平均水平太多。怎样算“太多”呢，因为standard deviation表示集合中元素到mean的平均偏移距离，因此最简单就是和它进行比较。在normal distribution（正态分布）中，99.73%的数据都在偏离mean 3个σ (standard deviation 标准差) 的范围内。如果某些datapoint到mean的距离超过这个范围，则认为是异常的。

z score

标准分，一个个体到集合mean的偏离，以标准差为单位，表达个体距mean相对“平均偏离水平（std dev表达）”的偏离程度，常用来比对来自不同集合的数据。

在模型中，z_score用来衡量窗口数据中，中间值的偏离程度。

算法流程：

1. 排除最后一个值
2. 求剩余序列的平均值
3. 全序列减去上面这个平均值
4. 求剩余序列的标准差
5. （ 中间三个数的平均值-全序列均值）/ 全序列标准差

Grubbs格拉斯测试

Grubbs测试是一种从样本中找出outlier的方法，所谓outlier，是指样本中偏离平均值过远的数据，他们有可能是极端情况下的正常数据，也有可能是测量过程中的错误数据。使用Grubbs测试需要总体是正态分布的。

算法流程：

1. 样本从小到大排序
2. 求样本的mean和std.dev
3. 计算min/max与mean的差距，更大的那个为可疑值
4. 求可疑值的z-score (standard score)，如果大于Grubbs临界值，那么就是outlier；

Grubbs临界值可以查表得到，它由两个值决定：检出水平$α$（越严格越小），样本数量$n$。排除outlier，对剩余序列循环做 1-4 步骤。由于这里需要的是异常判定，只需要判断tail_avg是否outlier即可。

moving average（移动平均）

给定一个时间序列和窗口长度$N$，moving average等于当前data point之前$N$个点（包括当前点）的平均值。不停地移动这个窗口，就得到移动平均曲线。

cumulative moving average（累加移动平均）

设$\{x_i:i \ge 1\}$是观察到的数据序列。 累积移动平均线是所有数据的未加权平均值。 如果若干天的值是$x_1,…,x_i$，那么

$$CMA_i=\frac{x_1+...+x_i}{i}$$

当有新的值$x_{i+1}$，那么累积移动平均为

$$\begin{align*} CMA_{i+1} &= \frac{x_1+...+x_i+x_{i+1}}{i+1} \\ & = \frac{x_{i+1}+n·CMA_i}{i+1}\\ & = CMA_i + \frac{x_{i+1}-CMA_i}{i+1} \end{align*}$$

weighted moving average（加权移动平均）

加权移动平均值是先前$w$个数据的加权平均值。 假设$\sum_{j=0}^{w-1}\text{weight}_j=1$，其中$\text{weight}_j>0$，则加权移动平均值为，

$$WMA_i=\sum_{j=0}^{w-1}\text{weight}_j · x_{i-j}$$

一般，

$$\text{weight}_j = \frac{w-j}{w+(w-1)+...+1}, for\ 0\le j \le w-1$$

所以，

$$WMA_i=\frac{wx_i+(w-1)x_{i-1}+...+2x_{i-w+2}+x_{i-w+1}}{w+(w-1)+...+1}$$

exponential weighted moving average（指数加权移动平均）

指数移动与移动平均有些不同：

1. 并没有时间窗口，用的是从时间序列第一个data point到当前data point之间的所有点。
2. 每个data point的权重不同，离当前时间点越近的点的权重越大，历史时间点的权重随着离当前时间点的距离呈指数衰减，从当前data point往前的data point，权重依次为$α$, $α(1-α)$, $α(1-α)^2$,…, $α(1-α)^n$。

该算法可以检测一个异常较短时间后发生另外一个异常的情况，异常持续一段时间后可能被判定为正常。

PS：https://www.jianshu.com/p/a2dbd47b3f1a

double exponential smoothing（双指数平滑）

假设$\{Y_t:t\ge1\}$是观测数据序列，有两个与双指数平滑相关的方程：

$$S_t = \alpha Y_t+(1-\alpha)(S_{t-1}+b_{t-1}),\\ b_t = \beta(S_t-S_{t-1})+(1-\beta)b_{t-1}$$

其中$\alpha \in [0,1]$是数据平滑因子，$\beta \in [0,1]$是趋势平滑因子。

这里，初始值为$S_1=Y_1$，$b_1$有三种可能：

$$b_1=Y_2-Y_1,\\ b_1 = \frac{(Y_2-Y_1)+(Y_3-Y_2)+(Y_4-Y_3)}{3}=\frac{Y_4-Y_1}{3},\\ b_1 = \frac{Y_n-Y_1}{n-1}$$

预测：

• 单个时刻：$F_{t+1}=S_t+b_t$
• $m$个时刻：$F_{t+m}=S_t+mb_t$

triple exponential smoothing（三指数平滑）（Holt-Winters）

• multiplicative seasonality

  设$\{Y_t:t \ge 1\}$是观察到的数据序列，则三指数平滑为

  $$S_t = \alpha \frac{Y_t}{c_{t-L}}+(1-\alpha)(S_{t-1}+b_{t-1}), \text{overall smoothing} \\ b_t = \beta(S_t-S_{t-1})+(1-\beta)b_{t-1}, \text{trend smoothing}\\ c_t = \gamma\frac{Y_t}{S_t}+(1-\gamma)c_{t-L}, \text{seasonal smoothing}$$

  其中$\alpha \in [0,1]$是数据平滑因子，$\beta \in [0,1]$是趋势平滑因子， $\gamma \in [0,1]$是季节变化平滑因子。

  预测$m$个时刻：$F_{t+m}=(S_t+mb_t)c_{(t-L+m) mod L}$

  初始值设定：

  $$S_1 = Y_1,\\ b_0 = \frac{(Y_{L+1}-Y_1)+(Y_{L+2}-Y_2)+...+(Y_{L+L}-Y_L)}{L},\\ c_i = \frac{1}{N}\sum_{j=1}^N\frac{Y_{L(j-1)+i}}{A_j}, \forall i \in \{1,...,L\}\\ A_j=\frac{\sum_{i=1}^LY_{L(j-1)+i}}{L}, \forall \in \{1,...,N\}$$

• additive seasonality

  设$\{Y_t:t \ge 1\}$是观察到的数据序列，则三指数平滑为

  $$S_t = \alpha(Y_t-c_{t-L})+(1-\alpha)(S_{t-1}+b_{t-1}), \text{overall smoothing} \\ b_t = \beta(S_t- S_{t-1})+(1-\beta)b_{t-1}, \text{trend smoothing} \\ c_t = \gamma(Y_t-S_{t-1}-b_{t-1})+(1-\gamma)c_{t-L}, \text{seasonal smoothing}$$

  其中$\alpha \in [0,1]$是数据平滑因子，$\beta \in [0,1]$是趋势平滑因子， $\gamma \in [0,1]$是季节变化平滑因子。

  预测$m$个时刻：$F_{t+m}=S_t+mb_t+c_{(t-L+m) mod L}$

stddev from average（平均值-标准差）

该算法类似于3sigma准则。当数据服从高斯分布时，数值分布在$(μ-3σ,μ+3σ)$区间内的概率为99.74。所以，可以这么认为，当数据分布区间超过这个区间时，即可认为是异常数据。该算法使用$(t-\text{mean})/\text{std}$ 作为特征，用于衡量中间三个值的平均值相对于$3σ$的距离。

该算法的特点是可以有效屏蔽“在一个点上突变到很大的异常值但在下一个点回落到正常水平”的情况，即屏蔽单点毛刺：因为它使用的是3个点的均值（有效缓和突变），和整个序列比较（均值可能被异常值拉大），导致判断正常。

算法流程：

1. 求窗口数据的平均值。 (mean)
2. 求窗口数据的标准差。 (std)
3. 求窗口数据中间3个值的平均值。（t）
4. 使用$(t-mean)/std$作为特征。

stddev from moving average（移动平均-标准差）

先求出最后一个点处的指数加权移动平均值，然后再用最新的点和三倍方差方法求异常。

stddev from ewma（指数加权移动平均-标准差）

类似于特征3，不过在计算$(t-\text{mean})/\text{std}$时，使用的mean，std分别为对窗口数据进行移动加权平均后的平均值以及方差。

histogram bins

该算法和以上都不同，它首先将timeseries划分成15个宽度相等的直方，然后判断tail_avg所在直方内的元素是否<=20，如果是，则异常。

直方的个数和元素个数判定需要根据自己的metrics调整，不然在数据量小的时候很容易就异常了。

median absolute deviation（中位数绝对偏差）

median：大部分情况下我们用mean来表达一个集合的平均水平（average），但是在某些情况下存在少数极大或极小的outlier，拉高或拉低了（skew）整体的mean，造成估计的不准确。此时可以用median（中位数）代替mean描述平均水平。Median的求法很简单，集合排序中间位置即是，如果集合总数为偶数，则取中间二者的平均值。

median of deviation（MAD）：同mean一样，对于median我们也需要类似standard deviation这样的指标来表达数据的紧凑/分散程度，即偏离average的平均距离，这就是MAD。MAD顾名思义，是deviation的median，而此时的deviation = abs( 个体 – median )，避免了少量outlier对结果的影响，更robust。

绝对中位差实际求法是用原数据减去中位数后得到的新数据的绝对值的中位数。

1. 原数据-中位值=新数据
2. 新数据的绝对值的中位数作为特征

mean subtraction cumulation（平均值减法累积？）

该特征类似于3-sigma准则。

算法流程

1. 排除全序列（暂称为all）最后一个值（last datapoint），求剩余序列（暂称为rest，0..length-2）的mean；
2. rest序列中每个元素减去rest的mean，再求标准差；
3. 求窗口数据中间点到rest mean的距离，即 abs(last datapoint – rest mean)；

first hour average（前若干-平均值）

和上述算法基本一致，只是比较对象不是整个序列，而是开始一个小时（其实这种这种思想可以推广，只要是时间序列刚开始的一段时间即可）的以内的数据，求出这段时间的均值和标准差和尾部数据（新产生的数据）用三倍方差的方法比较即可。

熵特征

为什么要研究时间序列的熵呢？请看下面两个时间序列：

时间序列（1）：(1,2,1,2,1,2,1,2,1,2,…)

时间序列（2）：(1,1,2,1,2,2,2,2,1,1,…)

在时间序列（1）中，1 和 2 是交替出现的，而在时间序列（2）中，1 和 2 是随机出现的。在这种情况下，时间序列（1）则更加确定，时间序列（2）则更加随机。并且在这种情况下，两个时间序列的统计特征，例如均值，方差，中位数等等则是几乎一致的，说明用之前的统计特征并不足以精准的区分这两种时间序列。

通常来说，要想描述一种确定性与不确定性，熵（entropy）是一种不错的指标。对于离散空间而言，一个系统的熵（entropy）可以这样来表示：

$$\text{entropy}(X) = -\sum_{i=1}^{\infty}P\{x=x_{i}\}\ln(P\{x=x_{i}\}).$$

如果一个系统的熵（entropy）越大，说明这个系统就越混乱；如果一个系统的熵越小，那么说明这个系统就更加确定。

提到时间序列的熵特征，一般来说有几个经典的例子，那就是 binned entropy，approximate entropy，sample entropy。下面来一一介绍时间序列中这几个经典的熵。

Binned Entropy

从熵的定义出发，可以考虑把时间序列 $X_{T}$的取值进行分桶的操作。例如，可以把$[\min(X_{T}), \max(X_{T})]$这个区间等分为十个小区间，那么时间序列的取值就会分散在这十个桶中。根据这个等距分桶的情况，就可以计算出这个概率分布的熵（entropy）。i.e. Binned Entropy 就可以定义为：

$$\text{binned entropy}(X) = -\sum_{k=0}^{\min(\text{maxbin}, \text{len}(X))} p_{k}\ln(p_{k})\cdot 1_{(p_{k}>0)},$$

其中 $p_{k}$ 表示时间序列 $X_{T}$ 的取值落在第 $k$ 个桶的比例（概率），$\text{maxbin}$ 表示桶的个数， $\text{len}(X_{T}) = T$ 表示时间序列 $X_{T}$ 的长度。

如果一个时间序列的 Binned Entropy 较大，说明这一段时间序列的取值是较为均匀的分布在 $[\min(X_{T}), \max(X_{T})]$之间的；如果一个时间序列的 Binned Entropy 较小，说明这一段时间序列的取值是集中在某一段上的。

Approximate Entropy

回到本节的问题，如何判断一个时间序列是否具备某种趋势还是随机出现呢？这就需要介绍 Approximate Entropy 的概念了，Approximate Entropy 的思想就是把一维空间的时间序列提升到高维空间中，通过高维空间的向量之间的距离或者相似度的判断，来推导出一维空间的时间序列是否存在某种趋势或者确定性。那么，我们现在可以假设时间序列 $X_{N}: \{x_{1},\cdots, x_{N}\}$ 的长度是 $N$ ，同时 Approximate Entropy 函数拥有两个参数， $m$ 与 $r$ ，下面来详细介绍 Approximate Entropy 的算法细节。

1. 固定两个参数，正整数 $m$ 和正数 $r$，正整数 $m$ 是为了把时间序列进行一个片段的提取，正数 $r$ 是表示时间序列距离的某个参数。i.e. 需要构造新的 $m$ 维向量如下：

$$X_{1}(m) = (x_{1},\cdots, x_{m})\in\mathbb{R}^{m},\\ X_{i}(m) = (x_{i},\cdots, x_{m+i-1})\in\mathbb{R}^{m},\\ X_{N-m+1}(m) = (x_{N-m+1},\cdots, x_{N})\in\mathbb{R}^{m}.$$

1. 通过新的向量$X_{1}(m),\cdots, X_{N-m+1}(m)$，可以计算出哪些向量与 $X_i$ 较为相似。即，

$$C_{i}^{m}(r) = (\text{number of }X_{j}(m)\text{ such that } d(X_{i}(m), X_{j}(m))\leq r)/(N-m+1),$$

​ 在这里，距离$d$选择$L^1$，$L^2$，$L^p$，$L^\infty$范数。在这个场景下，距离$d$通常选择为$L^\infty$范数。

1. 考虑函数

$$\Phi^{m}(r) = (N-m+1)^{-1}\cdot \sum_{i=1}^{N-m+1} \ln(C_{i}^{m}(r)),$$

1. Approximate Entropy 可以定义为：

$$\text{ApEn}(m,r) = \Phi^{m}(r)-\Phi^{m+1}(r).$$

备注：

1. 正整数 $m$ 一般可以取值为 $2$ 或者 $3$， $r>0$ 会基于具体的时间序列具体调整；
2. 如果某条时间序列具有很多重复的片段（repetitive pattern）或者自相似性（self-similarity pattern），那么它的 Approximate Entropy 就会相对小；反之，如果某条时间序列几乎是随机出现的，那么它的 Approximate Entropy 就会相对较大。

Sample Entropy

除了 Approximate Entropy，还有另外一个熵的指标可以衡量时间序列，那就是 Sample Entropy，通过自然对数的计算来表示时间序列是否具备某种自相似性。

按照以上 Approximate Entropy 的定义，可以基于$m$与$r$定义两个指标$A$与$B$，分别是

$$A = \#\{\text{vector pairs having } d(X_{i}(m+1),X_{j}(m+1))<r \text{ of length } m+1 \},\\ B = \#\{ \text{vector pairs having } d(X_{i}(m), X_{j}(m))<r \text{ of length } m\}.$$

其中，$#$ 表示集合的元素个数。根据度量 $d$ （无论是$L^1$，$L^2$，$L^\infty$ ）的定义可以知道 $A\le B$，因此 Sample Entropy 总是非负数，即

$$\text{SampEn} = -\ln(A/B) \geq 0.$$

备注：

1. Sample Entropy 总是非负数；
2. Sample Entropy 越小表示该时间序列具有越强的自相似性（self similarity）。
3. 通常来说，在 Sample Entropy 的参数选择中，可以选择$m=2,r=0.2 \times std$。

分段特征

即使时间序列有一定的自相似性（self-similarity），能否说明这两条时间序列就完全相似呢？其实答案是否定的，例如：两个长度都是 1000 的时间序列，

时间序列（1）： [1,2] * 500

时间序列（2）： [1,2,3,4,5,6,7,8,9,10] * 100

其中，时间序列（1）是 1 和 2 循环的，时间序列（2）是 1~10 这样循环的，它们从图像上看完全是不一样的曲线，并且它们的 Approximate Entropy 和 Sample Entropy 都是非常小的。那么问题来了，有没有办法提炼出信息，从而表示它们的不同点呢？答案是肯定的。

首先，我们可以回顾一下 Riemann 积分和 Lebesgue 积分的定义和不同之处。按照下面两幅图所示，Riemann 积分是为了算曲线下面所围成的面积，因此把横轴划分成一个又一个的小区间，按照长方形累加的算法来计算面积。而 Lebesgue 积分的算法恰好相反，它是把纵轴切分成一个又一个的小区间，然后也是按照长方形累加的算法来计算面积。

之前的 Binned Entropy 方案是根据值域来进行切分的，好比 Lebesgue 积分的计算方法。现在我们可以按照 Riemann 积分的计算方法来表示一个时间序列的特征，于是就有学者把时间序列按照横轴切分成很多段，每一段使用某个简单函数（线性函数等）来表示，于是就有了以下的方法：

1. 分段线性逼近（Piecewise Linear Approximation）
2. 分段聚合逼近（Piecewise Aggregate Approximation）
3. 分段常数逼近（Piecewise Constant Approximation）

说到这几种算法，其实最本质的思想就是进行数据降维的工作，用少数的数据来进行原始时间序列的表示（Representation）。用数学化的语言来描述时间序列的数据降维（Data Reduction）就是：把原始的时间序列 $\{x_1,…,x_N\}$ 用 $\{x_1’,…,x_D’\}$来表示，其中$D<N$。那么后者就是原始序列的一种表示（representation）。

分段聚合逼近（Piecewise Aggregate Approximation）—- 类似 Riemann 积分

在这种算法中，分段聚合逼近（Piecewise Aggregate Approximation）是一种非常经典的算法。假设原始的时间序列是 $C=\{x_1,…,x_N\}$，定义PAA的序列是：

$$\overline{C} = \{\overline{x}_{1},\cdots,\overline{x}_{w}\},$$

其中，

$$\overline{x}_{i} = \frac{w}{N} \cdot \sum_{j=\frac{N}{w}(i-1)+1}^{\frac{N}{w}i} x_{j}.$$

在这里$1\le i \le w$。用图像来表示那就是

至于分段线性逼近（Piecewise Linear Approximation）和分段常数逼近（Piecewise Constant Approximation），只需要在 $\overline{x}_{i}$ 的定义上稍作修改即可。

符号逼近（Symbolic Approximation）—- 类似 Riemann 积分

在推荐系统的特征工程里面，特征通常来说可以做归一化，二值化，离散化等操作。例如，用户的年龄特征，一般不会直接使用具体的年月日，而是划分为某个区间段，例如 0~6（婴幼儿时期），7~12（小学），13~17（中学），18~22（大学）等阶段。

其实在得到分段特征之后，分段特征在某种程度上来说依旧是某些连续值，能否把连续值划分为一些离散的值呢？于是就有学者使用一些符号来表示时间序列的关键特征，也就是所谓的符号表示法（Symbolic Representation）。下面来介绍经典的 SAX Representation。

如果我们希望使用 $\alpha$ 个符号，例如 $\{l_1,…,l_\alpha\}$来表示时间序列。同时考虑分布$N(0,1)$，用

$\{z_{1/\alpha},\cdots,z_{(\alpha-1)/\alpha}\}$来表示 Gauss 曲线下方的一些点，而这些点把 Gauss 曲线下方的面积等分成了 $\alpha$段。

算法流程：

1. 正规化（normalization）：也就是该时间序列被映射到均值为零，方差为一的区间内。
2. 分段表示（PAA）：$\{x_{1},\cdots, x_{N}\} \Rightarrow \{\overline{x}_{1},\cdots,\overline{x}_{w}\}.$
3. 符号表示（SAX）：如果$\overline{x}_{i}<z_{1/\alpha}$，那么$\hat{X}_{i}=l_{1}$；如果$z_{(j-1)/\alpha}\leq \overline{x}_{i}<z_{j/\alpha}$，那么$\hat{X}_{i} = l_{j}$；如果$\overline{x}_{i}\geq z_{(\alpha-1)/\alpha}$，那么$\hat{X}_{i} = l_{\alpha}$。

于是，我们就可以用$\{l_{1},\cdots,l_{\alpha}\}$这$\alpha$个字母来表示原始的时间序列了。

总特征

总的训练特征为：时间窗口特征 + OneHot特征 + Timestamp所属的星期作为特征，label选取时间窗口中间点的标签。

因为选择了窗口数据的中间位置的点作为label值，所以在窗口移动过程中，在数据起始点和终结点会丢失 datasize-(windowsize/2) 个数据。（datasize为数据的总量，windowsize为时间窗口的大小)

在最终提交的结果中，缺失的数据点的预测结果，用0补充。

分类方法

• 随机森林
• XGBoost
• DNN

长期环比（LS）

上面短期环比参考的是短期内的数据，而仅仅有短期内的数据是不够的，我们还需要参考更长时间内数据的总体走势。

通常使用一条曲线对该趋势进行拟合来反应曲线的走势，如果新的数据打破了这种趋势，使曲线变得不平滑，则该点就出现了异常。曲线拟合的方法有很多，比如回归、moving average……。在这里，我们使用EWMA，即指数权重移动平均方法来拟合曲线。在EWMA中，下一点的平均值是由上一点的平均值，加上当前点的实际值修正而来。对于每一个EWMA值，每个数据的权重是不一样的，最近的数据将拥有越高的权重。

有了平均值之后，我们就可以使用3-sigma理论来判断新的input是否超过了容忍范围。比较实际的值是否超出了这个范围就可以知道是否可以告警了。

自己理解：长期环比除了通过拟合、EWMA等方法，ARIMA和RNN也应该能适用

同比(chain)

很多监控项都具有一定的周期性，其中以一天为周期的情况比较常见，比如lvs流量在早上4点最低，而在晚上11点最高。为了将监控项的周期性考虑进去，我们选取了某个监控项过去14天的数据。对于某个时刻，将得到14个点可以作为参考值，我们记为$x_i$，其中$i=1,…,14$。

我们先考虑静态阈值的方法来判断$\text{input}$是否异常（突增和突减）。如果$\text{input}$比过去14天同一时刻的最小值乘以一个阈值还小，就会认为该输入为异常点（突减）；而如果$\text{input}$比过去14天同一时刻的最大值乘以一个阈值还大，就会认为该输入为异常点（突增）。

静态阈值的方法是根据历史经验得出的值，实际中如何给$\text{max threshold}$和$\text{min threshold}$是一个需要讨论的话题。根据目前动态阈值的经验规则来说，取平均值是一个比较好的思路。

自己的理解：同比的方法只能针对周期性的曲线使用，但是同比的方法个人认为只能用来辅助预测，更多的应该是用来检测异常。

同比振幅(CA)

同比的方法遇到有些情况就不能检测出异常。比如今天是11.18日，过去14天的历史曲线必然会比今天的曲线低很多。那么今天出了一个小故障，曲线下跌了，相对于过去14天的曲线仍然是高很多的。这样的故障使用方法二就检测不出来，那么我们将如何改进我们的方法呢？一个直觉的说法是，两个曲线虽然不一样高，但是“长得差不多”。那么怎么利用这种“长得差不多”呢？那就是振幅了。

怎么计算$t$时刻的振幅呢？ 我们使用$\frac{x(t) – x(t-1)}{x(t-1)}$来表示振幅。举个例子，例如$t$时刻的流量为900bit，$t-1$时刻的是1000bit，那么可以计算出掉线人数是10%。如果参考过去14天的数据，我们会得到14个振幅值。使用14个振幅的绝对值作为标准，如果$m$时刻的振幅$\frac{m(t) – m(t-1)}{m(t-1)}>\text{amplitude} \times \text{threshold}$并且$m$时刻的振幅大于0，则我们认为该时刻发生突增，而如果$m$时刻的振幅大于$\text{amplitude}\times \text{threshold}$并且$m$时刻的振幅小于0，则认为该时刻发生突减。其中，

$$\text{amplitude} = \text{max}[|\frac{x_i(t)-x_i(t-1)}{x_i(t-1)}|],x=1,2,...,14$$

算法组合

上面介绍了四种方法，这四种方法里面，SS和LS是针对非周期性数据的验证方法，而chain和CA是针对周期性数据的验证方法。那这四种方法应该如何选择和使用呢？下面我们介绍两种使用方法：

一、根据周期性的不同来选择合适的方法。这种方法需要首先验证序列是否具有周期性，如果具有周期性则进入左边分支的检测方法，如果没有周期性，则选择进入右分支的检测方法。

上面涉及到了如何检测数据周期的问题，可以使用差分的方法来检测数据是否具有周期性。比如取最近两天的数据来做差分，如果是周期数据，差分后就可以消除波动，然后结合方差阈值判断的判断方法来确定数据的周期性。当然，如果数据波动范围比较大，可以在差分之前先对数据进行归一化（比如z-score）。

二、不区分周期性，直接根据“少数服从多数”的方法来去检测，这种方法比较好理解，在此就不说明了。

其他方法

Smoothed z-score algorithm

主要思想：

1. 利用过去一段历史窗口针对下个节点值做预测（利用平均值，方差信息），若是其超过了一定的阈值，则是个异常点。
2. 对异常点的数值进行平滑，以便评估下下个点是否为异常点。因为不做平滑，由于当前是个异常点，对平均值、方差影响较大，若是下一个点仍是异常点，可能不会识别。

https://zhuanlan.zhihu.com/p/39453139

相空间重构

• 时间序列模型之相空间重构模型

论文集

针对周期型KPI的异常检测算法

• Time Series Decomposition: Yingying Chen, Ratul Mahajan, Baskar Sridharan, and Zhi-Li Zhang. A provider-side view of web search response time. In Proceedings of the ACM SIGCOMM 2013 confere nce on SIGCOMM, pages 243–254. ACM, 2013.
• Holtwinters: He Yan, Ashley Flavel, Zihui Ge, Alexandre Gerber, Daniel Massey, Christos Papadopoulos, Hiren Shah, and Jennifer Yates. Argus: End-to-end service anomaly detection and localization from an isp’s point of view. In INFOCOM, 2012 Proceedings IEEE, pages 2756–2760. IEEE, 2012.
• AutoEncoder（AE）：Anomaly Detection异常检测的几种方法
• Variational AutoEncoder（VAE）：AIOps探索：基于VAE模型的周期性KPI异常检测方法

针对稳定型KPI的异常检测算法

• 静态阈值: Amazon cloudwatch alarm.
• Moving Average: David R. Choffnes, Fabián E. Bustamante, and Zihui Ge. Crowdsourcing service-level network event monitoring. In Proceedings of the ACM SIGCOMM 2010 Conf.
• Weighted Moving Average: Balachander Krishnamurthy, Subhabrata Sen, Yin Zhang, and Yan Chen. Sketch-based change detection: methods, evaluation, and applications. In Proceedings of the 3rd ACM
  SIGCOMM conference on Internet measurement, pages 234–247. ACM, 2003.
• Exponentially Weighted Moving Average: Balachander Krishnamurthy, Subhabrata Sen, Yin Zhang, and Yan Chen. Sketch-based change detection: methods, evaluation, and applications. In Proceedings of the 3rd ACM SIGCOMM conference on Internet measurement, pages 234–247. ACM, 2003.
• ARIMA: Yin Zhang, Zihui Ge, Albert Greenberg, and Matthew Roughan. Network anomography. In Proceedings of the 5th ACM SIGCOMM Conference on Internet Measurement, IMC’05, pages 30–30, Berkeley, CA, USA, 2005. USENIX Association.
• 时间序列的自回归模型—从线性代数的角度来看

针对不稳定型KPI的异常检测算法

• Extreme Value Theory: Siffer A, Fouque P A, Termier A, et al. Anomaly Detection in Streams with Extreme Value Theory[C]//Proceedings of the 23rd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining. ACM, 2017: 1067-1075.
• Wavelet: Paul Barford, Jeffery Kline, David Plonka, and Amos Ron. A signal analysis of network traffic anomalies. In Proceedings of the 2nd ACM SIGCOMM Workshop on Internet measurment, pages 71–82. ACM, 2002.

针对异常数据量太少， 采用异常注入算法

• Fernando Silveira, Christophe Diot, Nina Taft, and Ramesh Govindan. Astute: Detecting a different class of traffic anomalies. In Proceedings of the ACM SIGCOMM 2010 Conference, SIGCOMM ’10, pages 267–278. ACM, 2010.
• Anukool Lakhina, Mark Crovella, and Christophe Diot. Mining anomalies using traffic feature distributions. In Proceedings of the 2005 Conference on Applications, Technologies, Architectures, and Protocols for Computer Communications, SIGCOMM ’05, pages 217–228. ACM, 2005.
• Anukool Lakhina, Mark Crovella, and Christophe Diot. Diagnosing network-wide traffic anomalies. In Proceedings of the 2004 Conference on Applications, Technologies, Architectures, and Protocols for Computer Communications, SIGCOMM ’04, pages 219–230. ACM, 2004

针对不同类型的KPI进行聚类

• Robust and Rapid Clustering of KPIs for Large-Scale Anomaly Detection

工具库

• sklearn.covariance.EllipticEnvelope（高斯分布的协方差估计）
• sklearn.ensemble.IsolationForest（随机森林）
• sklearn.neighbors.LocalOutlierFactor（LOF）

裴丹论文解读

Opprentice

• 智能运维系统（二）

Robust and Rapid Clustering of KPIs for Large-Scale Anomaly Detection

VAE

• AIOps探索：基于VAE模型的周期性KPI异常检测方法

KPI比赛分析

预处理

• 负样本稀少：负样本过采样
• 缺失值：一阶线性插值填充、均值填充
• 样本权重：增加头部异常样本的权重
• KPI三种形态：周期波动/稳定/不稳定（进行聚类：裴丹18论文）
• 异常值替换：格拉布斯准则剔除替换
• 归一化/正则化：z-score、min-max
• 长异常区间只取前8个时间点

特征提取

• 统计特征：均值、方差等
• 对比特征：差分、变分等
• 频域特征：频谱分析、小波分析等
• 拟合特征：移动平均、查分平均、权重平均等
• 原始特征
• 深度特征：AR+[LSTM、Seq2Seq]、AttentionLSTM+CNN
• 其他特征：Seasonal Trend Decomposition、PID（P：误差的值、I：误差的积分值、D：误差的微分值）、tsfresh抽取

模型训练

• DNN+Threshold
• LR+小波分析+随机森林+BiLSTM，加权投票
• xgb

改进：聚类分别训练模型+加入周期特征+加入无监督模型