SLS机器学习介绍（03）：时序异常检测建模-白红宇

文章系列链接

摘要与背景

虽然计算机软硬件的快速发展已经极大提高了应用程序的可靠性，但是在大型集群中仍然存在大量的软件错误和硬件故障。系统要求7x24小时不间断运行，因此，对这些系统进行持续监控至关重要。这就要求我们就被从系统中持续采集系统运行日志，业务运行日志的能力，并能快速的分析和监控当前状态曲线的异常，一旦发现异常，能第一时间将信息送到相关人员手中。因此，使用机器学习和数据挖掘的手段对系统自动化进行异常检测至关重要。

让我们来看时序数据的一些常见异常：

延时的毛刺：在对读写敏感的场景中，经常会有IO毛刺问题困扰，这些毛刺往往被平均难以通过肉眼或P99等监控方法发现

2.业务系统调用量异常：以云服务提供的API为例，因各种原因会有一个短暂的Burst和下降。对于重要接口可以在监控系统上发现，但当我们面向几百个接口的监控时也会力不从心

3.水位上升：新版本发布后，各个指标的形态与历史吻合，但是整体的平均水位有拔高

当业务组合(B)复杂，部署规模（N）变大后，依靠传统的人眼+依赖同比环比等绝对值算法来判断就捉襟见肘了，以下我们提供对这一类问题的分类与解法。

$$ Scale = B * N * NumberOfSeriesData $$

异常类型概述

通常说的异常大致分为异常值、波动值、异常时间序列等几种情况：

异常值（Outlier）
- 给定输入时间序列x，异常值是时间戳值对$(t,x_t)$，其中观测值$x_t$与该时间序列的期望值$(E(x_t))$不同。

波动点（Change Point）
- 给定输入时间序列x，波动点是指在某个时间t，其状态在这个时间序列上表现出与t前后的值不同的特性。

断层异常（Breakout）
- 时序系统中某一时刻的值比前一时刻的值陡增或者陡降很多。

异常时间序列（Anomalous Time Series）
- 给定一组时间序列$X={x_i}$，异常时间序列$x_i \in X$是在X上与大多数时间序列值不一致的部分。

PS：上述归纳的异常均是针对单条时序曲线的形态异常，在工作交流中，发现网络流量的同学不仅要关注流量数值的大小异常、增量异常，同时还要考虑在一定的窗口期内，针对单业务线的各个网络包成分的异常，这些特定场景的异常需要设计特定的算法进行判别。在此不展开讨论，做异常的同学可以私下钉钉我，多多交流！

开源异常检测系统

这里添加了一些我在调研中找到的异常检测算法包，有Java版本、R版本供各位同学参考使用

算法	具体描述
EGADS ExtremeLow 密度模型异常值	EGADS基于密度的异常检测
EGADS CP	EGADS中波动点检测
EGADS KSigmaModel 异常值	EGADS中经典的KSigma模型
Extreme R异常值	开源的异常值检测、阈值绝对值和残差检测异常
Twitter Outliter	基于广义的ESD方法的Twitter R语言异常检测库
BreakOut Twitter CP	Twitter一个基于ESD统计检测来实现变点检测的R语言库
Change Pt1-2-3 R	R语言版本的变点检测库

异常检测方法

时序预测模型概况.png

对时序曲线进行大致的归纳整理后，可以发现，单条曲线形态大致可以分成如下三个场景：

数据无规律波动，但是正常稳定情况下，时序数据大致在一个稳定的区间内波动，通常可以设计一组不同等级的阈值，进行阈值的异常报警;

数据长期波动幅度较大，但是正常情况下，短期内数据波动幅度较小，整体展现出较为平稳的上涨或者下降，这里对应的异常就是断层异常，是否在某个时刻数据陡增或者陡降;

数据有明显的规律性和周期性，通过学习周期变化、涨幅变化得到最终异常数据;

本文围绕着异常检测、变点检测、断层检测这三个主题进行展开，这一期主要利用基于预测的方法进行异常检测。这其中主要展开关于ARMA、Holt Winters和STL模型集中讨论。

