基于分段降維和路徑修正DTW的時序特征分類器設計

常炳國，臧虹穎

(湖南大學 信息科學與工程學院，長沙 410082)(*通信作者電子郵箱657475865@qq.com)

0 引言

當前，時間序列的分類和聚類研究正受到越來越多的關注，在圖像識別、信號處理、生物信息識別以及金融等領域得到廣泛應用。相似性度量作為時間序列挖掘工作的基礎步驟，其運算效率以及準確率直接影響時間序列挖掘的最終效果[1]。動態時間彎曲(Dynamic Time Warping， DTW)方法的提出最初是為了解決語音序列長短不一的模板匹配問題，逐漸被應用到不等長時間序列的相似性度量當中[2]。DTW度量允許數據點之間“一對多”的映射，通過動態規劃查找數據點之間的最佳匹配路徑，以實現數據在時間軸上的伸縮變化。由于其概念簡單、準確率高、魯棒性強等優點，已成為相似性度量中最常用的距離度量方法之一。針對DTW計算復雜度高、算法效率低等局限，Mei等[3]提出了基于Mahalanobis距離的DTW度量方法，該方法先利用馬氏距離建立了每個變量和類別之間的準確對應關系；然后再通過DTW對齊時間上不同步的序列。Sharabiani等[4]提出新的逼近方法來減小DTW輸入序列的長度，以此提高DTW搜索效率，他將該種降維方法命名為控制圖近似法(Control Chart Approximation， CCA)。Sun等[5]提出了一種全局約束度下的修剪動態規劃方法，從理論上證明了帶有全局約束的DTW匹配路徑對于度量修剪策略的有效性，并證明用它能得到序列之間近似最優解。李正欣等[6]在Keogh構造的下界距離，記為LB_Keogh[7]的基礎上，提出了一種支持DTW距離的多元時間序列的索引結構，并設計了早停機制以減少計算代價。李海林等[8]結合數值導數構造了新的特征序列，并設計了符合該特征序列的度量函數。綜合來看，現階段的研究主要從降低時間序列維度、制定全局約束條件以及設計新的下界函數三個方面展開[9]，以改善動態時間彎曲度量性能。

針對DTW方法僅關注路徑累積距離最短而出現過度彎曲現象，無法準確選擇最優彎曲路徑，進而影響時間序列分類和聚類準確率的問題，本文提出一種基于路徑修正的動態時間彎曲(Updated Dynamic Time Warping， UDTW)距離度量方法，通過給彎曲路徑設置懲罰系數，實現了對DTW彎曲變化率的動態調整，有利于從多條序列中選出形態相似的彎曲路徑。在進行UDTW距離度量之前，利用改進過的分段聚集近似(Piecewise Aggregate Approximation， PAA)方法提取原始時間序列的特征，以降低DTW計算代價，從而整體上提高時間序列相似性度量的效率和準確率。

1 時序特征提取

假定一條時間序列S的長度為m：S={s1,s2,…,si,…,sm}，另一條時間序列T的長度為n：T={t1,t2,…,tj,…,tn}。創建一個m×n的距離矩陣D，用d(ith,jth)表示(si,tj)兩點的路徑距離：d(i,j)=‖si-tj‖p，其中‖x‖p表示p-范數。計算DTW距離時，先構造一個累積距離矩陣R，R中對應元素γ(ith,jth)定義為：

γ(i,j)=d(i,j)+min {γ(i-1,j),

γ(i-1,j-1),γ(i,j-1)}

(1)

計算時，按照行(或列)的方向依次計算γ(i,j)的值，最后求得γ(m,n)即為序列S、T的DTW距離。由于對R中每一個元素進行了計算，DTW度量方法的時間復雜度為O(m×n)。由此可見，如果能夠縮短時間序列S、T的長度，便能大幅度減少DTW的計算代價。

