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

基于約束HMM的變點(diǎn)檢測(cè)算法①

2017-06-07 08:24:04何振峰
關(guān)鍵詞:檢測(cè)模型

莊 玉,何振峰

(福州大學(xué) 數(shù)學(xué)與計(jì)算機(jī)科學(xué)學(xué)院,福州 350108)

基于約束HMM的變點(diǎn)檢測(cè)算法①

莊 玉,何振峰

(福州大學(xué) 數(shù)學(xué)與計(jì)算機(jī)科學(xué)學(xué)院,福州 350108)

時(shí)間序列的變點(diǎn)分析在現(xiàn)今社會(huì)各個(gè)領(lǐng)域中都有著廣泛的應(yīng)用.針對(duì)時(shí)間序列進(jìn)行變點(diǎn)分析中要求變點(diǎn)狀態(tài)需要連續(xù)持續(xù)一定的時(shí)間的應(yīng)用背景,提出了一種結(jié)合狀態(tài)最短連續(xù)長(zhǎng)度約束的隱馬爾可夫模型.給出了約束Baum-Welch訓(xùn)練算法和約束Viterbi狀態(tài)提取算法.應(yīng)用在仿真數(shù)據(jù)和GNP數(shù)據(jù)集的實(shí)驗(yàn)表明,結(jié)合狀態(tài)最短連續(xù)長(zhǎng)度約束的HMM相比于一般HMM在時(shí)間序列變點(diǎn)檢測(cè)中效率較高.

變點(diǎn)檢測(cè);約束隱馬爾可夫模型;時(shí)間序列分割

時(shí)間序列是隨時(shí)間次序而變化的一系列數(shù)據(jù),是一類多維的復(fù)雜類型數(shù)據(jù)[1].當(dāng)所觀察的時(shí)間序列跨越時(shí)間越長(zhǎng)時(shí),形成的時(shí)間序列的隨機(jī)變量會(huì)由于某種條件的變化,從某時(shí)間點(diǎn)開始不再服從原來的分布,即出現(xiàn)了變點(diǎn).隨著變點(diǎn)問題應(yīng)用的愈加廣泛,在一些文獻(xiàn)中有許多同義詞,包括了結(jié)構(gòu)斷點(diǎn)[2],時(shí)間序列分割[3],和檢測(cè)“異常點(diǎn)”[4].

時(shí)間序列的變點(diǎn)檢測(cè)方法可以很好的用HMM來建模,其中時(shí)間序列數(shù)據(jù)就是HMM中的輸出符號(hào),可以通過時(shí)間序列數(shù)據(jù)來檢測(cè)所處的狀態(tài),判斷是否出現(xiàn)變點(diǎn).一般的HMM對(duì)隱狀態(tài)鏈沒有任何限制,即在隱狀態(tài)鏈中可以自由的從一個(gè)狀態(tài)轉(zhuǎn)移到其他任何一個(gè)狀態(tài).這樣的HMM稱為無約束HMM.在實(shí)際工作中,需要解決的問題一般有領(lǐng)域背景.希望訓(xùn)練出的HMM符合用戶的預(yù)期.例如,在經(jīng)濟(jì)學(xué)中,至少要有兩個(gè)連續(xù)的負(fù)增長(zhǎng)(收縮)狀態(tài),才能說這段時(shí)間處于經(jīng)濟(jì)衰退期[5];在遺傳學(xué)中,一個(gè)罕見的遺傳現(xiàn)象比如,至少要有幾百個(gè)堿基長(zhǎng),才能被認(rèn)為出現(xiàn)了CpG island(Aston and Martin,2007)[6].

針對(duì)于此,本文提出了一種結(jié)合狀態(tài)最短連續(xù)長(zhǎng)度約束的隱馬爾可夫模型,一般的HMM假設(shè)時(shí)間序列是由一系列隱狀態(tài)構(gòu)成,而系統(tǒng)的運(yùn)行本質(zhì)就是在不同的隱狀態(tài)間轉(zhuǎn)換,一般的HMM對(duì)隱狀態(tài)之間轉(zhuǎn)換沒有限制,本文的方法限制了狀態(tài)之間的轉(zhuǎn)換.首先擴(kuò)展?fàn)顟B(tài)數(shù)為原來的倍,并在訓(xùn)練模型時(shí)加約束,限制了狀態(tài)之間轉(zhuǎn)移,這樣學(xué)習(xí)出約束HMM就自帶限制最小變點(diǎn)連續(xù)長(zhǎng)度,滿足在一些特定的應(yīng)用情況下要求狀態(tài)持續(xù)的長(zhǎng)度限制.在解碼隱狀態(tài)序列時(shí),控制最后一個(gè)狀態(tài)一定是擴(kuò)展?fàn)顟B(tài)的最后一個(gè),就能保證只要序列中出現(xiàn)狀態(tài)變化即出現(xiàn)變點(diǎn),并且狀態(tài)已經(jīng)持續(xù)至少h.并將該模型應(yīng)用在兩組已控制變點(diǎn)位置的模擬數(shù)據(jù)和GNP數(shù)據(jù)集,檢測(cè)其變點(diǎn).

