基于相似性連接的時間序列Shapelets提取

張振國 王 超 溫延龍 袁曉潔

1(延邊大學計算機科學與技術系 吉林延吉 133002)2(南開大學計算機學院 天津 300350)3 (南開大學網絡空間安全學院 天津 300350) (zhangzhenguo@dbis.nankai.edu.cn)

近年來，對于時間序列數據(time series，簡稱時序數據)的研究與分析受到越來越多的關注，已成為數據挖掘領域的一個重要的研究課題[1-2]，其原因在于時序數據廣泛地存在于我們日常生活的諸多領域，如醫療、金融、運動軌跡以及天文星系觀測等[3].一般來說，時序數據是由某一觀測量隨時間變化產生的，每一數據值之間有嚴格的時間先后順序，這是時序數據區別于其他類型數據的本質特征.廣義上，任何具有嚴格先后順序的數值構成的序列，都可以稱為時序數據.作為一項基礎工作，時序數據的分類是被研究最多的問題之一，具有重要的應用價值，如在心電圖信號分析中，對其進行分類能夠快速鑒別信號的正常與異常，有助于輔助醫療；在工業生產中，對設備性能數據的采集與分類，能夠及時了解設備的運行狀況.相比于其他數據的分類，在提取時序數據特征，構建分類器的過程中，數據之間的先后順序是必須要考慮的關鍵因素.因此，以時序數據整體作為處理對象，利用數據間相似性的最近鄰分類器(nearest neighbor classifier, NN)被首先用于時序數據的分類問題[4-5].

基于NN分類器的時序分類算法的核心在于如何計算時序數據間的相似性.由于噪音、采集時間上的對齊等因素，時序數據在獲取過程中可能會存在一些輕微的幅值偏移，這就使得在計算時序數據間的相似性時，相似性算法必須能夠處理這種時間軸上的數據錯位現象.有2種類型的相似性計算方式被廣泛應用：1)基于動態時間規整(dynamic time warping, DTW)的距離計算[6]；2)基于編輯距離(edit distance)的相似度計算[7].DTW及其相關改進算法，如導數DTW(derivative DTW)[8]、權重DTW(weighted DTW)[9]，通過尋找一條有效的規整路徑來最小化對應的時序數據點間的總距離并借此計算2條時序數據的相似性.DTW已被認為是時序數據相似性計算的基準[5].同DTW一樣，編輯距離也通過尋找數據間的最佳匹配來計算相似度，所不同的是，基于編輯距離的方法，如時間規整編輯(time warp edit, TWE)[10]、移動-分裂-融合(move-split-merge, MSM)[11]，允許數據間存在跨越和改變.這2種類型的時序數據相似性計算被稱為彈性距離度量(elastic distance measures).在很長一段時間內，使用這些彈性距離度量的簡單NN算法在時序數據分類問題上占據主要地位且很難被超過[12].然而，由于NN算法需要存儲和搜索整個時序數據集，而這些彈性距離度量的計算又比較耗時(通常情況下，時間復雜度為O(m2)，m為時序數據長度)，這就導致在處理較大時序數據集時，時間和空間的消耗都比較大，限制了其應用.另一方面，這種方法是將時序數據作為整體來考慮的，然而更多的時候，對于時序數據類別區分起關鍵作用的是局部信息而非全局信息，同時，局部信息更能體現出時序數據間的本質差異.

為了解決上述2個問題，Ye等人[13]提出使用子序列(subsequence)作為區分不同類別的依據，并將具有高類別可區分的子序列稱為Shapelets.相比于使用NN分類器的算法，Shapelets提供了良好的可解釋性，易于體現不同類別時序數據之間的關鍵差異.因此，Shapelets一經提出便受到研究者的關注.如何快速提取高質量的Shapelets已成為最近幾年時序數據處理領域的一個熱點問題.然而，由于任意長度的時序子序列都可能是Shapelets，而子序列數量龐大(例如長度為60，時序數量為600的數據集，子序列數量達1.098×106)，所以逐一判斷子序列的類別可分能力是一項非常耗時的工作.針對該問題，Ye等人[13]提出使用剪枝策略加速子序列分類能力的判斷.這些策略對于小數據集是可行的，但對于大數據集仍難以處理.為了提高Shapelets提取效率，許多使用加速技術的算法被提出，如使用空間換時間的Logical Shapelets算法[14]、借助時序數據的離散化SAX(symbolic aggregate approximation)[15]表示的Fast Shapelets算法[16]等.這些方法雖然能夠在一定程度上加速Shapelets的提取，但仍需要大的時間消耗，且會損失一定的預測準確率.除此之外，Hills等人[17]更改子序列分類能力度量方法，使用易于計算的Kruskal-Wallis度量和F-statistic作為衡量標準，但加速效果并不明顯，并且降低了分類準確率.

除了上述直接使用提取的Shapelets構建分類器外，Lines等人[18]提出了一種使用Shapelets構建對應時序數據特征向量的方法，可以使用當前比較成熟的分類器，如SVM，而不必拘束于Ye等人[13]使用的決策樹分類器.使用時序數據的特征向量表示，一方面可以保證Shapelets的可解釋特性，另一方面可以充分利用成熟分類器的優點，提高分類準確率.該方法的關鍵問題仍然是如何提取出高質量的Shapelets，然而作者并沒有改變辨別子序列分類能力的方式，故其方法仍然非常耗時.

