999精品在线视频,手机成人午夜在线视频,久久不卡国产精品无码,中日无码在线观看,成人av手机在线观看,日韩精品亚洲一区中文字幕,亚洲av无码人妻,四虎国产在线观看 ?

一種時間序列鑒別性特征字典構建算法?

2020-01-02 03:45:46
軟件學報 2020年10期
關鍵詞:單詞分類特征

(北京交通大學 計算機與信息技術學院,北京 100044)

時間序列是一系列按時間進行排序的實值數據組成的集合.在許多研究領域或實際應用領域之中存在著大量的時間序列數據,例如惡意軟件檢測、風能預測、工業自動化、電壓穩定評估、移動設備追蹤等領域[1-3].目前,獲得的時間序列通常具有如下兩個特點:時間序列數據集的規模很大,同時,每條時間序列數據的維度都很高.

時間序列分類(time series classification,簡稱TSC)技術的研究涉及到許多方面的技術問題,這些問題可能包括時間序列特征的發現或生成以及如何存儲或壓縮時間序列等.Esling 等人對時間序列的研究領域給出了詳細介紹[4].Bagnall 等人對各種時間序列分類算法進行了詳細分析[5].一般TSC 算法可大致分為兩類:基于完整時間序列的方法和基于局部特征的方法[6].前者基于全局相似性進行分類預測,主要研究相似性的不同度量方式和使用方法[7-11];而后者基于時間序列的局部特征進行分類預測,通常利用不同的特征生成和轉換技術來發現局部特征,然后基于建立的特征集合直接構建分類模型或對時間序列數據進行轉換[12-15].然而,目前多數的TSC算法無法以足夠的精度和速度充分處理大規模的數據.一些精度較高的分類算法,例如,基于局部鑒別性特征Shapelet 的轉換算法的時間復雜度為O(N2×n4)[16],其中,N為時間序列數據集中的實例數,n為時間序列的長度.基于時間序列轉換的集成分類算法(the collective of transformation-based ensembles,簡稱COTE)[17]的分類精度顯著地優于常見的時間序列分類算法,但該算法中包含多個時間復雜度非常高的基分類器.

為了提高TSC 算法的分類效率,時間序列的快速轉換表示方法成為當前的一個研究熱點.本文對基于特征包(bag-of-pattern,簡稱BOP)的分類模型進行研究,這類模型具有分類精度高和運行速度快的特點.BOP 模型將時間序列分成一系列子結構,并將這些子結構作為離散化特征構建特征字典,最后將基于特征頻數的向量作為模型訓練和分類的基礎.最早的特征包模型(bag-of-patterns using symbolic aggregate approximation,簡稱BOPSAX)[18]使用固定長度的滑動窗口和符號聚合近似技術(SAX)將每個窗口中的子序列轉換成離散化特征,然后使用基于頻數向量的歐幾里德距離作為相似性度量方式,最后通過1NN 分類器進行分類預測.Senin 等人利用SAX 技術和向量空間模型構建了一種新的時間序列分類模型(symbolic aggregate approximation and vector space model,簡稱SAX-VSM)[19],該模型基于詞頻-逆向文檔頻率(term frequency-inverse document frequency,簡稱tf-idf)對符號化特征進行加權,同時使用余弦距離代替歐式距離進行相似性度量;此外,它只為每個類構建一個特征向量,而不是每個樣本一個向量,這極大地縮短了模型的運行時間.Neuyen 等人[14]對基于SAX 的轉換方法的可變長度的單詞進行了研究,并將序列學習算法用于轉換后的時間序列分類問題.SAX 技術本質上仍然是在時域空間對時間序列進行處理,然而,實際中的一些問題在時域空間進行處理時會比較困難,而在頻域空間更容易取得好的處理效果.例如,通常可在頻域空間將代表噪聲的頻率成分去除來實現降噪.此外,頻域中包含一些時域難以發現的鑒別性信息.

離散傅里葉變換(discrete Fourier transform,簡稱DFT)[20,21]可將時域空間的時間序列轉換為一組頻域空間的不同頻率的正弦波.因此,基于離散傅里葉變換的時間序列符號化表示方法受到各國研究人員的重視.Schafer等人[22]提出用符號傅里葉近似(symbolic Fourier approximation,簡稱SFA)來代替SAX 對時間序列進行離散化表示.接著,Schafer 又提出了一種基于SFA 的特征包算法(bag-of-SFA-symbol,簡稱BOSS)[23],但該算法對時間序列離散化過程中只簡單挑選前l個傅里葉值,未考慮傅里葉值的鑒別性.為此,Schafer 等人進一步提出一種用于時間序列的單詞抽取方法(word extraction for time series classification,簡稱WEASEL)[24],該方法用鑒別性較強的top-l個傅里葉值代替前l個傅里葉值對時間序列進行符號化,但是該方法存在一個明顯的缺陷,其將所有長度子序列轉化得到的單詞設為定長,即,所有子序列經離散傅里葉變換后保留的傅里葉值的個數相同.然而,實際中不同周期的時間序列子序列所含有的鑒別性頻域信息量可能不同,單一的固定長度會導致損失大量鑒別性信息或包含冗余信息.為此,針對目前基于SFA 的時間序列離散化表示方法存在的問題,本文首先提出一種用于時間序列分類的基于SFA 的可變長度單詞抽取算法(variable-length word extraction algorithm,簡稱VLWEA).該算法為每個窗口長度學習性能最優的單詞長度.圖1 展示了本文提出的變長單詞抽取算法相對于定長單詞抽取算法所具有的優勢.圖中,w表示滑動窗口長度,lw表示長度為w的滑動窗口中提取的單詞長度,c1和c2表示類別.

圖1(b)中給出了長度為6 和8 的滑動窗口獲得的4 個子序列經過SFA 轉換得來的完整字符序列.圖1(a)所示為用定長單詞抽取算法得到的特征(假定共同單詞長度為4),其中,“6aabc”表示一個長度為6 的滑動窗口子序列抽取出的長為4 的字符序列,即特征.圖1(c)中給出了變長單詞抽取算法獲得的兩種長度滑動窗口對應的特征,變長單詞抽取算法旨在獲得具有更強鑒別性的子序列.從圖1 可以看出,定長單詞“8abab”無法區分長為8的兩類子序列,而變長單詞“8ababc”和“8ababd”可區分兩類子序列.此外,定長單詞“6aabc”和“6acbc”雖可區分長度為6 的兩類子序列,但與變長單詞“6aa”“6ac”相比,包含了更多的冗余成分.

