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

基于條件熵擴維的多變量混沌時間序列相空間重構*

2011-09-28 07:06:38張春濤馬千里姜友誼
物理學報 2011年2期
關鍵詞:系統

張春濤 馬千里 彭 宏 姜友誼

1)(重慶三峽學院數學與計算機科學學院,萬州 404000)

2)(華南理工大學計算機科學與工程學院,廣州 510006)

基于條件熵擴維的多變量混沌時間序列相空間重構*

張春濤馬千里2)彭 宏2)姜友誼1)

1)(重慶三峽學院數學與計算機科學學院,萬州 404000)

2)(華南理工大學計算機科學與工程學院,廣州 510006)

(2010年4月12日收到;2010年5月22日收到修改稿)

提出一種多變量混沌時間序列相空間重構的條件熵擴維方法.首先使用互信息法求解每個變量的時間延遲,其次按條件熵最大原則逐步擴展相空間的嵌入維數,使得重構坐標從低維到高維的轉換保持較強的獨立性,最終的重構相空間具有較低的冗余度,為多變量時間序列的預測構造了有效的模型輸入向量.通過對幾個經典多變量混沌時間序列進行數值實驗,結果表明該方法比單變量預測和已有多變量預測方法具有更好的預測效果,說明了該重構方法的有效性.

多變量混沌時間序列,相空間重構,條件熵,神經網絡預測

PACS:05.45.Pq,05.45.Tp

1.引 言

復雜系統中任一分量的演化都可由系統的其他分量所決定,因此任一分量的發展過程都包含其他分量的信息,即通過該分量發展的時間序列重構相空間能充分反映原系統的動力學特征,這就是單變量混沌時間序列相空間重構[1—8].在現實狀況下由于無法從單變量的時間序列獲全系統信息,以及在分析生物學系統時,由于其結構和動力系統的復雜性,很難由一維時間序列通過延遲嵌入法重構出相空間[9],因此不能保證實際系統中任何給定的單變量時間序列足以重構原系統[10,11].而多變量時間序列由于包含了更豐富完整的系統信息,故能重構出更為準確的相空間,因此多變量混沌時間序列的研究得到越來越多的關注[12—18].

目前,對多變量時間序列已取得了許多研究成果.王海燕等[14]提出由單變量推廣的虛假鄰近點法計算多變量時間序列的相空間嵌入維數,但其缺少每次擴維選擇標準,使重構相空間具有不定性. Broomhead等[15]提出奇異值分解法求取多變量時間序列的嵌入維數.Cao等[10]對多變量的相空間重構問題進行理論和應用研究,并論證了其優越性;樊重?。?6]討論了多觀測變量狀態空間重構技術,為研究多變量時間序列的相空間重構提供了理論基礎.

信息熵是系統整體不確定性的一種量度,在混沌時間序列相空間重構中有重要應用.田玉楚等[19]系統地論述了在不變測度意義下,將混沌系統視為一類特殊的隨機系統加以處理;肖方紅等[6]對一維時間序列使用條件熵求取相空間的嵌入維數;我們也研究了兩個重構參數的信息熵模型[8].這些研究都集中在單變量混沌時間序列上,對多變量混沌系統,主要用互信息法求解單個變量的延遲時間.

本文沿用單變量時間序列建模的思想,結合信息熵和神經網絡對多變量時間序列進行相空間重構和預測,提出了多變量混沌時間序列的條件熵擴維方法(CEED).首先使用互信息法求解單個變量的時間延遲;然后使用條件熵擴維方法,按條件熵最大原則選擇重構分量,把低維相空間逐漸擴展到高維相空間,使得最終的相空間具有較少的冗余度.為了檢驗該重構方法的有效性,最后使用神經網絡對2個典型的多變量混沌系統進行仿真實驗,數值實驗表明在該重構方法上進行的多變量混沌時間序列預測能取得比單變量和已有的多變量預測更好的效果,說明CEED方法重構了高質量的相空間.

2.多變量混沌時間序列相空間重構

在一定程度上單變量時間序列是多變量時間序列的特殊情況,因此可將單變量的延遲重構應用到多變量中.對單變量時間序列x(n)(n=1,2,…,N),在確定嵌入維數 m和時間延遲 τ后,重構相空間

