李國榮，冶繼民，甄遠(yuǎn)婷

（西安電子科技大學(xué)數(shù)學(xué)與統(tǒng)計學(xué)院，西安 710126）

（*通信作者電子郵箱jmye@mail.xidian.edu.cn）

0 引言

隨著計算機(jī)技術(shù)的發(fā)展，時間序列數(shù)據(jù)作為一種重要的數(shù)據(jù)形式存在于生活各個領(lǐng)域，例如：醫(yī)學(xué)、金融股票、天氣數(shù)據(jù)、圖像分割等［1］。面對這些龐大的數(shù)據(jù)集，如何對其進(jìn)行有效處理，挖掘出有效信息已經(jīng)成為了當(dāng)今研究的熱點話題。時間序列聚類是其中最典型的一種方法，其目的是將未標(biāo)記的時間序列數(shù)據(jù)劃分到不同的組，并確保組內(nèi)序列的相似程度最大、組間序列相似程度最小，從而識別出序列集合的結(jié)構(gòu)，提取需要的重要信息［2］。由于在收集這些數(shù)據(jù)的過程中，會遇到很多不穩(wěn)定因素，導(dǎo)致數(shù)據(jù)集中存在一些異常值點。這些異常值會直接影響到最終的聚類結(jié)果，從而影響后續(xù)對數(shù)據(jù)進(jìn)行分析、預(yù)測的結(jié)果，最終可能做出錯誤的決策。因此，采取適當(dāng)?shù)姆椒ㄏ惓Ｖ祵垲惤Y(jié)果的影響非常有必要。

目前已有一些處理異常值的方法，例如：彭柳青等［3］將樣本加權(quán)思想與子空間聚類算法相結(jié)合，為每一個樣本分配反映離群程度的尺度參數(shù)，提出了一種魯棒的子空間聚類算法；曹晨曦等［4］將統(tǒng)計的方法應(yīng)用到時間序列數(shù)據(jù)中，提出了一種基于統(tǒng)計方法的異常值檢測與校正方法；Grané 等［5］提出了一種基于小波的通用的異常值檢測與校正的方法；Lafuente-Rego 等［6］通過分別考慮度量、噪聲、分段方法，并基于比較樣本分位數(shù)自協(xié)方差，提出了三種模糊C-mesoids 模型，三種方法從不同的角度對異常值進(jìn)行了處理。但是這些工作基本都需要對數(shù)據(jù)進(jìn)行處理，而處理過程中涉及的特征提取、參數(shù)選擇等工作都很容易造成數(shù)據(jù)中有用信息的丟失，從而導(dǎo)致聚類結(jié)果不準(zhǔn)確；另外在對異常值進(jìn)行修正時也可能會產(chǎn)生較大誤差，影響最終的聚類結(jié)果。

為了避免上述方法的缺點，本文通過增強(qiáng)算法過程中使用的相似性度量的魯棒性，提出了直接基于原數(shù)據(jù)的魯棒時間序列聚類方法。傳統(tǒng)的衡量時間序列間相似性的度量，如：歐氏距離、馬氏距離、閔可夫斯基距離、Pearson 相關(guān)系數(shù)、動態(tài)時間規(guī)整（Dynamic Time Warping，DTW）距離等，除DTW 距離外，其余的相似性度量一般都沒有很強(qiáng)的魯棒性，當(dāng)時間序列數(shù)據(jù)中存在異常值點時這些度量不再能準(zhǔn)確衡量出序列間的相關(guān)關(guān)系。DTW 距離雖然有較強(qiáng)的魯棒性，但其計算復(fù)雜度非常高，導(dǎo)致聚類算法效率降低。本文在文獻(xiàn)［7］中提出的廣義互相關(guān)（Generalized Cross-Correlation，GCC）度量的基礎(chǔ)上進(jìn)行改進(jìn)，提出了一種魯棒廣義互相關(guān)度量。該度量是非參的且不需要對數(shù)據(jù)進(jìn)行任何處理過程，在一定程度上可以減少對數(shù)據(jù)進(jìn)行處理與修正過程帶來的誤差。

1 廣義互相關(guān)度量

GCC是衡量兩個平穩(wěn)時間序列之間線性相關(guān)性的一般度量［7］。設(shè)xt、yt為兩個平穩(wěn)時間序列，Ryx，k是(yt，yt-1，…，yt-k，xt，xt-1，…，xt-k)的協(xié)方差矩陣（k為時間滯后），其對角線部分為Ryy，k和Rxx，k，上三角部分為，下三角部分為Cxy，k。這里的Ryy，k為(yt，yt-1，…，yt-k)的協(xié)方差矩陣、Rxx，k為(xt，xt-1，…，xt-k)的協(xié)方差矩陣、Cxy，k是由兩個時間序列間的互相關(guān)系數(shù)構(gòu)成的矩陣，并且有|Ryx，k|=|Rxx，k||Ryy，k-則兩序列間的廣義互相關(guān)度量為：