與此同時,針對基于VLWEA 建立的特征字典中存在大量冗余特征的問題,本文基于tf-idf 的基本思想定義了一種新的特征鑒別性強弱度量方式來對鑒別性特征進行選擇,并能根據不同周期生成的特征的整體分類性能動態設定各周期生成特征的選擇閾值.本文的主要貢獻概括如下.

(1)本文提出了一種基于網格搜索的滑動窗口單詞長度學習算法.該算法為不同窗口長度學習分類性能最優的單詞長度,不同窗口長度可能生成不同長度的單詞.變長單詞可有效減少鑒別性信息的損失,在此基礎上可以獲得更優的分類效果.與此同時,我們分析了Bigrams 語法模型對基于變長單詞特征字典的模型分類性能的影響.

(2)本文基于tf-idf 的基本思想定義了一種新的特征鑒別性評價統計量用于特征選擇.

(3)針對采用固定閾值進行特征選擇忽略了不同周期轉換得來的特征存在性能差異的問題,我們提出一種根據不同滑動窗口生成特征的整體分類性能以動態地設定各窗口生成特征選擇閾值的機制,該機制可有效縮小時間序列特征空間的規模并提高模型的分類性能.

本文第1 節介紹基本概念和理論基礎.第2 節介紹本文使用的相關算法.第3 節給出實驗方式和實驗結果.第4 節總結全文.

1 定義與理論基礎

本節介紹時間序列的基本概念和基于SFA 的時間序列轉換表示方法的理論基礎.

1.1 時間序列

下面我們首先介紹本文用到的一些基本概念.

定義1(時間序列).時間序列T是由n個有序的實際觀測值t0,t1,…,tn-1組成的一個實值序列,即T={t0,t1,…,tn-1},其中,ti∈R.

用D={T0,T2,…,TN-1}表示包含N條時間序列的數據集合,其中,每條時間序列Ti={ti1,ti2,…,tin}.

定義2(時間序列子序列).給定一條完整的時間序列T={t0,t1,…,tn-1},該序列T中長為ω的窗口S之中包含的ω個連續值組成的序列稱為時間序列T的子序列,若其在時間序列T上的起始位置為a,則S(a,ω)=(ta,ta+1,…,ta+ω–1),其中,0≤a≤n-ω.

定義3(特征字典).特征字典是用于表示時間序列數據的特征集合.

本文構建的特征字典中的每個元素為一個指定長度滑動窗口生成的字符序列.表1 給出文中涉及的常用符號.

1.2 時間序列離散傅里葉變換

本節我們介紹時間序列的離散傅里葉變換過程[20,21].

給定一條由n個數值組成的離散序列T={t0,t1,…,tn-1},其離散傅里葉變換公式為

其中,j為虛數單位,即,W 為

經公式(1)可將離散數值序列T轉換為序列x,轉換過程可表示為

上式可簡單記作:

其中,F稱為離散傅里葉變換矩陣.

矩陣F中第i行表示第i個正弦波,xi表示時間序列在這個正弦波上的投影,即,時間序列T中包含的第i個頻率正弦波成分的多少.它反映了時間序列與該頻率正弦波的相關性.

數據x={x0,x1,…,xn-1}為經離散傅里葉變換得到的數據,稱為頻域向量.時間序列T可通過逆變換來恢復,即,

其中,FH表示正交矩陣F的共軛轉置.

定理1(Parseval theorem)[25].若x為序列T經離散傅里葉變換得到的序列,則有

上述定理說明序列T在時域空間的能量與在頻域空間的能量相同.

時間序列的每個元素可以表示為

復數xf通常稱作傅里葉系數(Fourier coefficient).

通過歐拉公式可以將公式(1)表示為

傅里葉系數xf可分為realf實數部分和imagf虛數部分:

經過傅里葉變換,我們得到的傅里葉系數值序列可以表示為

由于realn-f=realf,imagn–f=-imagf,imag0=0,所以上述2n維數組可用n維數組表示.

當n為偶數時,imagn/2=0,x可以表示為

當n為奇數時,x可以表示為

考慮到傅里葉值的對稱性,為了減少計算量和存儲空間,一般用n維數組對傅里葉值進行存儲.由于x中的傅里葉值代表了固定周期內的時間序列的頻率成分,因此,符號傅里葉的近似過程就是對這些傅里葉值進行選擇和符號化轉換.下面我們介紹基于SFA 的特征生成過程.

1.3 基于SFA的特征生成

基于SFA 的特征生成過程分為兩步:首先,從子序列轉換得到的傅里葉值數組中選擇鑒別性最強的top-l個傅里葉值;其次,利用分箱技術將選擇的傅里葉值依序轉換為符號組成單詞,即,特征.

傅里葉值的鑒別性度量本質上是根據一系列標定類屬性的實值來判斷這組值對應的頻率成分的鑒別性強弱.Schafer 等人[24]提出利用F統計量對傅里葉值的鑒別性強弱進行度量,然后選擇鑒別性最強的前l個傅里葉值.這樣處理比簡單挑選top-l個傅里葉值生成的特征更具有鑒別性.本文我們也使用F統計量進行傅里葉值鑒別性度量,與Schafer 所提方法不同的是,我們為每個窗口長度選擇分類性能最優的傅里葉值個數進行特征生成,而Schafer 等人使用的方法所有滑動窗口都選擇鑒別性top-l個傅里葉值.計算reali或imagj對應的實數集上的F統計量的類內方差MSW和類間方差MSB的計算公式[24]為

其中,Q為reali或imagj對應的傅里葉值集合,q為集合Q中實數的個數,k為集合Q中的類屬性取值個數,x為集合Q中的任意實數.Qc為集合Q中類屬性值為c的傅里葉數值集合,qc為集合Qc中實數的個數,為集合Q中所有實數的均值,為Qc中所有實數的均值.

F值的計算公式為

當F值用于傅里葉值選擇時,我們希望找到最大的F值,其等價于不同類均值之間的最大差,此時,傅里葉值的鑒別性最強.

例如,給定一個長為10 的子序列,假設我們選擇的前4 個傅里葉值為“0.11,0.23,0.02,0.63”,利用分箱技術我們可將該實值序列轉換成為一個字符序列,假設為“abcd”,由于不同長度的滑動窗口代表不同周期的頻率成分,它們生成的特征不同,為此,每個單詞都對應固定的滑動窗口長度,上述生成的單詞通常表示為“10abcd”.Schafer在文獻[23]中對SFA 的具體符號化過程進行了詳細介紹,這里不再贅述.

