基于系數矩陣弧微分的時間序列相似度量

王智博,林 意,曹洋洋

(江南大學 數字媒體學院,江蘇 無錫 214122)

0 概述

將某一個統計指標的各個數值按時間先后順序排列便構成了時間序列。從金融領域到科學工程,從天文氣象到社會學[1-2],時間序列無處不在。由于實際應用中的時間序列往往具有高維、規模巨大、易受噪聲干擾等特點[3-4],直接在原始時間序列上進行數據分析、處理和挖掘變得非常困難,因此在對時間序列挖掘之前進行有效的預處理成為解決上述問題的關鍵。這其中時間序列特征表示和相似度量是預處理的關鍵[5-6]。

相似度量是時間序列挖掘中一項重要的基礎任務,主流的度量算法通常自定義一個距離函數,選取的自變量為離散序列點坐標及其變形,序列之間的距離越小則序列越相似。常見的算法有:歐氏距離[7](Euclidean Distance,ED),設定單一距離閾值,容易理解且算法簡單;動態時間規整[8](Dynamic Time Warping,DTW),借鑒語音數據處理的思路并運用動態規劃思想,通過彎曲時間軸來實現相似性度量;符號化距離[9],將時間序列預處理為字符串,利用查詢等概率劃分的正態分布完成相似度量;基于條件復雜性距離[10],嵌入信息論和計算理論,關注算法運行過程中的連接和壓縮操作,借助壓縮率來反映數據之間的相似性。

本文引入數理統計和回歸分析中的最小二乘思想,通過若干離散序列點并利用偏微分構建系數矩陣方程,從而求得擬合多項式的不定參數。該參數以向量形式存在且刻畫曲線的基本形態,被稱為向量基。由于實現了離散序列點的連續化,因此可利用幾何連續的性質來繼續研究問題。本文給出最小相似點和微分三角形的概念,分析對比連續曲線的弧微分與曲率半徑微分的關系,發現當以最小相似點為端點構成的微分三角形相似時,可以使得微分三角形對應的弧微分與曲率微分成等比關系,從而得出弧微分相似判定等式,最后根據分治、遞歸思想,判斷若在某一連續區間內2條曲線所有最小相似點都在判定等式的合理閾值范圍內,可以得出這2個序列相似的情況。當候選的2條序列時間粒度不相同時,本文算法具有無需人工干預也能完成序列的相似度量的優點;當候選的2條序列長度不相同時,該算法也能彌補傳統算法不能實現形態相似度量的不足。

1 時間序列相關定義與問題描述

1.1 時間序列相關定義

定義1(時間序列) 時間序列是由記錄時間和記錄值組成的有序集合。對于給定的有限時間集T、非空狀態屬性集A=〈A1,A2,…,Am〉及其對應值域DAj,時間序列X表示如下:

X=〈X1,X2,…,Xn〉

(1)

定義2(時間序列的模式表示)[12]時間序列模式指時間序列的某種變化特征,通過提取時間序列的模式將其變換到模式空間,即得到時間序列的模式表示。設有時間序列X=〈x1,x2,…,xn〉,其模式表示為:

X(t)=f(w)+e(t)

(2)

其中,f(w)是時間序列的模式表示,e(t)是時間序列與其模式表示之間的誤差。

定義3(時間序列的分段線性表示) 設有時間序列X=〈x1,x2,…,xn〉,則其分段線性表示為:

(3)

其中,fi(t,wi)表示連接時間序列分段點的線性函數,ei(t)是時間段內時間序列與其分段線性表示之間的誤差。

定義4(閔氏度量) 設有2條長度為n的時間序列Q=〈Q1,Q2,…,Qn〉和C=〈C1,C2,…,Cn〉,則它們之間的閔氏度量為:

(4)

其中,p為可變參數。當p=2時,閔氏度量即為使用最為廣泛的歐式距離[7]。

1.2 問題描述