為此，本文在Yeh等人[19]所提出的時間序列matrix profile基礎上，提出了一種基于相似性連接的時間序列數據Shapelets提取方法.該方法舍棄了現有工作中先生成子序列，再對子序列進行區分能力判斷的方案，而是通過差異向量體現不同類別時序數據本身的差異，利用相似性連接策略快速計算出隱含高分類能力Shapelets的候選矩陣，然后通過對比候選矩陣中的數值差異，篩選出Shapelets集合.以子序列為單位的時序數據間的相似性可以通過快速傅里葉變換計算得到，大大減少了距離計算帶來的時間消耗.在獲得Shapelets集合之后，利用這些Shapelets將原時序數據集轉換成相應的特征向量表示，使用SVM分類器完成分類工作.本文在多個領域的真實數據集上進行一系列實驗，對比了目前主要Shapelets提取算法.實驗結果表明：本文提出的Shapelets提取算法具有很好的時間效率，同時，又能保證很高的分類準確率，與其他算法相比具有更好的實用性.

1 相關工作

目前的研究工作中對于Shapelets的提取思路主要有2種：基于搜索的策略[13-14,16-17,20-23]和基于學習的策略[24-25].

1.1 基于搜索策略的Shapelets提取

在Shapelets提取最初的研究中，其基本思路是考慮訓練數據中的所有子序列，然后使用信息增益(information gain, IG)來評估每一個子序列對于類別的區分能力，通過循環獲取當前訓練數據中具有最優信息增益度量的子序列來構建決策樹分類器，實現時序數據的分類[13].類似于信息增益度量，一些新的度量方式，如F-Statistic，Kruskall-Wallis等，也被應用于衡量子序列的區分能力，但從效果上來看，信息增益要優于這些度量方式[17].

由于時序數據子序列數目眾多，逐一進行判斷的brute-force算法需要大量的運行時間，其時間復雜度為O(n2m4)(n為訓練集中時序數據的條數)，因此，相應的剪枝方法和加速技術是研究的重點.基于最小距離的提前終止計算和基于信息增益的熵剪枝是最早提出的策略[13]，這些策略雖然在一定程度上可以加速搜索過程，但仍然十分耗時.Mueen等人[14]提出以空間換時間的方法，使用5個充分統計量使得在距離計算過程中可以部分重用之前的結果.另外，使用新的度量方式(信息增益上界)來減少耗時的熵的計算.通過這些技術，算法的時間復雜度被減少至O(n2m3).在文獻[14]的基礎上，原繼東等人[20]提出利用Shapelets之間的邏輯組合關系來提高分類準確率，但并沒有提升Shapelets提取效率.為進一步壓縮搜索空間，Rakthanmanon等人[16]借助于時序數據符號化表示，將時序數據離散化，在低維空間中通過隨機映射方法快速過濾掉類別區分能力低的子序列.對于篩選出的有可能成為Shapelets的子序列，則繼續使用信息增益加以判斷.該方法能夠大大減小搜索空間，使得時間復雜度降低至O(nm2).

除了使用加速技術改進算法外，其他方式的策略也被用于加快Shapelets提取.Chang等人[21]借助GPU運算對不同的子序列使用并行化方法來判斷子序列分類能力.然而，這一方面需要對距離計算方式做較大的修改，另一方面，其Shapelets的搜索時間非常依賴于硬件性能.Wistuba等人[22]基于有區分能力的Shapelets應當在時序數據中經常出現的假設，使用隨機采樣的方式構成Shapelets.這種策略確實可以減少運行時間，但分類準確率相對較低.Zhang等人[23]提出沒有數值變化的子序列成為Shapelets的可能性很小，只利用有數值變化的子序列構成Shapelets候選空間.這種方法大大減小搜索空間，有很好的加速效果，但本質上并沒有改變算法時間復雜度.

1.2 基于學習策略的Shapelets提取

與搜索整個子序列空間不同，基于學習的策略試圖尋找針對Shapelets提取的數學表示作為時序數據分類的目標函數，通過優化算法從時序數據集中學習到Shapelets.以這種方式得到的Shapelets有可能不是時序數據中的子序列，但確能將不同類別的時序更好地區分.Grabocka等人[24]使用logistic回歸分類模型構建Shapelets提取的目標函數，稱為LTS(learning time series)，以隨機梯度的方式不斷地優化目標函數，使損失將至最低，進而構造出Shapelets.該方法得到的Shapelets往往不是某條時序數據中的子序列，而是能夠對數據進行劃分的特征.目前，LTS仍然是在不使用集成策略前提下的分類效果最好的算法[26].但是這種策略的最大問題是時間和空間消耗都很大，在硬件資源受限情況下無法得到分類結果.另一種基于學習的方法是由Hou等人[25]提出的FLAG(fused lasso generalized eigenvector method)算法，與LTS直接優化產生Shapelets不同，FLAG通過構建目標函數來尋找具有好的分類能力的數據所在的位置信息，然后提取這些位置上的子序列形成Shapelets.這種方法的優點是速度快，但是分類準確率有一定的降低.

2 相關定義與符號表示

定義1. 時序數據.一條時序數據是一組實數值變量組成的向量，表示為T=(t1,t2,…,tm)，其中m為時序數據的長度.

定義2. 子序列.給定長度為m時序數據T，子序列S是T中的一個連續片段，記為S=(tp,tp+1,…,tp+l-1)，其中p為子序列在T中的起始位置，l為子序列長度，1≤p≤m-l+1.

長度為m的時序數據T中包含m-l+1條長度為l的子序列，可以使用滑動窗口逐一取出.

定義3. 全子序列集合.時序數據T的全子序列集合是指由長度為l的滑動窗口在T上生成的子序列集合，記為A={S1,l,S2,l,…,Sm-l+1,l}.為了表示方便，本文使用A[i]表示Si,l.

時序數據間的相似性通常是通過距離大小來度量的，長度相等的2條時序數據T與R之間的距離

(1)

式(1)只能度量長度相等的時序數據，而本文主要的處理對象是子序列，也就是在不等長情況下度量序列間的相似性，為此，需要定義子序列與時序數據間的距離表示.

