基于子序列相似性的時間序列語義挖掘算法

陸怡，王鵬，汪衛

（1.復旦大學 軟件學院，上海 201203；2.復旦大學 計算機科學技術學院，上海 201203）

0 概述

時間序列數據是在一段時間內以固定的時間間隔采集的數據點序列，用于描述現象隨時間變化的情況，已成為日常生活中重要的信息記錄形式。隨著大數據時代的到來以及大規模計算能力的提升，時間序列數據分析與挖掘已成為研究熱點，廣泛應用于醫療、交通、金融等領域。在現實場景中，時序數據雖然存在不同的演化規律，例如記錄服務器CPU 使用情況的監控數據、反映患者健康狀況的心電數據以及表征市場行情的商品銷售數據，但這些領域各異的數據的內核通常是一致的，會根據觀測對象或系統的實際物理狀態，呈現出不同的波動形式。因此，挖掘時間序列中潛在的語義信息，實際上是識別被監測系統的物理狀態，通過將時間序列轉換為狀態序列，實現數據壓縮或者異常檢測。

時間序列具有數據量大、復雜度高、干擾信息多等特性，一般采用特征提取手段對時間序列進行加工處理。針對時間序列的語義挖掘，主要分為基于領域特征、基于全局特征、基于子序列相似性三類。基于領域特征的時間序列語義挖掘算法除了時間序列數據之外，還引入了領域知識或者帶標簽的數據。文獻［1］針對金融時間序列的自相關特征，提出一種基于ARMA 模型［2］的時間序列分割算法。文獻［3］根據心電時間序列的特點，提出一種基于殘差平衡及邊界約束的分段線性回歸方法以得到時間序列的分段表示。文獻［4］通過從訓練集中構建特征的高斯概率分布模型，實現有監督的狀態挖掘。基于全局特征的時間序列語義挖掘算法使用全局特征對時間序列進行概括。全局特征又可細分為斜率和均值兩種類型。一種以pHMM［5］為代表，從斜率的角度將時間序列簡化成線段序列，并利用隱馬爾科夫模型（Hidden Markov Model，HMM）［6］指導線段的聚類，最終獲得狀態信息。另一種以AutoPlait［7］和StreamScope［8］為代表，從均值的角度描述不同的狀態，得到多層的隱馬爾科夫模型，并通過最小描述長度（Minimum Description Length，MDL）［9］對可選的模型進行評估，從而得到時間序列的最優表示。基于均值的算法還可進一步應用于多維或高階時間序列，例如DBSE［10］以及在AutoPlait 基礎上提出的CUBEMARKER［11］。在基于子序列相似性的時間序列語義挖掘算法中，具有代表性的子序列會頻繁出現。實際上，大量對于時間序列的分析是從子序列入手的，例如文獻［12-13］實現的基序列挖掘，文獻［14］提出的基于子序列的時間序列分類方法。但是，基于子序列相似性的語義挖掘算法仍在初步研究階段。FLUSS［15-16］和ESPRESSO［17］作為典型代表，僅能識別出狀態的轉變點。

目前，多數時間序列語義挖掘算法對于時序數據特征有一定的約束條件，因此難以處理海量且形式各異的時序數據。本文提出一種基于子序列相似性的時間序列語義挖掘算法SEDIS，通過計算子序列的相似性，將時間序列分割成片段序列進行兩級聚類，識別出時間序列中潛在的物理狀態，使其具備在不同應用場景下的通用處理能力。

1 問題描述

定義1（時間序列）時間序列是對某個事物或系統進行連續同間隔的測量得到的數值序列，可表示為T=(t1，t2，…，tL)，其中，L為時間序列長度。

定義2（子序列）子序列是由時間序列中一段連續的值組成的數值序列。Ti，l=(ti，ti+1，…，ti+l-1)表示從第i個時間點開始，長度為l的子序列。

