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

常炳國，臧虹穎

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

0 引言

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

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

1 時序特征提取

假定一條時間序列S的長度為m：S={s1,s2,…,si,…,sm}，另一條時間序列T的長度為n：T={t1,t2,…,tj,…,tn}。創(chuàng)建一個m×n的距離矩陣D，用d(ith,jth)表示(si,tj)兩點的路徑距離：d(i,j)=‖si-tj‖p，其中‖x‖p表示p-范數(shù)。計算DTW距離時，先構(gòu)造一個累積距離矩陣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]，利用最小二乘法求得分段序列的最佳線性擬合曲線，該方法能較好地還原原始序列的形態(tài)，但對于長度為m，分段數(shù)目為L的序列來說，時間復雜度高達O(m2L)。類似的方法還有分段多項式表示(Piecewise Polynomial Representation， PPR)、分段回歸近似(Piecewise Regression Approximation， PRA)和自適應分段常量近似(Adaptive Piecewise Constant Approximation， APCA)等，雖然這些方法都對序列作了分段處理以達到降維目的，但這些方法本身時間復雜度較高，算法效率較低[11]。

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

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

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

l=1,2,…,L

(2)

其中AVG(Sl′)表示第l個窗口內(nèi)數(shù)據(jù)子集Sl′的均值，窗口長度w=m/L(或m/L+1)。通過將每段的均值數(shù)據(jù)相連接形成新序列SPAA，從而實現(xiàn)數(shù)據(jù)降維的目的。PAA具有概念簡單、參數(shù)少(僅只有一個參數(shù)L)、時間復雜度低等優(yōu)點，是一種有效的特征提取方法，但它僅反映了序列整體的變化趨勢，而忽略了時間序列的形態(tài)特征，并且，PAA對均值平穩(wěn)的獨立噪聲數(shù)據(jù)不敏感，存在重要數(shù)據(jù)點或異常數(shù)據(jù)點信息丟失的不足[12]。綜合考慮PAA的優(yōu)缺點，本文采用窗口最大值代替平均值的數(shù)據(jù)提取策略，然后用分段平均值對最大值進行適度平滑處理以消除噪聲影響，本文將此種方法稱為分段局部最大值平滑法(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個窗口最大值與平均值的差值，β為平滑常數(shù)，用于調(diào)整數(shù)據(jù)的變化范圍，由于前期對序列作了Z-score規(guī)范化處理，β取1即可。圖1給出了長度為48的序列S分別經(jīng)過PLM和PAA方法提取特征值后所形成的特征序列SPLM、SPAA，此處窗口長度w=6。

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

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

2 修正算法

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

2.1 懲罰函數(shù)

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

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

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

圖2 彎曲路徑與點對點匹配結(jié)果

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

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

(4)

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

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

2.2 修正動態(tài)彎曲距離

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

步驟1 初始化計數(shù)矩陣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行數(shù)據(jù)γ(s1,T)以及第1列數(shù)據(jù)γ(S,t1)。

步驟4 按照行(或列)的方向依次求解其他γ(i,j)的值，同時，當路徑更新到γ(i,j)時，記錄下當前點si和tj被使用的次數(shù)，并更新計數(shù)矩陣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 重復執(zhí)行步驟4，直至求得修正動態(tài)彎曲距離UDTW(S,T)=γ(m,n)。

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

算法1 GENMATRIX(S,T)。

輸入：時間序列S、T。

輸出：計數(shù)矩陣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中給出了求解修正動態(tài)彎曲距離的偽代碼。

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

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

輸出：修正動態(tài)彎曲距離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的基礎(chǔ)上使用UDTW度量方法，因此可簡記為PLM-UDTW(Piecewise Local Max-smoothing-Updated Time Dynamic Warping)。首先通過PLM方法對序列進行特征提取，將長度為m的序列壓縮成長度為L的短序列，此處時間復雜度為O(m)。對降維后的序列作UDTW度量時增加了一個更新計數(shù)矩陣的操作，但該操作并沒有增加其時間復雜度，因此，對兩條長度為m的序列作UDTW度量的時間復雜度為O(2m+L2)，當壓縮率(m/L)為10%，m=128時，PLM-UDTW的時間復雜度降低了97.4%。

3 實驗研究與分析

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

3.1 分段特征提取

(7)

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

(8)

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

表1 時間序列數(shù)據(jù)集信息

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

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

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

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

3.2 UDTW度量效果

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

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

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

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

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

(9)

表3 PLM-UDTW與傳統(tǒng)度量方法的錯誤率對比

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

3.3 計算時間開銷

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

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

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

4 結(jié)語

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