定義4. 子序列與時序數據的距離.給定時序數據T(全子序列集合為A)及子序列S，兩者之間的距離定義為子序列S在T中最優匹配時的歐氏距離：

sd(T,S)=min(d(A[1],S),d(A[2],S),…,
d(A[m-l+1],S)).

(2)

圖1直觀地描述了距離sd(T,S)中最小值的含義.

Fig.1 Illustration of best matching in T for subsequence S圖1 子序列S在T中的最優匹配示例

盡管式(2)易于理解，但由于子序列數量龐大，其計算需要大量的運行時間，Yeh等人[19]提出借助均值、方差等充分統計量和點積運算來降低計算的復雜度，即S與T中每一條子序列之間的計算為

(3)

其中，μA[i]，μS，σA[i]，σS分別為子序列A[i]，S的均值和方差.通常情況下，對于每一條子序列，均值和方差的計算時間復雜度為O(l)，但可以通過存儲T中數據值和數據值平方的累積和向量，使得在每一次運算中都可以直接用來計算任意長度的子序列的均值和方差[27].這樣，在式(3)中，A[i]與S的點積運算就成為制約計算效率的關鍵因素，3.2.1節將給出快速計算該點積操作的算法.

定義5. 最近鄰函數.給定2條時序數據的全子序列集合A和B，對于任意一子序列對A[i],B[j]，如果A[i]在B中的最優匹配是B[j]，則最近鄰函數θ值為1，否則為0，即:

(4)

其中,k∈[1,m-l+1]且k≠j.使用最近鄰函數，可以得到2條時序數據的相似性連接向量.

定義6. 相似性連接向量.給定2條時序數據的全子序列集合A和B，對于A中的每一條子序列A[i]，使用最近鄰函數在B中尋找其最優匹配B[j]，每一個匹配對之間的距離所構成的向量稱為A與B之間的相似性連接向量，記為PA B.

需要注意的是，PA B是有序的，在計算過程中嚴格按照第1條時序的全子序列集合A中子序列的先后順序尋找在B中的最近鄰，也就是說A中的子序列一定會在θ=1的最優匹配對中出現，而B中的子序列不一定全部在θ=1的匹配對中，因此，PA B并不是對稱的，即PA B≠PBA.

定義7. 自相似性連接向量.在計算PA B時，若B=A，這樣構成的相似性連接向量稱為自相似性連接向量，記為PA A.

PA A可看作表示一條時序數據的元數據.由于每一條子序列都來自該時序數據，所以總能找到距離為0的最優匹配以及在最優匹配附近近似等于0的匹配子序列，這樣的匹配對于描述數據特性是沒有意義的，通常稱之為平凡匹配(trivial matching)[28].因此，在計算過程中，需要避免這樣的匹配.2條時序數據之間的相似性可以通過PA A和PA B表現出來，為此，我們定義差異向量來描述時序數據之間的差異.

定義8. 差異向量.對于2條不同時序數據的全子序A和B，則兩者的差異向量DiffA B定義為以A為基準的相似性連接向量PA B和A的自相似性連接向量PA A之間差的絕對值，形式化表示為：若

(5)

在時序數據集中，如果同一類的時序數據之間在DiffA B的某些分量上值比較小，而在不同類的數據上值比較大，則可認為相應的子序列是具有較高的類別可分性，是比較好的Shapelets.在第3節中，將詳細描述如何使用這些定義提取高質量的Shapelets.

3 Shapelets提取

本文所提出的基于相似性連接的時序數據Shapelets提取主要由3部分構成：1)數據的預處理過程，主要完成對時序數據集的精簡工作，只保留對Shapelets提取起關鍵作用的訓練數據以加速后續處理；2)通過計算相似性連接向量、自相似性連接向量及差異向量來生成Shapelets候選矩陣，然后從候選矩陣中提取相對應的子序列，并對這些子序列進行合并操作實現Shapelets提取，這是本文的核心；3)完成時序數據的轉換工作，借助已提取的Shapelets將時序數據轉換為由距離值構成的特征向量表示.對于時序數據的特征向量表示，使用現有成熟的分類器，如SVM，完成分類工作.本文方法的整體框架及流程如圖2所示：

Fig. 2 The overall framework and whole process of the proposed method圖2 本文方法的整體框架及流程

3.1 關鍵時序數據的選擇

在時序數據集中，往往存在著一些高度相似的數據，其全子序列集合也非常相近.由于使用相似性連接向量計算差異向量的過程相對復雜，因此，對原始時序數據集進行精簡，只選擇相似度比較小的關鍵時序數據進行計算是十分必要的.Wistuba等人[22]基于Shapelets在時序數據集中出現的頻率高的假設所做的工作也驗證了只選用關鍵時序數據提取高質量的Shapelets是可行的.

由于在選擇關鍵時序數據時不僅要求整體上有極大的相似度，還需要子序列之間也要高度相似，因此本文使用不具有任何偏移處理的歐氏距離(式(1))計算時序數據間的相似程度.通常情況下，同類別中存在高相似性的數據的可能性大，而不同類之間，時序數據的差異較大.基于此，本文按類別分別進行關鍵時序數據的選擇.對于每一類中的時序數據，使用隊列作為核心數據，計算隊首時序與其他時序的距離，若距離小于預先設定的閾值，則認為兩者有高度相似性，然后提取隊首時序作為關鍵時序數據，而與隊首時序相似的時序就可以舍棄.具體算法如算法1所示：

算法1. 關鍵時序選取算法.

輸入：時序數據集dataSet、相似度閾值mthres；

輸出：關鍵時序數據集nDataSet.

①tempDataSet=partitionByClasses(dataSet)；

②nDataSet=[]；

③ for eachclassDataintempDataSet

④queue=index(classData)；*使用隊列存儲當前類時序的索引*

⑤ whilequeueis not empty