n0元語法模型(n0-gram)是自然語言處理中很重要的統計語言模型.該模型在實際應用中通常假設某個詞出現的概率只與它前面的一個或幾個詞有關,即,馬爾可夫假設.當n0=i時(i為正整數),稱為i元語法,也稱為i階馬爾可夫鏈,此時第j個詞出現的概率只與其前i-1 個詞有關.n0取1、2 和3 時,基于n0元語法表示獲得的詞組稱為一元單詞(unigram word)、二元單詞(bigram word)和三元單詞(trigram word).此外,若特征字典有m個特征,則基于n0元語法要考慮的詞組合可能有個.因此,為了縮小特征字典的規模,通常只考慮一元和二元詞組組成的特征.為了彌補所有周期抽取定長單詞所導致的鑒別性信息的損失,Schafer 等人[24]將n0元語法用于時間序列的特征構建過程中,他們將連續出現的兩個特征組成一個新的特征加入到特征字典.本文將一元和二元語法模型分別記為Unigram 和Bigrams,并對Bigrams 語法模型對變長單詞特征字典的性能影響進行了實驗分析.

1.4 特征的鑒別性評價

tf-idf 是一種用于尋找文本中關鍵詞的統計方法[26,27].它通常用來評估一個詞與一篇文檔的主題的相關程度.詞對一篇文檔的重要性與其在該文檔中出現的頻率成正比關系,同時,與其在語料庫中出現的頻率成反比.如果某個詞在所有的文檔中都出現,則意味著與主題并不相關.本文基于tf-idf 的基本思想設計了一種新的tfidf 統計量對特征的鑒別性進行評價.新定義的鑒別性評價指標主要從兩個方面對各類特征的鑒別性進行評價.

(1)某類特征和該類的相關程度.主要通過tf 值來度量,我們用某類特征在該類所有實例中的出現頻率來度量該特征和所屬類別的相關性.

(2)某類特征對該類的鑒別性強弱.實際中,某類中高頻特征也可能是其他類中的常見特征.此時,該特征不能有效區分不同類別.為此,我們用idf 度量該特征對所屬類別的鑒別性強弱.

下面首先對本文定義的tf-idf 計算公式進行介紹.我們用特征在某類中出現的相對頻率代替其在某個實例中出現的頻率以衡量它與某類實例的相關程度.我們用f(t,c)表示特征t在其類屬性c對應的實例集合中出現的次數,即,類屬性為c的特征t在類屬性為c的所有實例中出現的次數總和.特征t在其類屬性c對應的特征字典中的頻率tf(t,c)的計算公式為

其中,dict(c)表示類屬性c對應的特征字典,p表示任意特征.

我們在自然對數尺度下對特征頻率進行比較分析,對頻率公式tf(t,c)進行如下處理:

如果tf(t,c)值越大,則特征t在類屬性c對應的實例中出現的頻率越高,意味著它與類屬性越相關,但不能說明它對于類屬性c的鑒別性越強.為了準確度量特征的鑒別性,我們提出一個新的idf 計算公式.在定義的idf 公式的分子中用實例總數減去類屬性為c的實例數,即,只考慮類屬性不為c的實例中包含特征t的情況;idf 公式分母中只計數類屬性不為c且包含特征t的實例數,這樣可以直接反映特征t在其他各類實例中的出現頻率.我們使用的類屬性為c的特征t的逆文檔頻率idf(t,c)的計算公式為

其中,N為實例總數,Nc為類屬性為c的實例數,為數據集中包含特征t同時類屬性不為c的實例數.idf(t,c)值越大,說明類屬性為c的特征t在其他類中出現的頻率越低.基于tf-idf 的基本思想,我們定義的類屬性為c的特征t的鑒別性度量值d(t,c)的計算公式為

上式說明:類屬性為c的特征t在類屬性為c的實例中出現的相對頻率越高(tf(t,c)值越大),而在其他類實例中出現的次數越少(idf(t,c)值越大),則該特征對類屬性c的鑒別性越強.即,類屬性為c的特征t的d(t,c)值越大,其對類屬性c的鑒別性越強.為了給出我們定義的特征的鑒別性度量公式的直觀解釋,下面我們通過一個計算實例來與其他研究人員常使用的tf-idf 公式進行對比.

Table 2 A calculation example of formula tf-idf表2 tf-idf 公式計算實例

表2 中ti表示第i個特征,其所在列中數字表示包含該特征的各類別實例數,ci表示第i類,Ni表示數據集中第i類實例的個數.這里,為了計算方便,我們假定每個特征在每個包含該特征的實例中的出現頻數為1.

Senin 等人[19]和Schafer[28]使用的tf-idf 計算公式雖有小的不同,但本質上是一致的.這里,我們使用Schafer在文獻[28]中的公式進行說明.基于文獻[28]中的公式,特征t1對類屬性c1的tf 值為該特征在c1類中頻數取對數加 1,即,tf(t1,c1)=1+ln(5)=2.6;idf 值為數據類別總數和包含該特征的類別數的比取對數加 1,即,idf(t1,c1)=1+ln(2/2)=1;特征t1對類屬性c1的鑒別性評價值為df(t1,c1)=2.6×1=2.6.

同理可得:

tf(t1,c2)=1+ln(50)=4.9;idf(t1,c2)=1+ln(2/2)=1;df(t1,c2)=4.9×1=4.9,

tf(t3,c1)=1 +ln(90)=5.5;idf(t3,c1)=1 +ln(2/2)=1;df(t3,c1)=5.5×1=5.5,

tf(t3,c2)=1 +ln(2)=1.7;idf(t3,c2)=1 +ln(2/2)=1;df(t3,c2)=1.7×1=1.7.

上面的計算公式主要存在兩個問題.

(1)在tf 值的計算過程中直接對頻數取對數,極大地縮小了頻數間的差異.例如,c1類中特征t1和c2類中特征t1的頻數比為1:10,而計算tf 值后比例變為1:1.9.

(2)idf 值不能準確反映特征的鑒別性.例如,特征t1和t3對各類別的idf 值都為1,這導致實際的鑒別性評價值僅為tf 值,不能準確反映這兩個特征對不同類別的鑒別性.這是由于Schafer 在文獻[28]中的idf 計算公式使用類屬性取值個數和包含某特征的類的個數的比值來衡量特征對于某一類的鑒別性.然而,實際上,由于類內變異的存在,某類中的少數實例中可能具有其他類的特征,這會導致idf 度量出現偏差.

下面我們介紹基于本文定義的tf-idf 計算公式得到的特征鑒別性評價值.

首先,由于我們假定每個包含特征t的實例中特征t出現的頻數都為1,因此,某類特征的最大頻數不超過該類實例數.為了簡便起見,在計算過程中,我們假定各類特征的最大頻數為該類實例數.同時,根據表2,我們可知類別為c1的特征t1的頻數為5,則其頻率為5/100=0.05,因此tf 值為tf(t1,c1)=ln(1+0.05)=0.05.