時間序列作為事物或系統某個剖面的觀測結果，按照時間序列反映的物理意義將其片段組成若干分段，并根據分段的內在聯系形成狀態序列，稱為分段表示與狀態表示。

定義3（分段表示）時間序列T的分段表示可定義為SEG(T)={s1，s2，…，sn}，其中si(1 ≤i≤n)是T的一個子序列。

定義4（狀態表示）時間序列T的狀態表示在分段表示的基礎上，引入了分段的物理狀態，可定義為STATE(T)={＜s1，c1＞，＜s2，c2＞，…，＜sn，cn＞}，其中ci表示分段si的狀態。假設m表示物理狀態個數，則ci是一個介于1和m之間的正整數，且1 ≤m≤n。

根據上文描述，SEDIS 的目標任務是尋找一個符合用戶期望的狀態表示STATE(T)以描述或壓縮原時間序列T，其中包含的物理狀態個數m由用戶定義。實際上，如果將時間序列中的每個時間點視作一個對象，該任務也可理解為聚類任務，目標是將反映同一狀態的觀測值聚集在同一類簇中。

2 時間序列語義挖掘算法

SEDIS 基于以下基本假設：在同一狀態下，由于具有代表性的子序列會頻繁出現，因此屬于同一狀態的分段之間必然存在多個相似的子序列。SEDIS分為兩個階段：1）基于子序列相似性結合基于密度的聚類方法構建組成分段的片段；2）采用貪心策略對片段進行預聚類，再通過k-means 聚類識別出片段對應的物理狀態。

2.1 片段分割

為避免兩兩計算子序列相似性帶來的巨大計算開銷，SEDIS 采用基于概率的抽樣方法，每輪選取一個子序列作為參考子序列。如果兩個子序列相似，則它們與參考子序列的距離相近。

2.1.1 子序列相似性的優化計算

子序列相似性通過標準化歐式距離進行度量。在歐式距離的基礎上，預先對每個子序列進行標準化操作，使得子序列各個時間點對應的數值的均值等于0、標準差等于1，從而消除了幅度縮放、基線漂移等對波形相似性的影響。具體而言，兩個長度等于l的子序列X=(x1，x2，…，xl)和Y=(y1，y2，…，yl)的標準化歐式距離可通過式（1）計算：

其中：μ表示均值；σ表示標準差。

假設時間序列的長度是L，參考子序列的長度是l，需要計算參考子序列與其余L-l+1 個長度同樣等于l的子序列的標準化歐式距離。由于時間序列的數據量通常較大，因此SEDIS 采用文獻［18］提出的基于快速傅里葉變換［19］的算法，加速子序列的相似性計算。首先對時間序列和參考子序列使用快速傅里葉變換等操作得到兩者的滑動點積，然后以滑窗的形式增量地記錄子序列的平均值和標準差，從而基于O（LlogaL）的時間復雜度求出L-l+1 個子序列對的標準化歐式距離。

由于相似是一個主觀的評判標準，難以通過量化后的數值完成二分類，因此為得到多個相似子序列以支持分段，通過最大堆篩選出與參考子序列距離最近的k個子序列完成初過濾。

2.1.2 基于密度的優化聚類

由于參數k的取值將極大影響k個子序列的相似情況，但時間序列的平滑特性以及全自動的應用要求，難以通過類似手肘法的方案提供一個k的建議值，因此SEDIS 采用基于密度的聚類方法［20］對k個子序列進行再過濾，其核心思想為屬于同一分段的子序列在時間軸上會緊密地聚集在一起。基于密度的聚類方法能夠自動地將鄰近的子序列分到一個類中，進而識別出潛在噪聲。具體地，聚類對象是k個子序列，聚類相似性標準取決于子序列的起始時間點。子序列對Ti，l和Tj，l的距離可簡單表示為|i-j|，即兩個子序列在時間軸上越接近，越有可能聚成一類。