⑥pos=distance(classData,queue,mthres);

⑦queue=remove(queue,pos);

⑧nDataSet=[nDataSet;classData(queue[1])];

⑨ end while

⑩ end for

算法1所消耗的運行時間是很少的.首先，算法在選取關鍵時序數據時是以類別為單位的，不同類的時序數據間不進行距離計算，這使得參與相似度計算的時序數據條數大大減少.其次，在每次迭代過程中，算法在只保留1條關鍵時序數據的同時，從隊列中刪除與之相似的數據，減少了下次迭代過程中參與計算的時序數量.具體地說，若時序數據集(大小為n)的第i類中含有ni條時序，則在最好情況下，所有時序均相似，每一條時序只參與一次距離計算(ni-1次)；在最壞情況下，所有時序均不相似，需進行ni(ni-1)2次距離計算，故算法1在處理第i類時序數據時距離計算的次數介于ni-1至ni(ni-1)2之間.因此，算法1的距離計算次數約為n～n(ni-1)2.通常情況下，ni要比n小得多且數據集中存在大量的相似時序數據，所以算法1具有很高的時間效率.

3.2 候選矩陣

候選矩陣(candidate matrix)是本文基于相似性連接策略提取時序數據Shapelets的核心，候選矩陣中隱含著可用于表示類別差異的Shapelets，其值表示相連接的時序在相應子序列上的差異，值越大，對應子序列區分能力越強，因此，可從候選矩陣中快速定位并提取出高類別可分性的子序列.候選矩陣生成的時間消耗是影響算法整體效率的關鍵因素，如何快速的生成候選矩陣是本文的重點研究內容.

計算候選矩陣的關鍵是相似性連接，而相似性連接屬于全匹配問題，即需要計算出時序數據全子序列集合中每一個子序列在其他集合中的最近鄰子序列.候選矩陣的計算分為3個部分：1)計算距離向量，即子序列與時序數據全部子序列間的歐氏距離；2)根據相似性連接策略計算兩條時序數據全子序列集合間的最優匹配，得到兩者之間的相似性連接向量；3)利用2條時序數據的相似性連接向量和其中1條時序數據的自相似性連接向量構造差異向量，將由同一條時序數據的自相似性連接向量得到的差異向量合并，即可組成候選矩陣.

3.2.1 距離向量的快速計算

距離向量是某一子序列S與一條時序數據T全部子序列之間的歐氏距離構成的向量，其核心是2個子序列距離的計算問題.然而，由于子序列數量龐大，逐個計算子序列間的距離會消耗大量的運行時間，降低算法效率.為此，本文根據Rakthanmanon等人[27]提出的方法，使用式(3)的計算方式，實現子序列與時序數據全子序列集合之間距離的快速計算.式(3)中最關鍵部分是子序列之間內積的運算，Mueen等人[29]已經證明，借助快速傅里葉變換，可以使得子序列與時序數據全部子序列之間的內積在1次計算中完成.距離向量的具體實現如下：

算法2. 距離向量計算算法distanceVector.

輸入：時序數據T及長度m、子序列S及長度l；

輸出：S與T的距離向量D.

① [μT,μS,σT,σS]=meanStd(S,T);*根據文獻0計算均值方差*

②S′=reverse(S)；*將S翻轉*

③Ta=append(T,m)；*對T擴充，填m個0*

⑥Taf=FFT(Ta);

⑧ST=inverseFFT(STt);

⑨D=distance(l,ST,μT,μS,σT,σS).*根據式(3)計算距離并組成距離向量*

算法行⑧得到的ST是子序列S與T中每一個子序列點積所構成的有序向量，即ST[i]表示S與T中第i條子序列A[i]的點積.因此，在使用式(3)計算距離向量(行⑨)時，只需要從中取出相應的值，而不用每次都計算內積.該算法的一個優勢在于算法的時間復雜度與子序列長度無關，這也就不存在最好情況與最壞情況的差異，其主要時間消耗在于傅里葉變換(O(mlbm))，故算法的運行效率高.

3.2.2 提取相似性連接向量

對于2條時序數據T，R(其全子序列集合分別為A，B)，可以迭代使用算法2計算A中每一條子序列與R的距離向量.T和R的相似性連接向量PA B就可以通過查找這些距離向量的最小值得到.在計算過程中，本文采用在每次計算距離向量之后更新PA B的方式，逐步得到與T中每一條子序列最優匹配的R中的子序列，記錄其距離和位置，具體算法如下：

算法3.PA B提取算法similarityJoinVector.

輸入：時序數據T，R及長度m、子序列長度l；

輸出：PA B.

①B=allSubseqSet(R);

② initializePA B;

③index=1;*B中子序列索引變量*

④ for eachS∈B

⑤D=distanceVector(T,m,S,l);

⑥PA B=updating(PA B,D,index);

⑦ end for

在算法3中，PA B是以T為基準計算的，每一次的迭代中，算法都更新T相應位置上的距離值，直至B中所有子序列計算完成.這種更新方式的好處是不必開辟內存空間去存儲所有子序列對之間的距離，減少了空間消耗.如果算法3的輸入中R替換為T，則算法得到的是T的自相似性連接向量PA A，但必須要處理平凡匹配導致的距離為0的問題.本文采用Yeh等人[19]的處理方式，直接跳過相關區域的計算，且不更新PA A中對應位置上的值.圖3直觀描述了算法3生成PA A的過程.當PA A初始化為Inf后，按順序計算每一條子序列與T的距離向量并更新PA A，為了方便展示，本文只選取第i條和第j條子序列作為示例來更新PA A.在得到A[i]與T的距離向量(以Di表示)之后，按照逐元素取最小值的策略更新PA A對應位置上的值(圖3中第1行和第2行)，得到更新后的PA A.當遇到平凡匹配的情況時，即距離向量中值為0的情況，則保留PA A原來的值.對于A[j]，思路相同，所不同的是生成的距離向量(以Dj表示)比較的對象是上次更新后的PA A(圖3中第3行和第4行).