常用的特征提取的方法有分段線性近似(Piecewise Linear Approximation， PLA)[10]，利用最小二乘法求得分段序列的最佳線性擬合曲線，該方法能較好地還原原始序列的形態，但對于長度為m，分段數目為L的序列來說，時間復雜度高達O(m2L)。類似的方法還有分段多項式表示(Piecewise Polynomial Representation， PPR)、分段回歸近似(Piecewise Regression Approximation， PRA)和自適應分段常量近似(Adaptive Piecewise Constant Approximation， APCA)等，雖然這些方法都對序列作了分段處理以達到降維目的，但這些方法本身時間復雜度較高，算法效率較低[11]。

分段聚合近似(PAA)用等長度窗口分割時間序列，每個窗口內序列特征用窗口的平均值來表示。

定義1 PAA。將時間序列S={s1,s2,…,si,…,sm}分為L段的PAA模型表示為：

SPAA={AVG(S1′),AVG(S2′),…,AVG(Sl′)};

l=1,2,…,L

(2)

其中AVG(Sl′)表示第l個窗口內數據子集Sl′的均值，窗口長度w=m/L(或m/L+1)。通過將每段的均值數據相連接形成新序列SPAA，從而實現數據降維的目的。PAA具有概念簡單、參數少(僅只有一個參數L)、時間復雜度低等優點，是一種有效的特征提取方法，但它僅反映了序列整體的變化趨勢，而忽略了時間序列的形態特征，并且，PAA對均值平穩的獨立噪聲數據不敏感，存在重要數據點或異常數據點信息丟失的不足[12]。綜合考慮PAA的優缺點，本文采用窗口最大值代替平均值的數據提取策略，然后用分段平均值對最大值進行適度平滑處理以消除噪聲影響，本文將此種方法稱為分段局部最大值平滑法(Piecewise Local Max-smoothing， PLM)。

定義2 PLM。將時間序列S={s1,s2,…,si,…,sm}分為L段的PLM模型表示為：

l=1,2,…,L

(3)

用αl=MAX(Sl′)-AVG(Sl′)表示第l個窗口最大值與平均值的差值，β為平滑常數，用于調整數據的變化范圍，由于前期對序列作了Z-score規范化處理，β取1即可。圖1給出了長度為48的序列S分別經過PLM和PAA方法提取特征值后所形成的特征序列SPLM、SPAA，此處窗口長度w=6。

圖1 時間序列S及其特征表示舉例

PLM方法只需一次遍歷序列即可完成特征提取，時間復雜度為O(m)，與PAA相比，PLM能更好地提取數據的形態特征，保留了主要的趨勢轉折點數據，同時適度的平滑處理消除了部分異常噪聲數據的影響，增強了對異常數據的容錯能力。

2 修正算法

與歐氏距離(Euclidean)度量方法[13]相比，DTW方法解決了時間點不對齊、形態之間伸縮、擴展的問題，但是DTW在選擇最優路徑時，為了達到累積距離最短，可能會出現多個點對應到同一點的現象，最終選擇出一條在垂直方向或水平方向過度復制的彎曲路徑。圖2(a)描述了序列A與序列B經過DTW度量后產生的最短彎曲路徑距離為5.09；圖2(b)給出了序列A與序列B對應的點對點匹配關系；圖2(c)描述了序列A與序列C經過DTW度量后產生的最短彎曲路徑距離為5.12；圖2(d)給出了序列A與序列C對應的點對點匹配關系。在進行相似性度量時，由于序列B與序列A的DTW距離更短，往往會把序列B匹配到序列A所在類別，盡管從形態上看，序列C的變化趨勢更符合序列A的特征。在圖像檢索、模式識別等領域內，序列模式形態的相似性是匹配成功與否的一個重要標準，而DTW度量方法僅僅只基于彎曲距離最短這一限制，極有可能錯失正確的分類。

2.1 懲罰函數

兩個時間序列之間最理想的對應關系應該是相同特征之間的一一對應，比如一條時間序列上的波峰應與另一條序列的波峰對應，而不是波谷。DTW度量在選擇最優彎曲路徑時，僅僅只考慮了累積距離最短，沒有對其彎曲的步數進行一定的限制，從而忽略了形態相似性的重要性，因此，給過度彎曲的路徑設置一定的懲罰系數是有必要的，對于在路徑的垂直和水平方向連續轉移的步數，要以一定的增長率逐漸增加其懲罰系數，達到約束其彎曲路徑的目的。理想的懲罰函數q(x)應滿足以下三個條件：