針對該距離度量，本文提出一種基于密度的優化聚類方法。在確定每個子序列Ti，l是否為核心點時，需要檢索其EEps鄰域內存在的子序列數量。由于相似性僅考慮起始時間點，因此問題可簡化成找到滿足{Tj，l|i-EEps≤j≤i+EEps}的子序列集合。假設子序列已經按照起始時間點升序排序，則可以通過兩次二分搜索快速定位第一個滿足j≥i-EEps的子序列和最后一個滿足j≤i+EEps的子序列。兩者之間包含的子序列數量即Ti，l的EEps鄰域內存在的子序列數量。通過二分搜索，該步驟的時間復雜度可降至對數級。

經過基于密度的聚類方法過濾掉一部分的噪聲子序列后，剩余的子序列間可能存在重疊。為便于后續處理，在每一個形成的類簇中，找出其中起始時間點最小的子序列Ti，l和起始時間點最大的子序列Tj，l，以此構建一個新的子序列Ti，j+l-i，該子序列即為一個候選分段。

2.1.3 基于概率的迭代模式

由單個參考子序列得到的候選分段僅能表征一個狀態，因此需要引入多個不同的參考子序列使其涵蓋所有的物理狀態。本文提出一種基于概率的迭代模式，根據候選分段的情況動態地調整子序列被選為參考子序列的概率，設計思想為對于沒有被包含在候選分段中的子序列給予更多的機會。具體地，對于所有的子序列，賦予相同的初始概率，以模擬均勻分布。在每一輪完成后，將包含在候選分段中的子序列的概率減半，從而直接影響下一輪參考子序列的選取情況。

至此，SEDIS 產生若干候選分段，一組候選分段對應一個參考子序列，多組候選分段間可能存在重疊。候選分段的起始點和終止點表示候選的狀態轉變點。根據這些候選點對所有候選分段進行分割，產生的子序列稱為片段。一個分段至少產生一個片段。值得注意的是，片段不一定足以覆蓋整個時間序列，對于漏檢片段的處理，將在下文進行具體介紹。

2.2 狀態識別

狀態識別是根據片段反映的物理狀態對其進行聚類，同一類簇的片段對應同一狀態。因此，狀態識別的核心為片段在物理狀態層面上的相似性定義。

2.2.1 基于貪心策略的預聚類

在聚類前，時間序列通常是平滑演進的。排除監測設備的突發故障，時間序列在狀態轉變點附近的波動通常較小，因此會導致附近的子序列被分割成多個片段。為提高效率，本文對所有生成的片段采用貪心策略進行預聚類，使用棧模擬貪心策略的執行，具體步驟如下：

1）將所有片段按照起始點的降序依次推入棧中。

2）從棧中推出棧頂片段Ti，p，找出所有包含該片段的候選分段C={Tj，q|j≤i＜i+p≤j+q}。

3）計算集合C中所有候選分段的終止點的平均值A=，其中|C|表示集合C的元素個數。

4）從C中選擇終止點最接近A的候選分段。將Ti，p的起始點i和選中的候選分段的終止點j+q-1 組合形成新的片段Ti，j+q-i。

5）從棧中推出新的棧頂片段，如果被Ti，j+q-i包含，則略過；否則，重復步驟2～步驟5 直到棧中無剩余元素。

基于上述步驟，在保留片段原有語義信息的同時，通過片段間的合并減少了片段數量，提高了算法運行效率。

2.2.2 基于k-means 的再聚類

根據聚類的任務描述：將沒有分類標簽的數據集分為若干個簇［21］，再結合狀態識別目標，假設給定一個在物理狀態層面上的相似性度量，通過聚類可將片段組織成多個類簇，每個類簇對應一個狀態。

SEDIS 的基本假設為屬于同一狀態的片段之間存在多個相似的子序列，因此基于片段分割的過程，提出一種新的相似性度量函數。為便于描述，給定一個片段Pi，輔助向量Oi是一個r維的向量，其中r是參考子序列的總個數。Oi的每一維是0 或者1，表示Pi是否被由該參考子序列生成的候選分段包含，進而定義兩個片段Pi和Pj的相似性，如式（2）所示：