1 相關(guān)工作

變點(diǎn)問題的研究始于Page在Biometrika上發(fā)表的一篇關(guān)于連續(xù)抽樣檢驗(yàn)的文章[7].近二十年來,變點(diǎn)問題的研究在理論和世紀(jì)應(yīng)用等方面有了快速的發(fā)展.理論上已有了一系列較為成熟的結(jié)果,國(guó)外的變點(diǎn)的應(yīng)用研究主要涉及金融股票市場(chǎng)檢測(cè)[8]、環(huán)境檢測(cè)[9]、醫(yī)藥與生物工程[10]、系統(tǒng)維護(hù)[11]等.我國(guó)學(xué)者對(duì)變點(diǎn)問題的研究始于上世紀(jì)80年代,由陳希孺教授和繆柏其教授利用“局部法”開始了變點(diǎn)的研究[12].在我國(guó)對(duì)于變點(diǎn)問題的應(yīng)用還是很少的,只在少數(shù)領(lǐng)域中獲得了應(yīng)用,如張學(xué)新等[13]研究了最小二乘法檢測(cè)多個(gè)變點(diǎn)的性能,并成功檢測(cè)出1952-2003年中國(guó)主要經(jīng)濟(jì)部門GNP的變點(diǎn).

對(duì)實(shí)際應(yīng)用背景下,Chib[14]和 Luong[15]提出了一種帶約束的HMM即隱狀態(tài)鏈中的狀態(tài)轉(zhuǎn)移是有一定限制的.Nam[16]提出了基于HMM的限制最小連續(xù)長(zhǎng)度變點(diǎn)的定義,描述了一個(gè)基于該定義的約束HMM,但是他沒有給出明確的約束模型的定義,沒有說明狀態(tài)遷移矩陣以及訓(xùn)練該模型的算法.他是用一般的訓(xùn)練算法來學(xué)習(xí)HMM,對(duì)于限制的最小變點(diǎn)連續(xù)長(zhǎng)度是在分析隱序列時(shí)來約束,而且Nam提出的模型是用來做變點(diǎn)的風(fēng)險(xiǎn)分析.而本文提出的結(jié)合狀態(tài)最短連續(xù)長(zhǎng)度的約束HMM,描述其狀態(tài)之間的遷移限制,給出了約束HMM的訓(xùn)練算法以及狀態(tài)提取算法,用來解決實(shí)際問題中對(duì)于最小狀態(tài)連續(xù)長(zhǎng)度的限制的問題.

2 結(jié)合狀態(tài)最短連續(xù)長(zhǎng)度約束的HMM

2.1 隱馬爾可夫模型

HMM是一種雙重隨機(jī)過程,一個(gè)是隱含的有限狀態(tài)馬爾可夫鏈即xt,它描述狀態(tài)的轉(zhuǎn)移,另一個(gè)描述狀態(tài)與觀察值之間的統(tǒng)計(jì)對(duì)應(yīng)關(guān)系.在實(shí)際問題中我們只能看到觀察值,而不能直接看到隱狀態(tài).只能通過研究觀察值序列yt去推斷隱狀態(tài)序列xt[17].隱馬爾可夫模型為一五元組其中:

狀態(tài)轉(zhuǎn)移概率矩陣:

A=(aij)N′N其中即由狀態(tài)si轉(zhuǎn)移到狀態(tài)sj的概率且

2.2 變點(diǎn)

使用HMM為時(shí)間序列建模,在已知隱狀態(tài)序列的情況下,時(shí)間序列變點(diǎn)的一般定義如定義1.

定義1.在t時(shí)刻隱狀態(tài)鏈的狀態(tài)改變,即t時(shí)刻出現(xiàn)變點(diǎn),即:

然而,對(duì)于一些實(shí)際的應(yīng)用中,確定一個(gè)變點(diǎn),要求這個(gè)狀態(tài)變化要持續(xù)一定的時(shí)間,比如第四部分的分析GNP數(shù)據(jù)的經(jīng)濟(jì)周期等,而用一般的HMM分析時(shí),無法限制.在這些應(yīng)用情況下,Nam定義一個(gè)持續(xù)的變點(diǎn)如定義2.

定義2的兩個(gè)重要的特性:第一,與其他基于HMM 的變點(diǎn)方法相似,變點(diǎn)的分析是從隱狀態(tài)序列推斷的;第二,相比于分析隱馬爾科夫鏈的狀態(tài)變化,更重要的是分析在隱馬爾科夫鏈中最短持續(xù)長(zhǎng)度h的狀態(tài)變化.第二個(gè)點(diǎn)就是本文提出結(jié)合最短連續(xù)長(zhǎng)度約束HMM的主要思想.

2.3 結(jié)合狀態(tài)最短連續(xù)長(zhǎng)度約束的HMM確定在時(shí)間t出現(xiàn)變點(diǎn),即:

針對(duì)Nam定義的持續(xù)的變點(diǎn),本文提出結(jié)合狀態(tài)最短連續(xù)長(zhǎng)度約束的HMM,通過控制了狀態(tài)之間的轉(zhuǎn)移,來控制狀態(tài)的最短連續(xù)長(zhǎng)度,一個(gè)約束HMM為一新的五元組

與一般HMM不同的是,約束HMM在學(xué)習(xí)狀態(tài)轉(zhuǎn)移矩陣時(shí),對(duì)狀態(tài)之間的轉(zhuǎn)移加上了限制.如圖1為狀態(tài)數(shù)為2,h=2時(shí)的狀態(tài)轉(zhuǎn)移圖.初始狀態(tài)只能從擴(kuò)展?fàn)顟B(tài)的第一個(gè)狀態(tài)開始,如圖中的s11,s21,當(dāng)處于同一個(gè)狀態(tài)時(shí),s11只能轉(zhuǎn)移到s12,而處于s12時(shí),說明s1這個(gè)狀態(tài)已經(jīng)持續(xù)了h=2,s12可達(dá)的狀態(tài)可以是自身s12另一個(gè)狀態(tài)的第一個(gè)擴(kuò)展?fàn)顟B(tài)s21.在這樣的限制下,就可以保證一個(gè)狀態(tài)可以持續(xù)至少h=2.

圖 1 2-狀態(tài),h=2狀態(tài)轉(zhuǎn)移圖

當(dāng)h=1時(shí),約束HMM就是一般HMM模型.

2.4 約束HMM的學(xué)習(xí)算法

Baum-Welch算法是隱馬爾科夫模型學(xué)習(xí)問題的一種常用解決方法.約束HMM的學(xué)習(xí)算法是在Baum-Welch算法基礎(chǔ)上加入狀態(tài)最短連續(xù)長(zhǎng)度的約束.訓(xùn)練約束HMM的過程如下:

2)計(jì)算輔助變量:

3)參數(shù)更新,重估公式:

與一般HMM不同,由于隱狀態(tài)擴(kuò)展為原來的h倍,初始狀態(tài)只能從開始,即m11時(shí),.

在每次循環(huán)的參數(shù)重估結(jié)束后,都要修改轉(zhuǎn)移矩陣 A:當(dāng)i=j時(shí),當(dāng)如果其他情況下當(dāng)時(shí),當(dāng)m=1,n=h或者時(shí),其他情況下控制狀態(tài)之間的轉(zhuǎn)移.

本文的方法與一般的Baum-Welch算法不同的地方就在于在參數(shù)重估步驟中,對(duì)初始概率,輸出矩陣和狀態(tài)轉(zhuǎn)移矩陣的修改.

2.5 約束HMM的狀態(tài)提取算法

解決給定模型及觀測(cè)序列O,求生成此觀測(cè)序列的最大可能的隱狀態(tài)序列,一般采用的是Viterbi算法,解碼最大可能的隱狀態(tài)序列.

在分析隱狀態(tài)序列變點(diǎn)時(shí),有兩種情況一種是隱狀態(tài)序列最后一個(gè)狀態(tài)允許是一半狀態(tài),另一種是狀態(tài)必須滿足h個(gè),本文的方法是用到后一種情況來分析Viterbi序列.限制了隱狀態(tài)序列最后一個(gè)狀態(tài)一定是擴(kuò)展?fàn)顟B(tài)的最后一個(gè)狀態(tài),保證了最后一個(gè)狀態(tài)的持續(xù)時(shí)間也不小于h.算法1為約束HMM的狀態(tài)提取算法.