傳統基于點距離的時間序列相似度量,如歐式距離[7]是利用式(4),通過計算2條序列一一對應點之間的距離得到最終度量。該算法的實現依賴2個充分條件:1)候選的2條序列等長;2)2條序列一一對應的點坐標在時間軸上的投影重合。

如圖1所示,有3條時間序列A、B、C,根據歐式距離公式,代入對應的序列點坐標值得出D2(A,C)

圖1 時間序列示意圖

圖2所示為2條時間跨度不相等的時間序列。根據歐式距離的充分條件,當候選序列時間跨度不相同時,歐式距離算法失效。解決辦法是應用動態時間軸彎曲距離算法(DTW)[8]。其中的2個核心步驟是:1)動態時軸彎曲或動態時間規整;2)距離測度計算。DTW算法的本質是尋找一個合適的函數j=w(i),將序列A的時間軸非線性地映射到序列B的時間軸上,使得A的第i個序列點與B的第j個序列點對齊,并且使每組對齊點達到距離最小,如圖3所示。但該算法時間效率較低,不利于大量較長時間序列的相似度量。圖3為DTW算法示意圖。

圖2 時間跨度不相同的情況

圖3 DTW算法示意圖

設A、B時間序列是某一超市同類的2種商品1天內銷售額的序列。它們跨度相同,都為12個月。但序列A統計的是月銷售額,序列B統計的是季度銷售額。2條序列時間跨度相同描述間單位不同,即刻畫序列單位的粒度不同。為能應用經典的基于點對點距離的算法,需要人工干預,使用時刻對等使得待比較的2條序列具有相同的粒度,即算法對時間粒度的敏感性不強,如圖4所示。

圖4 時刻對等示意圖

綜上可以發現,造成這些問題的原因是所選的序列之間距離度量函數的自變量為單一離散點在坐標軸中絕對位置坐標,沒有很好捕捉到因各種原因造成的坐標偏移,使算法只能局限于序列微觀上的相似度量(距離相近)缺少了對于序列宏觀的相似度量的魯棒性(形態相近);并且距離函數自變量的選取忽視了對序列的形態識別能力,造成度量結果的不合理;與此同時,離散的思維使得不能應用更成熟的連續幾何性質去做繼續的研究,不能挖掘到每一條時間序列的本質權值,造成不能對規模巨大的待比較序列根據權值做分類處理,不利于數據挖掘的后續工作,例如時間序列的相似性搜索[13]。

2 基于系數矩陣弧微分的相似度量理論分析

2.1 系數矩陣

為了能夠應用連續幾何的性質,首先需要解決的問題就是離散點的連續化。本文給出方法是利用統計學[14-15]中回歸模型——最小二乘法,其一般形式為y=f[x|θ]+ε,其由參數θ決定的回歸函數,ε是不可觀測的隨機誤差。目標是使得觀測點和估計點的距離平方達到最小,從而誤差達到最小。

設待擬合曲線的函數為:

y=a0+a1x+a2x2+…+akxk

(5)

假設時間序列的已知觀測點個數為n,根據最小二乘理論可知各點到這條曲線的距離和為:

(6)

求使得Q(θ)最小的a0,a1,…,ak值,對每一個a求偏導:

令上式的偏導都為0,化簡得:

將這組等式表示成矩陣形式:

即XA=Y,解此矩陣方程求出A,即可得到最佳的擬合曲線。其中X即為系數矩陣,由已知的觀測點和曲線方程的最高階數所決定,構造系數矩陣的目的就是把研究的關注點從離散轉換成連續。

2.2 向量基

在矩陣方程XA=Y中,求得的A向量是擬合曲線未知數的系數,顯然,這些系數刻畫著曲線的形態。形態不隨著曲線在坐標系中絕對位置的改變而改變;每一條曲線都有自己的固有形態,又因為形態由A向量所決定,所以,本文將A向量稱之為向量基。

2.3 微分三角形

