基于新的鲁棒相似性度量的时间序列聚类

李国荣，冶继民，甄远婷

（西安电子科技大学数学与统计学院，西安 710126）

（*通信作者电子邮箱jmye@mail.xidian.edu.cn）

0 引言

随着计算机技术的发展，时间序列数据作为一种重要的数据形式存在于生活各个领域，例如：医学、金融股票、天气数据、图像分割等［1］。面对这些庞大的数据集，如何对其进行有效处理，挖掘出有效信息已经成为了当今研究的热点话题。时间序列聚类是其中最典型的一种方法，其目的是将未标记的时间序列数据划分到不同的组，并确保组内序列的相似程度最大、组间序列相似程度最小，从而识别出序列集合的结构，提取需要的重要信息［2］。由于在收集这些数据的过程中，会遇到很多不稳定因素，导致数据集中存在一些异常值点。这些异常值会直接影响到最终的聚类结果，从而影响后续对数据进行分析、预测的结果，最终可能做出错误的决策。因此，采取适当的方法消除异常值对聚类结果的影响非常有必要。

目前已有一些处理异常值的方法，例如：彭柳青等［3］将样本加权思想与子空间聚类算法相结合，为每一个样本分配反映离群程度的尺度参数，提出了一种鲁棒的子空间聚类算法；曹晨曦等［4］将统计的方法应用到时间序列数据中，提出了一种基于统计方法的异常值检测与校正方法；Grané 等［5］提出了一种基于小波的通用的异常值检测与校正的方法；Lafuente-Rego 等［6］通过分别考虑度量、噪声、分段方法，并基于比较样本分位数自协方差，提出了三种模糊C-mesoids 模型，三种方法从不同的角度对异常值进行了处理。但是这些工作基本都需要对数据进行处理，而处理过程中涉及的特征提取、参数选择等工作都很容易造成数据中有用信息的丢失，从而导致聚类结果不准确；另外在对异常值进行修正时也可能会产生较大误差，影响最终的聚类结果。

为了避免上述方法的缺点，本文通过增强算法过程中使用的相似性度量的鲁棒性，提出了直接基于原数据的鲁棒时间序列聚类方法。传统的衡量时间序列间相似性的度量，如：欧氏距离、马氏距离、闵可夫斯基距离、Pearson 相关系数、动态时间规整（Dynamic Time Warping，DTW）距离等，除DTW 距离外，其余的相似性度量一般都没有很强的鲁棒性，当时间序列数据中存在异常值点时这些度量不再能准确衡量出序列间的相关关系。DTW 距离虽然有较强的鲁棒性，但其计算复杂度非常高，导致聚类算法效率降低。本文在文献［7］中提出的广义互相关（Generalized Cross-Correlation，GCC）度量的基础上进行改进，提出了一种鲁棒广义互相关度量。该度量是非参的且不需要对数据进行任何处理过程，在一定程度上可以减少对数据进行处理与修正过程带来的误差。

1 广义互相关度量

GCC是衡量两个平稳时间序列之间线性相关性的一般度量［7］。设xt、yt为两个平稳时间序列，Ryx，k是(yt，yt-1，…，yt-k，xt，xt-1，…，xt-k)的协方差矩阵（k为时间滞后），其对角线部分为Ryy，k和Rxx，k，上三角部分为，下三角部分为Cxy，k。这里的Ryy，k为(yt，yt-1，…，yt-k)的协方差矩阵、Rxx，k为(xt，xt-1，…，xt-k)的协方差矩阵、Cxy，k是由两个时间序列间的互相关系数构成的矩阵，并且有|Ryx，k|=|Rxx，k||Ryy，k-则两序列间的广义互相关度量为：

2 鲁棒广义互相关度量

上述的协方差矩阵Ryx，k中的元素是利用的Pearson 相关系数计算出来的，Pearson 相关系数对异常值点敏感，从而导致利用它构造出来的序列间的相似性度量GCC 也对异常值敏感，而生活中大多数的时间序列数据中都存在异常值点，所以选择适当的方法增强聚类算法中使用的相似性度量的鲁棒性很有必要。

一种可行的提高GCC 鲁棒性的方法就是采用随机变量之间的鲁棒相关系数代替Pearson 相关系数构造GCC。该方法从另一个角度即相关系数的鲁棒估计来增加时间序列间相似性度量的鲁棒性，可以解决现有的时间序列聚类算法不适用于存在异常值的数据集的问题，而现有的方法中很少有从此角度进行改进的。

文献［8-9］中提出了一种基于Qn统计量的鲁棒相关系数，Qn统计量的表达式为：