與一般HMM不同的是,約束HMM要限制隱狀態(tài)序列最后一個(gè)狀態(tài)一定是擴(kuò)展?fàn)顟B(tài)的最后一個(gè)狀態(tài),即一定是sih,本文狀態(tài)提取算法與Viterbi算法不同之處是把其他擴(kuò)展?fàn)顟B(tài)的設(shè)置為 -1.

3 實(shí)驗(yàn)

3.1 實(shí)驗(yàn)設(shè)計(jì)與數(shù)據(jù)準(zhǔn)備

為了驗(yàn)證本文提出的約束HMM比一般的HMM在檢測(cè)實(shí)際應(yīng)用上的便捷性,我們將約束HMM應(yīng)用到實(shí)際數(shù)據(jù)的處理與分析工作之中.

實(shí)驗(yàn)中用到的數(shù)據(jù)集:

①已控制變點(diǎn)位置的模擬數(shù)據(jù),變點(diǎn)出現(xiàn)在t=10, t=20處.

第一組序列:

第二組序列:

根據(jù)經(jīng)濟(jì)周期理論,一個(gè)完整的經(jīng)濟(jì)周期應(yīng)包括復(fù)蘇、高漲、衰退、蕭條等四個(gè)階段,故將GNP數(shù)據(jù)即觀測(cè)序列離散化為四個(gè)觀測(cè)值,離散化的結(jié)果如圖2.

圖 2 離散化后的GNP數(shù)據(jù)

基于約束HMM的變點(diǎn)檢測(cè)方法步驟:

① 數(shù)據(jù)集,整理數(shù)據(jù)

②初始化約束HMM

③ 用約束Baum-Welch算法訓(xùn)練約束HMM

④ 用約束Viterbi算法解碼隱狀態(tài)序列

⑤ 根據(jù)實(shí)際應(yīng)用背景分析隱狀態(tài)序列變點(diǎn)

3.2 實(shí)驗(yàn)結(jié)果及分析

3.2.1 分別用約束HMM和HMM分析模擬數(shù)據(jù)

第一組序列:

用HMM為觀察序列建模,得出隱狀態(tài)鏈得到的是:

設(shè)定h=3,約束HMM得到的隱狀態(tài)鏈為:

第二組序列:

用HMM得到隱狀態(tài)鏈得到的是:

用約束HMM得到的隱狀態(tài)鏈為:

對(duì)比序列1和2的結(jié)果,一般HMM容易檢測(cè)出小波動(dòng)的變點(diǎn),即狀態(tài)改變只有一兩個(gè)時(shí)間點(diǎn),容易檢測(cè)出額外的變點(diǎn)如t=1,t=29,而使用約束HMM就能夠準(zhǔn)確地檢測(cè)到變點(diǎn)的位置,并保證了狀態(tài)至少持續(xù)了h=3時(shí)間.

3.2.2 分別用約束HMM和HMM分析GNP數(shù)據(jù)的經(jīng)濟(jì)周期

根據(jù)離散化的GNP數(shù)據(jù),觀測(cè)值有四個(gè)為(“F”,“G”,“S”,“X”)分別代表一個(gè)完整經(jīng)濟(jì)周期的四個(gè)階段:復(fù)蘇、高漲、衰退和蕭條.隱狀態(tài)為(“K”,“S”)分別代表經(jīng)濟(jì)活動(dòng)的擴(kuò)張和收縮狀態(tài).在要求限定下,設(shè)置h=2.約束HMM隱狀態(tài)擴(kuò)展為(“K1”,“K2”,“S1”,“S2”).

用一般的HMM檢測(cè)到的變點(diǎn)分布如圖3所示.

圖 3 一般HMM檢測(cè)到的變點(diǎn)分布

從圖3可觀察到檢測(cè)出了12個(gè)變點(diǎn),分別分布在1951/2~1952/2,1953/3~1954/2,1955/4~1958/2,1959/3~1 959/4,1960/2~1960/4,1962/4,1968/3~1970/4,1971/2~197 1/4,1973/2~1975/1,1977/4~1978/1,1978/3~1980/3,1981/ 2~1982/4,1984/3~1984/4

基于約束HMM檢測(cè)到的變點(diǎn)如圖4.

圖 4 約束HMM檢測(cè)到的變點(diǎn)分布

