0 引言

时间序列作为数据挖掘问题中一种重要的数据形式，是聚类领域的重点研究对象之一。从变量数目的角度看，时间序列可以分为单维时间序列和多维时间序列(Multivariate Time Series, MTS)。MTS是多维变量按照时间顺序所记录的一系列观察值的集合，广泛存在于商业、金融业、航空业、社会学、生物学等众多应用背景中，MTS数据集有其重要的数据特点：1)MTS的各维之间可能存在一些固有的时序关系，简单地拆分和降维容易造成关键信息的丢失；2)MTS样本之间长度往往不等，故与样本长度相关的预处理算法不利于后期的聚类分析；3)每个MTS样本由多个单维时间序列组成，这些时间序列的数据类型可能不同。以上特点增加了MTS数据集的聚类难度，使得面向MTS数据集的聚类算法有别于一般的方法。

现有的时间序列聚类算法主要包含三类^[1]：基于原始数据的聚类、基于特征的聚类和基于模型的聚类。然而在这些聚类算法中，用于MTS数据集的聚类算法相对较少。文献[2]利用单维时间序列的聚类思想，给多维时间序列的各个维度赋予特定的权值，每个行向量作为一个时间点。由于MTS样本长度不等，样本之间的相似度使用动态时间弯曲(Dynamic Time Warping, DTW)度量，最佳匹配路径上每一对时间点的多维向量之间的距离利用闵可夫斯基参数模型计算。该算法需要领域知识为各个变量赋予权值，且DTW距离度量方法的计算量较大。文献[3]提出基于变量相关性的MTS特征表示方法，通过协方差反映系统中各个参数的相关关系，将MTS样本转化为协方差矩阵；MTS集所有的协方差矩阵拼接为综合协方差矩阵，对该协方差矩阵进行主成分分析得到各MTS的特征矩阵。该方法可以将数值型不等长MTS数据集转变为大小相同的特征矩阵集合，处理结果可用于聚类分析。文献[4]提出了一种基于参数交互关系的MTS聚类方法，指出MTS中的任一维变量都可以被其他解释变量近似线性组合表示，且将一维线性关系纳入了考虑范畴，假定这些变量间的线性相关关系可以用来进行聚类，其不足之处在于模型计算时间会随着样本数量变大而增加，也不能处理非数值型变量。文献[5]将每一维时间序列转化为一个统计特征数组，MTS样本由各维变量统计特征数组拼接成的向量来表示。该算法可以处理不等长时间序列，但要求各维选取的统计特征必须一致导致其在处理混合型MTS数据集时会遇到困难。文献[6]针对MTS数据集存在的样本之间不等长、数据类型多样和噪声等问题，提出了一种基于协方差矩阵与测地线距离(geodesic-based distance)的MTS聚类算法。该算法首先将MTS样本转化为协方差矩阵；然后将协方差矩阵从黎曼空间映射到欧氏空间；最后对矩阵集进行聚类。如果使用基于距离的聚类算法，上述映射过程可以省略，协方差矩阵之间的距离度量方法使用测地线距离。Zhou等^[7]提出了一种基于模型的多维时间序列聚类算法——MUTSCA〈LR〉(Multivariate Time Series Clustering Algorithm 〈Lift Ratio〉)，该聚类算法假设目标数据集由一系列概率分布模型系统生成，不同的系统将生成相异的多维时间序列。该算法先将连续型数值符号化；然后在符号化样本上计算由LR(Lift Ratio)向量表示的时序模式，将时序模式累加生成用来表示MTS样本的模型向量；最后对模型向量集进行聚类。它不需要特定的领域知识，同时可以处理包含数值和非数值型变量的混合型MTS数据集。但实际应用中的MTS数据集，各样本的长度往往不等，该算法在处理不等长多维时间序列时需要利用DTW^[8]或滑动窗口^[9]来度量模型向量之间的相似性，造成算法的时间开销较大。