Fig. 3 The process of PA A updating圖3 PA A更新過程

3.2.3 候選矩陣生成

由3.2.2節可知，PA B是以子序列為單位描述兩時序數據的，其相似性可以通過PA B與PA A的差值來體現，即差異向量DiffA B.然而由于Shapelets是具有類別可分性的子序列，體現的是類別差異，因此，需要考慮數據集中的全部時序數據，將得到的差異向量進行綜合分析，為此，本文將所得到的差異向量進行組合，形成候選矩陣，用于后續處理.

① http://www.cs.ucr.edu/～eamonn/time_series_data/

一個關鍵的問題是并不是所有時序數據之間的差異向量對于提取Shapelets都是有用的.如果2條時序數據具有相同的類別標識，其差異向量描述的是兩者在數據生成過程中的差異，并不是體現不同類別時序數據的本質區別.因此本文以類別為依據，在計算差異向量時，從不同類別中選取時序數據，而同類別時序數據間不進行計算.若某一時序數據集中包含3個類別，分別為Ⅰ，Ⅱ，Ⅲ類，差異向量的計算過程為：從Ⅰ類中選出1條時序數據T，以T為基準依次計算其與Ⅱ類和Ⅲ類中所有時序數據的差異向量.對T與每一類時序數據得到差異向量進行逐點求平均值，得到平均差異向量，將其合并到候選矩陣中，之后對Ⅰ類中的其他數據進行相同的操作.然后再以Ⅱ類中時序數據為基準計算其與Ⅲ類中時序數據的平均差異向量.本質上，平均差異向量從統計意義上表示了T與某一類別時序數據間的不同.候選矩陣生成的具體算法過程如下：

算法4. 候選矩陣生成算法.

輸入：關鍵時序數據集nDataSet、子序列長度l；

輸出：候選矩陣M.

①tempDataSet=partitionByClasses(nDataSet)；*按類別劃分關鍵時序數據集*

② initializeM；*初始化候選矩陣*

③ for eachcurrClassintempDataSet

④ for eachT∈currClass

⑤PA A=similarityJoinVector(T,T,l);

⑥ initializeZ；*初始化差異矩陣*

⑦ for eachotherclassaftercurrClass

⑧ for eachR∈otherclass

⑨PA B=similarityJoinVector(T,R,l);

⑩Z=[Z;abs(PA B-PA A)];

在候選矩陣生成過程中，需要對數據集中所有的時序對計算PA B和PA A，算法4本身的復雜度較高，但由于本文在預處理階段借助算法1大大減小了數據集規模，僅使用關鍵時序數據用于計算，因此，在實驗中不會消耗很多的運行時間.候選矩陣中的每一行代表著1條時序數據與另一類時序數據的平均差異，使用平均值的目的在于剔除這一類時序數據之間的差異，只使用共性的特征來提取Shapelets.

3.3 基于候選矩陣的Shapelets提取

從3.2節可知，候選矩陣中的每一個值表示1條子序列與1類時序數據中最優匹配的差異，其值的大小代表該子序列區分能力，因此，Shapelets可以通過提取候選矩陣中較大的值對應的子序列得到.為了更形象地描述Shapelets提取過程，本文首先以2條時序數據為例，展示其相似性連接向量、自相似性連接向量及差異向量的計算過程，提取出能夠表示兩者之間差異的子序列，然后將其擴展到整個時序數據集上.在差異向量上選取可能成為Shapelets的子序列時，閾值的選擇對于Shapelets的數量有很大影響.閾值設置過大，會導致僅有少數子序列構成Shapelets，不能夠全面體現時序間的差異；閾值過小，則會使得大量可區分性不明顯的子序列也被加入到Shapelets中，這一方面會導致后續處理計算量的增加，另一方面會對分類過程形成干擾.在實驗中，本文通過交叉驗證的方式來確定合適的閾值.時序數據的Shapelets提取地過程如圖4所示.

2條示例時序數據T和R選自于UCR archives①的UWaveGestureLibraryAll時序數據集，且屬于不同的類別，如圖4(a)所示.假設T和R的全子序列集合分別為A和B，則由算法2得到的PA B與PA A如圖4(b)所示.從圖4(b)中可以看出，PA B與PA A之間在某些位置上存在著明顯的差異，這些差異所對應的子序列就可以將T和R進行區分.圖4(c)給出了兩者的差異向量，根據其值的大小可以很容易得到有區分能力的子序列所在的位置，如圖4(c)中的點P，然后在原始時序數據T中即可將相應的子序列A[P]提取出來.在構成Shapelets時，若子序列位置之間的距離小于子序列長度(如圖4(c)中以P和Q開始的子序列)，本文將其合并形成Shapelets，如圖4(d)所示.

如果由算法1所提取的關鍵時序數據集中只包含2類數據且每類中僅有1條時序數據，按照圖4的流程就是Shapelets提取的全過程，差異向量DiffA B即為候選矩陣.但在一般情況下，數據集中會包含多類，而每類中含有若干時序數據，因此在提取Shapelets時，需要逐行掃描候選矩陣，在每一條時序上均提取子序列，實現算法如下：

Fig. 4 The process of Shapelets extraction between T and R圖4 時序數據T和R的Shapelets提取過程

算法5. Shapelets提取算法.

輸入：候選矩陣M、子序列長度l、閾值th、時序數據集DataSet；

輸出：Shapelets集合SL.

① initializeSL；*初始化Shapelets集合*

② for eachrowinM

③ initializetempSL;

④ fori=1:length(row)

⑤ ifrow[i]>th