大量的仿真实验发现，原始的Qn统计量存在一定的小样本偏差。小样本偏差可以通过给Qn统计量乘一个小样本校正因子［11］校正，即校正后的Qn统计量为（为方便，仍采用原记号）：

其中调节因子dn的具体取值如下。

1）当n≤9时，调节因子的取值如表1所示。

表1 调节因子dn的取值Tab.1 Value of adjustment factor dn

2）当n＞9时：

注：当样本量n≥10 000 时，Qn将不需要再乘小样本校正因子。

基于校正后的Qn统计量构造鲁棒相关系数

式中：med(x)为随机变量x的中位数；med(y)为随机变量y的中位数。

计算向量Zt，2(k+1)的协方差矩阵：

得到两个时间序列间的鲁棒广义互相关（Robust Generalized Cross-Correlation，RGCC）度量：

采用AIC 准则（Akaike Information Criterion）和BIC 准则（Bayesian Information Criterion）选取时间滞后k的值［7］。自回归模型中，k一般取1 或2；动态因子模型中，k一般在1～3 中取值。

另外，可以验证RGCC(xt，yt)也满足距离度量的基本性质，即：

1）对称性：RGCC(xt，yt)=RGCC(yt，xt)；

2）0 ≤RGCC(xt，yt) ≤1；

3）当且仅当序列之间有精确的线性关系时，RGCC(xt，yt)=1；

4）当且仅当所有的互相关系数为零时，RGCC(xt，yt)=0。

证明 1）、2）易见，这里仅给出3）、4）的证明。

对于3），已知当且仅当存在一个线性组合aTZt=0时，有|=0，即这两个序列是线性相关的［7］。而RGCC(xt，yt)=1则一定有|=0，所以性质（3）成立。

下面证4），首先对RGCC(xt，yt)变形

根据Hadamarsd’s 不等式，右边的行列式小于等于1，且等式成立等价于矩阵为对角矩阵，即=0xy，k。 证毕。

3 仿真实验

这一部分本文用两个常见的时间序列模型生成仿真数据集，给模型生成的数据集中随机加入异常值点，分别用GCC与RGCC 计算其距离矩阵，作为层次聚类算法的输入对数据进行聚类，并对聚类的结果进行分析。

本实验用3 个常见的聚类评价指标比较基于两种度量聚类得到的聚类结果。假设U 是真实的聚类结果，V 是某次实验得到的聚类结果。

1）调整的兰德系数（Adjusted Rand Index，ARI）［13］，ARI∈[ -1，1]且ARI 的值越大代表聚类结果越好。具体计算表达式为：

其中RI（Rand Index）为兰德系数，计算公式为：

这里的a等于在U 中属于同一组，在V 中也属于同一组的对象的对数；b等于在U 中属于同一组，在V 中属于不同组的对象的对数；c等于在U 中不属于同一组，在V 中属于同一组的对象的对数；d等于在U 中不属于同一组，在V 中也不属于同一组的对象的对数。E(RI)是兰德系数的期望。

2）F-measure，F-measure 是精确率P与召回率R的调和平均［14］，F∈[0，1]，且F的值越大代表聚类结果越好。具体的计算表达式为：

3）Jaccard 系数（Jaccard Coefficient，JC），JC 是用来比较两个有限集合间的相似性与差异性的度量，JC∈[0，1]且JC 的值越大，集合间的相似性越高，代表聚类的结果越好。具体的计算表达式为：

用F-measure 与Jaccard 系数评价聚类的结果时，可以得到与ARI指标一样的结论，但由于ARI指标的变化更加直观，所以本文实验结果中仅列出ARI 指标的值。另外，本文各表中给出的实验结果均为实验100次取平均的结果。

本文的实验环境为Windows 7 系统下8 GB 内存的Intel 3.40 GHz PC，Matlab2017a软件。

3.1 自回归模型生成的仿真数据

用简单的AR（1）模型（AutorRegressive model）生成一组线性、平稳的时间序列数据：

i=1，2，…，5 时，ai=0.9；i=6，7，…，15 时，ai=0.2。通过对误差项εi的限制引入时间序列间的交叉依赖关系，r(i，j)=E(εi，t，εj，t)。

情形1

情形2

情形3

情形4

情形5

从原始AR（1）模型看，该模型生成的时间序列数据存在明显的两大类，即前5个序列相关，后10个序列相关。但当引入交叉依赖项εi后，序列间的相关关系发生了变化。情形1、3中显然有6大类，即前10个序列相关，后5个序列独立；情形2、4、5中显然有两大类，即前10个序列相关，后5个序列相关。

表2 为数据不存在异常值的情况下，T=100、200（T为每个时间序列的长度）时，GCC(xt，yt)与RGCC(xt，yt)的实验结果对比。