1)有界性。q(0)=0，qmax為常數，即最小懲罰系數為0(不懲罰)，最大的懲罰系數可控。

2)單調性。n-m≥0 ?q(n)-q(m)≥0。隨著連續彎曲步數的增大，懲罰系數也要逐漸增大。

圖2 彎曲路徑與點對點匹配結果

為了滿足以上3個條件，本文定義懲罰函數為：

q(x)=-qmax(cos(πx/2μ)-1);x∈(0,μ)

(4)

其中：qmax為懲罰函數的上界值，可以根據需要進行調整；μ為達到上界值對應的最大彎曲步數，也是算法允許在同一方向連續轉移的最大彎曲步數，q(μ)=qmax。圖3給出了上界值qmax=10時不同μ值對應懲罰函數曲線，μ值越大，曲線越平緩，則度量時允許連續彎曲的步數越多。最大彎曲步數μ是一個經驗值，通過設定μ值的大小來控制懲罰函數的傾斜度，實現對懲罰系數的動態調整。

圖3 不同μ值對應的懲罰函數q(x)

2.2 修正動態彎曲距離

假定時間序列S的長度為m：S={s1,s2,…,si,…,sm}，時間序列T的長度為n：T={t1,t2,…,tj,…,tn}。基于懲罰函數計算動態彎曲距離時，采用計數矩陣A記錄時間序列S上每個數據點在水平和垂直方向連續轉移的步數，計數矩陣B記錄時間序列T的數據點連續轉移的步數。UDTW的計算步驟如下：

步驟1 初始化計數矩陣A=(ai)1×m和B=(bj)1×n，A與B初始化為零矩陣。

步驟2 計算時間序列S與T之間的距離矩陣D=(dij)m×n，其中dij=d(i,j)=|si-tj|。

步驟3 計算累積距離矩陣R=(γij)m×n中第1行數據γ(s1,T)以及第1列數據γ(S,t1)。

步驟4 按照行(或列)的方向依次求解其他γ(i,j)的值，同時，當路徑更新到γ(i,j)時，記錄下當前點si和tj被使用的次數，并更新計數矩陣ai和bj的值。在計算γ(i,j)時，按照如下公式求解：

γ(i,j)=min {c(bj)×d(i,j)+γ(i-1,j),d(i,j)+

γ(i-1,j-1),c(aj)×d(i,j)+γ(i,j-1)}

(5)

其中

(6)

步驟5 重復執行步驟4，直至求得修正動態彎曲距離UDTW(S,T)=γ(m,n)。

算法1中給出了UDTW中生成計數矩陣和距離矩陣的偽代碼。

算法1 GENMATRIX(S,T)。

輸入：時間序列S、T。

輸出：計數矩陣A、B，距離矩陣D。

m=len(S),n=len(T)

A=zeros(m),B=zeros(n)

fori=1:m

forj=1:n

d(i,j)=|si-tj|

returnA、B、D

算法2中給出了求解修正動態彎曲距離的偽代碼。

算法2 GENUDTW(A,B,D,q(x))。

輸入：計數矩陣A、B，距離矩陣D，懲罰函數q(x)；

輸出：修正動態彎曲距離UDTW(S,T)。

fori=1:m

R(i,1)=c(B[1])×d(i,1)+R(i-1,1)

B[1]=B[1]+1,A[i]=A[i]+1

forj=1:n

R(1,j)=c(A[1])×d(1,j)+R(1,j-1)

A[1]=A[1]+1,B[j]=B[j]+1

forj=2:n

fori=2:m

T1=c(B[j])×d(i,j)+R(i-1,j)

T1=d(i,j)+R(i-1,j-1)

T3=c(A[i])×d(i,j)+R(i,j-1)

Index,R(i,j)=MIN{T1,T2,T3}

