基于三次樣條趨勢濾波和知識 指導(dǎo)的大規(guī)模不等長時(shí)序族聚類

冀東麗

( 山東特殊教育職業(yè)學(xué)院，250022，濟(jì)南 )

1 引 言

時(shí)間序列是指在一定的時(shí)間間隔內(nèi)所記錄的一系列具有時(shí)間意義的序列數(shù)據(jù)，在金融、經(jīng)濟(jì)、氣象、通信工程、醫(yī)療衛(wèi)生等領(lǐng)域有著廣泛的應(yīng)用.從大量時(shí)序數(shù)據(jù)中挖掘潛在的信息和模式，是數(shù)據(jù)挖掘領(lǐng)域中最具挑戰(zhàn)性的大問題之一[1]，而聚類分析能有效挖掘數(shù)據(jù)中的重要信息，是數(shù)據(jù)挖掘中的重要任務(wù)之一.在實(shí)際問題中，由于時(shí)間段的數(shù)據(jù)丟失及數(shù)據(jù)保密性等，會導(dǎo)致遇到的時(shí)序數(shù)據(jù)族并不等長，因此開展大規(guī)模不等長時(shí)序數(shù)據(jù)族的智能聚類研究更具有現(xiàn)實(shí)意義.

目前智能聚類研究主要包括兩個(gè)方面：一是由于數(shù)據(jù)族中各數(shù)據(jù)長度不一致，即包含的數(shù)據(jù)點(diǎn)的個(gè)數(shù)不相同，因此如何度量任意兩數(shù)據(jù)之間的相似關(guān)系便成為首要解決的問題.常用的歐氏距離能夠快速實(shí)現(xiàn)對兩個(gè)等長時(shí)間序列的距離度量，但對于不等長時(shí)間序列卻無能為力，因此，常通過動態(tài)時(shí)間扭曲方法(dynamic time warping, DTW)[2-4]實(shí)現(xiàn)不等長序列的最佳匹配.由于DTW算法時(shí)間較長，大量研究主要集中在提高DTW的時(shí)間性能上[5-9].這些方法雖然能提高時(shí)間效率，但往往以犧牲時(shí)間特征為代價(jià)，聚類結(jié)果不夠理想.另一方面，涌現(xiàn)了大量關(guān)于時(shí)間序列形態(tài)的研究，通過各種方式度量時(shí)間序列之間的相似性，比如Keogh的分段線性表示[10]，Hung的分段聚合近似[11]，Lin的符號化表示[12]，Chan的離散小波[13]，Agrawa的離散傅里葉變換[14]及Fuchs的多項(xiàng)式回歸分析[15]等.這些方法的共同之處是雖然注意了時(shí)間序列的形態(tài)特征，但都是只側(cè)重于數(shù)值軸的信息，并沒有完全反映時(shí)間軸的信息.