表2 无异常值时，GCC、RGCC性能对比Tab.2 Performance comparison of GCC and RGCC without outliers

不存在异常值的情况下，从表2 可以看出两个相似性度量得到的聚类结果没有明显的差异。T=100时，单链情况下GCC 的结果更精确一些，全链情况下RGCC 的结果更精确一些，均链接的情况下两个度量的结果几乎一致；T=200时，基于两个相似性度量得到的结果均更加精确，ARI 指标的值相差不大。

给模型在五种场景下生成的仿真数据中随机加入一些异常值点，再分别利用GCC 与RGCC 对存在异常值点的数据进行聚类。仿真实验结果见表3。图1 为在情形1 中，使用单链接方式用两种度量进行聚类的层次聚类图。

表3 存在异常值时，GCC、RGCC性能对比Tab.3 Performance comparison of GCC and RGCC with outliers

从表3 可以看出，在时间序列数据中存在异常值时，三种链接方式下的五种情形中，RGCC 得到的聚类结果均比GCC的聚类结果更精确。

情形1 中，理想的聚类结果应该是第1～10 个时间序列为一组，第11、12、13、14、15个序列各自单独成一组，从图1（a）、图1（b）可以看到，由于异常值点的干扰，GCC 不能识别出时间序列间的某些相关关系，导致第9、10、12、14 个时间序列被错分。

图1 同一数据集上，GCC、RGCC进行聚类得到的树形图对比Fig.1 Comparison of dendrogram obtained by clustering with GCC or RGCC on the same dataset

3.2 动态因子模型生成的仿真数据集

这部分本文使用文献［15］中提出的三种场景，即带有组因子结构的动态因子模型（Dynamic Factor Model，DFM）作为数据生成的过程。考虑三个组，每个组有3 个特定的组因子r1=r2=r3=3 和相同的数量的序列N1=N2=N3=100。DFM可以表示为：

这里的N=300，T=100。G={g1，g2，g3}表示每组中的成员，Nj(j=1，2，3)表示每组中序列的个数，xit是一个80×1 的可观测变量，其每个元素均是由在区间[-2，2]上的均匀分布产生，β=(1，2，3，0，…，0)1×80T。是一个rj× 1维的不可观测的组特异因子向量，其每个元素由均值为j、方差为1 的高斯分布生成的。的每个元素是由高斯分布N(0，j)生成。通过对εi，t的性质做不同的假设，产生三种情形。

情形1 误差项εi，t服从标准高斯分布。

情形2 εi，t=其中，如果t 为奇数，则δt=1；否则，δt=0；并且 ε1=与 ε2=服从均值向量为0、协方差矩阵为Σ=(σi，j)，σi，j=0.3|i-j| 的多元高斯分布。

情形3 εi，t=0.2εi，t-1+ ei，j，这里et=(e1，t，e2，t，…，eN，t)′服从均值向量为0、协方差矩阵为Σ=(σi，j)，σi，j=0.3|i- j|的多元高斯分布。

给模型生成的三种情形的数据中随机加入一些异常值点，分别用GCC与RGCC计算距离矩阵，作为层次聚类算法的输入进行聚类。表4 为在该模型下的两种度量的聚类结果对比。观察表4 中得到的ARI 指标的值很直观可以看到，对动态因子模型生成的带有异常值点的时间序列数据进行聚类，RGCC仍可以得到比GCC更精确的聚类结果。

表4 动态因子模型中，GCC、RGCC性能对比Tab.4 Performance comparison of GCC and RGCC in dynamic factor model

通过对两个仿真数据集的聚类结果分析可知，本文提出的序列间的相似性度量RGCC 不仅在数据无异常值时，可以发现序列间的相关关系，而且在数据中存在异常值点时仍有效。

4 结语

本文提出了一种基于相关系数鲁棒估计的时间序列间的鲁棒相似性度量——RGCC。仿真实验结果显示，针对存在异常值点的时间序列数据，基于RGCC 进行聚类得到的结果，明显比基于GCC 进行聚类得到的结果更接近真实的聚类结果。该度量具有很强的鲁棒性，更加适用于现实生活中的时间序列数据集。此外，基于RGCC 对存在异常值点的时间序列数据聚类时，不需要异常值检测与处理的过程，可以有效避免查找异常值点过多或者过少，或者对异常值点进行修正时产生的误差；同时，该度量的计算不需要提前设置任何参数，避免了参数选取不当对聚类结果的影响。

由RGCC 的结构可知，在小样本数据情况下，选取不同的小样本校正因子会导致相关矩阵的变化，从而影响到距离矩阵的计算，最终造成聚类结果的不准确。因此，未来的工作可以从优化小样本校正因子选取的角度对本文提出的度量进行改进。