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

分段聚合近似和數值導數的動態時間彎曲方法

2016-05-24 12:01:38李海林梁葉
智能系統學報 2016年2期
關鍵詞:分類方法

李海林,梁葉

(華僑大學 信息管理系, 福建 泉州 362021)

?

分段聚合近似和數值導數的動態時間彎曲方法

李海林,梁葉

(華僑大學 信息管理系, 福建 泉州 362021)

摘要:針對動態彎曲方法對時間序列數據相似性度量的質量和效率的局限性,本文提出一種基于分段聚合近似和數值導數的動態時間彎曲方法。該方法通過分段聚合近似將時間序列數據進行有效地降維,再結合數值導數對降維后的特征序列構建新特征序列,并且設計符合該特征序列相似性度量方法。實驗結果分析表明,與傳統動態彎曲方法相比,新方法具有較好的度量質量,能在時間序列數據挖掘中得到較好的分類效果,且在低維空間具有較高的分類效率,具有一定的優越性。

關鍵詞:動態時間彎曲;時間序列;分段聚合近似;數值導數;相似性度量;分類;數據降維;特征表示

中文引用格式:李海林,梁葉. 分段聚合近似和數值導數的動態時間彎曲方法[J]. 智能系統學報, 2016, 11(2): 249-256.

英文引用格式:LI Hailin, LIANG Ye. Dynamic time warping based on piecewise aggregate approximation and data derivatives[J]. CAAI transactions on intelligent systems, 2016, 11(2): 249-256.

近年來,時間序列成為了數據挖掘研究領域中的熱點,其研究成果已被廣泛地應用在工業、金融、氣象、航天、生物醫學等各種領域[1-5]。然而,相似性度量方法成為時間序列數據挖掘領域中基礎而又極其關鍵的工作之一,很多挖掘任務的結果易受其相似性度量質量和效率的影響。目前,動態時間彎曲(dynamic time warping, DTW)是一種魯棒性較強的方法之一,其最早應用于語音識別[6],如今已被廣泛應用于時間序列相似性的度量[7-9]。歐氏距離雖能快速度量時間序列之間的相似性,但其限于相同時間點上的數據匹配,而且對異常數據點較為敏感。DTW不僅可以彎曲度量不等長時間序列之間的相似性,其度量質量能較好地反映數據之間的異步關系[10]。然而,由于DTW需要在代價矩陣中計算最優彎曲路徑,使其具有較高的時間復雜度或較好的度量精度[11-13]。

時間序列降維表示是通過某種方法將原始序列進行轉換或者特征提取,以達到將原始時間序列在低維度表示的目的。由于直接使用和計算原始時間序列可能需要付出較大的代價,如消耗較多存儲空間,運行效率低等問題,加上原始時間序列包含很多“噪音點”,容易導致相似性度量結果不準確。因此,為簡化數據模型和算法的復雜性,提高數據挖掘技術的性能,有必要對原始時間序列進行降維處理,達到效率和準確性的平衡。目前,針對時間序列降維的方法已有很多且較為成熟,如分段聚合近似(piecewise aggregate approximation, PAA)[14],該方法通過將時間序列進行平均分段,并以分段均值反映分段信息,最終達到數據降維的目的;分段線性近似(piecewise liner approximation, PLA)[15]利用線性模型對時間序列進行分割,根據不同的分割策略可以得到不同的時間序列降維效果;基于域變換的離散傅里葉變換(discrete fourier transform, DFT)[16]和離散小波變換(discrete wavelet transform, DWT)[17],利用變換后的系數對時間序列進行特征表示;奇異值分解(singular value decomposition, SVD)[18],利用數值計算將高維數據轉換到低維空間,不僅應用在時間序列數據降維及索引,而且也被廣泛地應用于模式識別、圖像壓縮等等。由于降維方法對整個研究過程起著十分關鍵的作用,不僅要求算法能夠盡可能簡單快速,以降低時間消耗,也要盡量充分反映時間序列的信息。因此研究過程中應選擇合適的降維方法來解決問題。