檢測(cè)出 8個(gè)變點(diǎn)分別分布在 1953/2~1954/2, 1957/2~1958/1,1960/2~1960/1,1969/4~1970/2,1971/2~ 1971/4,1973/3~1975/1,1979/1~1980/2,1981/2~1982/3.

由一般HMM檢測(cè)到的變點(diǎn)和約束HMM檢測(cè)到的變點(diǎn)對(duì)比分析可得到,一般HMM檢測(cè)到的變點(diǎn),無法滿足經(jīng)濟(jì)周期拐點(diǎn)的定義,即無法滿足經(jīng)濟(jì)周期的衰退期至少持續(xù)兩個(gè)季度,比如變點(diǎn)1962/4,該變點(diǎn)只有一個(gè)季度,無法判定為經(jīng)濟(jì)周期,并且變點(diǎn)與變點(diǎn)之間相距時(shí)間太短,如 1977/4~1978/1和1978/3~1980/3之間,約束HMM因?yàn)橄拗屏薶故不會(huì)出現(xiàn)這種情況.還可以觀察到一般HMM會(huì)把數(shù)據(jù)開始和結(jié)束當(dāng)做變點(diǎn),這在約束HMM結(jié)果中不會(huì)出現(xiàn).

4 結(jié)語

基于一般HMM的變點(diǎn)檢測(cè)方法沒有考慮到實(shí)際應(yīng)用背景中要求變點(diǎn)狀態(tài)持續(xù)一定時(shí)間,比如本文實(shí)驗(yàn)中檢測(cè)經(jīng)濟(jì)周期時(shí),需要衰退期至少持續(xù)兩個(gè)季度.本文提出結(jié)合狀態(tài)最短連續(xù)長(zhǎng)度約束HMM通過控制狀態(tài)之間的轉(zhuǎn)移,能夠有效的檢測(cè)出時(shí)間序列的變點(diǎn),并且滿足相關(guān)應(yīng)用背景對(duì)狀態(tài)持續(xù)一定時(shí)間的要求.實(shí)驗(yàn)結(jié)果表明,較之一般HMM,約束HMM在保證檢測(cè)變點(diǎn)準(zhǔn)確度前提下,滿足了狀態(tài)最短持續(xù)時(shí)間,提高了變點(diǎn)檢測(cè)效率.

1 Janacek G.Time series analysis forecasting and control. Journal of Time SeriesAnalysis,2010,31(4):303.

2 Davis RA,Lee TCM,Rodriguez-Yam G A.Structural break estimation for nonstationary time series models.Journal of theAmerican StatisticalAssociation,2006,101(473):223–239.

3 Cheong SA,Fornia RP,Lee GHT,et al.The Japanese economy in crises:A time series segmentation study. Economics:The Open-Access,Open-Assessment E-Journal, 2012,6(2012-5):1–81.

4 Zaccarelli N,Li BL,Petrosillo I,et al.Order and disorder in ecologicaltime-series:Introducing normalized spectral entropy.Ecological Indicators,2013,28(5):22–30.

5 Stock JH,Watson M W.Has the business cycle changed and why?Nber MacroeconomicsAnnual,2002,17(1):159–218.

6 Aston JAD,Martin DEK.Distributions associated with general runs and patterns in hidden Markov models.Annals ofApplied Statistics,2007,1(2):585–611.

7 Page ES.Continuous inspection schemes.Biometrika,1954, 41(1/2):100–115.

8 Tseng YH,Durbin P,Tzeng GH.Using a fuzzy piecewise regression analysis to predict the nonlinear time-series of turbulent flows with automatic change-point detection.Flow Turbulence&Combustion,2001,67(2):81–106.

9 Leonte D,Nott DJ,Dunsmuir WTM.Smoothing and change point detection for gamma ray count data.Mathematical Geology,2003,35(2):175–194.

10 Patra K,Dey DK.A general class of change point and change curve modeling for life time data.Annals of the Institute of Statistical Mathematics,2002,54(3):517–530.

11 Moltchanov D.State description of wireless channels using change-point statistical tests. Wired/wireless Internet Communications,2006,(3970):275–286.

12陳希孺.只有一個(gè)轉(zhuǎn)變點(diǎn)的模型的假設(shè)檢驗(yàn)和區(qū)間估計(jì).中國(guó)科學(xué),1988,31(8):817–827.

13張學(xué)新,段志霞.最小二乘法對(duì)多變點(diǎn)檢驗(yàn)的性能研究.河南師范大學(xué)學(xué)報(bào):自然科學(xué)版,2009,37(6):7–10.

14 Chib S.Estimation and comparison of multiple change-point models.Journal of Econometrics,1998,86(2):221–241.