擺在我們面前的一個(gè)問題就是如何引入能同時(shí)體現(xiàn)時(shí)間軸和數(shù)值軸信息的新度量方法，讓時(shí)間序列動態(tài)特性得到充分體現(xiàn)的同時(shí)，還能準(zhǔn)確度量距離，在時(shí)間效率上也能得到保障.目前這方面的研究較少，文獻(xiàn)[16]中首次引入基于線段面積的度量方式，但該方法對于由光滑曲線產(chǎn)生的時(shí)間序列，聚類效果較差.基于以上原因，本文提出基于三次樣條規(guī)則化的改進(jìn)的動態(tài)時(shí)間規(guī)整算法(SPDTW)，并將該算法與知識指導(dǎo)相結(jié)合，應(yīng)用于大規(guī)模不等長時(shí)間數(shù)據(jù)族的聚類問題中.全文共分為四部分：第一部分借助三次樣條規(guī)則化濾波，將原始時(shí)間序列用三次樣條近似表示，找出反映時(shí)間序列趨勢的關(guān)鍵點(diǎn)；第二部分定義任意兩時(shí)間序列之間的距離，這是對于大規(guī)模不等長時(shí)間數(shù)據(jù)聚類的關(guān)鍵.首先考慮任意兩曲線段的趨勢差異，利用函數(shù)的導(dǎo)數(shù)構(gòu)造曲線段間的距離，從而得到一個(gè)距離矩陣，利用動態(tài)時(shí)間規(guī)整算法在距離矩陣上尋找最優(yōu)路徑，得到最終的累積距離作為兩個(gè)不等長時(shí)間序列之間的距離；第三部分基于SPDTW和知識指導(dǎo)的聚類過程.將大規(guī)模不等長時(shí)間序列分成若干部分，通過SPDTW得到各部分的距離矩陣，利用最大樹法得到各部分子序列族的聚類結(jié)果，統(tǒng)計(jì)任意兩個(gè)時(shí)間序列出現(xiàn)在同一類中的次數(shù)，得到一個(gè)對稱矩陣，通過歸一化處理，得到相似矩陣，再次利用最大樹法得到最終聚類結(jié)果；最后一部分是以仿真實(shí)驗(yàn)來驗(yàn)證算法的有效性及時(shí)間效率.

2 三次樣條規(guī)則化趨勢濾波

目前有很多關(guān)于趨勢濾波的算法，例如Hodrick-Precott(H-P)濾波[17],l1趨勢濾波[18],l1規(guī)則化趨勢濾波[19]，l2規(guī)則化趨勢濾波[19]，移動平均濾波[20]等，這些方法在數(shù)據(jù)挖掘中起了重要的作用，但每種濾波方式都有一定局限性.例如：l1和H-P趨勢濾波對原始數(shù)據(jù)呈線性且比較平滑的情況，濾波后得到的趨勢與原數(shù)據(jù)相似性較好，但對于一些特殊異常值的處理不如謝邦昌等提出的l1和l2規(guī)則化趨勢濾波[19]效果好.但是如果原始數(shù)據(jù)對應(yīng)的函數(shù)是非常光滑的，謝提出的方法得到的趨勢整體誤差較大，且在連接點(diǎn)處的光滑性難以保證.指數(shù)濾波解決了光滑性的問題[21]，但在很多情況下，曲線并不是指數(shù)形的，而是具有波動性.Christian H. Reinsch分別在文獻(xiàn)[22]和文獻(xiàn)[23]中討論了三次和高次樣條函數(shù)的光滑性，為我們提供了一種濾波的算法，但他們均未對異常值進(jìn)行處理.受Christian H. Reinsch和Huber[24]的啟發(fā)，鑒于三次樣條比k(k>3)次樣條計(jì)算量小這一原因，本文給出三次樣條規(guī)則化趨勢濾波算法，該方法既能保證挖掘的數(shù)據(jù)關(guān)鍵點(diǎn)與原始數(shù)據(jù)在趨勢上有較高的吻合度，又可以解決最小二乘損失函數(shù)對異常值的干擾，可以得到穩(wěn)健的趨勢估計(jì).

假設(shè)時(shí)間序列y=(y1,…,yi,…,yn)包含一個(gè)基本趨勢x=(x1,…,xi,…,xn)和隨機(jī)成分z=(z1,…,zi,…,zn)，其中yi=y(ti),xi=x(ti),zi=z(ti)，y(t),x(t),z(t)都是時(shí)間t的函數(shù)，ti是時(shí)間采樣點(diǎn).本文假設(shè)時(shí)間序列y有有限個(gè)異常值，它是由幾乎處處光滑的函數(shù)生成的.不失一般性，不妨假設(shè)只有一個(gè)異常值，記為yi0=y(ti0)，且異常值不出現(xiàn)在兩個(gè)端點(diǎn)處，否則在眾多的數(shù)據(jù)點(diǎn)中刪除該異常值點(diǎn).若有多個(gè)異常值，則將全體數(shù)據(jù)分割成若干部分，使得每個(gè)部分包含一個(gè)異常值，然后按相同的方式進(jìn)行處理.