鑒于傳統DTW方法在度量時間序列相似性時具有計算時間效率較低和缺少考慮序列的形態(即凹凸性)等問題,并且從使用較為簡單的數據降維方法以簡化模型的角度出發,本文提出通過利用分段聚合近似(PAA)方法對時間序列進行數據降維,獲得保持原始序列主要特征的數據序列,再構造基于數據點的一階導數的新特征序列,結合動態時間彎曲提出一種新的距離度量方法,即近似導數動態時間彎曲方法。數值實驗分類結果和效率分析表明,該方法能從形態視角出發,較好地從數據集中找出相似序列,提高數據挖掘領域中分類算法的質量,具有一定的優越性。

1相關理論基礎

1.1分段聚合近似(PAA)

分段聚合近似是一種有效的數據降維方法,其在減小數據存儲和提高計算效率方面起到了很大的作用。通過對長度為n的序列S=(s1,s2,…,sn)轉化為另一條長度為m的序列Q=(q1,q2,…,qm),實現時間序列的數據降維和特征表示,其中n>m,且令k=n/m。新序列中任意元素qi滿足:

(1)

對序列S進行分段后求出每段的均值,每段均值形成序列Q,進而達到序列S數據降維的目的。如圖1,子圖(a)顯示了長度為60的時間序列S通過PAA被平均分成10段,每段均值代表相應片段的特征;子圖(b)顯示由10個均值組成的新序列Q。

圖1 基于PAA的時間序列數據降維和特征表示Fig.1 Dimensionality reduction and feature representation of time series based on PAA

分段聚合近似方法通過均值來表示序列片段的特征,容易忽略數據的局部形態變化情況。然而,在實際運用中,時間序列的整體形態通常是關注和研究的重點。PAA得到的序列特征不僅能夠很好地反映較長時間序列數據形態的整體變化趨勢,而且還能對時間序列數據進行數據降維,起到提高相關數據挖掘算法效率的作用。

1.2動態時間彎曲

動態時間彎曲最初被應用于語音識別中,常被運用到比較2條時間序列的相似性。針對2條時間序列S=(s1,s2,…,sn)和Q=(q1,q2,…,qm),對任意兩點之間的距離構建一個n×m的距離矩陣D,其中d(i,j)表示時間序列數據點si和qj之間的距離,即d(i,j)=(si-qj)2。DTW的基本思想就是從距離矩陣中尋找一條使得2條序列之間的累計距離最小的路徑,其最小累積距離值為

(2)

彎曲路徑是一條具有連續K個距離矩陣元素的集合W=(w1,w2,…,wK),第k個元素為wk=(i,j)k且max(n,m)≤K≤n+m-1。與此同時,彎曲路徑通常要遵循著3個條件:

1)邊界性:w1=(1,1),wk=(n,m);

2)連續性:wk=(a,b),wK-1=(a′,b′),則a-a′≤1,b-b′≤1;

3)單調性:wk=(a,b),wK-1=(a′,b′),則a-a′≥0,b-b′≥0;

最優彎曲路徑的求解可以通過動態規劃來找出最小累積距離,即

(3)

式(3)說明當前元素的累積距離(i,j)是當前元素距離值與鄰近3個元素累積距離的最小值之和,且γ(n,m)表示DTW度量時間序列S和Q的最小距離 DTW(S,Q)=γ(n,m)。

如圖2(a)所示,動態時間彎曲方法可以很好地匹配不同時間點但具有相似形態的序列特征,使得時間序列的相似性特征能夠較好地得到體現。另外,通過累積矩陣的構建可以獲取其最優彎曲路徑和最小度量距離,如圖2(b)所示。

圖2 動態時間彎曲及其最優路徑Fig.2 Dynamic time warping and the best path