設復雜系統的N維多觀測變量序列為 X1,X2,…,XK,其中 K表示觀測變量個數,Xi=[xi(1),xi(2),…,xi(N)]表示第i個變量的時間序列.如果每個變量都選擇恰當的時間延遲τi和嵌入維數mi(i=1,2,…,K),則多變量時間序列重構的相空間

其中τi和mi是重構相空間成敗的關鍵,目前還沒有一個統一的選擇標準,下面使用信息論知識來選擇這K對重構參數.

2.1.多變量時間序列延遲時間的選擇

設兩個離散型隨機變量X∈{x1,x2,…,xn}和Y∈{y1,y2,…,yk},其 先 驗 概 率 分 別 為 {p (xi)}i=1,2,…n和{p(yj)}j=1,2,…,k.信息熵可定義為

它描述了隨機變量X的不定性.類似地可以定義聯合熵

其中p(xi,yj)為聯合概率.

變量X和Y之間的互信息定義為

根據(7)式可推出時間序列和其延遲時間序列的互信息

根據互信息法[5],計算出互信息函數 I(τi)在不同延遲τi下的值,取I(τi)的第一個局部極小值點對應的時間τi為第 i個變量時間序列 Xi的延遲時間.

2.2.多變量時間序列嵌入維數的確定

一般認為,如果重構相空間的維數足夠大,相空間就能夠刻畫出原系統的吸引子.對多變量時間序列,一個很自然的想法就是對每個變量的時間序列都分別求一個嵌入維數mi,然后將所有變量的重構相空間進行簡單迭加形成多變量時間序列的相空間如(2)式.由于混沌系統任一變量的時間序列包含其他變量的信息,即使單個變量的重構相空間無冗余信息,經過簡單迭加而成的多變量時間序列的重構相空間也包含大量的冗余信息,不利于預測.如何減少這些冗余信息,是確定多變量時間序列嵌入維數的關鍵.主成分分析(PCA)方法[17,18]為我們提供了一種解決途徑,但是PCA方法只能消除各重構分量間的線性關系,對非線性的時間序列作用有限.下面根據條件熵能反映變量之間的獨立性的特征,通過條件熵消除各重構分量間的非線性關系,使重構相空間包含較少冗余信息,以達到降低相空間維數的目的.

記多變量重構相空間(2)式的延遲坐標向量的各分量為

相空間重構也就是用重構分量 X11,X12,…,XK1,XK2,…,XKmK近似表征原系統的各變量,同時X11,X12,…,XK1,XK2,…,XKmK也可以看成是重構相空間的坐標,因此原系統表現出來的混沌吸引子也將在重構分量X11,X12,…,XK1,XK2,…,XKmK的重構相空間中表現出來.這樣隨著重構分量的逐漸增多,重構向量就能越來越完整地表征原系統的動力學特征.當重構的分量個數增加到一定值時,這時重構的分量包含了原系統充分多的動力學信息,進一步增加重構分量不會增加原系統的動力學信息,只會攜帶冗余信息.因此可以設想分量X11,X21,…,XK1的熵H(X11,X21,…,XK1)和各分量之間的條件熵的值隨著重構分量個數的增加由大到小逐漸減小,直至重構分量包含充分多的原系統的動力學信息,進一步增加重構分量個數,這時條件熵將出現一個較明顯的轉折點,隨后的條件熵將趨于平穩.如果取條件熵出現較為明顯的轉折且趨于平穩處的重構分量個數作為多變量時間序列重構相空間的嵌入維數,可以看出此時重構向量已包含原系統充分多的動力學信息,應該是一個合適的嵌入維數[6].

多變量時間序列的重構向量中增加一個重構分量有諸多選擇,即在已知重構分量為 X11,X12,X21的條件下再增加一個分量,有 K種不同選擇:X13,X22,X31,X41,…,XK1.對這 K個條件熵做怎樣的取舍?由上所述,多變量時間序列重構相空間應具有較少的冗余信息,即重構分量之間應最大限度地相互獨立.因此增加一個重構分量的最佳選擇應是K個條件熵的最大值所對應的分量.因為條件熵最大,表明增加的這個分量和原有分量的獨立性最大,這樣不斷增加重構分量到整個相空間就能夠保證有最大的獨立性.下面給出從低維到高維的擴維(增加相空間重構分量個數)算法步驟:

步驟1 對每個變量的時間序列使用互信息法求時間延遲τi,i=1,2,…,K;