针对上述算法在处理MTS数据集时的不足，本文将在MUTSCA〈LR〉的基础上进行改进，提出一种基于LR分量提取的多维时间序列聚类算法(MUltivariate Time Series Clustering Algorithm based on LR Component Extraction，MUTSCA〈LRCE〉)。该算法使用等频离散化方法对MTS数据集进行符号化；在符号化后的MTS样本每一维时间序列之间计算时序模式向量LR，对LR向量进行排序并提取向量两端固定数目不同数值的分量，同一MTS样本内的提取分量拼接成一个模型向量(Model Vector, MV)，MTS数据集由MV向量的集合表示；最后对模型向量集使用k-means算法进行聚类分析。实验结果表明，改进算法可以快速地完成不等长MTS数据集的聚类分析。

1 相关知识背景 1.1 符号表示与相关定义

定义1  多维时间序列集。一个包含N个多维时间序列样本的多维时间序列集可表示为S={S¹, S², …, S^N}。若每个样本Sⁱ包含m个变量，每一维时间序列长度为n，则Sⁱ可表示为Sⁱ ={S₁ⁱ, S₂ⁱ, …, S_mⁱ}，其中第k维时间序列为S_kⁱ={S_{k, 1}ⁱ, S_{k, 2}ⁱ, …, S_{k, n}ⁱ}。文中符号化后的多维时间序列集则表示为S^*={S^1*, S^2*, …, S^N*}。

定义2  时序模式。用LR向量表示的时序模式的计算方法^[6]如下：

$ \begin{array}{l} Lift(S_{j, t-\tau }^{i{\rm{*}}} \Rightarrow S_{j, t}^{i{\rm{*}}}) = \frac{{{\rm{Pr}}(S_{j, t}^{i*}|S_{j, t-\tau }^{i*})}}{{P(S_{j, t}^{i*})}} = \\ \;\;\;\;\;\;\frac{{freq(S_{j, t-\tau }^{i*} \cap S_{j, t}^{i*})}}{{freq(S_{j, t}^{i*})*freq(S_{j, t - \tau }^{i*})}};1 \le \tau < t \end{array} $

(1)

其中：τ是常数，表示时间偏移量; ${S_{j, t}^{i*}} $与${S_{j, t-\tau }^{i*}} $分别表示第i个样本的第j维时间序列在t和t-τ时刻所对应的符号值；$ {P(S_{j, t}^{i*})}$表示$ {S_{j, t}^{i*}}$的发生概率；${{\rm{Pr}}(S_{j, t}^{i*}|S_{j, t-\tau }^{i*})} $表示当${S_{j, t-\tau }^{i*}} $发生时${S_{j, t}^{i*}} $发生的条件概率；${freq(S_{j, t}^{i*})} $表示${S_{j, t}^{i*}} $在时间序列S_j^i*中出现的次数；t是时间变量，与MTS长度n有关，因此该公式将生成一个长度为n-τ的向量，称为Lift Ratio(LR)。它反映了每个时刻t-τ的$ {S_{j, t-\tau }^{i*}}$对${S_{j, t}^{i*}} $的依赖程度，以此反映样本S^i*的变量j的时间序列S_j^i*内部的时序模式。

式(1)反映的是同一维时间序列S_j^i*内部的LR计算公式，当式(1)中的$ {S_{j, t}^{i*}}$换成$ S_{k, t}^{i*}\left( {k \ne j} \right)$时，则表示不同维的时间序列S_j^i*与S_k^i*之间的LR计算公式。同一变量序列内部计算的LR称为intra-pattern，不同变量时间序列之间的LR称为inter-pattern。

1.2 MUTSCA〈LR〉算法