DTW在很多領域的應用中都有著很好的效果,然而較高的時間復雜度成為了該方法的瓶頸問題。目前也有部分學者對其進行了改進,使其具有較好的度量效率。另外,由于DTW有時也存在過度“變態”彎曲的狀況,即一條時間序列的一個值對應著另一條序列的較長子序列,進而造成計算結果不精確的問題。針對該問題,Keogh和Pazzani提出了導數動態時間彎曲(DDTW),有效地克服了DTW對時間軸的過度彎曲時造成的影響。

(4)

2近似導數動態時間彎曲方法

若直接使用動態時間彎曲方法對時間序列進行相似性度量,則時間消耗太大,并且也會造成時間過度彎曲等問題,不利于大規模的長度較長的時間序列數據挖掘。為此,在使用動態時間彎曲度量相似性之前,通常先對時間序列數據進行數據降維,不僅可以得到反映時間序列大部分信息的低維特征表示,而且還能近似表示原始時間序列的整體性信息。除此之外,在分類過程中,針對低頻時間序列數據缺少具體局部信息的狀況,序列被歸為哪一類通常可能取決于其整體形狀,而不是局限于嚴格的數值比較。因此,相似性度量的結果不僅需要充分考慮序列的具體數值,也要考慮序列整體形態(凹凸性)。

傳統研究方法是對時間序列進行數據降維,然后使用傳統距離公式直接進行相似性度量。由于數據降維會損耗時間序列的部分信息,實驗效果也就相當于以準確性來換取時間效率。鑒于時間序列具體數據值的重要性、形態的特殊性、度量準確性以及高運行效率的必要性,結合分段聚合近似以及數值導數,提出一種基于分段聚合近似和數值導數的動態時間彎曲方法,稱為分段聚合導數距離(piecewiseapproximationderivativedistance,PADD)。該方法主要綜合考慮時間序列的具體數值信息和形態特征,以獲得更好的度量效果,提高運行效率,減少異常點對度量的影響為目標,同時從時間序列的整體數據匹配和形態的凸凹性出發,更為完善地度量時間序列異步相關性。

(5)

式中:α、β是通過

(6)

求解得到。如圖3所示,參數α、β可以通過三角公式求解得到,兩者取值范圍均在[0,1],且式(6)確保了兩者之和為1,使其特征序列的動態時間彎曲距離是線性組合關系。為方便起見,系數α、β可由參數θ來確定。

圖3 參數α、β選取原理Fig.3 The principle of selecting parametersαandβ

結合上述思想,本文提出的分段聚合導數距離(PADD)算法步驟如下:

輸入時間序列S=(s1,s2,…,sn)和Q=(q1,q2,…,qm),分段數目ws和wq以及參數θ;

輸出度量距離DPADD(S,Q)。

1)分別對時間序列進行平均分段,每個子序列長度分別為ks=n/ws和kq=m/wq。

在上述算法中,1)實現了時間序列的分段,使用等劃分的方法實現數據劃分并獲得相應的子序列片段;2)中方法對子序列段的均值特征表示,而且大量研究結果表明[19],對于高維時間序列來說,PAA方法可以保留足夠的數據信息來反映原序列的波動形態;3)實現均值特征序列的導數表示,以便更好地反映特征序列的形態變化情況。最后,利用式(6)同時從均值特征序列和導數特征序列來反映原始時間序列的數據信息,以確保距離度量函數能夠較好地通過特征序列之間的數值關系來反映原始時間序列之間的數值和形態變化關系。

另外,從算法的整個描述過程,容易分析得到新方法的計算時間效率。長度為n的時間序列進行分段聚合近似時,其時間復雜度為O(n),計算一階導數的時間復雜度也為O(n)。另外,對PAA特征序列和導數序列進行DTW距離度量所需要的時間復雜度分別為O(wswq)和O((ws-2)(wq-2)),故PADD整個算法過程的時間復雜度近似為O(2wswq+4n)。由于分段數目ws和wq通常遠小于原始時間序列的長度,即ws?n和wq?n,故PADD的時間復雜度O(2wswq+4n)要低于與傳統DTW和DDTW的時間復雜度O(n2)。特別地,數值實驗的分類結果和時間計算效率分析表明,PADD能夠較好地度量時間序列的相似性,反映數據形態之間的變化關系,進而提高時間序列數據挖掘算法的分類結果,還能在低維特征空間下獲得較高的時間效率。