2 魯棒廣義互相關(guān)度量

上述的協(xié)方差矩陣Ryx，k中的元素是利用的Pearson 相關(guān)系數(shù)計算出來的，Pearson 相關(guān)系數(shù)對異常值點敏感，從而導(dǎo)致利用它構(gòu)造出來的序列間的相似性度量GCC 也對異常值敏感，而生活中大多數(shù)的時間序列數(shù)據(jù)中都存在異常值點，所以選擇適當(dāng)?shù)姆椒ㄔ鰪?qiáng)聚類算法中使用的相似性度量的魯棒性很有必要。

一種可行的提高GCC 魯棒性的方法就是采用隨機(jī)變量之間的魯棒相關(guān)系數(shù)代替Pearson 相關(guān)系數(shù)構(gòu)造GCC。該方法從另一個角度即相關(guān)系數(shù)的魯棒估計來增加時間序列間相似性度量的魯棒性，可以解決現(xiàn)有的時間序列聚類算法不適用于存在異常值的數(shù)據(jù)集的問題，而現(xiàn)有的方法中很少有從此角度進(jìn)行改進(jìn)的。

文獻(xiàn)［8-9］中提出了一種基于Qn統(tǒng)計量的魯棒相關(guān)系數(shù)，Qn統(tǒng)計量的表達(dá)式為：

大量的仿真實驗發(fā)現(xiàn)，原始的Qn統(tǒng)計量存在一定的小樣本偏差。小樣本偏差可以通過給Qn統(tǒng)計量乘一個小樣本校正因子［11］校正，即校正后的Qn統(tǒng)計量為（為方便，仍采用原記號）：

其中調(diào)節(jié)因子dn的具體取值如下。

1）當(dāng)n≤9時，調(diào)節(jié)因子的取值如表1所示。

表1 調(diào)節(jié)因子dn的取值Tab.1 Value of adjustment factor dn

2）當(dāng)n＞9時：

注：當(dāng)樣本量n≥10 000 時，Qn將不需要再乘小樣本校正因子。

基于校正后的Qn統(tǒng)計量構(gòu)造魯棒相關(guān)系數(shù)

式中：med(x)為隨機(jī)變量x的中位數(shù)；med(y)為隨機(jī)變量y的中位數(shù)。

計算向量Zt，2(k+1)的協(xié)方差矩陣：

得到兩個時間序列間的魯棒廣義互相關(guān)（Robust Generalized Cross-Correlation，RGCC）度量：

采用AIC 準(zhǔn)則（Akaike Information Criterion）和BIC 準(zhǔn)則（Bayesian Information Criterion）選取時間滯后k的值［7］。自回歸模型中，k一般取1 或2；動態(tài)因子模型中，k一般在1～3 中取值。

另外，可以驗證RGCC(xt，yt)也滿足距離度量的基本性質(zhì)，即：

1）對稱性：RGCC(xt，yt)=RGCC(yt，xt)；

2）0 ≤RGCC(xt，yt) ≤1；

3）當(dāng)且僅當(dāng)序列之間有精確的線性關(guān)系時，RGCC(xt，yt)=1；

4）當(dāng)且僅當(dāng)所有的互相關(guān)系數(shù)為零時，RGCC(xt，yt)=0。

證明 1）、2）易見，這里僅給出3）、4）的證明。

對于3），已知當(dāng)且僅當(dāng)存在一個線性組合aTZt=0時，有|=0，即這兩個序列是線性相關(guān)的［7］。而RGCC(xt，yt)=1則一定有|=0，所以性質(zhì)（3）成立。

下面證4），首先對RGCC(xt，yt)變形

根據(jù)Hadamarsd’s 不等式，右邊的行列式小于等于1，且等式成立等價于矩陣為對角矩陣，即=0xy，k。 證畢。

3 仿真實驗

這一部分本文用兩個常見的時間序列模型生成仿真數(shù)據(jù)集，給模型生成的數(shù)據(jù)集中隨機(jī)加入異常值點，分別用GCC與RGCC 計算其距離矩陣，作為層次聚類算法的輸入對數(shù)據(jù)進(jìn)行聚類，并對聚類的結(jié)果進(jìn)行分析。