MUTSCA〈LR〉算法首先使用MDD(Mode-Driven Discretization)算法^[10]对MTS样本进行符号化，MDD符号化过程如下：MDD算法将MTS作为由若干个向量组成的矩阵，行为向量，列为变量维；首先选取MTS中独立冗余度(Multiple interdependence Redundancy, MR)值最大的变量维作为标签属性，标签维每一行的符号值为该行向量的分类标签；然后利用有监督的离散化算法OCDD(Optimal Class-Dependent Discretization)^[11]对其他数值型变量维进行符号化。在符号化后的MTS样本上利用式(1)计算时序模式LR并累加得到模型向量。最后，利用k-means算法对模型向量集进行聚类分析。MUTSCA〈LR〉算法模型向量计算方法如下：

输入  S^*={S^1*, S^2*, …, S^N*}，其中S^i*={ S₁^i*, S₂^i*, …, S_m^i*}；

输出  N个模型向量组成的向量集finalLR。

For每个多维时间序列S^i*

  For每个变量S_j^i*

    利用式(1)计算intra-pattern;

    For每个变量S_k^i*(k≠j)

      利用式(1)计算inter-pattern;

      finalLR_i+=inter-pattern;

    End

    finalLR_i+=intra-pattern;

  End

  finalLR=finalLR∪{finalLR_i};

End

MUTSCA〈LR〉算法存在两个问题：

1) MDD算法中采用的OCDD利用动态规划的思想选取分割点集，若处理长时间序列需要较多的离散化符号，则要求候选分割点集包含较多的元素，导致OCDD的开销过高，离散化执行效率低下。

2) 由于j, k∈{1, 2, …, m}，每个样本计算得到m²个时序模式向量LR，多维时间序列S^i*最终生成的时序模型finalLR_i由m²个LR向量累加求和得到。如果多维时间序列集S^*中的各个样本S^i*长度不等，则MUTSCA〈LR〉生成的模型向量集finalLR中的各个finalLR_i(i∈{1, 2, …, N})长度也不相等。上述问题增加了相似性度量的难度，造成该算法在聚类过程的耗时较长。

2 MUTSCA〈LRCE〉算法 2.1 多维时间序列的等频离散化

针对MDD算法效率较低的问题，文中采用等频离散化(Equal Frequency Discretization, EFD)算法进行符号化。EFD是一种简单的离散化方法，它需要事先给定一个参数来决定离散化后最终的离散符号个数，记为num_bin。在没有领域知识的情况下，各个样本的num_bin往往难以确定。传统的等频离散化方法需要用户随机给出num_bin的取值，如此会导致MTS聚类分析的效果不稳定。num_bin的取值同MTS样本长度有关，但MTS数据集各样本长度差异过大会导致num_bin的值域范围较宽。鉴于以上原因，本文基于样本长度和变异系数提出一个用来为样本Sⁱ(i∈{1, 2, …, N})选取num_bin值的计算方法：

$ num\_bin = ({\mathop{\rm int}} )C*\sqrt n ;num\_bin \le n/2 $

(2)

其中：C为由MTS样本集S确定的唯一常数，n是Sⁱ的时间序列长度。输入样本集S={S¹, S², …, S^N}，统计S中各MTS长度的均值M和标准差Ve，计算变异系数C_V=Ve/M，C为变异系数的倒数1/C_V。从式(2)可知，当MTS样本集各样本长度的C_V取值越大，参数C的取值越小，从而缓解各MTS样本离散符号数目差异大的问题。文中的等频离散化方法具体描述如下。

输入  样本集S={S¹, S², …, S^N}，其中Sⁱ ={S₁ⁱ, S₂ⁱ, …, S_mⁱ}；

输出  符号化后的样本集S^*={S^1*, S^2*, …, S^N*}。

统计样本集各个时间序列的均值和方差，并计算系数C；

For样本集每一个多维时间序列Sⁱ

  利用式(2)计算Sⁱ的num_bin;

  For所有变量序列S_jⁱ(j∈{1, 2, …, m});

    根据num_bin的取值，对该S_jⁱ进行等频离散化，得到S_j^i*;

    S^i* =S^i*∪S_j^i*;

  END

  S^*=S^*∪S^i*;