3數值實驗

為了更好地理解和檢驗PADD方法的性能,實驗通過對時間序列數據集進行分類實驗,與傳統DTW和DDTW方法進行比較,驗證新方法在不同時間序列數據集中的距離度量質量和時間計算效率。

3.1分類實驗

實驗一共進行兩個階段,1)訓練模型,使用“留一法”(leave-one-out)求出合適的參數θ構建度量距離;2)測試模型,得到分類效果。由于時間序列取值范圍差距較大,度量相似性結果可能會出現很大偏差,因此有必要對所有時間序列數據進行歸一化處理,消除量綱。

使用訓練集訓練模型過程中,給定時間序列分段數目w以及范圍為[0,π/2]且步長為0.01的θ,利用分段聚合近似以及一階求導分別得到新的特征序列。訓練階段結合“枚舉法”測試在w情況下每個θ的分類情況,選擇錯誤率最小的θ作為該w情況下的最優參數。一旦求得參數立即構建度量距離公式(5),進入2)利用最近鄰(1-NN)分類方法檢測該距離公式分類效果。

圖4給出了數據集ECG在不同w取值下的最佳參數θ。根據式(5)與式(6)易知,θ取值越大,時間序列相似性度量形態匹配所占比重越大;反之,則形態匹配的比重就越小。圖4也說明了時間序列數據集在大部分情況下更傾向于形態匹配而不是具體數值上的比較。

圖4 不同w下的最優參數θFig.4 The best parameterθfor different w

本次實驗使用Keogh教授提供的數據集[20],并隨機抽取20個數據集進行分類實驗。如表1所示,包含了各個數據集的基本詳細信息,包括數據集名稱、類別個數、訓練集個數、測試集個數、時間序列長度。表2包含了DTW、DDTW、PADD的平均錯誤率,并且給出了PADD的最小錯誤率,其中括號表示能夠取得PADD的最小分類錯誤率的w,例如當CBF數據被降維到原始序列長度的70%時,可以取得最小錯誤率0。

表1 時間序列數據集信息

注:表頭NO.C、S.Tr、S.Tst、Length分別表示類別個數、訓練集個數、測試集個數、序列長度。數據集名D.S.Z.、E.F.D.、I.P.D.、M.I.、M.S.、S.A.R.S.、S.A.R.S II、S.C.分別表示Diatom Size Reduction、ECG Five Days、Italy Power Demand、Medical Images、Mote Strain、Sony AIBO Robot Surface、Sony AIBO Robot Surface II、Synthetic Control.

表2 實驗分類錯誤率

從實驗結果來看, PADD的平均分類效果與DTW和DDTW相比有著較大的優勢,并且從最小分類錯誤率的均值及標準差上來看也取得了不錯的效果,進而驗證了PADD在時間序列數據挖掘中距離度量的有效性。實驗結果反映了時間序列分段聚合后的特征序列長度大多數為原始序列長度的50%~100%,即w大多數取值為5~10時,PADD將會取到比較好的效果。這也說明在進行序列轉換過程中,將時間序列的整體信息以及具體細節盡量完整保存下來會取得更好的效果。

為了更為清楚地剖析分類結果在不同參數下的情況,給出了Sony AIBO Robot Surface和ECG這兩個數據集的分類情況,如圖5所示。從圖5 中可以看出,分類錯誤率呈高低波動狀態。對于某些數據集如Sony AIBO Robot Surface來說,在不同w情況下并非所有的PADD錯誤率都是最小的,但是在最小錯誤率上卻能夠優于DTW和DDTW的錯誤率。然而,在有些數據集如ECG,不管w取何值,新方法的分類錯誤率都是低于另外兩者的錯誤率。