⑥subSequence=DataSet(row,i:i+l-1);*提取子序列*

⑦tempSL=[tempSL;[i,subSequence]];

⑧ end if

⑨ end for

⑩tSL=mergeSubsequence(tempSL);

算法5與圖4所描述過程本質上是一樣的，唯一的區別在于候選矩陣值的計算方式有所不同.假設時序數據集中包含3類：Ⅰ，Ⅱ，Ⅲ類，則多類情況下的Shapelets提取基本過程如圖5所示.圖5中第Ⅰ類所提取的Shapelets，可以將Ⅰ類時序數據從數據集中區分開來，但并不能將Ⅱ類和Ⅲ類的時序數據進行區分，因此還需要以Ⅱ類時序數據為基準從Ⅱ類中提取Shapelets，此時，Ⅰ類中的時序數據不再參與計算以減少時間消耗.

Fig. 5 The process of shapelets extraction圖5 Shapelets 提取過程

Fig. 6 The space transformation of time series圖 6 時序數據的空間變換示意圖

3.4 算法分析

從算法2～5描述過程可以看出，本文所提出的Shapelets提取算法以不同類別時序的差異為基本出發點，從時序數據中提取能夠區分不同類別的子序列是一種最直觀的方式.與現有Shapelets提取算法不同的是本文直接從子序列的角度描述時序數據間的相似性，其優勢在于易于發現時序間差異所對應的子序列.研究的核心在于如何以子序列為基準快速地計算不同類別時序間的差異.本文通過2個方面來加速計算：1)通過選取關鍵時序數據來減少后續參與計算的時序數量，通過實驗驗證(5.2節)該步驟是合理和有效的；2)借助傅里葉變換快速完成時序間的相似性連接向量和自相似性連接向量.本質上，本文是通過兩兩比較來構成Shapelets集合的，這樣得到的每一條子序列至少能區分來自不同類別的2條時序，由其構成的集合也能夠很好的將不同類別的時序區分開來，因此，算法是合理可行的.

4 時間序列空間變換

在以Shapelets為特征的時序數據分類算法中，決策樹是使用最為廣泛的分類器，其優點是容易構造，對分類結果有直觀的解釋，但缺點也很明顯.決策樹以Shapelets為節點，將數據集不斷地進行劃分，直至每個子數據集不能再分或者達到某些終止條件.這會導致2方面的問題：如果只使用少數Shapelets構建淺層二叉樹，其分類準確率受到限制，并且很多Shapelets具有非常相近的分類能力，不恰當的選擇也會導致分類準確率降低；而如果使用的Shapelets太多，樹的結構復雜，又會導致運行時間加大.Bagnall等人[30]指出將時序數據映射到特征易于被檢測的空間是解決時序數據分類問題的一個好方法，本文借鑒Lines等人[18]的策略，將時序數據映射到Shapelets所構成的空間，形成原始時序數據的特征表示.

圖6描述了基于距離運算的空間轉換的方法：

① http://www.uea.ac.uk/computing/machine-learning/shapelets/shapelet-data

假設提取的Shapelets數目為k，對于每一條時序數據Ti，依次計算其與k條Shapelets(SL1,SL2,…,SLk)之間的距離(式(2))，將這些距離組合起來形成1個距離向量，該距離向量即為原時序數據在Shapelets空間中的表示.這種表示一方面保持了Shapelets在時序數據分類問題中的優勢，因為向量中每一個值都是以Shapelets為單位得到的；另一方面，向量中的每個值都蘊含著時序數據中數據的先后次序關系，因此，向量中各元素之間具有相互獨立性，可以使用現有成熟的分類器(本文使用SVM)進行分類.

5 實驗與分析

本文通過對時序數據進行分類來驗證所提出的Shapelets提取方法的有效性.本節將首先介紹相關的實驗設置，包括實驗所使用的時序數據集、評價指標以及用于比較的Shapelets提取算法，之后給出實驗結果并進行分析和討論.

5.1 實驗設置

5.1.1 數據集

本文使用UCR archives和UEA時序數據集①中的26個數據集來評估方法的性能，這些數據集是目前Shapelets相關研究中普遍使用的數據集[16,25].數據來自于多個領域，包括傳感器數據、圖像輪廓信息、人體心電圖以及動作數據等，數據集的基本信息如表1所示:

Table 1 Datasets Used in the Experiments表1 實驗中所使用的數據集

Continued (Table 1)

從表1可以看出，所使用的數據集類型多樣：從二分類到多分類，且多分類問題較多，最多可達37類(Adiac)；長度不一，最短為24(ItalyPower.)，最長為512(Otoliths)；訓練集和測試集的數據差異也比較大，因此，可以全面的度量算法的性能.為了便于性能比較，本文采用默認的訓練集和測試集劃分，實驗中相關的參數均通過交叉驗證的方式獲得.

5.1.2 方法對比

本文算法(similarity join for Shapelets extraction, SJS)與目前主要的7個Shapelets提取算法進行性能比較：

1) 標準基于Shapelets時序數據分類算法使用不同的度量指標：信息增益(IG)[13]、Kruskal-Wallis檢驗(KW)[17]、F-統計量(FS)[17].

2) FSH(fast-Shapelets)[16]算法.將時序數據轉換為離散的低維表示，在低維空間中過濾掉大部分子序列以加速搜索的算法.

3) LTS[24]算法.使用梯度下降法從數據中學習出最優Shapelets和分類決策函數的算法.

4) IGSVM[18]算法.首次對時序數據進行空間變換的算法，以信息增益作為Shapelets區分能力度量，使用線性SVM作為分類器完成分類工作.

5) FLAG[25]算法.采用學習優化的策略，使用廣義特征向量的方法對時序數據進行降維，進而提取相應子序列作為Shapelets的算法.