END

2.2 模型向量集的生成方法及其聚类

首先在样本S^i*每一对变量的时间序列(S_v1^i*, S_v2^i*)(v1, v2∈{1, 2, …, m})上计算时序模式LR向量；然后将LR向量进行排序并选取排序后向量两端K个不同数值的分量，数值较大的分量反映了时序模式的主要特征，数值较小的分量反映了时序模式的特殊状态；最后S^i*的全部提取分量将组成模型向量MV_i。根据v1、v2的取值，样本S^i*利用式(1)计算出m²个LR向量，每个LR向量中提取K个分量，因此MUTSCA〈LRCE〉生成的模型向量MV_i长度L仅与分量个数K、参数个数m相关，L=m²*K。模型向量集生成方法描述如下。

输入  离散化后的样本集S^*={S^1*, S^2*, …, S^N*}，S^i* = {S₁^i*, S₂^i*, …, S_m^i*}，分量个数K;

输出  模型向量集MV={MV₁, MV₂, …, MV_N}。

For每一个多维时间序列S^i*

  For多维时间序列的所有参数v1(1≤v1≤m)

    For多维时间序列的所有参数v2(1≤v2≤m)

      利用式(1)计算LR向量;

      对LR进行排序，从中提取首尾K个数值不同的分量，并拼接到MV_i;

    End

  End

  MV=MV∪MV_i;

End

文中采用k-means算法对模型向量集MV进行聚类分析，首先随机选取k个模型向量作为初始簇中心，计算模型向量与各簇中心的欧氏距离，并将该向量分配给最相似的簇，MV分类完毕后更新簇中心。重复上述步骤，直至分簇结果不再变化，输出分簇结果向量KM=[km₁, km₂, …, km_N]，其中km_i(km_i∈{1, 2, …, k})表示样本Sⁱ的分簇编号。聚类过程如下。

输入  模型向量集MV={MV₁, MV₂, …, MV_N}，分簇个数k；

输出  分簇结果向量KM=[km₁, km₂, …, km_N]。

从MV中随机选取k个向量作为初始簇中心集合Core={core₁, core₂, …, core_k};

构建一个临时簇中心集合Core′ ={core′₁, core′₂, …, core′_k};

构建一个数组Num[k];

初始化KM=[0, 0, …, 0];

Do

初始化Core′为0向量集，Num[k]为0向量;

For模型向量MV_i(1≤i≤N)

  计算每个模型向量与Core的k个簇中心的相似度，并将它分配到最相似的簇θ，其中θ∈{1, 2, …, k};

  km_i=θ;

  core′_kmi +=MV_i;

  Num[km_i-1]++;

End

利用Core′与Num[k]更新簇中心，结果返还给Core; Until KM no change

3 实验与结果分析

算法使用Java进行编程实现，实验在一台配备Intel四核3.80 GHz处理器、4 GB内存、装有Window 7系统的PC上进行。

3.1 实验准备与数据集介绍

选用4个来自UCI的MTS数据集：EMGPAD(EMG Physical Action Data set)、EMGLL(EMG dataset in Lower Limb)、AReM(Activity Recognition system based on Multisensor data fusion)、DSAD(Daily and Sports Activities Data set)，其中AReM与DSAD为等长MTS数据集，详见表 1。

1) EMGPAD。样本数目80，包括3位男性和1位女性实验者。每个实验者做20个动作(20个样本)，包括10个攻击性动作和10个一般动作，样本长度大多在10000左右。实验使用动作的性质为标签，即攻击性和非攻击性。

2) EMGLL。样本数目66，包括11位膝关节患者和11位正常人。每个实验者做3种运动，各样本长度波动较大。实验使用自然人的分类作为标签，即患者和正常人。