基于预测的异常检测

针对单条时序数据，可以利用这篇文章中的方法进行建模预测，根据预测出来的时序曲线和历史数据求出时序的残差，对残差序列建模，利用KSigma或者分位数等方法进行异常检测。

数据平滑

以网络流量曲线为例，该曲线局部都懂的很厉害，若对原始数据建模，则会吸收很很大的噪音，影响模型的效果，可以先用平滑函数对数据进行平滑处理。下图中，蓝色曲线是原始的网络流量数据，紫色是经过长度为5的矩形窗口平滑之后的结果。

* | select ts_smooth_fir( stamp, val, 'rectangle', 5, 1, 'avg') limit 1000

数据窗口平滑.png

为了让大家有对比的看清楚平滑的效果，上图中两条曲线是画在不同的Y轴坐标系下的。先将两条曲线绘制在同一个Y轴下面，计算对应的点的残差序列，可以较好的找到异常的点。如下图中的红色圆圈标记出来的结果。

阈值检测

为了避免单点抖动产生的误报，需要将求取累积的窗口均值进行阈值判别

具体的累积就是通过窗口进行操作

$$ \hat{x}(t) = \frac{x_t + x_{t-1} + ... + x_{t-w+1}}{w} $$

同时业务同学可以根据自己的接受范围对阈值进行分级设置，不同的阈值对应不同的操作

KSigma异常检测

假设数据集由一个正太分布产生，该分布可以用 $N(\mu,\sigma)$ 表示，其中$\mu$是序列的均值，$\sigma$是序列的标准差，数据落在$(\mu-3\sigma, \mu+3\sigma)$之外的概率仅有0.27%，落在$(\mu-4\sigma, \mu+4\sigma)$之外的区域的概率仅有0.01%，可以根据对业务的理解和时序曲线，找到合适的K值用来作为不同级别的异常报警。

分位数(Box Plot)异常检测

箱型图，是一种用作显示一组数据分散情况资料的统计图。主要用于反映原始数据分布的特征，还可以进行多组数据分布特征的比较，其绘制方法是：先找出一组数据的最大值、最小值、中位数和上下两个四分位数。通过不同分位数来划分异常值和疑似异常值。

boxplot

平台实验结果

* | select ts_predicate_simple(stamp, val, 6, 1, 'avg') from log limit 500* | select ts_predicate_arma(stamp, val, 5, 2, 6, 1, 'avg') from log limit 500

基于统计的异常检测

变点的描述

使用数学的方法描述上述状况，给定时序数据，$z_1, z_2,...,z_n$，如果$\exists \tau$，使得$z_1, z_2, ..., z_{\tau}$在某统计特性上区别于$z_{\tau+1}, z_{\tau+2}, ..., z_n$，则$\tau$就是我们需要寻找的Change Ponit。

以均值为统计特性，去分析存在变点的条件：

$$ z_t = \begin{cases} \mu_1 & if 1 \leq t \leq \tau_1 \\ \mu_2 & if \tau_2 < \leq \tau_2 \\ ... & ... \\ \mu_{k+1} & if \tau_k < t \leq \tau_{k+1} = n \end{cases} $$

1. 查找单个变点的模型

假设，时间序列的分布遵循如下模型：

$$ Z_t|\theta_t \sim N(\theta_t, 1) $$

其中，$\theta_t$可以使用一个分段函数进行描述，设计相关的算法找到分段函数的分段点。我们使用如下模型进行求解：

$$ LR = \max \limits_{\tau}\{ l(z_{1:\tau}) + l(z_{\tau+1:n}) - l(z_{1:n}) \} $$

给定一个惩罚系数，当$LR > \lambda$时，去求取如下参数：

$$ \tau = argmax \{ l(z_{1:\tau}) + l(z_{\tau+1:n}) - l(z_{1:n}) \} $$

2. 查找多个变点的模型

假设，时序数据中存在多个变点，这里假设为k个，$\tau = (\tau_0, \tau_1, ... , \tau_{k+1})$，这里$\tau_0 = 0$，$\tau_{k+1}=n$，可以得到如下定义：