步驟2 初始化:令Y=[X11,X21,…,XK1]為初始的重構向量,每個變量的初始嵌入維數mi=1,i= 1,2,…,K;

步驟3 計算條件熵并擴維:

3)擴維mk=mk+1,此時的重構向量 Y=[Y,Xkmk];

4)重復1)—3)直至最大條件熵趨于平穩.

步驟4 循環結束后的重構向量空間就是多變量時間序列的重構相空間,此時的mi即為第i觀測變量的嵌入維數,重構相空間的嵌入維數m=m1+ m2+…+mK.

對條件熵按(6)式轉化為聯合熵來計算,本文采用矩陣標識法[6]計算聯合熵.

3.多變量混沌時間序列仿真實驗

為了驗證 CEED的有效性,我們做了 Henon,Lorenz系統的單步和多步預測實驗,并與現有的單變量和多變量預測方法進行比較.

3.1.預測模型的建立

假設重構的狀態空間如(2)式,相空間中預測的基本思想是:根據歷史狀態數據集建立第 i時刻狀態Y(i)到第i+h時刻狀態Y(i+h)的映射Fh,然后由映射Fh和當前n時刻狀態Y(n)去推斷未來n+h時刻狀態Y(n+h)的估計值Z(n+h),其中h為預測步長.根據嵌入定理,狀態映射Fh是存在的,即在m維的重構相空間存在光滑函數fh:Rm→R.對多變量混沌時間序列,因此有預測模型

其中,zi(n+h)為模型對未來(n+h)時刻第i個變量狀態xi(n+h)的預測值(i=1,2,…,K),h為預測步長.由于混沌系統的高度非線性,要確定函數fih的解析形式是困難的,我們采用神經網絡的學習原理來逼近這些映射函數.結合相空間重構和神經網絡建立預測模型的結構如圖1所示.

圖1 多變量時間序列預測模型結構

本文采用的神經網絡是 Elman遞歸神經網絡[20],它能夠通過內部的遞歸連接反映非線性動力系統的特征,在混沌時間序列的預測中有重要應用[2,3].

3.2.仿真實驗

為了評價預測效果采用均方誤差(MSE)、均方根誤差(RMSE)、標準化均方誤差(NMSE)以及絕對誤差e(t)作為誤差評價標準:

3.2.1.Henon混沌時間序列

法國天文學家Henon于1976年提出一個簡單的二維平方映像如下:

通過計算得出變量x,y的3000個點的數據集,使用互信息法分別求解x和y的延遲時間如圖2和圖3所示,可知τ1=15,τ2=14,然后通過條件熵擴維如圖4所示.

由圖4可知,在嵌入維數為5時條件熵出現一個明顯的轉折,當嵌入維數增加到6后條件熵幾乎不發生變化,所以我們取Henon重構相空間的嵌入維數為6,再根據每次的重構分量選擇m=m1+m2=4+2=6.

通過Elman神經網絡按(9)式對變量x進行預測.前2000個數據點進行訓練,后1000個數據進行預測,其中單步預測結果與單變量預測方法:在線最小二乘支持向量機(LS-SVM)[21],最小二乘支持向量(LS-SVD)[22],模糊樹Ⅰ(FTⅠ)[23],FTⅡ[23],模糊聚類[24],基于自適應網絡的模糊推理系統(ANFIS)[25],局域離散余弦變換領域Vloterra(LDCTDV)[26],模糊邊界模塊化神經網絡 FBMNN[2]以及多變量預測方法的比較如表1所示,其中非線性獨立分量分析(NICA)為文獻[27]所取的多變量相空間結構:m1=2,τ1=1;m2=2,τ2=1,然后使用Elman神經網絡的預測.

圖2 x的互信息

圖3 y的互信息

圖4 條件熵與嵌入維數的關系

從表1可以看出CEED取得比單變量預測和多變量預測更好的預測效果.單步預測效果和預測誤差分別如圖5和圖6所示.

從圖6可以看出,對Henon混沌時間序列的單步預測誤差達到10-5數量級,比文獻[28]的10-2高3個數量級,說明使用多變量時間序列預測 Henon是高效.

表1 單步預測

圖5 單步預測效果

圖6 單步預測誤差

在進行的多步預測中,預測結果與正交最小二乘徑向基函數(RBF(OLS))[29],局域線性模型樹(LoliMoT)[29],FBMNN和NICA的比較如表2所示,其中節點數為隱含層節點數.

表2 Henon預測結果的NMSE比較