IfIndex=0:B[j]=B[j]+1,A[i]+1

elseifIndex=2:A[i]=A[i]+1,B[j]=1

elseA[i]=1,B[j]=1

returnR(m,n)

2.3 時間復雜度分析

本文所提出方法是在PLM的基礎上使用UDTW度量方法，因此可簡記為PLM-UDTW(Piecewise Local Max-smoothing-Updated Time Dynamic Warping)。首先通過PLM方法對序列進行特征提取，將長度為m的序列壓縮成長度為L的短序列，此處時間復雜度為O(m)。對降維后的序列作UDTW度量時增加了一個更新計數矩陣的操作，但該操作并沒有增加其時間復雜度，因此，對兩條長度為m的序列作UDTW度量的時間復雜度為O(2m+L2)，當壓縮率(m/L)為10%，m=128時，PLM-UDTW的時間復雜度降低了97.4%。

3 實驗研究與分析

采用來自UCR[14]時間序列數據集，經過Z-scores(均值為0、方差為1)規范化處理。運用本文提出的修正算法進行實驗分析并驗證其有效性。數據集分成訓練集和測試集兩部分，采用“1-近鄰”分類方法，運用訓練集學習生成分類器，運用測試集驗證分類器的準確率。本文算法均使用Python 3.6代碼實現。

3.1 分段特征提取

(7)

β=max{α2,α3,…,αi,…,α20}

(8)

此處wmax=20，αi是當窗口長度w=i時算法對應的分類準確率。表1給出了15個時間序列數據集的相關信息，其中I.P.D(ItalyPowerDemand)、M.I(MedicalImages)、T.L.ECG(TwoLeadECG)為對應數據集的縮寫模式。

表1 時間序列數據集信息

通過表2結果可以看出，PLM方法在Computers、ECG、Face(four)、Gun-Poin、tLightning- 7、M.I、OliveOil、T.L.ECG這8個數據集上的平均分類準確率和最高準確率均高于PAA方法，在Beef、Coffee、Lightning- 2、Trace這4個數據集上最高分類準確率與PAA方法相同，但平均分類準確率比PAA方法高。證明與PAA方法相比，PLM方法至少可以提高12個數據集上的分類準確率。為了進一步說明PAA方法和PLM方法在具體數據集上的表現，圖4給出了Gun-Point數據集在不同窗口長度下分類準確率的變化情況。通過觀察可得，隨著窗口長度的增大，經過PAA方法提取特征后的數據分類準確率逐漸降低，而PLM方法提取特征后的分類準確率在一定范圍內上下波動，當窗口長度w=20時(壓縮率為5%)，PLM方法的分類準確率依然保持約90.67%，遠高于PAA方法的78.67%。實驗結果說明，PLM方法具有更好的魯棒性，其對窗口長度的依賴性更小。與PAA方法相比，在未增加時間復雜度的基礎上具有更好的特征提取效果。

表2 PAA與PLM方法分類準確率對比

圖4 Gun-Point數據集不同窗口長度w下的分類準確率對比

另外，通過對表2結果分析可得，最佳窗口長度大多集中在2～10，即當壓縮率范圍為50%～10%，分段特征提取方法能取得較好效果。

3.2 UDTW度量效果

UDTW中的懲罰函數控制了動態彎曲路徑的修正程度，有兩個重要參數對分類效果產生直接影響：一個是懲罰系數的上界值qmax，決定了最大懲罰系數；一個是達到上界值所對應的最大彎曲步數μ。由于這兩個參數都是對懲罰函數的傾斜度作調整，且上界值對所有彎曲路徑是同一標度下的，因此可將qmax固定為1。實驗通過調整μ值大小驗證UDTW算法的分類效果。圖5給出了μ值在[0,10]區間變化時，上述15個數據集分類準確率的波動情況。為提高UDTW計算效率，實驗先采用PLM方法對時間序列進行特征提取，此處窗口長度設置為w=7。