設左邊曲線的方程為Y=Y(X),X∈D;相應地,右邊曲線的方程為y=y(x),x∈d,其中D、d為各自的定義域。設他們至少存在三階導數,且二階導數處處不等于0。OA、OB分別是曲線Y=Y(X)上A、B點的曲率中心,曲率圓半徑分別為RA、RB;相應地,rD、rE分別是曲線y=y(x)上點D、E對應的曲率圓半徑。過點A做割線AB的垂線,并截取AC=RA-RB;過點D做割線DE的垂線,并截取DE=RD-RE。于是,由點A、B、C和對應的點D、E,F組成一組對應的直角三角形。由已知得,2條曲線光滑且連續。

(7)

滿足上式的點A、D稱之為最小相似點,并稱2條曲線在A、D處最小相似,以點A和點D為直角端點構成的2個直角三角形稱之為微分三角形。

圖5 微分三角形

2.4 最小相似點判定方程

在連續幾何圖形中可知,曲率描述著曲線和的彎曲程度,由曲線弧的長度和切線夾角所決定,如圖6所示。

圖6 曲率定義示意圖

(8)

同理,有:

(9)

將式(8)和式(9)代入到最小相似點等式,得到最小相似點判定方程:

(10)

2.5 曲線相似的判定

首先引入2個幾何圖形相似的判定公理:如圖形D上點與圖形D′上兩對應點的線段之比,是恒定的非零常量,就認為圖形D與D′相似。

然后給出曲線完全相似的判定定理[16]:2條曲線在所考慮的區間內同向,對應的函數都存在至少三階導數且二階導數處處不為0。若滿足關系式(11),則2條曲線在給定的區間內完全相似。在式(11)中,C為非零常量。

(11)

最后,給出完全相似的證明過程。

圖7 曲線微小分量相似示意圖

當曲線上相鄰兩點Ak、Ak+1與對應的Bk、Bk+1表示兩點間距離最大者,且Ak→Ak+1、Bk→Bk+1時,得出曲線弧AiAi+1Ai+2…An～曲線弧BiBi+1Bi+2…Bn。再以An,Bn為新的對應頂點(An、Bn是最小相似點),重復以上步驟,又得新的相似弧段,且相似比仍為C。把各個相似弧段順序連接起來,得到2個邊長為微小量的相似多邊形。

綜上,將待比較的2條曲線分治為多個以最小相似點為頂點構成的多邊形,遞歸地論證個多邊形的相似比仍為C。根據公理,任意兩組相似點的之比即對應的對角線之比為常數,因此,2條曲線完全相似,證畢。

3 基于系數矩陣弧微分的相似度量算法

分析第2節各概念的推導過程可以發現:1)離散的原始序列點是通過最小二乘法構建的矩陣完成連續化的;2)微分三角形和最小相似點的定義借助了弧微分及其所建立的比例形式;3)最小相似點的判定方程是依靠弧微分導出的變量——曲率來完成的;4)曲線相似的定理及其證明過程也用到了弧微分的概念,不難發現,弧微分在本文理論中的核心地位,因此,將本文算法命名為基于系數矩陣弧微分的時間序列相似度量算法(CMAD)。

第2節給出的2條曲線相似的定理,是使用等式建立的,結果就是2條候選曲線完全相似或者不相似,無法對其進行相似性的其他微小度量,算法中不可以直接應用,因此,本節引入一個概念——互相關函數。此概念來自信號分析[17-18],描述了隨機信號x(t)、y(t)在任意2個不同時刻t1、t2之間的相關程度,它是在某一頻域內2個信號是否相關的一個判斷指標,定義為R(u)=x(t)*y(-t),其中*表示卷積,其結果大小所表示的意義在統計學界通常有如表1所示的結論。

表1 相關程度

根據第3節式(9),當有如式(12)所示關系時,2條曲線完全相似。

(12)

構建2個互相關函數:

其中:

可做以下分析:當函數A(u)的結果處于相關程度較高的區間占整個結果的比率越高,說明a(Y)與b(y)越相關,越滿足式(12)的第1個等式;同理,當函數B(v)的結果在滿足上一個條件的同時,并且在給定的區間內結果穩定,說明c(Y)與d(y)越相關,越滿足式(12)的第2個等式。當2個條件都滿足時,說明2條曲線越相似,又曲線是原始時間序列根據刻畫的形態向量基連續化而來,因此兩條時間序列也就越相似,達到了相似性度量的目的。本文算法描述如下:

算法基于系數矩陣弧微分的時間序列相似度量算法(CMAD)。

輸入原始時間序列X=,Y=

輸出互相關函數結果序列。

1)判斷2條原始序列的長度是否為大量,如果超出一定規模,使用線性分段表示對序列進行壓縮、降維,否則轉步驟2)。

2)對2條序列的離散點應用系數矩陣進行連續化。

3)依據式(11)求出各等式要素。

4)使用各要素構建互相關函數。

5)在給定的區間內,對互相關函數結果進行檢測(當m=n時,只需在原始離散點所在同一區間進行檢測;當m≠n時,為了說明程序的魯棒性,需要在原始離散點的2個不同區間分別檢測)。

6)根據檢測結果,對時間序列X、Y的相似性進行綜合判定。

在上述算法中,步驟1)判斷原始序列可通過線性掃描,在對序列壓縮可選用線性分段表示,這兩小步的時間復雜度均為O(n)。步驟2)和步驟3)通過最小二乘思想求系數矩陣的元素,計算(n+1)個偏導數方程:T1(n)=O((n+1)×f(n))=O(f(n))=O(n)。步驟4)～步驟5)所構建互相關函數的最高次冪為2,所以,T2(n)=O(f2(n))=O(n2)。綜上所述,CMAD算法時間復雜度為O(n2)。

4 實驗結果及分析

4.1 實驗環境

本文實驗運行環境為:CPU2.0 GHz、內存8 GB、500 GB硬盤,Windows7系統上實現。開發工具為Matlab2014a。

4.2 實驗方法

為驗證CMAD算法的可行性和優越性,本文將做3個實驗。實驗1的關注點是:待比較序列是等跨度且時間粒度不同;實驗2關注點是:對原始序列進行分段線性表示后,再對2條等跨度、粒度不同的序列進行相似度量;實驗3的關注點是:2條時間序列跨度不相同。CMVB算法針對擬合函數,y=a0+a1x+a2x2+…+akxk要求連續化后的函數至少存在三階導數且二階導數處處不為零,為了算法簡便,3個實驗統一把函數未知數x的最高次冪定為3。

4.2.1 實驗1

實驗數據來自http://www.bundesbank.de/網站,是證券市場版塊下的Time series WU0053:Gross sales of domestic debt securities at nominal value。2個時間序列數據都選自1985年—1989年,時間序列A的取樣周期為每個季度一次,共20個數據;時間序列B的取樣周期為每個月一次,共60個數據。2個序列如圖8所示。

圖8 不同粒度的序列

時間序列A、B是同一debt securities的不同時間粒度的2組觀測值,本質上是同一事物,如果想要應用傳統的基于離散點距離函數的算法,首先需要人工時刻對等;其次結果往往是不能很好判定2條序列相似的。應用CMVB算法,互相關函數在同一區間結果如表2所示,穩定性分析結果如圖9所示。

表2 相關程度分布 %

圖9 穩定性分析結果

根據表2可知:1)函數A(u)的結果處于相關程度高的區間的比率達58.3%;2)函數B(v)的結果較為穩定且在高相關的區間比率高。由上述結論可知時間序列A、B相似,符合真實數據本來結果,驗證了CMVB算法的可行性。

如果用傳統的時間序列相似度量方法,需要增加以下時刻對等步驟:

1)將tA序列和tB序列進行歸并,序列值合并并去除重復值,得到對等后的標準時刻序列:t={t1,t2,…,tw}。