本實驗用3 個常見的聚類評價指標(biāo)比較基于兩種度量聚類得到的聚類結(jié)果。假設(shè)U 是真實的聚類結(jié)果，V 是某次實驗得到的聚類結(jié)果。

1）調(diào)整的蘭德系數(shù)（Adjusted Rand Index，ARI）［13］，ARI∈[ -1，1]且ARI 的值越大代表聚類結(jié)果越好。具體計算表達(dá)式為：

其中RI（Rand Index）為蘭德系數(shù)，計算公式為：

這里的a等于在U 中屬于同一組，在V 中也屬于同一組的對象的對數(shù)；b等于在U 中屬于同一組，在V 中屬于不同組的對象的對數(shù)；c等于在U 中不屬于同一組，在V 中屬于同一組的對象的對數(shù)；d等于在U 中不屬于同一組，在V 中也不屬于同一組的對象的對數(shù)。E(RI)是蘭德系數(shù)的期望。

2）F-measure，F(xiàn)-measure 是精確率P與召回率R的調(diào)和平均［14］，F(xiàn)∈[0，1]，且F的值越大代表聚類結(jié)果越好。具體的計算表達(dá)式為：

3）Jaccard 系數(shù)（Jaccard Coefficient，JC），JC 是用來比較兩個有限集合間的相似性與差異性的度量，JC∈[0，1]且JC 的值越大，集合間的相似性越高，代表聚類的結(jié)果越好。具體的計算表達(dá)式為：

用F-measure 與Jaccard 系數(shù)評價聚類的結(jié)果時，可以得到與ARI指標(biāo)一樣的結(jié)論，但由于ARI指標(biāo)的變化更加直觀，所以本文實驗結(jié)果中僅列出ARI 指標(biāo)的值。另外，本文各表中給出的實驗結(jié)果均為實驗100次取平均的結(jié)果。

本文的實驗環(huán)境為Windows 7 系統(tǒng)下8 GB 內(nèi)存的Intel 3.40 GHz PC，Matlab2017a軟件。

3.1 自回歸模型生成的仿真數(shù)據(jù)

用簡單的AR（1）模型（AutorRegressive model）生成一組線性、平穩(wěn)的時間序列數(shù)據(jù)：

i=1，2，…，5 時，ai=0.9；i=6，7，…，15 時，ai=0.2。通過對誤差項εi的限制引入時間序列間的交叉依賴關(guān)系，r(i，j)=E(εi，t，εj，t)。

情形1

情形2

情形3

情形4

情形5

從原始AR（1）模型看，該模型生成的時間序列數(shù)據(jù)存在明顯的兩大類，即前5個序列相關(guān)，后10個序列相關(guān)。但當(dāng)引入交叉依賴項εi后，序列間的相關(guān)關(guān)系發(fā)生了變化。情形1、3中顯然有6大類，即前10個序列相關(guān)，后5個序列獨立；情形2、4、5中顯然有兩大類，即前10個序列相關(guān)，后5個序列相關(guān)。

表2 為數(shù)據(jù)不存在異常值的情況下，T=100、200（T為每個時間序列的長度）時，GCC(xt，yt)與RGCC(xt，yt)的實驗結(jié)果對比。

表2 無異常值時，GCC、RGCC性能對比Tab.2 Performance comparison of GCC and RGCC without outliers

不存在異常值的情況下，從表2 可以看出兩個相似性度量得到的聚類結(jié)果沒有明顯的差異。T=100時，單鏈情況下GCC 的結(jié)果更精確一些，全鏈情況下RGCC 的結(jié)果更精確一些，均鏈接的情況下兩個度量的結(jié)果幾乎一致；T=200時，基于兩個相似性度量得到的結(jié)果均更加精確，ARI 指標(biāo)的值相差不大。

給模型在五種場景下生成的仿真數(shù)據(jù)中隨機(jī)加入一些異常值點，再分別利用GCC 與RGCC 對存在異常值點的數(shù)據(jù)進(jìn)行聚類。仿真實驗結(jié)果見表3。圖1 為在情形1 中，使用單鏈接方式用兩種度量進(jìn)行聚類的層次聚類圖。

表3 存在異常值時，GCC、RGCC性能對比Tab.3 Performance comparison of GCC and RGCC with outliers

從表3 可以看出，在時間序列數(shù)據(jù)中存在異常值時，三種鏈接方式下的五種情形中，RGCC 得到的聚類結(jié)果均比GCC的聚類結(jié)果更精確。