由于在其他類實例(即,c2類,該類實例數為100)中包含類c1的特征t1的實例數為50,因此,根據公式(9),idf值為idf(t1,c1)=ln[(200-100+1)/51]=0.7.進而可得df(t1,c1)=0.05×0.7=0.035.

同理可得:

tf(t1,c2)=0.4;idf(t1,c2)=2.8;df(t1,c2)=0.4 ×2.8=1.12,

tf(t3,c1)=0.6;idf(t3,c1)=3.5;df(t3,c1)=0.6 ×3.5=2.1,

tf(t3,c2)=0.02;idf(t3,c2)=0.1;df(t3,c2)=0.02 ×0.1=0.002.

從上面的計算結果可以看出,基于本文的計算公式,特征t1對類c1和c2的tf 值比為1:8,與實際的頻數比1:10接近,更準確地反映了特征和各類別的相關性.另一方面,特征t3的對類別c1和c2的idf 值分別為3.5 和0.1.從中可以看出,特征t3對類別c1具有更強的鑒別性.這樣就有效解決了由簡單地通過特征所屬的類的個數進行鑒別性評價所導致的誤差.因此,從計算得到的df(t,c)值更容易區分同一特征對不同類別的鑒別性.例如,基于Schafer等人給出的公式,特征t3對類別c1和c2的df(t,c)值比為3.2:1,而基于本文定義的公式比值為1050:1.綜上所述,利用本文定義的鑒別性度量公式,更有利于區分各類別的鑒別性特征.

1.5 特征選擇閾值設定

通過固定閾值進行特征選擇忽略了不同長度滑動窗口生成的特征的分類性能,這會降低生成特征字典的的質量.為此,我們提出了一種動態閾值設定機制,該機制利用單一長度滑動窗口轉換得到的特征字典對應的分類性能為其生成的特征設定閾值.首先,我們用交叉驗證在單一長度滑動窗口轉換得到的訓練集上進行分類預測,假設獲得的精度為a;然后,基于該精度和最大精度的差值為其生成特征設定選擇閾值.滑動窗口對應精度越高,則其特征選擇閾值越小.函數f(a)將某窗口長度精度映射為該窗口生成特征的閾值因子,定義如下:

其中,amax為單一滑動窗口獲得的最大交叉驗證精度.

下面我們給出不同長度窗口對應的閾值計算公式:

其中,θ為加權因子(θ>0),al為長度為l的滑動窗口對應的交叉驗證精度.

加權因子θ用于校正閾值因子存在的偏差.滑動窗口對應的交叉驗證精度和最優值差距越大,則對應窗口生成特征的選擇閾值越大,從而根據整體分類性能對特征進行選擇可有效提升所建立的特征字典的有效性.

2 算法設計

本節對我們提出的特征字典構建算法進行詳細介紹.

2.1 滑動窗口最優單詞長度學習算法

各種特征包方法(BOP)的不同之處在于將實值序列轉換為單詞的具體過程.本文我們使用的數值序列離散化方法是SFA[22].本文提出為每個窗口長度在指定范圍內動態學習性能最優的單詞長度.下面首先介紹本文提出的為每個窗口長度學習最優單詞長度的算法.

算法1.computeBestWindowsF(D,min,max,minF,maxF,k).

算法1 給出了本文從指定區間為每個窗口長度尋找最優傅里葉值個數(即,最優單詞長度)的過程.首先,我們使用監督SFA 對訓練集D進行轉換得到數組AF(第2 行),其中,AF第1 維對應不同長度滑動窗口的序號,維度為max-min+1;第2 維對應訓練實例的序號,維度為N;第3 維為轉換后的實例,維度為對應滑動窗口從該實例上獲得的子序列的個數,且實例中的每個單詞長度為maxF(若滑動窗口長度小于maxF,則單詞長度為窗口長度).由文獻[22]可知,使用單一滑動窗口對N個實例進行監督傅里葉符號化轉換的時間復雜度為O(|Wi|×wilogwi),其中,wi為第i個滑動窗口長度,Wi為該滑動窗口在數據集D上生成的子序列集合,|Wi|表示集合Wi中元素個數.由于|Wi|包含O(N×n)個子序列,wi≤n,因此,算法1 第2 步的算法復雜度為O(N×n3logn).然后,我們在轉換后的數據AF的基礎上從區間[minF,maxF]中學習各滑動窗口對應的性能最優的單詞長度(第3 行~第12 行).

在學習各窗口最優單詞長度的過程中,我們依次只使用單一滑動窗口長度對訓練實例進行轉換(第6 行),這一轉換過程只需要截取每個單詞序列(AF[i][x][y])的前j個字母即生成指定長度的單詞,這一過程的時間復雜度可忽略.然后在轉換得到的數據集上利用交叉驗證進行預測(第7 行),并將正確預測實例數最大值對應的傅里葉值的個數作為該窗口長度對應的最優單詞長度(第8 行~第11 行).整個單詞長度的學習過程需執行max×maxF×k次模型訓練和預測過程.本文我們使用的預測算法是一種適用于處理大規模稀疏表示數據的運行較快的邏輯回歸算法[29],其算法復雜度取決于使用的梯度下降算法的收斂速度,這不在本文討論范圍之內.我們假定使用的分類模型的時間復雜度為T(n),則算法1 的時間復雜度為O(maximum(N×n3logn,max×maxF×k×T(n))).實驗過程中,我們將交叉驗證重數統一設為10.

2.2 鑒別性特征生成

本節我們介紹如何對訓練集和測試集中的實例進行符號化轉換,這一過程是算法1 中數據符號化轉換的主要內容.轉換過程主要分為兩步.

(1)利用監督SFA 技術將時間序列或子序列轉換為一組傅里葉值序列;

(2)利用離散化技術將傅里葉值轉換為字母表中的字母.

傅里葉值的選擇關系到時間序列及其子序列轉換得來的特征的鑒別性強弱.不同位置的傅里葉值的鑒別性強弱不同.我們利用F統計量度量每個位置的傅里葉值的鑒別性,然后選擇整體分類性能最優時對應的傅里葉值個數作為該長度窗口生成單詞的長度,最后利用裝箱技術將時間序列子序列對應的各最優傅里葉值轉換為字母表中的字母,這些字母共同組成一個單詞,即,特征.下面介紹本文使用的傅里葉值的挑選算法.

算法2.SelectFourierCoefficient(wi,W,bestF).