圖5 分類錯誤率和時間序列分段聚合長度的關系Fig.5 The relationship of classification and the length of time series piecewise aggregation

從圖5中發現,不同數據集呈現出來的性能不一樣。由于分段聚合近似采用的是均值作為替代值,這并非最佳的降維方法,因此,若數據集中的序列數據波動較多、振幅較大,會對降維效果造成一定的影響。如圖6所示。

圖6 實驗數據集序列示例Fig.6 The example of datasets

圖6分別給出的是訓練集Sony AIBO Robot Surface和ECG的3條時間序列,可以觀察到子圖(a) 的序列在整體波動上比子圖(b)的多。因此,影響了序列降維的效果,給后續的相似性度量造成了一定的影響。由于數據集的特性會給本方法造成一定的影響,因此,本方法適合波動較為平緩的數據集,并會取得較好的效果。

圖7描述的是PADD的平均分類錯誤率以及最小分類錯誤率與DTW、DDTW的錯誤率的比較,為了便于直觀比較,數值均經過歸一化后取值范圍為[0,0.5],數值偏向方表示對應方法的錯誤率較大。

圖7 分類錯誤率的比較Fig.7 Comparison of classification error rates

圖7中子圖(a)、(b)分別描述PADD的平均錯誤率與DTW、DDTW的錯誤率比較,子圖(c)、(d)分別描述PADD最小錯誤率與DTW、DDTW錯誤率比較。結果分析表明,不管從平均錯誤率還是最小錯誤率角度來比較,PADD錯誤率所對應的坐標縱軸值相對較小,使其偏向于DTW和DDTW所代表的橫軸值較大,大多散點都偏向于DTW和DDTW,故說明PADD具有較小的分類錯誤率,進行驗證了本文新方法在時間序列度量中的有效性和優越性。

3.2時間效率分析

本實驗使用“留一法”(leave-one-out)求解參數,參數確定后立即構建度量距離公式。特征序列長度越長,計算所要消耗的時間越大,且成數量級增長。如圖8所示,描述了數據集Sony AIBO Robot Surface和ECG在不同w條件下PADD時間效率與DTW、DDTW的時間效率比較。

圖8 3種方法的時間效率Fig.8 Time efficiency of the three methods

4結束語

針對動態彎曲方法對時間序列數據相似性度量質量和效率的局限性,以及從簡化模型的角度考慮,本文提出一種基于分段聚合近似和數值導數的動態時間彎曲方法,結合分段聚合近似算法以及數值求導將時間序列進行轉換,得到符合要求的特征序列,構建新的距離公式對新的特征序列進行相似性度量。數值實驗結果表明,與傳統的動態時間彎曲度量方法相比,在分類效果和計算性能上具有一定的優勢。通過實驗發現,針對波動較為平緩的數據進行挖掘,將會取得更好的效果。特別地,在數據壓縮量較大的低維空間下,新方法在保證分類質量的前提下,能夠獲得較好的時間效率。另外,參數選擇是一個至關重要的步驟,不同參數對實驗結果具有一定的影響,將來需要對參數選擇做進一步探討和研究。

參考文獻:

[1]韓敏, 許美玲, 王新迎. 多元時間序列的子空間回聲狀態網絡預測模型[J]. 計算機學報, 2014, 37(11): 2268-2275.

HAN Min, XU Meiling, WANG Xinying. A multivariate time series prediction model based on subspace echo state network[J]. Chinese journal of computers, 2014, 37(11): 2268-2275.

[2]MARSZA EK A, BURCZY SKI T. Modeling and forecasting financial time series with ordered fuzzy candlesticks[J]. Information sciences, 2014, 273: 144-155.

[3]ZAMORA M, LAMBERT A, MONTERO G. Effect of some meteorological phenomena on the wind potential of Baja California[J]. Energy procedia, 2014, 57: 1327-1336.

[4]GRAVIO G D, MANCINI M, PATRIARCA R, et al. Overall safety performance of air traffic management system: forecasting and monitoring[J]. Safety science, 2015, 72: 351-362.