其中：&表示按位與；sum 表示求和；min 表示求最小值。

式（2）計算的是兩個片段在同一候選分段中的共現頻率，在一定程度上反映了兩個片段包含的相似子序列的數量。該值越大，表明兩者屬于同一狀態的可信度越大。

基于相似性定義，采用k-means 聚類［22］算法得到片段的隱含結構，即片段的物理狀態，具體步驟如下：

1）隨機選取m個片段作為初始的聚類中心，其中m是用戶指定的物理狀態個數。

2）將每個片段分配到相似性最高的聚類中心所在的類簇中。

3）對于每個類簇，選取其中與所有其他片段相似性之和最大的片段作為新的聚類中心。

4）重復步驟2 和步驟3，直到聚類中心不再變化。

至此，每個類簇對應一個物理狀態。

2.2.3 基于k 最近鄰的漏檢片段分類

需要指出的是，雖然SEDIS 在片段分割階段通過基于概率的迭代模式盡可能地保證參考子序列涵蓋全部的物理狀態，但是片段的漏檢現象仍可能存在，即在時間序列中，有部分子序列并未包含在任何候選分段中。對于這部分子序列，或稱為漏檢片段，SEDIS 采用k 最近鄰（k Nearest Neighbor，kNN）算法進行分類處理。在實際應用中，將每個漏檢片段作為參考子序列，考察在標準化歐式距離度量下與其最相似的k個子序列的物理狀態，取頻次最高的作為分類標簽。

2.3 時間復雜度分析

SEDIS 主要包括片段分割和狀態識別兩個階段。理論上，第一階段的時間復雜度與時間序列的長度強相關，第二階段的時間復雜度與片段的數量強相關，后者遠小于前者，同時在實際應用中證明，第一階段的時間開銷在總時間開銷中占比較大，因此本節將主要圍繞該階段展開分析。

片段分割包括對每個參考子序列執行相似性計算、相似子序列篩選、聚類3 個步驟。假設時間序列的長度等于n，參考子序列的個數等于r，每次只保留k個與參考子序列最相似的子序列，則該階段的時間復雜度為O（r×（nlogan+nlogak+klogak））。但需要指出的是，由執行快速傅里葉變換引入的O（nlogan）的時間復雜度在實際應用中是近似線性的，文獻［18］認為這可能是由于該算法已在工程領域中得到高度優化。此外，r和k作為輸入參數，與n不存在線性增長關系，這意味著對于長時間序列，r和k也可視作常量。因此，SEDIS 的總時間復雜度是近似線性的。

3 實驗與結果分析

實驗主要分為3 個部分：1）驗證SEDIS 所得的物理狀態的可解釋性；2）將SEDIS 在時間序列分割和狀態識別兩階段的準確率和運行效率與其他算法進行對比；3）驗證SEDIS 對于相關參數的魯棒性。

3.1 實驗設置

SEDIS 使用Matlab 實現。實驗運行于配置Intel Core i5 處理器、1.4 GHz主頻、帶 有8 GB 內存的MAC 筆記本電腦上。

實驗數據集主要細分為以下2 類：1）Barbet、Fetal、SP02 和ECG 數據集，這類數據集包括狀態變化更少且長度更短的時間序列，選自文獻［18］中的32 個數據集，其中涵蓋醫療、生物、工業等領域相關數據，以證明SEDIS 的通用性；2）PAMAP 數據集，這類數據集包括狀態更多且長度更長的時間序列，選自PAMAP 運動數據集［23-24］，是傳感器采集的監測數據，反映傳感器佩戴者進行的有氧運動種類，如慢走、跑步、踢足球等。5 個數據集的具體信息如表1所示。

表1 數據集信息統計Table 1 Data set information statistics