算法2 基于指定長度的子序列集合計算各個位置傅里葉值的鑒別性強弱.首先,對原始子序列進行離散傅里葉變換.由于單個長度為n的時間序列進行離散傅里葉變換的時間復雜度為O(nlogn),因此這一過程的時間復雜度為O(N×n2logn)(第2 行).矩陣A中每行對應一個子序列轉換得到的傅里葉值數組,A[i]表示所有子序列的第i個位置的傅里葉值組成的數組(即,矩陣A的第i列).算法2 中依次計算第i個位置傅里葉值的鑒別性強弱統計量F值,然后將對應的傅里葉值序號和F值組成一對元組添加到集合indexBestF中(第5 行~第7 行).最后,根據F值對集合indexBestF中的元素進行降序排列(第8 行),并返回其中前bestF個元組組成的集合(第9 行~第10 行).算法2 的時間復雜度為O(N×n2logn).

獲得子序列集合對應的最優傅里葉值序號數組后,我們對子序列進行符號化轉換,傅里葉值符號化的過程是將reali或imagj映射到符號空間的過程.下面給出將時間序列子序列轉換成特征的算法.

算法3(VLWEA).createFeature(indexBestF,S,alphabet).

算法3 給出了將時間序列子序列轉換成單詞的算法過程,首先,根據最優傅里葉值序號數組將轉換得到的子序列傅里葉值數組S中對應位置的傅里葉值依次利用裝箱技術映射為字母表中的字母(第3 行),并將獲得的字母依次組合起來構成一個單詞(第4 行),最后得到的單詞序列即為一個特征.這一符號化轉換過程需要在離散化區間中執行O(|S|×log|alphabet|)次運算[22].

2.3 特征選擇和鑒別性特征字典構建

本節我們給出基于變長單詞生成算法建立特征字典的過程,以及基于tf-idf 的特征選擇算法.首先,給出利用可變長度的單詞建立特征字典的算法過程.

算法4.createFeatureDictionary(min,max,windowBestF[],D,alphabet).

算法4 給出了基于訓練集建立特征字典的算法過程.首先,遍歷訓練集中的每個實例(第2 行);然后,利用長度從min 到max 的滑動窗口在時間序列上進行滑動獲得一系列子序列,并將這些子序列逐個轉換為單詞(第3行~第6 行);最后,將不重復的單詞作為特征加入到特征字典中(第7 行).算法4 中特征構建過程的主要部分為訓練集實例的SFA 轉換過程,因此,特征庫構建過程的時間復雜度也為O(N×n3logn).

基于特征包的方法轉換過程中會生成規模巨大的特征字典.為提高分類模型的效率,通常需要對輸入模型的特征進行選擇.本文提出了一種根據不同窗口生成特征的整體分類性能來動態設定特征選擇閾值的特征選算法.下面介紹本文提出的鑒別性特征字典構建算法.

算法5(TfIDfDynamicVLWEA).createTfIdfFeatureDict(dict,corrects[],θ).

算法5 介紹了本文使用的特征字典建立算法.首先計算各滑動窗口生成特征的選擇閾值(第2 行~第3 行).然后,遍歷字典dict中的每個特征(第4 行),計算特征的d(t,c)值(第6 行),然后將d(t,c)值大于等于其所對應的窗口閾值的特征加入到鑒別性特征字典TfIdfDict中(第9 行~第10 行).

綜上所述,算法1 最優參數學習和算法4 特征庫構建是本文模型的主要構成部分.因此,本文通過動態閾值構建變長特征字典模型的時間復雜度為O(maximum(N×n3logn,max×maxF×k×T(n))),其中,T(n)為參數學習過程中使用的分類模型算法復雜度.

2.4 基于特征字典的時間序列轉換表示

這節我們給出基于SFA 的BOP 模型中用于訓練集和測試集中實例轉換的算法[24].這一過程主要分為兩步:首先,在訓練集上利用某長度得到的所有不重疊子序列為每個窗口長度訓練用于時間序列轉換的監督符號傅里葉近似轉換對象,訓練轉換對象主要是為了計算滑動窗口子序列每個位置的F值(算法2)和用于將傅里葉值數組進行離散化的分箱區間(文獻[23]中有詳細介紹);然后,用訓練得到的轉換對象對時間序列進行轉換[24].下面我們給出時間序列的轉換算法.

算法6.SupervisedSFATransform(TfIdfDict,T,min,max,transferObjects[]).

算法6 給出了將一條時間序列轉換為符號序列頻數數組的算法.該算法用指定長度滑動窗口從給定時間序列的初始位置進行滑動依次獲得一系列可重疊的子序列(第2 行~第3 行),然后通過對應長度的監督符號傅里葉近似轉換對象將每個子時間序列依次轉換成一個單詞,即,時間序列包含的特征(第4 行).若轉換得到的特征在給定的鑒別性特征字典中,則將該特征加入到轉換后的實例中,作為轉換后實例的一個特征(第5 行~第6 行).此時,若轉換后的實例中沒有該特征,則為轉換實例加入該特征,并將該特征的取值設為1;若轉換后的實例中已有該特征,則該特征對應的值加1.最后返回獲得的特征頻數序列,即為轉換得到的時間序列.

邏輯回歸(logistic regression)是統計學中的經典分類方法.它是一個對數線性模型.本文我們使用基于L2正則化的邏輯回歸模型對轉換后的實例進行分類[29].此外,本文還利用邏輯回歸學習到的權重對特征字典進行分析.

3 實驗分析

本節對我們所提模型的相關實驗內容進行詳細介紹.主要包括4 個方面:模型參數分析、模型設計的方法分析、分類精度比較和模型的可解釋性分析.本文在UEA & UCR 時間序列知識庫提供的65 個長度小于750的時間序列數據集上對模型進行分析[30].表3 給出了各數據集的信息,其中,Train/Test 表示訓練集和測試集的實例數,n為時間序列長度,|C|為類屬性取值個數.我們的實驗環境CPU 是3.40GHz,內存為16G.

Tabel 3 Introduction to 65 time series datasets表3 65 個時間序列數據集介紹

3.1 模型參數實驗分析

本節我們對模型TfIDfDynamicVLWEA 的幾個重要參數進行實驗分析.

(1)訓練集規模和時間序列長度對模型運行時間的影響分析;

(2)模型精度和壓縮比(即,TfIdf 特征庫中特征數和初始特征庫中特征數的比值)對閾值加權因子θ的敏感性;

(3)模型精度對最大滑動窗口長度和最大單詞長度的敏感性.

首先,分析模型的運行時間分別與訓練集規模、時間序列長度的關系.我們使用一個生成的二分類數據集進行實驗.訓練集和測試集中各包含100 個長度為200 的時間序列,且每個數據集中兩類實例個數相同.在實例數遞增分析實驗中,我們設置初始訓練集和測試集各由原數據集中10%的實例構成,然后訓練集和測試集中實例個數以原訓練集實例個數的10%遞增,并始終保持數據集中類屬性分布和時間序列長度不變.我們分別統計模型在上述10 組數據上的運行時間.在時間序列長度遞增分析實驗中,我們將時間序列長度從原長度的10%開始按10%遞增,同時保持訓練集和測試集規模不變,這樣進行10 組實驗并統計模型運行時間.圖2 給出了模型TfIDfDynamicVLWEA 的運行時間的實驗結果.