[5]李海林, 郭韌 萬校基. 基于特征矩陣的多元時間序列最小距離度量方法[J]. 智能系統學報, 2015, 10(3): 442-447.

LI Hailin, GUO Ren, WAN Xiaoji. A minimum distance measurement method for multivariate time series based on the feature matrix[J]. CAAI transactions on intelligent systems, 2015, 10(3): 442-447.

[6]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.

[7]IZAKIAN H, PEDRYCZ W, JAMAL I. Fuzzy clustering of time series data using dynamic time warping distance[J]. Engineering applications of artificial intelligence, 2015, 39: 235-244.

[8]ZHANG Zheng, TANG Ping, DUAN Rubing. Dynamic time warping under pointwise shape context[J]. Information sciences, 2015, 315: 88-101.

[9]李正欣, 張鳳鳴, 李克武. 基于DTW的多元時間序列模式匹配方法[J]. 模式識別與人工智能, 2011, 24(3): 425-430.

LI Zhengxin, ZHANG Fengming, LI Kewu. DTW based pattern matching method for multivariate time series[J]. Pattern recognition and artificial intelligence, 2011, 24(3): 425-430.

[10]LI Hailin. Asynchronism-based principal component analysis for time series data mining[J]. Expert systems with applications, 2014, 41(6): 2842-2850.

[11]KEOGH E, PAZZANI M J. Derivative dynamic time warping[C]//Proceedings of the 1st SIAM International Conference on Data Mining. Chicago, IL, USA, 2001: 1-11.

[12]JEONG Y S, JAYARAMAN R. Support vector-based algorithms with weighted dynamic time warping kernel function for time series classification[J]. Knowledge-based systems, 2015, 75: 184-191.

[14]李海林, 楊麗彬. 時間序列數據降維和特征表示方法[J]. 控制與決策, 2013, 28(11): 1718-1722.

LI Hailin, YANG Libin. Method of dimensionality reduction and feature representation for time series[J]. Control and decision, 2013, 28(11): 1718-1722.

[16]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. Berlin Heidelberg: Springer, 1993, 730: 69-84.

[17]WANG Hongfa. Clustering of hydrological time series based on discrete wavelet transform[J]. Physics procedia, 2012, 25: 1966-1972.

[18]WENG Xiaoqing, SHEN Junyi. Classification of multivariate time series using two-dimensional singular value decomposition[J]. Knowledge-based systems, 2008, 21(7): 535-539.

[19]LIN J, KEOGH E, WEI Li, et al. Experiencing SAX: a novel symbolic representation of time series[J]. Data mining and knowledge discovery, 2008, 15(2): 107-144.

[20]KEOGH E, XI X, WEI L, et al. The UCR time series classification/clustering homepage[EB/OL]. 2007. http://www.cs.ucr.edu/~eamonn/time_series_data/.

李海林,男,1982年生,副教授,博士,主要研究方向為數據挖掘與決策支持,主持國家自然科學基金、省部級基金多項,發表學術論文30余篇,其中被SCI檢索11篇,EI檢索10余篇。

梁葉,女,1992年生,碩士研究生,主要研究方向為數據挖掘與金融數據分析。

Dynamic time warping based on piecewise aggregate approximation and data derivatives

LI Hailin, LIANG Ye

(Department of Information Management, Huaqiao University, Quanzhou 362021, China)

Abstract:Dynamic time warping (DTW) is often used to measure the similarity of time series data; however, it has efficiency and quality limitations. In this study, a novel DTW method combining piecewise aggregate approximation (PAA) and derivatives is proposed to measure the similarity of time series. The dimensionality of the time series data was effectively reduced by PAA, and the feature sequence was transformed into new sequences by combining the numerical derivatives after the dimensionality reduction. Furthermore, the DTW design corresponded to the similarity measurement method of the feature sequence. The experimental results demonstrate that the proposed method is superior because it has better measurement quality, obtains a better classification effect in time series data mining, and has high efficiency in lower dimensional spaces.