情形1 中，理想的聚類結(jié)果應(yīng)該是第1～10 個時間序列為一組，第11、12、13、14、15個序列各自單獨成一組，從圖1（a）、圖1（b）可以看到，由于異常值點的干擾，GCC 不能識別出時間序列間的某些相關(guān)關(guān)系，導(dǎo)致第9、10、12、14 個時間序列被錯分。

圖1 同一數(shù)據(jù)集上，GCC、RGCC進(jìn)行聚類得到的樹形圖對比Fig.1 Comparison of dendrogram obtained by clustering with GCC or RGCC on the same dataset

3.2 動態(tài)因子模型生成的仿真數(shù)據(jù)集

這部分本文使用文獻(xiàn)［15］中提出的三種場景，即帶有組因子結(jié)構(gòu)的動態(tài)因子模型（Dynamic Factor Model，DFM）作為數(shù)據(jù)生成的過程?？紤]三個組，每個組有3 個特定的組因子r1=r2=r3=3 和相同的數(shù)量的序列N1=N2=N3=100。DFM可以表示為：

這里的N=300，T=100。G={g1，g2，g3}表示每組中的成員，Nj(j=1，2，3)表示每組中序列的個數(shù)，xit是一個80×1 的可觀測變量，其每個元素均是由在區(qū)間[-2，2]上的均勻分布產(chǎn)生，β=(1，2，3，0，…，0)1×80T。是一個rj× 1維的不可觀測的組特異因子向量，其每個元素由均值為j、方差為1 的高斯分布生成的。的每個元素是由高斯分布N(0，j)生成。通過對εi，t的性質(zhì)做不同的假設(shè)，產(chǎn)生三種情形。

情形1 誤差項εi，t服從標(biāo)準(zhǔn)高斯分布。

情形2 εi，t=其中，如果t 為奇數(shù)，則δt=1；否則，δt=0；并且 ε1=與 ε2=服從均值向量為0、協(xié)方差矩陣為Σ=(σi，j)，σi，j=0.3|i-j| 的多元高斯分布。

情形3 εi，t=0.2εi，t-1+ ei，j，這里et=(e1，t，e2，t，…，eN，t)′服從均值向量為0、協(xié)方差矩陣為Σ=(σi，j)，σi，j=0.3|i- j|的多元高斯分布。

給模型生成的三種情形的數(shù)據(jù)中隨機(jī)加入一些異常值點，分別用GCC與RGCC計算距離矩陣，作為層次聚類算法的輸入進(jìn)行聚類。表4 為在該模型下的兩種度量的聚類結(jié)果對比。觀察表4 中得到的ARI 指標(biāo)的值很直觀可以看到，對動態(tài)因子模型生成的帶有異常值點的時間序列數(shù)據(jù)進(jìn)行聚類，RGCC仍可以得到比GCC更精確的聚類結(jié)果。

表4 動態(tài)因子模型中，GCC、RGCC性能對比Tab.4 Performance comparison of GCC and RGCC in dynamic factor model

通過對兩個仿真數(shù)據(jù)集的聚類結(jié)果分析可知，本文提出的序列間的相似性度量RGCC 不僅在數(shù)據(jù)無異常值時，可以發(fā)現(xiàn)序列間的相關(guān)關(guān)系，而且在數(shù)據(jù)中存在異常值點時仍有效。

4 結(jié)語

本文提出了一種基于相關(guān)系數(shù)魯棒估計的時間序列間的魯棒相似性度量——RGCC。仿真實驗結(jié)果顯示，針對存在異常值點的時間序列數(shù)據(jù)，基于RGCC 進(jìn)行聚類得到的結(jié)果，明顯比基于GCC 進(jìn)行聚類得到的結(jié)果更接近真實的聚類結(jié)果。該度量具有很強(qiáng)的魯棒性，更加適用于現(xiàn)實生活中的時間序列數(shù)據(jù)集。此外，基于RGCC 對存在異常值點的時間序列數(shù)據(jù)聚類時，不需要異常值檢測與處理的過程，可以有效避免查找異常值點過多或者過少，或者對異常值點進(jìn)行修正時產(chǎn)生的誤差；同時，該度量的計算不需要提前設(shè)置任何參數(shù)，避免了參數(shù)選取不當(dāng)對聚類結(jié)果的影響。

由RGCC 的結(jié)構(gòu)可知，在小樣本數(shù)據(jù)情況下，選取不同的小樣本校正因子會導(dǎo)致相關(guān)矩陣的變化，從而影響到距離矩陣的計算，最終造成聚類結(jié)果的不準(zhǔn)確。因此，未來的工作可以從優(yōu)化小樣本校正因子選取的角度對本文提出的度量進(jìn)行改進(jìn)。