现如今，大型软硬件系统为了确保能够安全、可靠地运行，需要有专业的运维人员来进行系统的部署、运行和维护。尤其是随着互联网产业的发展，互联网公司业务的稳定运行和公司的利益密切相关，服务的不稳定甚至崩溃往往会导致巨大的经济损失。

从运维工作发展的历史中来看，运维工作是在朝着自动化、高效化方向不断发展的。尤其随着 AIOps的出现 ，让运维管理模式从ITIL向DevOps演化。相比而言，AIOps作为一个全新的技术发展和应用方向，更适合运维，而且更有格局，甚至有人认为AIOps 是运维的归宿。

本期“智汇华云”专栏为您带来“AIOps之异常检测——SSA模型详解”。

下面，就给大家讲解一下我们用到的SSA算法，在检测指标异常中扮演重要角色。

SSA的全称是Singular Spectrum Analysis，中文是奇异频谱分析，是一种时间序列分解技术。SSA可以把时间序列分解成一组成分，这些成分组合起来就是趋势trend、周期periodicity和噪音noise。趋势性和周期性成分就可以视为指标的主要成分，计算这些成分和指标之间的残差，再利用Grubbs检验就可以判定出异常值。

为了展示SSA的魔力，我们先定义一个玩具时间序列，包含趋势、周期和噪音成分：

这是一个抛物线趋势函数，两个周期函数和一个噪音函数组成的时间序列。

SSA的第一步就是将时间序列F映射为多维滞后向量序列。取L作为窗口长度，2<=L<=N/2。我们沿着时间序列滑动这个窗口，得到一组列向量Xi：

这些列向量叠到一起，就形成了L-轨迹矩阵：

可以看到，反对角线上的元素就是相等的，这种矩阵也叫Hankel矩阵。

对于我们的玩具矩阵，我们将窗长设为70，我们后面还会讨论如何选择合适的窗长。K=N-L+1表示轨迹矩阵的列数。我们把X的列向量叫做L滞后向量，行向量叫做K滞后向量。

SSA的第二步就是使用奇异值分解singular value decomposition (SVD)来分解轨迹矩阵，

U是一个L*L的酉矩阵，包含了X的左奇异向量正交集合作为列；

Σ是一个L*K的矩形对角矩阵，包含了X的L个奇异值，从大到小排列；

V是一个K*K的酉矩阵，包含了X的右奇异向量正交集合作为列。

轨迹矩阵的SVD便可以重写成：

σi是第i个奇异值，Ui和Vi是U和V的第i列向量，d<=L是轨迹矩阵的秩，Xi是X的第i个初等矩阵。d是时间序列轨迹空间的本征维数，我们可以把前r<d个初等矩阵加起来，得到一个X的低维近似。

上图便是X的前12个初等矩阵，可以看出，X0和X1的L滞后向量和K滞后向量在矩阵中变化较慢，暗示了它们可能是时间序列的总趋势。X2到X5均是棋盘形状，暗示了周期性。X6可能在周期性和趋势之间。X7之后看起来是在某些值之间快速变换，比较可能是原时间序列的噪音。

我们可以把前12个初等矩阵的对于轨迹矩阵的相对贡献和累积贡献画出来：

可以看到X0和X1对于X的展开分别贡献了52%和22%。前7个初等矩阵合计贡献了97%。当重建时间序列的时候，做出差不多贡献的初等矩阵要被集合起来。比如X2和X3，X4和X5应该被集合到一起。

在理想情况下，时间序列F的所有成分都是可分的，我们可以把X分解成j个初等矩阵，每个Xj都应该是Hankel结构的。但是在现实情况下，没有一个初等矩阵的反对角线的值是相等的。所以，我们需要把初等矩阵转换为Henkel矩阵，再转换为时间序列。

在这里，我们使用对角线平均，可以表示为：

左边是Hankel矩阵，右边是Hankel化运算符乘上初等矩阵。

Hankel矩阵中的元素如下，s=m+n

因为Hankel化运算符是线性运算符，所以对于轨迹矩阵X，

Hankel化之后的初等矩阵如下图所示：

图上可以很容易得出以下结论：

我们再把Hankel化的初等矩阵转换为时间序列：

我们从图上可以看出F0和F1看起来像趋势，F2和F3几乎重合，验证了我们想把它们集合起来的想法。我们把上面得到的集合结论画出来：

在我们集合了初等成分之后，可以看出SSA的分离原时间序列的效果很不错，特别是周期性成分。

我们上面是用眼睛观察来确定哪些成分需要被集合在一起的，对于复杂的时间序列，我们不能老用眼睛来判断，下面介绍一种集合判定的方法。

对于两组时间序列，长度N，窗长N，我们定义一个加权内积：

我们再定义一个d*d的加权相关矩阵Wcorr，测量了成分Fi和Fj的w正交偏差：

如果Fi和Fj相近，分子分母就接近，Wi,j就接近于1。如果Fi和Fj正交，Wi,j就为0。我们可以定义Wi,j>0.3，表示两个成分需要被集合在一起。

我们可以把成分0-6的Wcorr画出来：

可以看出F0和F1的W是0.4，说明它们应该被组合。F1和F6的W是0.39，说明F6应该和F0F1组合起来成为趋势成分。F6也和F4F5有一点相关，但是因为F4F5和F0F1没有相关性，所以我们还是只把F6和F0F1放在一起。

根据W相关性矩阵，我们可以得出这样的集合：F(周期性1)=F2+F3, F(周期性2)=F4+F5，因为W2,3=0.98，W4,5=0.98。

我们试着选择不同的L来探索一下窗长参数对时间序列分解的影响。

L=2

对于这么小的窗长，SSA开始从时间序列中分离出高频的噪音，给出了一个稍微去噪的F0。

L=5

我们可以看到F0是一个去噪很不错的版本了。F1是一个分解比较差的周期性成分，F2到F4只是噪音。

L=20

我们可以看到这时趋势和周期成分开始显现形状了。F0看起来像是抛物线趋势，F1+F2+F3算是比较好看的周期性。F3看上去比较讨厌，兼具周期性和噪音，这告诉我们需要增加窗长，看看能不能获得一个更好的周期与噪音的分离。

L=40

有趣的是，F0在L=40的时候恶化了，周期性在开始和结束的时候也恶化了，说明我们还需要进一步增加窗长。

L=60

现在我们已经接近了我们最开始的选择L=70。

W相关矩阵可以被分为两块：F0到F6和F7到F59。F0到F6是趋势和周期成分，F7到F59是噪音成分：

可以清楚的从图上看出，F1包含了和周期成分F2F3一样的震动频率，导致了明显的W相关性。当我们把窗长增加到L=70后，我们就能完全移除这些震动。当然我们也可以把窗长增加到L=N/2=100，但是超过窗长超过70后，前7个本征成分就对窗长不敏感了。

在SSA中，没有一个严格标准去选取窗长。但是需要更长的窗长（时间序列长度的30-45%）从趋势中充分分离出周期性成分。

总结一下SSA算法，主要步骤如下：

1. 对于一个时间序列F和窗长L，造一个轨迹矩阵X

2. 使用奇异值分解去分解X，

3. 得到d个初等矩阵

4. 对角线平均Xi去造出初等时间序列成分Fi

5. 对于Fi计算出加权相关性矩阵Wcorr

6. 聚合并分类Fi，最终得到时间序列的趋势性、周期性和噪音成分

https://www.kaggle.com/jdarcy/introducing-ssa-for-time-series-decomposition

https://medium.com/@jonathan_hui/machine-learning-singular-value-decomposition-svd-principal-component-analysis-pca-1d45e885e491