2)通過遍歷循環分別找到tA序列和tB序列的值在t序列中的位置序列:loc1序列和loc2序列。

3)遍歷A和B序列,對A和B序列進行插值補充,最終得到時刻對等后的序列Aplr和Bplr。

4.2.2 實驗2

實際應用中的時間序列數據往往是海量、高維、易受干擾的,直接對原始時間序列進行挖掘不僅時間和空間效率低下,而且算法的可靠性和準確性也容易受到影響,因此,在對時間序列挖掘之間往往需要壓縮等預處理,常見方法時分段線性表示方法,本文在此應用基于斜率提取邊緣點的時間序列分段線性表示方法PLR_SEEP。

設有時間序列X=,提取點集合為〈xi1,xi2,…,xik〉,且1≤i1≤i2≤…≤ik≤n。

根據式(3),時間序列的PLR_SEEP[19]表示為:

(11)

其中,L(x,y)表示連接趨勢點x和y之間的線性函數。公式可以簡單表示為:

XT=〈L(xi1,xi2),L(xi2,xi3),…,L(xik-1,xik)〉

本實驗的原始數據與實驗1相同,首先應用SEEP算法對序列B進行預處理,使其壓縮率為15%,保留45個數據;序列A點較少,不進行預處理,2條序列如圖10所示。

圖10 分段線性表示預處理

預處理后再對序列A,XT應用CMVB算法,互相關函數的結果如表3所示。由圖11和實驗1的分析可知,序列A、XT相似,從而序列A、B相似。

表3 分段線性表示后的相關程度分布 %

圖11 分段線性表示后的穩定性分析結果

4.2.3 實驗3

實驗數據來自我國某一海港港口,記錄了每個月的集裝箱月吞吐量。序列A記錄了2013年12個月的吞吐情況;序列B記錄了2013年—2014年24個月的吞吐情況,如圖12所示。

圖12 時間跨度不相同的2條序列

顯然待比較的2條序列時間跨度不相同,傳統的基于離散點距離函數的算法失效,基于動態時間彎曲的算法又過于復雜,在此應用CMVB算法。因為2條序列的區間不相同,為了驗證算法的可靠性,構造的互相關函數需要分別在2個區間上驗證,結果如表4、表5和圖13和圖14所示。可以看出,互相關函數的結果在2個區間都有較好反映,說明了2個區間原始序列較為相似,符合直觀認知。

表4 互相關函數12個月記錄相關程度分布 %

表5 互相關函數24個月記錄相關程度分布 %

圖13 CMAD算法12個月記錄穩定性分析結果

圖14 CMAD算法24個月記錄穩定性分析結果

綜合以上3個實驗,可以得出以下結論:CMAD算法對于跨度相同、時間粒度不同的候選序列不需要人工干預也可很好地完成相似性度量的任務,并且對于數據規模較大序列應用分段線性表示后,仍可以較好地進行相似性度量,算法有著較強的穩定性;對于時間跨度不相同的序列,可以完成2條序列宏觀意義上的相似性——形態相近,具有較強的魯棒性。

5 結束語

選擇一個合適的時間序列相似度量算法是時間序列數據挖掘的重要前提。本文提出的CMAD算法首先利用系數矩陣對離散序列點連續化,然后分析連續曲線的弧微分與曲率半徑微分的關系,找出相似性判定等式,最后通過互相關函數完成最終的時間序列相似性度量。實驗結果表明,該算法可同時完成距離相近度量和形狀相似度量,具有良好的適用性和可行性,利于后續數據挖掘的工作進程。下一步將在本文算法基礎上進行海量數據的時間序列相似性度量。

[1] ZHAI Yuanzheng,WAHG Jinsheng,TENG Yangguo,et al.Water Demand Forecasting of Beijing Using the Time Series Forecasting Method[J].Journal of Geographical Science,2012,22(5):919-932.