SEDIS 主要涉及參考子序列長度l、參考子序列個數r以及與參考子序列最相似的子序列個數k這3 個參數。由于SEDIS 基于子序列相似性，因此l的推薦值是代表性子序列的長度。例如，對于心電圖（Electrocardiogram，ECG）序列，l可設置成一次心跳的時長。借助該推論提出一種自動化設置參考子序列長度的方法。具體地，利用快速傅里葉變換將時間序列從時域映射到頻域，此時能量最大的頻率f即主導頻率，頻率的倒數為代表性子序列的周期l=。根據此規則，對于5 個數據集，l分別取81、20、24、58、13，r統一取100，k分別取n/600、n/100、n/200、n/100、n/100，其中n表示時間序列長度。

實驗選取基于子序列相似性的FLUSS 算法［15］（僅支持時間序列的分割問題）、基于斜率的pHMM算法［5］、基于均值和標準差的AutoPlait算法［7］與SEDIS 進行性能對比。FLUSS 基于子序列相似性，需要輸入子序列長度，為保證一致性將其設置成SEDIS 使用的參考子序列的長度l。pHMM 通過εr和εc兩個參數限制使用線段擬合原始序列以及構建線段聚類結果產生的誤差，在實驗中通過調整這兩個參數來獲得最佳識別結果。AutoPlait 是全自動算法，不涉及參數設置。

3.2 算法可解釋性驗證

在SP02 數據集上對SEDIS 所得的物理狀態進行可解釋性驗證。SP02 數據集由傳感器采集的人體血氧飽和度數據組成，血壓在一定程度上會對該數據產生影響。為便于展示，對該時間序列進行降采樣處理，將長度縮減至350。由圖1（a）可以看出，時間序列中存在2 種不同的狀態，兩者對應的代表性子序列的長度和波形均不一致，2 種狀態分別表示傳感器佩戴者的血壓維持正常和血壓突然驟降。對比圖1（b）的原始標簽以及圖1（c）由SEDIS 得到的識別結果可知，SEDIS 識別出的物理狀態與真實狀態基本吻合，證明了通過子序列相似性區分狀態是可行的。pHMM 和AutoPlait的識別結果如圖1（d）和圖1（e）所示。由于FLUSS 僅支持序列分割，因此不在此進行展示。從選擇的特征進行分析，pHMM 基于斜率，將上升段和下降段分為2 種狀態，最終識別出3 種狀態，但無法從更宏觀的角度描述狀態，而AutoPlait基于均值和標準差，由于該案例中的2 種狀態在均值和標準差方面不具備區分性，因此AutoPlait 將整個時間序列視為同一狀態，識別效果不理想。

圖1 算法可解釋性驗證結果Fig.1 Verification results of algorithm interpretability

3.3 算法準確率分析

由于語義挖掘算法在識別出狀態的同時也監測出不同狀態間的轉變點，因此本節將從序列分割和狀態識別兩個維度對算法準確率進行量化評估。

3.3.1 序列分割準確率

FLUSS 僅能進行序列分割，而其他對比算法忽略了每個分段的狀態標識，僅考慮分割點的準確率。本文提出一種新的分割誤差指標sse，包括精確率指標和召回率指標，其中為所有預測分割點與其最接近的真值點的距離和，為所有不包含預測分割點的區域的長度和。將整個時間序列劃分成多個區域，區域的分界點由兩個相鄰真值點的中值決 定，越接近于0，表示分割點越準確。