$$ \min \limits_{k, \tau} \lgroup \sum_{i=1}^{k+1} \lbrack-l(z_{\tau_{i-1}:\tau_{i}}) \rbrack + \lambda f(k) \rgroup $$

针对上述问题，其解空间巨大：

对于含有N个点的时序数据，存在$2^{N-1}$中可能性

若k是已知的，这组合的空间也是巨大的$C_{n-1}^{k-1}$

常用的方法

局部峰值点检测策略

针对上述流程图的说明:

务必确保滤波之后的序列是平滑的，因为一阶差分对噪声数据很敏感，若噪声未能很好的过滤掉，会产生太多的候选点，影响查找效率

寻找局部峰值点

均值漂移策略：需要确定窗口的大小，具体的计算公式如下所示：

$$ \hat{x}(t) = \frac{x_t + x_{t-1} + ... + x_{t-w+1}}{x_{t-w} + x_{t-w-1} + ... + x_{t-2w+1}} $$

EDM(E-Divisive with Medians)

EDM算法是由Twitter在2014年提出来，具体细节介绍请见论文《Leveraging Cloud Data to Mitigate User Experience from ‘Breaking Bad’》。这里仅仅简单的描述下该算法的特点和原理。

算法特点
- EDM使用E-statstics统计方法来检测平均值的差异，通常，EDM算法也可以用于检测给定时间序列中的分布的变化情况。
- EDM使用较为鲁棒的统计指标，并通过组合测试的方法进行显著性的检验。
- EDM算法是非参数的，很多云端数据并不遵循简单意义上的正太分布，具有较好的适用性。

算法原理
- Parameters: $Z, \sigma, D$
- Let $T_{A}, T_{B}, T_{AB}$ be interval trees with $2^D$ leaf nodes
- // Initialize within distance trees
- for $i \leq i \leq \sigma$ do
  - for $i+1 \leq j \leq \sigma$ do
    - Insert $|Z_i - Z_j|$ to $T_A$
    - Insert $|Z_{i+\sigma} - Z_{j+\sigma}|$ to $T_B$
  - end
- end
- // Initialize between distance tree
- for $1 \leq i \leq \sigma$ do
  - for $1 \leq j \leq \sigma$ do
    - Insert $|Z_i - Z_{j+\sigma}|$ to $T_{AB}$
  - end
- end
- $ = approx - median $
- $bestStat = \frac{\tau(k-\tau)}{k}(2m_1 - m_2 - m_3)$
- $bestLoc = \sigma$
- $\tau = \sigma$
- forwardMove = 0
- // Update trees
- while $\tau \leq n - \sigma$ do
  - if forwardMove = 1 then
    - Perform ForwardUpdate
  - end
  - else
    - Perform BackwardUpdate
  - end
  - forwardMove = 1 - forwardMove
- end
- return bestLoc

其中关于BackwardUpdate和ForwardUpdate等具体操作请详见论文

涉及具体参数的选择等问题论文中有简单的说明，并通过实验不断调试理解

PELT(Pruned Exact Linear Time)

PELT算法是优化有了由Yao和Jackson提出的最优划分方法（The Optimal Partitioning Method），该方法的目标是最小化如下目标：

$$ \sum_{i=1}^{m+1}\lbrack C(y_{(\tau_{i-1}+1) : \tau_{i}}) + \beta) \rbrack $$

令$F(s)$等于上式，对于$y_{1:s}$序列中，存在$\Gamma_s = \{ \tau : 0 = \tau_0 < \tau_1 < \tau_1 < ... < \tau_m < \tau_{m+1} = s \}$是待检测出来的变点集合，令$F(0) = - \beta$，可以得到如下的表达式：

$$ F(s) = \min \limits_{\tau \in \Gamma_s} \{ \sum_{i=1}^{m+1} \lbrack C(y_{ (\tau_{i-1} + 1 : \tau_{i}) }) + \beta \rbrack \} $$