上述對比算法的實現均由提出該方法的原作者提供，其參數為作者推薦的設置.

5.1.3 評價指標

Fig. 8 Running time comparison of SJS and N-SJS algorithm圖8 SJS算法和N-SJS算法的運行時間對比

對于分類問題，分類準確率是評價算法性能最重要的標準，但對于多種算法在多個數據集上的結果比較，分類準確率難以直觀地度量各個算法的性能，因為在很多情況下，某一算法在一數據集上效果很好，但在另一數據集上，分類結果可能較差.因此，本文采用在機器學習領域廣泛使用的無參數Friedman測試作為評價算法性能的標準.本質上，Friedman測試是基于分類準確率計算的.在獲得K個算法在N個數據集上的分類結果后，對每一個數據集上的K個算法進行排序，分類準確率最高的算法標記為“1”，次高標記為“2”，以此類推，分類效果最差的算法標記為“K”，這樣就可以得到N×K的排序矩陣H，其中，hi j表示第i個數據集上第j個算法的標記值.然后，計算每個算法的平均標記值

(6)

(7)

可以用具有K-1個自由度的卡方分布近似.Dem?ar[31]通過分析前人的工作得出卡方分布是相當保守的近似，提出使用

(8)

作為度量.當實驗結果拒絕零假設時，Dem?a提出將所有算法按照F值分到不同的組，使得同一組內的算法在分類準確率上沒有顯著差異.這樣，不同算法的性能差異可以通過包含平均次序和無明顯差異算法組的關鍵差異圖(critical difference diagram)表示.

5.2 提取關鍵時序數據對實驗結果的影響

本節首先驗證僅僅使用關鍵時序數據所提取的Shapelets對于分類結果的影響.圖7顯示了在表1數據集上使用算法1和不使用算法1的情況下，所提取的Shapelets用于分類的結果對比，其中，N-SJS表示使用全部時序數據時的算法，SJS算法表示僅使用關鍵時序數據的算法.橫、縱軸均以分類準確率為單位.

Fig. 7 Classification accuracy comparison of SJS and N-SJS algorithm for datasets in Table 1圖1 SJS算法和N-SJS算法在表1數據集上的分類準確率比較

從表1數據集上的分類結果(圖7)可以看出，在絕大多數數據集上，SJS和N-SJS幾乎沒有區別，即分類結果位于中間對角線附近，說明通過算法1精簡數據集后，并沒有減弱本文方法提取高質量Shapelets的能力.同時，在部分數據集上，準確率略有提升，意味著冗余的Shapelets也可能導致分類器性能的下降.

另一方面，算法1的主要目的是減少參與計算的時序數據的數量，加速Shapelets的提取，圖8是兩者的運行時間對比圖(縱軸為對數坐標).在絕大多數的數據集上，使用算法1之后，加速效果明顯，如在Adiac數據集上的運行時間比為23.6 s∶626.1 s，Chlorine.數據集上為7.7 s∶957.3 s.僅有4個數據集沒有加速效果，這是因為這些數據集的訓練集本身非常小，且數據之間差異比較大，算法1將所有時序都作為關鍵數據.但這并沒有影響本文算法的效率，因為算法在這些數據集上使用全部時序數據的時間消耗非常小(Coffee:3.5 s，ItalyPower:0.7 s，Sony:0.1 s，TwoLeadECG:0.4 s).

5.3 實驗結果對比

本節從分類準確率和運行時間2方面來驗證本文算法SJS的性能.

5.3.1 分類準確率

表2列出了本文算法與所有對比算法在各個數據集上的分類準確率.其中，每個數據集上的最優算法以加粗的形式表示，需要說明的是IG，KW，FS這3個方法在部分數據集上運行時間過長，無法在合理的運行時間內結束(超過24 h)，標記為“*”.表2的最后一行表示各個算法在不同數據集上具有最好分類效果的次數.從該統計結果上來看，LTS具有一定的優勢，在13個數據集上取得最好的分類準確率，本文算法SJS僅次于LTS，在9個數據集上處于領先，而標準的使用IG和KW作為度量標準的Shapelets算法效果最差.

然而，僅從算法在各個數據集上獲勝的次數來評價算法性能是不全面的，為此，本文使用5.1.3節的方法來評價各個算法，其結果如圖9所示.從該關鍵差異圖可以看出，本文算法與最優的LTS算法差別很小(2.54∶2.35)，但與其他算法有比較明顯的差別，如其他算法中最好的FLAG算法其平均次序已到3.15.盡管IGSVM算法也在7個數據集上取得最好的準確率，但在平均次序上，與本文算法差距明顯(2.54∶3.58).

Table 2 Testing Accuracy of Different Methods on the Datasets表2 不同方法在各個數據集上的測試準確率 %

Note: The results highlighted in bold denote that the method gets the highest accuracy for this dataset and the “*” represents the corresponding method cannot finish in 24 hours.

Fig. 9 Critical difference diagram for SJS and other baselines圖9 SJS和對比算法的關鍵差異圖

綜合上述2方面的分析，可以看出本文所提出的SJS算法具有很好的分類準確率，要優于目前絕大多數Shapelets提取算法.除此之外，SJS的最大優勢在于算法的穩定性，在幾乎全部的數據集上，SJS都保持很高準確率，而LTS算法則有些波動，如在Adiac和Otoliths數據集上，LTS準確率不足60%.更顯著的是，當類別數增多時，SJS優勢明顯，如在37類的Adiac數據集和10類的MedicalImages數據集上，SJS都要優于LTS.

5.3.2 運行時間比較