由表2可以看出CEED使用較少的隱含層節點,從而降低了計算復雜性,卻能取得更好的預測效果.

從表1和表2可以看出,Henon映像是一個最簡單的二維混沌系統,使用單變量預測方法和已有的多變量相空間重構方法進行單步預測大多能取得較好的預測效果,但多步預測效果一般較差,而使用CEED的多變量預測方法在單步和多步預測上都取得了很好的預測效果.FBMNN,NICA和CEED都使用Elman進行預測,CEED上進行的預測取得了最好的預測效果,說明了 CEED重構方法對Henon映像的有效性.

3.2.2.洛倫茲方程

洛倫茲方程如下:

用四階Runge-Kutta方法求解方程(15),取積分步長h=0.02,得到變量x,y的兩個7500個點的數值解,去除前5000個點的暫態過程,最后考慮兩個觀察變量x,y的2500個點的數據集,將其歸一化到[0,1]區間.先通過互信息法分別求出x和y的延遲時間如圖7和圖8所示,由圖7,圖8可知,τ1=8,τ2=8.

然后通過條件熵擴維如圖9所示,由圖9可知,當嵌入維數增加到6后條件熵不再發生變化,再根據擴維分量選擇有m=m1+m2=2+4=6.

圖7 x的互信息

圖8 y的互信息

圖9 條件熵與嵌入維數的關系

通過Elman神經網絡按(9)式對變量x進行預測.前2000個數據點進行訓練,后500個數據進行預測,其中單步預測結果與單變量時間序列預測和多變量預測的比較如表3所示,其中IFNNMU為單變量推廣的虛假鄰近點法確定的相空間結構:m1= 2,τ1=8;m2=1,τ2=8,然后使用Elman神經網絡的預測.

從表3可以看出,對Lorenz系統CEED能取得比已有的單變量時間序列預測和多變量預測更好的單步預測效果.

為了與多變量預測方法[30]進行比較,將歸一化數據還原為原始數據,根據eRMSE的定義有,對Lorenz系統歸一化數據后的 eRMSE是原始數據的 1/40,即CEED原始數據的 eRMSE≤2×10-5,比文獻[30]的0.063高3個數量級.CEED的嵌入維數是6,而文獻[30]的嵌入維數為11,這樣CEED能簡化神經網絡結構,更利于實現.

單步預測效果如圖10所示,預測誤差如圖11所示.

圖10 單步預測效果

圖11 單步預測誤差

同時我們也進行了多步預測實驗,預測步長取為6和8,與共進化遞歸神經網絡(CERNN)[3]等相關預測方法的比較如表4和表5所示.

表4 6步預測

表5 8步預測

從表4和表5可以看出,對Lorenz系統,在CEED上進行的多步預測精度也高于單變量預測和多變量預測方法.從表 3—表 5可看出,FBMNN,CERNN,IFNNMU和CEED都使用了Elman網絡進行預測,在CEED重構相空間上進行的預測取得了更好的預測效果,說明CEED重構方法對Lorenz系統的有效性.

從表1—表5可以看出,對多維混沌系統,不管是與單變量時間序列預測方法,還是與已有的多變量預測方法相比,CEED重構相空間上進行單步和多步預測都能取得更好的預測效果,說明CEED重構的多維混沌系統的相空間是有效性.

4.結 論

本文提出一種多變量混沌時間序列相空間重構的新方法.該方法通過條件熵擴張相空間維數,按條件熵最大選擇重構分量,使得相空間從低維到高維的轉變過程中既保證重構坐標的獨立性又保持系統的動力學特征.因此通過CEED確定的重構相空間既能充分包含系統信息,又能有效避免冗余信息的產生.最后在CEED上進行多變量混沌時間序列的單步和多步預測,都取得了比使用單變量預測方法和已有的多變量預測方法更好的預測效果,表明CEED重構方法的有效性.

[1]Peng S G,Yu S M 2009 Chin.Phys.B 18 3758

[2]Ma Q L,Zheng Q L,Peng H,Qin J W 2009 Acta Phys.Sin.58 1410(in Chinese)[馬千里、鄭啟倫、彭 宏、覃姜維2009物理學報58 1410]

[3]Ma Q L,Zheng Q L,Peng H,Zhong T W,Qin J W 2008 Chin. Phys.B 17 536