從圖2 可以看出,模型TfIDfDynamicVLWEA 的運行時間與訓練集實例個數呈線性關系,而與時間序列長度呈多項式關系.這與第2 節分析獲得的模型算法復雜度相符.

下面我們在10 個數據集上對模型各參數的敏感性進行分析.圖3 和圖4 中加粗曲線表示10 個數據集上模型TfIDfDynamicVLWEA 的平均精度變化曲線,我們標出了每個點對應的平均精度.

首先,我們分析加權因子對模型TfIDfDynamicVLWEA 的精度和特征庫壓縮比的影響.從圖3(b)中可以看出,隨著θ的增大,壓縮比顯著減小,當θ大于5 后,減速變緩.同時,從圖3(a)中可以看出,隨著θ的增大,數據集精度呈現出不同的變化趨勢.例如,數據集Beef 和Ham 上的模型精度隨θ增大,先增大后減小;數據集BeetleFly、ShapeletSim 等的精度隨θ增大,變化不大;數據集Herring 的精度隨θ增大,呈現先減小后增大的趨勢.由于當θ為3時,模型的平均精度最大,壓縮比較優,因此綜合考慮,我們選擇將模型的參數θ設為3.

Fig.2 The running time of model TfIDfDynamicVLWEA圖2 模型TfIDfDynamicVLWEA 運行時間分析

Fig.3 Analysis of the influence of the parameter θ on model TfIDfDynamicVLWEA圖3 模型TfIDfDynamicVLWEA 參數θ的影響分析

圖4 給出了10 個數據集上模型TfIDfDynamicVLWEA 精度及平均精度隨最大滑動窗口長度和最大單詞長度的變化趨勢.

從圖4(a)和圖4(b)可以看出,不同數據集對不同參數的敏感性不同.從圖4(a)中可以看出,Lightning7 數據集上的精度隨最大滑動窗口長度的遞增總體呈增長趨勢,數據集 ShapeletSim、BeetleFly、Symbols、ToeSegmentation1 和Ham 對最大窗口長度不敏感,數據集Beef 和Herring 對最大滑動窗口長度較為敏感,隨窗口長度呈遞增趨勢,精度變化明顯.上述實驗結果表明,很難設定最優滑動窗口長度.由于最大滑動窗口長度取250 時,圖4(a)中平均精度值最大,同時也為了保證對比實驗的公平性,我們將模型TfIDfDynamicVLWEA 的最大窗口長度設為250.

從圖4(b)可以看出,數據集Beef 和Herring 對最大單詞長度較為敏感,數據集ShapeletSim、BeetleFly、Car、Symbols 和ToeSegmentation1 對最大窗口長度不敏感,數據集Ham、Lightning7 和ToeSegmentation2 隨最大單詞長度發生變化,精度有波動.在區間[8,18]中10 個數據集上的平均精度在15 時最大,因此,我們將模型中TfIDf DynamicVLWEA 的最大單詞長度設為15.

此外,由于目前基于SFA 對時間序列進行符號轉換的研究[19,23,24]表明字母表大小設置為4 時,BOP 模型具有更強的魯棒性.因此,本文將字母表大小固定為4.

Fig.4 Sensitivity analysis of model TfIDfDynamicVLWEA accuracy to parameters max and maxF圖4 模型TfIDfDynamicVLWEA 精度對參數max 和maxF 的敏感性分析

3.2 設計方法分析

我們在表3 中的65 個數據集上對本文所提模型中設計的新方法進行分析.主要包括如下3 個方面.

(1)可變單詞長度和固定單詞長度比較;

(2)Bigrams 特征對可變長度特征字典性能的影響;

(3)基于tf-idf 統計量的特征選擇算法和基于卡方統計量特征選擇算法的比較.

我們使用符號VLWEA 表示采用可變長度單詞生成算法的特征字典建立模型,FLWEA 表示采用固定長度單詞生成算法的特征字典建立模型.VLWEA_U、FLWEA_U 分別表示對應模型在特征字典建立過程中只將單個單詞作為特征,且不對特征進行選擇.FLWEA_B、VLWEA_B 分別表示模型在特征字典建立過程中使用Bigrams 語法模型進行特征生成,即,將連續的兩個單詞組成一個新的單詞作為特征.符號Chiδ和TfIdfδ分別表示模型中使用單一閾值的特征選擇算法,其中,δ表示使用的閾值.例如,Chi3VLWEA_U 表示模型VLWEA_U 中使用閾值為3 的卡方統計量進行特征選擇,即,保留卡方統計值大于等于3 的特征.TfIDfDynamicVLWEA 表示本文提出的使用動態閾值設定的VLWEA 模型,其中加權因子θ統一取為3.我們使用的對比模型包括:不利用Bigrams 語法進行特征生成,只利用可變長度單詞或固定長度單詞組成的特征字典建立模型VLWEA_U 和FLWEA_U,結合不同特征選擇算法和不同閾值的模型Chi3FLWEA_U、Chi3VLWEA_U、TfIdf0.3VLWEA_U、TfIdf0.3VLWEA_B 和TfIDfDynamicVLWEA_U.實驗中,上述所有VLWEA 和FLWEA 模型的最小和最大滑動窗口長度都分別設為4 和250,字母表的大小統一設定為4,其他處理方式與定長單詞生成模型WEASEL 相同.同時,使用相同參數設置的分類模型對轉換后的測試集進行分類預測.

接下來,我們利用Demsar 提出的模型分類性能顯著性和平均排名比較方法在65 個數據集上對本文新提出的方法進行分析[31].表4 中給出了多種條件下建立的特征字典對應的邏輯回歸模型在65 個數據集上的分類精度、平均精度和模型在65 個數據集上精度最高的數據集個數.表4 中給出的模型實驗結果都是通過在每個數據集上進行5 次實驗取均值獲得的.在對特征選擇算法的性能對比過程中,為了對比實驗的公平性,我們通過選取適當的閾值,使得不同特征選擇算法在65 個數據集上的平均壓縮比相近,即,65 個數據集上選擇后的特征字典大小和原特征字典大小的比的平均值相近.模型TfIdf0.3VLWEA_U、TfIdf0.3VLWEA_B、Chi3VLWEA_U和Chi3FLWEA_U 的平均壓縮比分別為27.3%、24.6%、32.6%和36.3%.