[2] FCUHS E,GRUBER T,NITSCHKE J,et al.On-line Segmentation of Time Series Based on Polynomial Least-squares Approximation[J].IEEE Transactions on Pattern Analysis and Machine Intelligence,2010,32(12):2232-2245.

[3] MACIEJ K,GRAZYNA S.An Approach to Dimensionality Reduction in Time Series[J].Information Science,2014,26(6):15-36.

[4] GUERRERO J L,BERLANGA A,GARCIA J,et al.Piecewise Linear Representation Segmentation as a Multiobjective Optimization Problem[M]//JANUSZ K.Advances in Intelligent and Soft Computing.Berlin,Germany:Springer,2010:267-274.

[5] MUEEN A,DING H,TRAJCEVSKI G,et al.Experimental Comparison of Representation Method and Distance Measures for Time Series Data[J].Data Mining and Knowledge Discovery,2012,26(2):275-309.

[6] EHMKE J F,MEISEL S,MATTFELD D C.Floating Car Based Travel Times for City Logistics[J].Trans-portation Research,Part C:Emerging Technologies,2012,21(1):338-352.

[7] AGRAWAL R,FALOUTSOS C,SWAMI A.Efficient Similarity Search in Sequence Databases[C]//Proceedings of the 4th International Conference on Foundations of Data Organization and Algorithms.Washington D.C.,USA:IEEE Computer Society,1993:69-84.

[8] KEOGH E,PAZZANI M.Derivative Dynamic Time Warping[C]//Proceedings of the 1st SIAM Inter-national Conference on Data Mining.Chicago,USA:SIAM,2001:1-11.

[9] LIN J,KEOGH E,LONARDI S,et al.A Symbolic Representation of Time Series,with Implications for Streaming Algorithms[C]//Proceedings of the 8th ACM SIGMOD Workshop on Research Issues in Data Mining and Knowledge Discovery.New York,USA:ACM Press,2003:2-11.

[10] KEOGH E,LONARDI S,RATANAMAMATANA C A,et al.Compression-based Data Mining of Sequential Data[J].Data Mining and Knowledge Discovery,2007,14(1):99-129.

[11] 潘 定,沈鈞毅.時態數據挖掘的相似性發現技術[J].軟件學報,2007,18(2):246-258.

[12] NOPIAH Z M,KHAIRIR M I,ABDULLAH S,et al.Peakvalley Segmentation Algorithm for Kurtosis Analysis and Classification of Fatigue Time Series Data[J].European Journal of Scientific Research,2009,29(1):113-125.

[13] SAKOE H,CHIBA S.Dynamic Programming Algorithm Optimization for Spoken Word Recognition[J].IEEE Transactions on Acoustics,Speech,and Signal Processing,1978,26(1):43-49.

[14] GUESTRINT C,BODIKZ P,THIBAUXT R,et al.Distributed Regression:An Efficient Framework for Modeling Sensor Network Data[C]//Proceedings of ACM International Conference on Sensor Networks.New York,USA:ACM Press,2004:1-10.

[15] DELIGIANNAKIS A,KOTIDIS Y,ROUSSOPOULOS N.Compressing Historical Information in Sensor Networks[C]//Proceedings of ACM SIGMOD International Conference on Management of Data.New York,USA:ACM Press,2004:527-538.

[16] 張智廣,趙學敏.平面曲線相似性初探[J].天津師范大學學報,1998,18(2):65-72.

[17] BECKOUCHE S,MA Jianwei.Simultaneous Dictionary Learning and Denoising for Seismic Data[J].Geophysics,2014,79(3):27-31.

[18] SONG Jun,LIU Yu,WANG Xudong.Improved Denoising Algorithm for Narrow-band Signal and Its Application[J].Journal of Vibration and Shock,2013,32(16):59-62.

[19] 詹艷艷,徐榮聰,陳曉云.基于斜率提取邊緣點的時間序列分段線性表示方法[J].計算機科學,2006,33(11):139-142.