對于算法的評價，除了準確率之外，算法的運行時間也是一個重要的評價指標.對于實際應用，如金融、醫療等數據的處理，算法在可接受的時間內得到結果顯得尤為重要.本文算法的核心在于時序數據間的相似性連接，子序列與1條時序數據的距離計算借助快速傅里葉變換完成(算法2)，其算法復雜度為O(mlbm)(m為時序數據長度)，因此，計算2條時序數據相似性連接的時間復雜度為O((m-l+1)mlbm)，即近似于O(m2lbm).對于時序數據集，如果要處理所有時序數據對的相似性，算法復雜度會很高，但本文通過算法1提取每類中的關鍵數據，且相似性連接計算只在不同類的時序數據間進行，使得參與計算的時序條數大大減少，降低了運行時間消耗.

本文算法與其他對比算法的運行時間如表3所示(所有結果均在使用Intel i7 CPU和16 GB內存的計算機上得到)，括號中的數字表示算法運行時間在每一個數據集上的排序次序.同表2一樣，算法如果無法在合理的時間內完成，標記為“*”，且次序用“(8)”表示.為了方便對比，本文按照式(6)計算每一算法的平均次序，表3中的最后一行Average Rank給出了相應的結果.該統計結果表明:FLAG算法具有最好的時間效率(平均次序為1.2)，SJS次之(平均次序為2)，標準使用IG作為度量的Shapelets提取算法的時間消耗最多.

Table 3 Running Time of Different Methods on the Datasets表3 不同方法在各個數據集上的運行時間

Continued (Table 3)

Note: The results highlighted in bold denote that the method gets the highest accuracy for this dataset and the “*” represents the corresponding method cannot finish in 24 hours.

雖然運行時間的平均次序能夠體現出算法的時間效率，但沒有體現具體運行時間差異的大小.從表3中可以看出，盡管FLAG要優于SJS，但在大部分數據集上，兩者幾乎在同一個數量級范圍內，均在幾秒甚至更短的時間內得出結果.相反，其他算法則需要很長的時間，特別是準確率最高的LTS算法，稍大一點的數據集都需要幾千甚至上萬秒的時間消耗.

本節的實驗表明：LTS具有最好的分類效果，但運行時間長，尤其在類別數多時，時間消耗大，如在Adiac數據集上，差不多需要24 h的時間；FLAG的時間效率最好，但分類準確率不高，而本文算法SJS不僅在分類準確率上接近LTS算法，而且在運行時間上與最優的FLAG具有同等的數量級.因此，綜合分類準確率和運行時間2方面，SJS的實用性更好.

5.4 算法可擴展性

本文以UCR archives中的最大時序數據集之一StarLightCurves來驗證算法的可擴展性.該數據集包含3類的“星光”數據，共計9 236條，其中Eclipsed Binaries有2 580條，Cepheids有1 329條，RR Lyrae Variables有5 236條.實驗中采用默認的訓練集和數據集劃分.Shapelets的提取時間主要受2方面的影響：1)時序的長度;2)訓練集的大小.本文從這兩方面進行實驗.由于對比算法IG，KW，FS，IGSVM在StarLightCurves數據集上的運行時間太長，難以獲得實驗結果，因此，本文舍棄與這些算法的比較.

Fig.10 Experimental results for varying the length of time series圖10 時序數據長度變化時的實驗結果

首先驗證時序長度的變化對于分類準確率和運行時間的影響，具體設置為：固定訓練集中時序的數量，將時序的長度從100變化至1 024，以100為間隔單元.實驗中，LTS算法受限于時間和空間消耗，在長度大于300時難以得到結果，因而只有3個實驗數據.實驗結果如圖10所示，圖10(a)是隨著時序長度的增加，各對比算法在測試集上分類準確率的變化，圖10(b)是對應運行時間的變化.從圖10(a)中可以看出，本文算法SJS受時序數據長度變化影響最小且相對穩定，而其他算法在長度為100和200時，分類準確率很低且在長度增大過程中分類結果有波動，這說明，SJS在提取Shapelets過程中能夠很好的發現當前處理的數據集中最具類別可分性的子序列.運行時間對比的結果與表3基本一致，即FLAG具有最優的運行時間，SJS次之，但長度在100～600之間時，SJS要優于FLAG，這充分驗證了SJS的主要時間消耗來自時序數據間的相似性連接.

Fig.11 Experimental results for varying the number of training time series圖11 訓練集時序數據數量變化時的實驗結果

其次，固定時序數據的長度為1 024，將訓練集中時序數據數量從100增加至1 000，使用默認的測試集.由于LTS算法無法在合理的運行時間內得到結果，所以本實驗中僅有SJS，FSH，FLAG這3個算法的結果，如圖11所示.與FSH和FLAG在訓練集增大時，分類準確率有起伏不同，SJS一直處于上升的趨勢，如圖11(a)，說明當訓練集增大時，FSH和FLAG所提取的Shapelets可能會不同，而SJS則會在原有提取出的Shapelets基礎上，對新增加的訓練集數據提取Shapelets，進而提高分類準確率.這一特性對于優化算法的運行效率更為重要，如圖11(b)所示，盡管在同等規模的數據集上FLAG要優于SJS(如圖11(b)所示)，但SJS在應對數據集變化時可以重用之前結果，能夠減少后續計算量，如圖11(b)中的SJS-R曲線，因此，SJS有更好的增量擴展性.

6 結 論

本文針對時序數據中的Shapelets提取問題進行研究，提出了基于相似性連接策略的Shapelets提取算法SJS.該方法以子序列為基本單元描述時序數據之間的相似性，通過不同類別的時序數據之間的差異向量構建候選矩陣，然后提取候選矩陣中大差異值對應的子序列構成Shapelets.在大量真實時序數據集上進行分類的實驗表明，本文方法所提取的Shapelets不僅能獲得良好的分類準確率，還具有很高的時間效率，實用性強.同時，可擴展性實驗表明:在訓練數據集增大時，本文算法能夠增量提取Shapelets，可以減少大量重復計算.