我們的目的是挖掘趨勢信息x，一方面希望x和y盡可能有相同的趨勢，另一方面希望隨機(jī)成分(殘差或者波動項(xiàng))盡可能小，同時(shí)減少異常值的干擾.提出如下優(yōu)化問題：

(1)

其中積分項(xiàng)反映的是趨勢的光滑性，第二項(xiàng)求和反映的是原始時(shí)間序列與基本趨勢之間的誤差，p是一個(gè)非負(fù)參數(shù)，用來控制光滑性和殘差在優(yōu)化問題中的權(quán)重.

下面分別對x(t)和HM(u)進(jìn)行說明：

1)y1=x1,yn=xn.

(2)

2) 對于任意的時(shí)間點(diǎn)ti(i=2,…,n-1)，

x(k)(ti)-=x(k)(ti)+,k=0,1,2，

(3)

(4)

4) 規(guī)定x(4)(t)=0.

從以上假設(shè)可知x(t)是[t1,tn]上的三次樣條函數(shù)[22].

(5)

6)M是一個(gè)給定的參數(shù)，用來控制數(shù)據(jù)中異常值對趨勢估計(jì)的影響，它描述了損失函數(shù)HM(u)從二次函數(shù)變?yōu)榫€性函數(shù)的過渡點(diǎn).Hubert建議當(dāng)殘差服從正態(tài)分布時(shí)，M取1.345可以在保證95%統(tǒng)計(jì)效率的前提下，得到盡量多的穩(wěn)健性[24].

7)δyi>0(i=1,…,n)是殘差的分布參數(shù)，通常可取原始時(shí)間序列的標(biāo)準(zhǔn)方差.

(6)

由于x(t)是三次樣條函數(shù)，設(shè)第i段的三次多項(xiàng)式的形式為

xi(t)=ai+bi(t-ti)+ci(t-ti)2+di(t-ti)3(i=1,…,n-1).

(7)

(8)

根據(jù)(2)、(4)及(8)式，I可寫成如下形式：

(9)

令H=H4H3-H2H1，根據(jù)上述四式容易計(jì)算J≈xAxT.將J≈xAxT及(9)式代入(1)式，得到一個(gè) (1)式近似的優(yōu)化問題：

(10)

(11)

(12)

由最小二乘法知(12)式目標(biāo)函數(shù)中的最優(yōu)解形式為

(13)

從(13)式可看出，當(dāng)p→時(shí)，x→y.通過解優(yōu)化問題，可以得到原始時(shí)間序列的關(guān)鍵點(diǎn)xi(i=1,…,n),利用Hermite插值或者三彎矩法可以得到三次樣條函數(shù)xi(t)，(t∈[ti,ti+1] ，i=1,…,n-1).將xi(t)的4個(gè)系數(shù)(ai,bi,ci,di)作為列向量存儲在矩陣B4×(n-1)中.

3 兩時(shí)間序列之間的距離

3.1曲線段間的距離度量設(shè)f(t)和g(t)是區(qū)間[t′,t″]上的兩條函數(shù).現(xiàn)給出它們間的距離度量規(guī)則.

1)f(t)和g(t)信息完整且趨勢相同(如圖1所示).

圖1 f(t)和g(t)信息完整且趨勢相同

圖2 f(t)和g(t)信息完整但趨勢不同

2)f(t)和g(t)信息完整但趨勢不同(如圖2所示).

不存在常數(shù)c，使得f(t)=g(t)+c，此時(shí)各點(diǎn)之間的趨勢差異仍可用變化率來表示，因此對于這種情況定義

(14)

可見1)是2)的特殊情況.若f(t)和g(t)是三次多項(xiàng)式，即