[4]Zhang S Q,Jia J,Gao M,Han X 2010 Acta Phys.Sin.59 1576 (in Chinese)[張淑清、賈 健、高 敏、韓 敘 2010物理學報59 1576]

[5]Fraser A M,Swinney H L 1986 Phys.Rev.A 33 11340

[6]Xiao F H,Yan G R,Han Y H 2005 Acta Phys.Sin.54 550(in Chinese)[肖方紅、閻桂榮、韓宇航2005物理學報54 550]

[7]Yan H,Wei P,Xiao X C 2009 Chin.Phys.B 18 3287

[8]Zhang C T,Ma Q L,Peng H 2010 Acta Phys.Sin.7623(in Chinese)[張春濤、馬千里、彭 宏2010物理學報7623]

[9]Kantz H,Schreiber T 1997 Nonlinear Time Series Analysis (Cambridge:Cambridge University Press)p7

[10]Cao L Y,Mees A,Judd K 1998 Physica D 121 75

[11]Rombowts S A R B,Keunen R WM,Stom C J 1995 Phys.Lett. A 202 352

[12]Lu S,Wang H Y 2006 Acta Phys.Sin.55 572(in Chinese)[盧 山、王海燕2006物理學報55 572]

[13]Zhang Y,Guan W 2009 Acta Phys.Sin.58 756(in Chinese)[張 勇、關 偉2009物理學報58 756]

[14]Wang H Y,Sheng Z H,Zhang J 2003 J.Southeast Univ. (Natural Science Edition)33 115(in Chinese)[王海燕、盛昭瀚、張 進2003東南大學學報(自然科學版)33 115]

[15]Broomhead D S,King G P 1986 Physica D 20 217

[16]Fan C J 2009 J. UniversityofShanghaiforScienceand Technology 31 273(in Chinese)[樊重俊2009上海理工大學學報31 273]

[17]Han M,Fan M M,Xi J H 2005 Lecture Notes in Computer Science 3497 618

[18]Han M,Wang Y 2009 Expert Systems with Applications 36 1280

[19]Tian Y C,Fu X T,Lu Y Z 1993 Control and Decision 8 345(in Chinese)[田玉楚、符雪桐、呂勇哉1993控制與決策8 345]

[20]Elman J L 1990 Cognitive Science 14 179

[21]Ye M Y,Wang X D,Zhang H R 2005 Acta Phys.Sin.54 2568 (in Chinese)[葉美盈、汪曉東、張浩然 2005物理學報 54 2568]

[22]Ren R,Xu J,Zhu S H 2006 Acta Phys.Sin.55 555(in Chinese)[任 韌、徐 進、朱世華2006物理學報55 555]

[23]Mao J Q,Yao J,Ding H S 2009 Acta Phys.Sin.58 2220(in Chinese)[毛劍琴、姚 健、丁海山2009物理學報58 2220]

[24]Chiu S 1994 J.Intell.Fuzzy Syst.2 267

[25]Jang J S R 1993 IEEE Trans.Syst.Man 23 665

[26]Zhang J S,Li H C,Xiao X C 2005 Chin.Phys.14 49

[27]Jan S,Guatavo D 1997 Physics D 108 335

[28]Zhang J F,Hu S S 2008 Acta Phys.Sin.57 2708(in Chinese)[張軍峰、胡壽松2008物理學報57 2708]

[29]Gholpour A,Araabi B N,Lucas C 2006 Neural Proc.Lett.24 217

[30]Han M,Wei R 2008 14th International Conference of Neural Information Processing Kitakyushu,Japan,November 13—16,2007 p415

PACS:05.45.Pq,05.45.Tp

Multivariate chaotic time series phase space reconstruction based on extending dimension by conditional entropy*

Zhang Chun-TaoMa Qian-Li2)Peng Hong2)Jiang You-Yi1)
1)(College of Mathematic and Computer Science,Chongqing Three Gorges University,Wanzhou 404000,China)
2)(College of Computer and Engineering,South China University of Technology,Guangzhou 510006,China)

12 April 2010;revised manuscript

22 May 2010)

For multivariate chaotic time series,a method of conditional entropy extending dimension(CEED)in the reconstructed phase space is proposed.First,the delay time of any variable time series is selected by mutual information method,and then the embedding dimension of phase space is extended by the conditional entropy.This method can ensure the independence of reconstructed coordinates from low space to high space and eliminate the redundancy of phase space,because the largest condition entropy is choosen.The effective input vector for the prediction of multivariate time series is given.Simulations of the Lorenz system and Henon system show that the neural network predictions of multivariate time series are much better than the prediction of univariate and existing multivariate.Therefore,CEED is effective for multivariate chaotic systems.