$$ F(s) = \min \limits_{t} \{ \min \limits_{\tau \in \Gamma_{t}} \sum_{i=1}^{m} \lbrack C(y_{(\tau_{i-1} + 1):\tau_i}) + \beta + C(y_{(t+1):n}) + \beta \} $$

$$ F(s) = \min \limits_{t} \{ F(t) + C(y_{(t+1):n}) + \beta \} $$

根据上述公式推导，得到如下的算法描述(Optimal Partitioning)
- Input
  - A set of data of the form, $(y_1,y_2,...,y_n)$ where $y_i \in R$
  - A measure of fit $C(.)$ dependent on the data
  - A penalty constant $\beta$ which does not depend on the number or location of changepoints
- Initialise:
  - Let n=length of data and set $F(0)=-\beta$, $cp(0)=NULL$
- Interate for $\tau^{*} = 1, 2, ..., n$
  - 1. Calculate $F(\tau^{*}) = \min \limits_{0 \leq \tau < \tau^{*}} \lbrack F(\tau) + C(y_{(\tau + 1) : \tau^{*}}) + \beta \rbrack$
  - 1. Let $\tau^{'} = arg\{ \min \limits_{0 \leq \tau < \tau^{*}} \lbrack F(\tau) + C(y_{(\tau+1):\tau^{*}}) + \beta \rbrack \}$
  - 1. Set $cp(\tau^{*}) = (cp(\tau^{'}, \tau^{'}))$
- Output the change points recorded in cp(n)

PELT是将上述算法修改称线性时间，其中设计到一些推论证明，具体的请参照论文《Optimal detection of changepoints with a linear computational cost》

折点描述

在物理学中，折点是不可导点，这表明该店没有切线，所以也不会有速度，在到折点之前，必定要经过一个减速阶段，将速度变成零，然后在向另一个方向加速。

在时序描述中，要发现数据中的尖刺点，同时也要发现时序中的断层点（或者解释成突增、突降点检测）。

1. 峰值分位数检测

将原始时序数据$x(t)$经过平滑后，求序列的一阶差分序列$x_{diff}(t)$，利用如下方法寻找异常值：

使用KSigma方法找到时序中的异常点

使用Quantile方法找到时序中的峰值点

2. 均值漂移检测

给定检测窗口长度$w$，通过下式可以计算$t \in \{2w, N\}$之前每个点的均值漂移值。

$$ r(t) = \frac{x_{t} + x_{t-1} + ... + x_{t-w+1}}{x_{t-w} + x_{t-w-1} + ... + x_{t-2w+1}} $$

根据检测出来$r(t)$值，利用KSigma异常检测方法，可以较好的得到时序中的断层点信息。

平台实验结果

调用命令

* | select ts_cp_detect(stamp, val, 1.0, 1, 'avg') limit 500* | select ts_breakout_detect(stamp, val, 5, 1, 'avg') limit 500

针对CPU用户态使用的5min时序数据两种方法的检测结果如下

针对15min内网络流量时序数据两种方法的检测结果如下

说明
- 平台针对不同方法提供了不同调用形式的接口，用户可以使用默认参数，平台也将重要的参数开放出来供大家调参使用。

硬广时间

日志进阶

阿里云日志服务针对日志提供了完整的解决方案，以下相关功能是日志进阶的必备良药：

机器学习语法与函数:

日志上下文查询：

快速查询：

实时分析：

快速分析：

基于日志设置告警：

配置大盘：

更多日志进阶内容可以参考：。

联系我们

纠错或者帮助文档以及最佳实践贡献，请联系：悟冥

问题咨询请加钉钉群：

文章系列链接

摘要与背景

异常类型概述

开源异常检测系统

异常检测方法

基于预测的异常检测

数据平滑

阈值检测

KSigma异常检测

分位数(Box Plot)异常检测

平台实验结果

基于统计的异常检测

变点的描述

1. 查找单个变点的模型

2. 查找多个变点的模型

常用的方法

局部峰值点检测策略

EDM(E-Divisive with Medians)

PELT(Pruned Exact Linear Time)

折点描述

1. 峰值分位数检测

2. 均值漂移检测

平台实验结果

硬广时间

日志进阶

联系我们

参考文献