f(t)=af+bf(t-t')+cf(t-t')2+df(t-t')3，

g(t)=ag+bg(t-t')+cg(t-t′)2+dg(t-t′)3，

(15)

3)f(t)和g(t)信息不完整且趨勢不同(如圖3-圖6所示).

圖3 f(t)信息完整，g(t)后續(xù)信息不完整且趨勢不同

圖4 f(t)信息完整，g(t)前部信息不完整且趨勢不同

圖5 f(t)信息完整，g(t)前后信息均不完整且趨勢不同

圖6 f(t)和g(t)信息均不完整且趨勢不同

圖3-圖5以及兩條曲線段不在同一時(shí)間段的情形均可利用圖6的處理方式來定義距離.因此本文主要考慮圖6的定義規(guī)則.對于區(qū)間[s',s'']中的趨勢差異仍可用(14)的公式進(jìn)行度量，即

(16)

對于[t',s']∪[s″,t″]這部分，只有一個(gè)函數(shù)的信息可用.在這部分，兩個(gè)函數(shù)之間的差異既有時(shí)間和數(shù)值上的差異，又有趨勢變化的差異.綜合考慮，可用這三種信息的加權(quán)平均來反映兩條函數(shù)的差異.定義

(17)

來表示信息不全的部分的距離，其中αi(i=1,2,3,4)是各項(xiàng)的權(quán)重.到此為止，任意兩個(gè)曲線段之間的距離度量公式可表示為

d(f,g)=dcom(f,g)+dincom(f,g).

(18)

3.2基于三次樣條規(guī)則化的動態(tài)時(shí)間規(guī)整算法(SPDTW)

步驟1預(yù)處理(Spline規(guī)則化趨勢濾波).

步驟2求三次樣條函數(shù).

利用所獲關(guān)鍵點(diǎn)求出三次樣條函數(shù)，用分段三次多項(xiàng)式近似表示原始時(shí)間序列.

步驟3計(jì)算任意兩個(gè)曲線段間的距離.

利用(17)式計(jì)算任意時(shí)間段里兩曲線段之間的距離，得到一個(gè)K1×K2距離矩陣.

步驟4計(jì)算原始時(shí)間序列X和Y之間的距離.

運(yùn)用動態(tài)時(shí)間規(guī)整算法,在距離矩陣上尋找最優(yōu)路徑

W=c1,…,ck,…,cK，max(K1,K2)≤K

k表示路徑的長度，并得到最終的累積距離作為不等長時(shí)間序列X和Y之間的距離.

動態(tài)時(shí)間規(guī)整算法搜索的約束條件為:

1) 端點(diǎn)約束：搜索路徑的起始點(diǎn)和終止點(diǎn)正好位于兩個(gè)時(shí)間序列的兩個(gè)端點(diǎn)，即

c1=(1,1)，ck=(K1,K2)；

2) 局部連續(xù)約束：路徑W中的點(diǎn)在時(shí)間上不能后退，且每一個(gè)彎折步必須彼此相鄰(包括對角線元素)，即ck=(a,b),ck-1=(a',b') 0≤a-a'≤1,0≤b-b'≤1.動態(tài)時(shí)間規(guī)整采用遞推和迭代的方法發(fā)現(xiàn)最短距離

(19)

下面通過動態(tài)規(guī)劃的方法找到最短路徑，即

(20)

最終時(shí)間序列彎曲路徑最小累加值W(m,n)為時(shí)間序列X和Y的動態(tài)彎曲距離.

4 基于SPDTW和知識指導(dǎo)的大規(guī)模不等長時(shí)序族的聚類

4.1問題描述N條不等長的時(shí)間序列X1,…,XN組成時(shí)間序列族，時(shí)間序列Xi可表示為

4.2聚類過程對每個(gè)子序列族Partq，應(yīng)用SPDTW算法，得到一個(gè)Nq×Nq對稱距離矩陣

統(tǒng)計(jì)任意Xi和Xj(i,j=1,…,N)在每個(gè)Partq中出現(xiàn)在同一類的次數(shù)并求和，記