multivariate chaotic time series,phase space reconstruction,conditional entropy,neural network prediction

*廣東省自然科學基金(批準號:9451064101003233)、華南理工大學中央高校基本科研業務費專項資金(批準號:2009ZM0125,2009ZM0189,2009ZM0255)和重慶三峽學院重點項目(批準號:10ZD-16)資助的課題.

*Project supported by the Natural Science Foundation of Guangdong Province,China(Grant No.9451064101003233),the Fundamental Research Funds for the Central Universities,South China University of Technology(Grant Nos.2009ZM0125,2009ZM0189,2009ZM0255)and the Key Program of the Chongqing Three Gorges University,China(Grant No.10ZD-16).

猜你喜歡
系統
Smartflower POP 一體式光伏系統
工業設計(2022年8期)2022-09-09 07:43:20
WJ-700無人機系統
ZC系列無人機遙感系統
北京測繪(2020年12期)2020-12-29 01:33:58
基于PowerPC+FPGA顯示系統
基于UG的發射箱自動化虛擬裝配系統開發
半沸制皂系統(下)
FAO系統特有功能分析及互聯互通探討
連通與提升系統的最后一塊拼圖 Audiolab 傲立 M-DAC mini
一德系統 德行天下
PLC在多段調速系統中的應用
主站蜘蛛池模板: 91精品福利自产拍在线观看| 18禁黄无遮挡免费动漫网站| 久久天天躁狠狠躁夜夜躁| 亚洲首页在线观看| 国产特级毛片aaaaaaa高清| 国产欧美日韩18| 色婷婷在线播放| 天天综合网色中文字幕| 国产人碰人摸人爱免费视频| 欧美无遮挡国产欧美另类| 亚洲大学生视频在线播放| 日韩欧美国产精品| 午夜日本永久乱码免费播放片| 美女视频黄频a免费高清不卡| 在线观看亚洲天堂| 国产色偷丝袜婷婷无码麻豆制服| 国内精品九九久久久精品| 日a本亚洲中文在线观看| 亚洲另类第一页| 国产欧美日韩va另类在线播放| 永久免费无码成人网站| 亚洲人成影院午夜网站| 欧美日韩综合网| 亚洲三级a| 中文字幕免费在线视频| 国产噜噜噜视频在线观看| 性欧美久久| 很黄的网站在线观看| 久久99热66这里只有精品一| 国内精品一区二区在线观看| 精品人妻一区二区三区蜜桃AⅤ| 日本人妻丰满熟妇区| 国产区精品高清在线观看| 69av在线| 成人免费视频一区二区三区| 国产欧美精品一区aⅴ影院| 就去色综合| 国产视频 第一页| 国产电话自拍伊人| 丝袜美女被出水视频一区| 免费亚洲成人| a毛片免费在线观看| 国产无码网站在线观看| 亚洲天堂视频在线播放| 色窝窝免费一区二区三区| 99久久精品视香蕉蕉| 国产毛片不卡| 日韩av手机在线| 久久成人国产精品免费软件| 亚洲精品自拍区在线观看| 国内精品九九久久久精品| 91在线日韩在线播放| 免费Aⅴ片在线观看蜜芽Tⅴ| 国产夜色视频| av尤物免费在线观看| 日本国产一区在线观看| 国产成人欧美| 国产欧美日韩91| 国产男女免费完整版视频| 亚洲中文字幕在线精品一区| 国产成在线观看免费视频| 国产欧美视频综合二区 | 国产美女精品一区二区| 国产在线观看第二页| 国产在线观看91精品亚瑟| 制服丝袜无码每日更新| 国产99久久亚洲综合精品西瓜tv| 免费A级毛片无码免费视频| 大香网伊人久久综合网2020| 国产精品任我爽爆在线播放6080| 国产精品亚洲天堂| 无码丝袜人妻| 福利在线不卡| 精品欧美视频| 激情影院内射美女| 亚洲九九视频| 一区二区三区国产精品视频| 伊人AV天堂| 色噜噜狠狠狠综合曰曰曰| 亚洲性日韩精品一区二区| 欧美成一级| 亚洲精品少妇熟女|