從圖5 中我們可以看出:使用閾值為3 的卡方統計量進行特征選擇的定長單詞生成模型Chi3FLWEA_U 的性能排名最低,這說明可變長度單詞生成算法有助于提升模型分類性能.使用本文定義的鑒別性特征評價統計量進行特征選擇的變長單詞生成模型TfIdf0.3VLWEA_U 在65 個數據集上特征字典的平均壓縮比優于基于卡方統計量的特征選擇算法模型(Chi3VLWEA_U 和Chi3FLWEA_U)壓縮比的條件下,分類性能排名更優.這說明,本文提出的基于tf-idf 的特征選擇算法與基于卡方統計量的特征選擇算法相比,能夠更有效地進行特征選擇.與此同時,從圖5 中我們還可以看出:在相同閾值條件下,結合Bigrams 語法的特征生成模型TfIdf0.3VLWEA_B 的排名低于模型TfIdf0.3VLWEA_U,這說明,Bigrams 語法模型不能有效地提升變長特征生成算法的分類性能.因此,綜合考慮生成特征字典規模和分類模型的運行效率,在接下來的對比實驗中,我們的模型中不再采用Bigrams 語法進行特征生成.此外,基于動態閾值的特征選擇模型TfIDfDynamicVLWEA_U 在65 個數據集上的性能平均排名第1,這顯示了本文提出的動態閾值設定算法的有效性.從表4 中我們還可以看出,TfIDfDynamic VLWEA_U 在65 個數據集上的平均精度最高和精度最高的數據集個數最多.因此,在接下來的實驗分析中,我們使用TfIDfDynamicVLWEA_U 模型與其他各類模型進行對比,并將其簡記為VLWEA.

Fig.5 The classification performance significance analysis and the average rank of the dictionaries built under 6 conditions on 65 datasets圖5 6 種條件下構建的字典在65 個數據集上的分類性能顯著性分析和模型平均排名

3.3 分類精度比較

本節我們將本文所提鑒別性特征字典建立模型VLWEA 分別與基于特征包的模型、1NN 模型、基于shapelet 的模型和集成分類模型的分類精度進行對比分析.兩個常用的基準1NN 分類模型為:基于歐式距離的1NN(ED1NN)和基于動態時間規整的1NN(DTW1NN);兩個基于Shapelet 的非集成分類模型包括:快速Shapelet分類(fast Shapelet,簡稱FS)[32]算法和Grabocka 等人提出的基于啟發式Shapelet 搜索算法的分類(learning Shapelets,簡稱LS)模型[33];6 種特征包模型為:SAXVSM、BOSS、TSBF、使用Bigrams 語法進行特征生成的模型WEASEL_B、使用Unigram 語法進行特征生成的WEASEL 模型WEASEL_U 以及Kate 等人提出的基于DTW 距離的特征生成算法(DTW features,簡稱DTW_F)[34];3 種集成分類模型包括:Deng 等人提出的基于11 種近鄰分類模型的集成分類算法(elastic ensemble,簡稱EE)[35]、基于Bostrom 等人提出的Shapelet 轉換表示方法建立的集成分類模型(ST)[16]和COTE.模型WEASEL_B 和WEASEL_U 的參數設置采用Schafer 等人提供的代碼中的默認設置,模型預測精度我們取5 次實驗的平均值(見表4,其中,C2FB 代表WEASEL_B,C2FU 代表WEASEL_U),其他模型采用Bagnall 等人[5]提供的實驗結果.下面我們將本文模型和各類模型在65 個數據集上分別進行對比,并給出了VLWEA 和對比模型相比在65 個數據集中的精度高、平、低的數據集個數,例如,“30/20/15”表示和對比模型相比,VLWEA 有30 個更準確,20 個相同,15 個更差.

首先,我們對VLWEA 和6 個BOP 模型的比較結果進行分析.從圖6 和具體的“高/平/低”對比結果統計可以看出,VLWEA 模型與模型 WEASEL_B(35/9/21)、WEASEL_U(41/6/18)、TSBF(48/5/12)、BOSS(35/2/28)、SAXVSM(55/2/8)及DTW_F(51/2/12)相比,分類精度高的數據集個數都是最多的.VLWEA 在不使用Bigrams 語法進行特征生成的條件下,比WEASEL_B 在更多的數據集上取得更好的分類效果.這再次說明本文可變長度單詞生成算法在特征生成上的有效性.

從圖7 中和具體的“高/平/低”對比結果統計可以看出,VLWEA 與模型ED1NN(60/1/4)和DTW1NN (53/3/9)相比顯著更好,分類效果具有絕對優勢.

圖8 和對比結果統計說明,VLWEA 和2 個不使用集成分類算法的基于Shapelet 的分類模型LS(49/6/10)、FS(60/3/2)對比同樣具有顯著優勢.

圖9 和具體的精度“高/平/低”對比結果統計說明,VLWEA 和3 個集成分類模型COTE(27/7/31)、EE(45/3/17)、ST(40/6/19)相比,比EE 和ST 更好,但比COTE 略差.

Tabel 4 Accuracies of the feature dictionaries built under eight different conditions (%)表4 8 種不同條件下建立的特征字典對應的模型分類精度(%)

Tabel 4 Accuracies of the feature dictionaries built under eight different conditions (%)(Continued)表4 8 種不同條件下建立的特征字典對應的模型分類精度(%)(續)

上面表4 中我們將各特征字典構建模型符號表示中的TfIdf 簡記為T,Chi 簡記為C,Dynamic 簡記為D,VLWEA_X 和FLWEA_X 分別簡記為VX 和FX,X 表示使用的n0元特征生成模型,U 表示Unigram,B 表示Bigrams.例如,模型TfIDfDynamicVLWEA_U 記作TDVU.表4 中每行加粗數值表示對應行的最優值.

Fig.6 Comparison of accuracy between VLWEA and 6 BOP models on 65 datasets圖6 VLWEA 和6 個BOP 模型在65 個數據集上的分類精度比較

Fig.7 Comparison of accuracy between VLWEA and 2 1NN models on 65 datasets圖7 VLWEA 和兩個1NN 模型在65 個數據集上的分類精度比較

Fig.8 Comparison of accuracy between VLWEA and 2 shapelet models on 65 datasets圖8 VLWEA 和兩個Shapelet 分類模型在65 個數據集上的分類精度比較

Fig.9 Comparison of accuracy between VLWEA and 3 ensemble models on 65 datasets圖9 VLWEA 和3 個集成分類模型在65 個數據集上的分類精度比較