15 Luong TM,Rozenholc Y,Nuel G.Fast estimation of posterior probabilities in change-point analysis through a constrained hidden Markov model.Computational Statistics &DataAnalysis,2013,(68):129–140.

16 Nam CFH,Aston JAD,Johansen AM.Quantifying the uncertainty in change points.Journal of Time Series Analysis,2012,33(5):807–823.

17 Concha OP,Xu RYD,Moghaddam Z,et al.HMM-MIO:An enhanced hidden Markov model for action recognition. CVPR Workshops,2011:62–69.

Change Point Detection Based on Constrained Hidden Markov Model

ZHUANG Yu,HE Zhen-Feng

(School of Mathematics and Computer Science,Fuzhou University,Fuzhou 350108,China)

The change point detection of time series is widely applied in various fields.In some applications,a minimum period is required before a state change.Motivated by such applications,a constrained Hidden Markov Model,which combines with the shortest state continuous length constraint,is proposed in this study.Moreover,a constrained Baum-Welch training algorithm and a constrained Viterbi state extraction algorithm are also given.And experimental results based on the simulation data and GNP data sets indicate that the constrained HMM has higher performance than the general HMM.

change point detection;constrained Hidden Markov Model;time series segmentation

2016-08-06;收到修改稿時(shí)間:2016-09-20

10.15888/j.cnki.csa.005712

猜你喜歡
檢測(cè)模型
一半模型
“不等式”檢測(cè)題
“一元一次不等式”檢測(cè)題
“一元一次不等式組”檢測(cè)題
“幾何圖形”檢測(cè)題
“角”檢測(cè)題
重要模型『一線三等角』
重尾非線性自回歸模型自加權(quán)M-估計(jì)的漸近分布
3D打印中的模型分割與打包
小波變換在PCB缺陷檢測(cè)中的應(yīng)用
主站蜘蛛池模板: 99久久精品免费观看国产| 国产精品3p视频| 日韩精品久久久久久久电影蜜臀| 日韩在线成年视频人网站观看| 午夜福利免费视频| 国产精品99一区不卡| 91精品国产丝袜| 麻豆精品在线| 区国产精品搜索视频| 三上悠亚精品二区在线观看| 国产精品三区四区| 不卡无码网| 精品久久久久无码| A级全黄试看30分钟小视频| 欧美a在线视频| 成人亚洲天堂| 亚洲第一天堂无码专区| 精品久久蜜桃| 欧洲高清无码在线| 不卡国产视频第一页| 久久久久久国产精品mv| 欧美一级爱操视频| 国产呦精品一区二区三区网站| 日本一区二区三区精品视频| 亚洲无码A视频在线| 国产成人夜色91| 久久男人视频| 这里只有精品免费视频| 亚洲成年人片| 亚洲欧美极品| 欧美成人日韩| 男人天堂亚洲天堂| 中文字幕 日韩 欧美| 国产a v无码专区亚洲av| 亚洲一级毛片免费看| 男人天堂亚洲天堂| 久久青草免费91观看| 欧美一级高清视频在线播放| 中国国产A一级毛片| 亚洲天堂自拍| 日韩视频免费| 日本黄网在线观看| 欧美精品一区二区三区中文字幕| 日韩精品亚洲精品第一页| 热久久国产| 日韩a在线观看免费观看| 一区二区无码在线视频| 高清久久精品亚洲日韩Av| 国产成人AV男人的天堂| 欧美日韩va| 欧美一区精品| 久久影院一区二区h| 国产农村妇女精品一二区| 高h视频在线| 曰韩免费无码AV一区二区| 亚洲欧美日韩动漫| 毛片免费视频| 中文字幕天无码久久精品视频免费| 国产成人艳妇AA视频在线| 国产成人夜色91| 婷婷亚洲视频| 亚洲精品中文字幕午夜| 天天做天天爱天天爽综合区| 毛片一级在线| 2021精品国产自在现线看| 欧美国产精品不卡在线观看| 无码国内精品人妻少妇蜜桃视频 | 亚洲色图欧美| 91区国产福利在线观看午夜| 亚洲一级色| 国内熟女少妇一线天| 国产精品白浆在线播放| 综合色天天| 精品91在线| 91最新精品视频发布页| 在线观看无码av免费不卡网站| 成人国产三级在线播放| 国产精品主播| 亚洲色精品国产一区二区三区| 国产波多野结衣中文在线播放| 一级毛片网| 亚洲天堂网2014|