參數μ的調整實際是在歐氏距離和DTW距離之間找到一個最優平衡，當μ=0時，UDTW距離向歐氏距離靠近，當μ足夠大時，UDTW距離無限逼近DTW距離。針對不同數據集，達到最佳分類準確率的μ值有所不同，說明不同數據集對時間序列形態特征的關注度有所不同。例如數據集ECG的最優參數為μ=0，在該參數下達到最高分類準確率91%，這說明對ECG數據集的分類更關注數據的形態特征，因此傾向于用歐氏距離進行相似性度量；文獻[15]也證明了這一觀點，即對于ECG數據集來說，歐氏距離作為度量距離能取得更好的分類效果。隨著μ值增大，同一數據集的分類準確率會呈現一定的波動，但對于I.P.D、M.I等數據集，分類準確率沒有明顯變化，經過分析發現原因在于原始序列長度較短，在分段特征提取之后保留的形態特征過少，因此對于這些數據集來說，一個解決的方法是減小特征提取的窗口長度，或者不作分段降維處理。

圖5 15個數據集的不同μ值下的分類準確率對比

進一步，為了說明PLM-UDTW的有效性，選擇歐氏距離、DTW、導數動態彎曲距離(Derivative Dynamic Time Warping， DDTW)[15]、PLM-UDTW四種距離度量方法，采用1-近鄰分類方法直接對上述15個數據集進行分類，表3給出了4種算法在每一數據集下的分類錯誤率，由于PLM-UDTW對序列進行了分段處理，因此給出達到該分類錯誤率下的一種窗口長度w和對應的最大彎曲步數μ。錯誤率的計算公式如下：

Errorrate=測試集中分類錯誤的個數/測試集大小

(9)

表3 PLM-UDTW與傳統度量方法的錯誤率對比

從表3中可以看出：PLM-UDTW方法在CBF、Coffee、Trace這3個數據集上錯誤率為0，即達到了100%分類正確；在Beef、Computers、ECG等14個數據集上分類準確率均為4種方法中最高；在Lightning- 2數據集上準確率僅次于DTW方法。與DDTW度量方法相比，PLM-UDTW在CBF、Computers、Ham三個數據集上準確率分別提高了71.8%、62.6%和47.07%，在15個數據集上準確率平均提高了27.7%。與歐氏距離度量方法和傳統DTW度量方法相比，PLM-UDTW的分類準確率也有顯著提高，最高分別提升了65.7%和109.4%；在15個數據集上準確率平均提高了21%和17.7%，進而證明了本文方法在時間序列分類中的有效性和優越性。

3.3 計算時間開銷

根據2.3節描述，分段降維方法能有效降低UDTW的時間復雜度。將時間序列的分類算法分成三個操作過程：1)分段降維過程；2)測試集、訓練集距離度量過程；3)1-近鄰(1-Nearest Neighbor, 1-NN)分類過程。表4給出ECG數據集在UDTW、DTW、DDTW三種度量方法下的時間開銷。

表4 3種算法的平均運行時間對比 s

根據表4分析，DDTW運行時間最長，這是因為DDTW方法比DTW方法增加了一個求數值導數的操作。運行時間最短的是UDTW，盡管與傳統DTW和DDTW方法相比新增加了一個分段降維過程，但是算法消耗的總時間減少了很大比例，這個比例主要依賴于窗口長度w。當w=7時(壓縮率約為14%)，耗時約減少99%。綜上所述，PLM方法能有效降低時序特征分類算法的時間復雜度，與DTW、DDTW方法相比，PLM-UDTW在不影響分類準確率的前提下極大地提高了計算效率。

4 結語

本文提出了一種基于路徑修正的DTW(UDTW)度量方法，解決了DTW易陷入過度彎曲而忽略時間序列形態相似性的問題。通過對原始序列進行分段特征提取降低了DTW度量的計算代價，提高了算法的整體效率；同時，對PAA分段降維方法作了改進，使其能更好地提取時間序列的曲線形態特征。實驗結果表明，PLM-UDTW算法可以提高15個時間序列數據集中大部分數據集分類準確率，并明顯提高分類速度。