最后,我們對本文提出的模型VLWEA 分別與其同類型和異類型模型的性能顯著性和平均排名進行對比分析.從后文所示圖10 可以看出,VLWEA 與其他6 個BOP 模型在65 個數據集上的性能相比沒有顯著差異,但是VLWEA 的分類精度排名最高.與非BOP 模型相比,VLWEA 的排名只比COTE 差,比ST、EE、LS、FS、ED1NN和DTW1NN 的排名都更好.與當前最先進的模型的對比結果說明了本文所提出的特征字典建立方法的有效性.

3.4 實例分析

本節對VLWEA 模型的可解釋性進行分析.我們選擇多分類數據集CBF 對模型學習到的最優單詞長度和生成特征的鑒別性進行分析.該數據集的訓練集實例共有3 類.我們用每類實例各特征的頻數平均值組成一個均值序列代表該類實例.圖11 給出了CBF 各窗口長度的最優單詞長度的箱型圖和9 個由原始序列生成的鑒別性子序列圖示.

從圖11(a)可以看出:在忽略極少數異常值的情況下(如圖11(a)中的長度15),25%的最優單詞長度落在[3,4]之間,50%的最優單詞長度落在區間[4,7]中,25%的最優單詞長度位于[7,11]之間,因此,固定長度單詞生成算法在單詞生成過程中會不可避免地損失鑒別性信息或攜帶大量冗余信息.從圖11(b)可以看出,本文提出的可變長度單詞生成算法可以有效地學習最優單詞長度.例如,圖中不同長度:24、27、50、16 和120 對應的特征分別有針對性地將原序列尾部的冗余信息“**…*”去除(符號“*”代表字母表中的任意字母),只保留具有鑒別性的部分.此外,對于圖11(b)所示長度為120 的4 個原始序列,若單詞長度為2,則只生成一個特征“cc”;單詞長度為3 時,生成特征“ccb”和“ccd”;單詞長度大于4 時,可以生成4 個特征,但會包含一定冗余信息.而本文算法可以有效學習單詞的最優長度,且不包含冗余信息.這也再次驗證了圖11 給出的示例.

Fig.10 The classification performance analysis and the average ranks of VLWEA and 13 models on 65 datasets圖10 VLWEA 和13 個分類模型在65 個數據集上的分類顯著性分析及平均排名

Fig.11 Optimal word lengths and 9 generation features obtained by VLWEA圖11 VLWEA 學習得到的最優單詞長度和9 個生成特征

由于建立的特征字典規模巨大,我們根據學習到的權重選擇top-10 個特征對數據集進行表示.圖12 中給出了3 類實例均值序列的直方圖.從圖12 中我們可以看出,這些特征具有明顯的鑒別性,例如,特征“74ccadd”只有類屬性為c0的實例具有,特征“29ccbbb”和“29acbbb”特征對于類屬性為c1的實例具有較強的鑒別性.具有特征“54ccbcd”和“74ccbdd”的實例類屬性為c1的可能性很小.

Fig.12 Top-10 discriminant features generated by VLWEA on dataset CBF圖12 數據集CBF 上VLWEA 生成的top-10 鑒別性特征

4 結 論

在進行時間序列數據挖掘之前,對時間序列數據重新表示是一個重要的研究課題,其目的是,一方面通過減少算法實際處理數據的量來提高算法的運行速度,另一方面,充分表達原始時間序列數據的本質內容以提高分類精度.針對目前基于SFA 的時間序列進行離散化表示方法存在的問題,本文提出了一種可變長度單詞抽取算法,該算法可以有效學習不同滑動窗口對應的最優單詞長度.與此同時,針對特征字典規模巨大的問題,本文定義了一種新的鑒別性特征選擇統計量,并設計了一種動態閾值設定機制來對生成的特征進行選擇,該方法在有效縮小特征字典規模的同時,可以獲得較高的分類精度.

猜你喜歡
單詞分類特征
分類算一算
單詞連一連
如何表達“特征”
不忠誠的四個特征
當代陜西(2019年10期)2019-06-03 10:12:04
分類討論求坐標
數據分析中的分類討論
看圖填單詞
教你一招:數的分類
抓住特征巧觀察
線性代數的應用特征
河南科技(2014年23期)2014-02-27 14:19:15
主站蜘蛛池模板: 高h视频在线| 青青草国产在线视频| 伊人久久大香线蕉aⅴ色| 欧美 亚洲 日韩 国产| 亚洲一区精品视频在线| 欧美性猛交一区二区三区| 伊人五月丁香综合AⅤ| 精品一区国产精品| 狠狠色婷婷丁香综合久久韩国| 欧美视频二区| 精品国产免费观看一区| 亚洲第一视频网| AV无码一区二区三区四区| 国产永久在线视频| 欧美成人a∨视频免费观看| 欧美不卡在线视频| 波多野结衣的av一区二区三区| 思思热精品在线8| 91视频青青草| 成人精品在线观看| 亚洲成a人片77777在线播放 | 亚洲二三区| 中文天堂在线视频| 欧美国产成人在线| 午夜毛片免费看| 亚洲动漫h| 精品国产成人高清在线| 国产91丝袜在线播放动漫 | 国产区人妖精品人妖精品视频| 青青草国产在线视频| 1769国产精品免费视频| 天天综合网在线| 乱色熟女综合一区二区| 在线无码九区| 久久无码免费束人妻| 91精品啪在线观看国产60岁| 国产人人射| 在线无码九区| 色老头综合网| a毛片基地免费大全| 国产精品第5页| 77777亚洲午夜久久多人| 99re免费视频| 色综合久久88色综合天天提莫| 精品国产一二三区| 国产欧美亚洲精品第3页在线| 日韩精品一区二区三区免费在线观看| 国产精品免费电影| 国产欧美高清| 天天激情综合| 国模粉嫩小泬视频在线观看| 毛片视频网| 黄色国产在线| 99热这里只有精品在线播放| 日韩无码黄色| 欧美h在线观看| 456亚洲人成高清在线| 国产日韩久久久久无码精品| 国产无人区一区二区三区| 国产亚洲精品精品精品| 日韩精品欧美国产在线| 亚洲AV无码不卡无码 | 国产成人永久免费视频| 国产亚洲精品97在线观看| 无码'专区第一页| 亚洲国产一成久久精品国产成人综合| 久久精品国产免费观看频道| 人人看人人鲁狠狠高清| av手机版在线播放| 精品国产美女福到在线不卡f| 亚洲综合专区| 亚洲国产精品成人久久综合影院| 欧美日韩免费在线视频| 欧美日韩导航| 欧美黄色a| 国产精品亚欧美一区二区三区| 亚洲欧美成人影院| 国产欧美高清| 人妻少妇久久久久久97人妻| 久久中文字幕不卡一二区| 久久亚洲黄色视频| 青青国产成人免费精品视频|