3) AReM。样本数目87，包括7个类型一的弯腰动作、5个类型二的弯腰动作、15个骑车动作、15个躺动作、15个坐动作、15个站立动作、15个走路动作，各样本长度均为480。实验使用动作的类型作为标签，共7种类别。该数据集是为了验证改进算法MUTSCA〈LRCE〉在等长MTS数据集上可以维持MUTSCA〈LR〉的聚类效果。

4) DSAD。样本数目9120，由8位实验参与者完成19种动作，每个动作包含60个样本。该数据集旨在检验算法MUTSCA〈LRCE〉在大样本集下的聚类效果。

实验过程中参数设置如下，LR计算的时间延时τ的值取5，k-means聚类算法的簇中心个数k取值为数据集的类别个数，提取排序后LR向量的首尾分量个数K取10。由于实验数据集有标签，故采用F-measure和信息熵作为实验中多维时间序列聚类算法的评价指标。

为了便于算法MUTSCA〈LR〉与MUTSCA〈LRCE〉对比，本实验使用两种方法对MUTSCA〈LR〉进行修改使其能够处理不等长MTS数据集：1)使用DTW计算不等长模型向量之间的距离，采用文献[12]提出的基于DTW的全局平均法进行k-means聚类中心点的更新；2)使用滑动窗口计算不等长模型向量之间的距离，簇中心点的更新方法如下：假设某次k-means迭代过程产生的簇中心集合为Core={core₁, core₂, …, core_k}，length(core_i)表示向量core_i的长度。设core_i所在簇中有n_i个模型向量，在簇中心core_i(1≤i≤k)的更新过程中，创建两个长度为length(core_i)的0向量，记为Sum和weight。对于core_i簇中的某模型向量MV_j(1≤j≤n_i)，若length(MV_j)≥length(core_i)，使用滑动窗口在MV_j上截取与core_i最相似的子序列累加至Sum、weight各分量数值加1；若length(MV_j) < length(core_i)，使用滑动窗口在core_i上找到与MV_j最相似子序列的起始位置st，将MV_j累加到Sum[st]至Sum[(st+length(MV_j))-1)]，且weight中相应位置的分量加1。更新后簇中心core_i′的计算公式为：core_i′=Sum[t]/weight[t](0≤t < length(core_i))。

3.2 本文的等频离散化方法评估

采用包含较长MTS样本的数据集EMGPAD和EMGLL对本文EFD与MDD算法进行评估。两种离散化方法的处理结果均使用上述基于滑动窗口的MUTSCA〈LR〉算法进行聚类分析。根据文献[10]的经验值，算法MDD对每个样本Sⁱ的分割点候选集元素个数取值为$\sqrt {length\left( {{S^i}} \right)/3} $，其中length(Sⁱ)表示样本Sⁱ的长度。本文测试了该参数与MDD符号化性能之间的关系，设实验选用的候选分割点集包含元素个数$ y = \theta *\sqrt {length\left( {{S^i}} \right)/3} $，计算结果进行“上取整”处理作为实际选取元素个数。这里用系数θ表示候选分割点集的大小，实验结果如图 1所示，图 1中横坐标表示候选分割点集大小θ。由图 1可以看出，随着候选分割点集元素个数增加，算法MDD得到的符号化MTS样本集聚类评价提高，但离散化执行时间会迅速增长。

表 2给出了算法EFD与MDD的测评结果，对于样本Sⁱ，MDD使用的候选分割点集元素个数取值为$ \sqrt {length\left( {{S^i}} \right)/3} $。从聚类评价指标值可以看出，样本集EMGPAD与EMGLL使用EFD离散化方法所得到的符号化样本集在进行聚类分析时效果较好，其F-measure值分别为0.896与0.904，而MDD算法相应的结果为0.870与0.742。其中，MDD在数据集EMGLL上符号化效果较差的现象与候选分割点集包含元素数目不足有关，但由图 1的实验结果可知，通过增加候选分割点集元素个数来提高聚类评价结果的方法时间代价过高。从离散化时间的角度看，本文EFD算法的执行时间分别为4065 ms和5348 ms，明显低于MDD算法的执行时间。