其中,

從而得到一個(gè)N×N的對稱矩陣R'，其主對角線上的元素均為N.對R'進(jìn)行歸一化處理，令

得到一個(gè)N×N的相似矩陣R，R中對角線元素全為1且R(i,j)∈[0,1].再次利用Prim方法構(gòu)造最大樹，即可得到不同閾值下的最終的清晰化的聚類結(jié)果.

5 仿真實(shí)驗(yàn)

5.1數(shù)據(jù)來源首先利用函數(shù)35sin(|t|)、25sin(t)、25cos(t)分別在[-400,300]、[-600,380]、[-200,500]以Δt=0.4為間隔生成三支不等長的時(shí)間序列X、Y、Z，長度分別為nX=1 751，nY=2 451，nZ=1 751.然后給出三組人造時(shí)序數(shù)據(jù)族，構(gòu)造方式如下：

經(jīng)過如此構(gòu)造，得到由X、Y、Z及人造時(shí)間序列XN、YN、ZN構(gòu)成的一族聚類對象.為方便起見，分別用{X1,X2,X3,X4,X5}，{X6,X7,X8,X9,X10}和{X11,X12,X13,X14,X15}標(biāo)記{XN1,XN2,NX3,NX4,X}，{YN1,YN2,YX3,YX4,Y}和{ZN1,ZN2,ZX3,ZX4,Z}.

5.2SPDTW和IDTW的比較為了更好展現(xiàn)SPDTW和IDTW的差別，先將IDTW的主要方法和步驟列出來.

1) 關(guān)鍵點(diǎn)的選取(l1趨勢濾波).

利用如下優(yōu)化問題

得到時(shí)間序列y=(y1,…,yi,…,yn)的基本趨勢x=(x1,…,xi,…,xn).

2) 利用分段直線近似表示原始時(shí)間序列.

4) 利用動態(tài)時(shí)間規(guī)整算法計(jì)算X和Y的最終距離.

5) 聚類過程.

由于每個(gè)Partq中所含數(shù)據(jù)的條數(shù)不相同，因此各個(gè)Partq對應(yīng)的相似矩陣Rq的階數(shù)也不相同，為了將所有數(shù)據(jù)最終用一個(gè)矩陣表示它們的相似關(guān)系，IDTW先找出一個(gè)階數(shù)最大的矩陣(為了方便敘述，假設(shè)最后一段階數(shù)最大)，然后利用

依次從后向前將每個(gè)Rq都補(bǔ)成與RQ相同階數(shù)的矩陣，再利用

得到所有數(shù)據(jù)族的相似矩陣，最后利用傳遞閉包進(jìn)行聚類.

下面通過仿真實(shí)驗(yàn)說明二者的差別.對于上述時(shí)間序列族，分別用本文給出的方法與基于IDTW和知識指導(dǎo)的聚類算法就聚類結(jié)果及時(shí)間耗費(fèi)進(jìn)行對比，對比結(jié)果分別列于表1和表2中.為了書寫方便，表1中數(shù)字1～15分別代表數(shù)據(jù)族中的時(shí)間序列X1～X15.

表1 SPDTW和IDTW聚類結(jié)果

表1結(jié)果表明：針對不等長時(shí)間序列族聚類問題，基于SPDTW和知識指導(dǎo)的聚類算法得到了預(yù)想的結(jié)果.5,10,15是數(shù)據(jù)族中的基本數(shù)列，不存在隨機(jī)性，而{1,2,34}，{6,7,8,9}和{11,12,13,14}分別是5,10,15通過加入隨機(jī)數(shù)并平移得到的，因此在趨勢上具有很大的相似性，因此在λ=0.994 5時(shí)，就相應(yīng)地聚為一類.而IDTW的聚類結(jié)果并不理想，在λ較小時(shí)，也沒能將趨勢相同的數(shù)據(jù)聚為一類.這是因?yàn)镮DTW雖然考慮了時(shí)間軸和數(shù)值軸的分布，但它是用折線段近似代替原始時(shí)間序列，這樣不能保證光滑性，且誤差大.此外在聚類過程中，文獻(xiàn)[16]中用最后一部分PartN的相似矩陣的元素?cái)U(kuò)充PartN-1對應(yīng)的相似矩陣中缺少的元素，按照這方式一直向前推，直到推到Part1對應(yīng)的相似矩陣.這種處理方式，主要的思想是用近期的數(shù)據(jù)推算之前的信息，作為短期的推測，得到的效果比較好，但若時(shí)間段很長，得到的趨勢差別就會很大.