表2 給出了4 種算法在5 個數據集上的分割誤差，其中，×表示該算法未識別出任何有效分割點，意味著將整個時間序列視作同一狀態。由表2可以看出，SEDIS在5 個數據集上均具有優異表現。FLUSS 同樣基于子序列相似性，較優于pHMM 和AutoPlait。但由于其難以處理狀態多次出現的情況，因此在部分數據集上表現不佳。pHMM 更適合區分斜率不同的波段而非復雜波形，因此在結果中出現了大量的分段碎片，導致召回率指標分值較低。AutoPlait難以區分均值和標準差較為一致的狀態，僅在兩個數據集上找到有效的分割點。由此可見，SEDIS 具備更強的通用性。但是，SEDIS 未能完全命中真值點的原因在于：時間序列是平滑演進的，分割點附近的子序列相對接近，導致分割點難以精準定位，然而相比于分割點的誤判和漏判，準確率方面的誤差是可接受的。

表2 4 種算法的分割誤差比較Table 2 Comparison of segmentation errors of four algorithms

3.3.2 狀態識別準確率

狀態識別問題實際上是一個特殊形式的聚類問題，假如將時間序列中的每個時間點視作一個聚類對象，那么目標就是將所有時間點劃分成多個不同的類簇，每個類簇反映一個物理狀態。本文引入調整蘭德系數（Adjusted Rnd Index，ARI）［25］作為狀態識別準確率的衡量標準。ARI的取值范圍在-1 到1 之間。ARI越大，表明聚類結果與實際分類情況越吻合。

表3給出了3種算法的狀態準確率比較結果，其中×表示該算法未在該數據集上進行狀態識別，即將整個時間序列視作一個狀態。由表3可以看出，SEDIS在5個數據集上均達到了90%以上的準確率，而pHMM 和AutoPlait 則表現一般，符合上文的理論分析。

表3 3 種算法的狀態識別準確率比較Table 3 Comparison of state recognition accuracy of three algorithms %

3.4 算法效率分析

基于PAMAP 數據集進行算法效率實驗，比較不同算法在不同時間序列長度下的運行時間，如圖2 所示。由圖2 可以看出，SEDIS 的運行效率是近似線性的，優于其他算法。當時間序列長度n等于694 380時，SEDIS的運行時間僅為112.73 s。

圖2 4 種算法的運行時間比較Fig.2 Comparison of the running time of four algorithms

3.5 算法參數設置對狀態識別結果的影響

通過實驗分析參考子序列長度l、參考子序列個數r以及與參考子序列最相似的子序列個數k這3 個參數在ECG 數據集上對于SEDIS 狀態識別結果的影響程度。假設將子序列長度表示為l′，實驗中固定其他參數，將子序列長度從0.5×l′增長至4×l′。圖3（a）給出了不同子序列長度對于狀態識別準確率的影響，可以看出SEDIS 對于子序列長度的依賴性較弱。圖3（b）給出了不同參考子序列個數對于狀態識別準確率的影響，可以看出由于SEDIS 中基于概率的迭代模式在參考子序列有限的情況下，仍舊具有較為穩定的識別能力。圖3（c）給出了與參考子序列最相似的子序列個數，實驗從50 個子序列增長至400 個子序列，可以看出盡管最相似的子序列個數相比其他兩個參數對于結果造成的影響更大，但是SEDIS 在根據這個指標篩選最相似子序列之后，還會通過基于密度的聚類算法對結果進行再過濾，因此子序列個數的決定性作用被弱化了，使得整體準確率仍保持在98%以上。

圖3 參數設置對狀態識別結果的影響Fig.3 Influence of parameter settings on state recognition results

4 結束語

針對現有時間序列語義挖掘算法缺乏通用性的問題，本文提出一種基于子序列相似性的時間序列語義挖掘算法。通過定義在物理狀態層面上的子序列相似性，并結合兩級聚類識別時間序列中的潛在狀態。引入基于概率的迭代模式，提高算法的識別準確率和運行效率。在真實數據集上的實驗結果證明了該算法的有效性和可解釋性，并且表明其具有較強的魯棒性。下一步將針對多維時間序列和流式數據進行擴展，以解決海量數據的實時分析問題，并從語義挖掘的角度出發，將所得狀態作為有監督學習樣本，進一步執行時間序列的異常檢測、預測和分類等流程，實現數據的高效利用。