3.3 MUTSCA〈LRCE〉算法评估

使用4个MTS数据集对改进算法MUTSCA〈LRCE〉进行评估，如表 3所示，其中MUTSCA〈LR〉使用MDD进行样本集符号化，MUTSCA〈LRCE〉则采用改进的EFD算法。聚类时间包括模型计算时间和k-means执行时间，由表 3可知，在不等长MTS数据集EMGPAD与EMGLL上，改进后的MUTSCA〈LRCE〉因包含向量排序过程，故时间略长，这里认为算法改进前后的模型计算时间基本一致。算法MUTSCA〈LRCE〉的k-means执行时间分别为35 ms和16 ms，而原算法k-means的执行时间明显较长。这是由于算法MUTSCA〈LRCE〉生成的模型向量长度仅取决于MTS参数个数m以及分量个数K，模型向量长度固定，使得算法聚类速度较快。此外，基于滑动窗口的MUTSCA〈LR〉在数据集EMGPAD上执行时间少于EMGLL的原因是EMGPAD中各样本长度波动较小，模型向量之间的相似性度量包含的窗口滑动次数较少。从聚类结果来看，MUTSCA〈LRCE〉与MUTSCA〈LR〉在EMGPAD上的聚类效果相当，而在数据集EMGLL上算法MUTSCA〈LR〉的F-measure值在0.75附近，由实验3.2节的分析可知其在该数据集上失效的原因在于该数据集存在长度较短的样本，它们的候选分割点集中元素较少，导致MDD离散化效果不佳，影响了后续的聚类分析。

数据集AReM是样本数目较少的等长多维时间序列数据集，由于样本长度较短，聚类时间相互之间区分度低，参考价值不大。从聚类结果的角度看，三种算法的熵和F-measure值基本一致。数据集DSAD为包含较多等长多维时间序列样本，该部分实验选取了4800个样本，由于本文算法生成的模型向量维度与变量个数有关，模型向量维度L大于样本长度，其中L=m²*K，样本长度为125。因此模型向量维度高，聚类时间较长。同时，较高维度的模型向量包含更多特征，故聚类效果优于算法MUTSCA〈LR〉。所以，算法MUTSCA〈LRCE〉在多标签的等长MTS数据集上仍然有效。

为了评估聚类簇个数k值对k-means聚类结果的影响，本文从数据集DSAD中选取两组样本子集进行实验，实验结果如表 4所示。第一组样本子集包含样本数目240个，分别来自4种不同类型动作的样本集，标签数目为4。第二组样本子集包含1920个样本，实验以同一动作类型的同一个行为主体的60个样本为一组，抽取32组，所以样本子集标签数目为32。

由表 4可以看出，当k值小于标签个数时聚类效果较差，当参数k略大于标签数时，聚类效果较好。这是因为k值较小时，不同标签的样本容易被合并到同一类簇中。当k值略大于样本标签数目时，k-means分簇更精细，噪声样本对聚类过程的影响降低。

4 结语

本文基于时序模式Lift Ratio向量的MTS表示方法提出了MUTSCA〈LRCE〉算法，该算法利用改进的等频离散化方法对MTS进行符号化，通过LR向量分量提取的方式将不等长MTS样本转化为等长的模型向量。实验结果表明本文算法可以更好地对不等长MTS数据集进行聚类分析。时间序列是数据之间具有严格上下文关系的一类特殊数据对象，LR向量所展现的MTS时序模式仅反映了时间点之间的时序关系，如何利用时间段之间的时序关系进行MTS聚类并减少时间段处理过程中造成的信息丢失有待进一步研究。