Keywords:dynamic time warping; time series; piecewise aggregate approximation; numerical derivative; similarity measure; classification;dimensionality reduction; feature representation

作者簡介:

中圖分類號:TP301

文獻標志碼:A

文章編號:1673-4785(2016)02-0249-08

通信作者:李海林. E-mail:hailin@hqu.edu.cn.

基金項目:國家自然科學基金項目(61300139);福建省中青年教育科研基金項目(JAS14024);華僑大學中青年教師科研提升計劃項目(ZQN-PY220).

收稿日期:2015-07-24. 網絡出版日期:2016-03-15.

DOI:10.11992/tis.201507064

網絡出版地址:http://www.cnki.net/kcms/detail/23.1538.TP.20160315.1239.012.html

猜你喜歡
分類方法
分類算一算
垃圾分類的困惑你有嗎
大眾健康(2021年6期)2021-06-08 19:30:06
學習方法
分類討論求坐標
數據分析中的分類討論
教你一招:數的分類
用對方法才能瘦
Coco薇(2016年2期)2016-03-22 02:42:52
四大方法 教你不再“坐以待病”!
Coco薇(2015年1期)2015-08-13 02:47:34
賺錢方法
給塑料分分類吧
主站蜘蛛池模板: 国产精品吹潮在线观看中文| 午夜福利网址| 国产99精品久久| vvvv98国产成人综合青青| 成人在线综合| 三区在线视频| 欧美日韩一区二区三| 三区在线视频| 国产精品9| 国产又黄又硬又粗| 国产毛片基地| 亚洲成人精品在线| 国产在线一区视频| 午夜激情婷婷| 97无码免费人妻超级碰碰碰| 正在播放久久| 沈阳少妇高潮在线| 九九九久久国产精品| 欧美色图第一页| 18禁黄无遮挡免费动漫网站| 成人韩免费网站| jizz亚洲高清在线观看| 亚洲精品色AV无码看| 波多野结衣第一页| 国产成人一区在线播放| 国产欧美日本在线观看| 美女毛片在线| 久久熟女AV| 毛片一区二区在线看| 91青青在线视频| 国模视频一区二区| 99热亚洲精品6码| 色综合激情网| 乱人伦视频中文字幕在线| 看国产毛片| 日本一区中文字幕最新在线| 欧美日韩国产在线人| 99热国产这里只有精品无卡顿"| 国产一区成人| 国产丝袜丝视频在线观看| 日韩高清无码免费| 亚洲欧洲日本在线| 国产一区二区免费播放| 欧美无专区| 在线高清亚洲精品二区| 伊人激情久久综合中文字幕| 久久亚洲中文字幕精品一区| 日本91在线| 青草视频免费在线观看| 青青青草国产| 国产在线精品网址你懂的| 国产二级毛片| 国产激情无码一区二区APP| 亚洲国产91人成在线| 91午夜福利在线观看精品| 亚洲一区二区约美女探花| 无遮挡一级毛片呦女视频| 中文字幕乱妇无码AV在线| 亚洲精品无码成人片在线观看| A级全黄试看30分钟小视频| 91九色国产在线| 岛国精品一区免费视频在线观看| 狠狠干综合| 一本色道久久88| 欧美亚洲国产精品久久蜜芽| 久久无码av三级| 国产美女无遮挡免费视频网站| 中文天堂在线视频| 欧美va亚洲va香蕉在线| 国产亚洲精品资源在线26u| 19国产精品麻豆免费观看| 国产91小视频| 久久91精品牛牛| 麻豆精品久久久久久久99蜜桃| 亚洲区一区| 国产男女XX00免费观看| 国产永久无码观看在线| 国产AV无码专区亚洲精品网站| 自拍偷拍欧美日韩| 视频国产精品丝袜第一页| 亚洲中文字幕无码爆乳| 亚洲日韩图片专区第1页|