本文沒有將各個(gè)部分的相似矩陣擴(kuò)充成階數(shù)相同的矩陣，而是直接在每個(gè)部分聚類，統(tǒng)計(jì)出任意兩個(gè)時(shí)間序列出現(xiàn)在同一類的次數(shù)，最后綜合考慮聚類結(jié)果，從而避免了對未出現(xiàn)的時(shí)間序列補(bǔ)充過多的不確定的信息，提高了聚類結(jié)果的準(zhǔn)確性.

5.3效率對比對于大規(guī)模時(shí)間序列的聚類算法，時(shí)間效率是一項(xiàng)重要的測試指標(biāo).基于SPDTW和知識指導(dǎo)的聚類算法包括三大部分：第一部分是每個(gè)Partq中三次樣條函數(shù)規(guī)則化趨勢濾波過程；第二部分是SPDTW距離度量過程；第三部分是知識指導(dǎo)聚類過程(統(tǒng)計(jì)頻數(shù)，生成最大樹，聚類).我們將基于SPDTW和知識指導(dǎo)的聚類算法及基于IDTW和知識指導(dǎo)的聚類算法所用時(shí)間進(jìn)行對比，結(jié)果見表2.“s”表示時(shí)間，單位：s.

表2 SPDTW和IDTW聚類時(shí)間消耗對比

表2表明：與基于IDTW和知識指導(dǎo)的聚類算法相比，基于SPDTW和知識指導(dǎo)的聚類算法在時(shí)間消耗上沒有太大差別，IDTW方法在數(shù)據(jù)擬合時(shí)，時(shí)間要優(yōu)于SPDTW，這是因?yàn)楹笳呱婕叭味囗?xiàng)式的四個(gè)系數(shù)存儲及運(yùn)算，而前者在計(jì)算距離時(shí)只用到關(guān)鍵點(diǎn)的坐標(biāo)，這在第一步濾波時(shí)就已得到.但在其他方面時(shí)間相差不多.在動態(tài)尋優(yōu)過程中SPDTW還略快于IDTW.計(jì)算整個(gè)過程的總時(shí)間，二者幾乎相等.

6 結(jié) 論

在文獻(xiàn)[16]中已經(jīng)驗(yàn)證IDTW與傳統(tǒng)的DTW相比，無論在聚類結(jié)果還是時(shí)間效率上都有很大的優(yōu)勢，而SPDTW對于光滑曲線產(chǎn)生的數(shù)據(jù)在聚類上優(yōu)于IDTW.由此看見本文提出的方法對于解決大規(guī)模不等長時(shí)間序列的聚類問題有一定的作用.但從表2中也可以看出，時(shí)間消耗主要集中在DWT尋找最短路徑上，因此本文提出的方法還有改進(jìn)的余地.主要可從以下兩方面考慮，一是改進(jìn)路徑的尋優(yōu)過程；二是給出兩時(shí)間序列動態(tài)彎曲距離的其他形式的計(jì)算公式.利用DTW尋找最短路徑的原因是兩個(gè)不等長的時(shí)間序列無法逐段比較趨勢變化，所以可從這點(diǎn)出發(fā)，將不等長時(shí)間序列改造成能夠逐段比較的兩個(gè)等價(jià)的對象.