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

路段通行能力不同的避難點選址模型及算法

2017-10-13 03:25:13劉克艷任佩瑜
中國管理科學 2017年9期

趙 容,劉克艷,任佩瑜

(四川大學商學院,四川 成都 610065)

路段通行能力不同的避難點選址模型及算法

趙 容,劉克艷,任佩瑜

(四川大學商學院,四川 成都 610065)

研究應對突發事件的避難點選址問題。假定一條直線型動態路徑網絡上有n個頂點,由n-1條邊相連,每個頂點有一個權重,每條邊有一個容量。邊的容量表示路段通行能力,是單位時間內允許進入該路段的最大聚集量。目標是在此網絡中選擇k個避難點,并為每個頂點指定一個避難點,使得所有頂點的權重到達各自避難點的最大時間最小。首先根據問題的性質,通過建立動態表結構,結合二分法的思想,在O(nlogn)時間內求解單個避難點選址問題。然后在此基礎上,針對k-避難點選址問題,通過更新動態表,結合動態規劃方法,設計了時間復雜度為O(knlogn)的遞歸算法求解。

道路通行能力;動態網絡;動態規劃;應急管理

1 引言

應急設施選址是應急管理中一項極其重要的內容,合理的應急設施選址能夠有效預防和降低突發事件的危害。Toregas等[1]在1971年首先提出應急設施選址問題,研究如何建立數量最少的應急救援點,在規定時間內給所有需要應急服務的需求點提供救援。此外還有集合覆蓋模型[2]和絕對中心點模型[3]等。這些模型都采用靜態網絡,只考慮距離因素,而在應急疏散過程中,由于道路通行能力的限制,在通往避難點的過程中,受災者可能會因為堵塞耗費大量時間,此時到達避難點的時間不能簡單地用權重與距離的乘積來表示,更應考慮道路通行能力等因素。

Cheng等[4]以2011年3月發生的東日本大地震為研究背景,將動態網絡應用到應急避難點選址問題中,相比靜態網絡,動態網絡中的邊有長度和容量屬性,其中容量表示單位時間允許通過該邊的最大量。動態網絡中的避難點選址問題可以描述為:將城市交通網絡抽象為一個動態網絡,網絡中頂點的權重對應初始時刻該點需要疏散的人數,如某棟建筑里的人數。網絡中邊的容量對應路段通行能力,表示單位時間內允許進入該路段的最大權重。

由于緊急避難時人群密度較大,可以認為人群行走速度一致。在各路段通行能力一致時,受災者到達避難點的時間由兩部分組成:一是受災者通過頂點進入路段時的擁堵時間;二是行走時間。問題的目標是在這樣的動態網絡中建立若干避難點,并為每個頂點指定一個避難點,使得在災害發生時所有頂點的權重能盡快到達各自的避難點。假設各頂點的權重在給定的區間內取值,以疏散完成時間的最大后悔值最小(Min-max regret)為目標,Higashikawa等[5]和Wang Haitao[6]改進了文獻[4]中的數據結構,求解算法復雜度由O(nlog2n)降為O(nlogn)。李紅梅[7]、倪冠群[8]、Arumugam[9]將其擴展到不確定型2-避難點和不確定型k-避難點選址問題,分別給出了自己的求解算法。隨后Bhattacharya[10]改進了前面文獻的算法,針對最小最大后悔值準則下的單個避難點選址問題給出了時間復雜度為O(n)的算法,針對2-避難點選址問題給出了時間復雜度為O(nlog4n)的算法,都是目前最優的結果。同時也有很多學者將動態路徑網絡上的避難點選址問題擴展到不同的網絡結構,如樹圖[11],環狀圖[12],方格網[13]和就近原則下的一般網絡結構[14]。假設各頂點的權重是確定的值,以疏散完成時間最小(Min-max)為目標,針對路徑動態網絡上的k-避難點選址問題,倪冠群等[15]基于“二分思想”提出了時間復雜度為O(nlogkn)的算法,適用于k值比較小的情形,Higashikawa等[16]基于動態規劃的方法設計了時間復雜度為O(knlogn)的算法,隨后進行優化將復雜度降為O(kn)[17]。針對樹圖上的單個避難點選址問題,Mamada等[18]設計了時間復雜度為O(nlog2n)的算法。

以上研究的基本假設是各路段通行能力一致。本文考慮各路段通行能力不同的情形。在路徑網絡上,避難點左右兩邊的權重互不影響,可以分別計算左右兩邊的疏散完成時間,取較大值即為所有權重的疏散完成時間。在向避難點疏散的過程中,會在左右兩邊分別產生一個最大擁堵點,左右兩邊的疏散完成時間將由這個最大擁堵點決定。當各路段通行能力相同時,將避難點向左移動(未越過該點)左邊最大擁堵點不變,這意味著移動避難點后可以在O(1)內求得新的疏散完成時間;而當各路段通行能力不同時,將避難點的向左(右)移動時,兩邊的最大擁堵點的位置也會隨之變化。這使得上述研究[4-17]中采用的平衡二叉樹的結構不再適用。本文根據新問題中擁堵點和疏散完成時間之間的關系,設計由擁堵點序列、與之對應的疏散時間以及該擁堵點和避難點間路段的最小通行能力構成的動態表結構,通過查找動態表求得左右兩邊的疏散完成時間,然后根據左右兩邊疏散完成時間的單調性,采用二分查找算法求解單個避難點選址問題?;趯蝹€避難點選址問題的分析,針對k-避難點選址問題,通過更新動態表,采用動態規劃算法求解。

2 問題描述與基本模型

給定路徑網絡P=(V,E),設路徑上依次排列有頂點v1,v2,…,vn構成頂點集V,頂點vi和vi+1由路段ei連接,所有路段的集合為E={e1,e2,…,en-1}。對于vi∈V,定義權重ωi為該頂點的權重,對于ei∈E,定義容量ci為該路段的通行能力,表示單位時間內允許流入路段ei的最大權重。定義τ為權重在各路段上流動單位距離所需的時間。設區間[vi,vj]中有一個避難點x,各權重在向避難點疏散的過程中不允許出現交叉流。定義區間[vi,x]中的疏散點屬于避難點x的左邊,區間[x,vj]中的疏散點屬于避難點x的右邊。特別地,當x=vl時,vl處的權重視為立即到達避難點x。

根據以上定義,初始時刻,權重集中在各頂點處。疏散過程開始時,所有頂點處的權重同時向各自指定的避難點流動。這時首先有了頂點處的權重進入路段的時間,頂點v處的權重ω進入容量為c的路段的時間為ω/c。它們全部進入路段后即變為流量為c,總量為ω的權重,這里的流量表示單位時間內通過路段橫截面的權重。由于各路段通行能力不同,權重到達避難點的時間除了從頂點處進入路段的時間和行走時間外,還有第三部分時間,是權重從通行能力高的路段進入通行能力低的路段時的擁堵時間。根據以上討論,流量為c的權重ω通過容量為ce的路段e時,可能出現以下2種情形:

(1)c≤ce。如下圖1中a圖所示,流量為c的權重ω通過路段e時不會遇到堵塞,即沒有等待時間,所有權重始終以速度1/τ流動,它們通過路段e的時間為leτ。

圖1 流量為c的權重ω通過路段e的兩種情形

由以上討論可知,流量為c的權重ω通過長度為le、容量為ce的路段e的時間te可以表示為:

(1)

圖2 權重ω1到達避難點x的過程中各路段流量

根據以上定義,對于區間[vi,vj],用Li(x)表示x左邊的疏散完成的時間,Rj(x)表示x右邊的疏散完成的時間,Θi,j(x)表示區間[vi,vj]中所有的點的疏散完成時間。

(2)

(3)

Θi,j(x)=max{Li(x),Rj(x)}

(4)

用x=(x1,…,xk)表示k個避難點的位置坐標,用d=(d1,…,dk-1)表示對應劃分點的角標,如下圖3所示,劃分點vdi-1和vdi指定了區間[vdi-1+1,vdi]中的權重全部疏散到避難點xi。問題的目標是在路徑P上選擇k個避難點,同時確定k-1個劃分點,將頂點集V劃分為k個子集,使得k個子集中的所有疏散點到達各自避難點的最大時間最小。

圖3 避難點及其劃分點

用Θi(x,d)表示區間[vdi-1+1,vdi]中的權重全部到達避難點xi的時間,路徑P上所有權重的疏散完成的時間則可以表示為

Θ(x,d)=max{Θi(x,d)|i=1,…,k}。路徑P上的k-避難點選址問題可以表示為:

Pk:minmax{Θ(x,d)|x∈Pk,d∈{1,...,n}k-1}

3 單個避難點選址問題

本節主要介紹如何在區間[vi,vj]內求解單個避難點的選址問題。由(4)式可知,此問題可以表示為:

P1:min{Θi,j(x)|x∈[vi,vj]}。

根據第1節的描述,首先有如下引理1。

圖4 路段通行能力相同與不同時圖形比較

引理2Li(x)隨x的增加而增加,Rj(x)隨x的增加而減小。

引理3 問題P1存在唯一解。

證明:由(4)式和引理2可知Θi,j(x)在區間[vi,vj]內是關于x的凸函數,所以存在唯一xopt=argmin{Θi,j(x)|x∈[vi,vj]}。

引理4 設x∈[vi,vj],若Li(x)≤Rj(x),則xopt≥x;若Li(x)≥Rj(x),則xopt≤x。

證明:根據(4)式和引理2可知,Θi,j(x)是關于x的凸函數,當Li(x)≥Rj(x)時,若將x向右移動,則Li(x)將繼續增加,由(4)式可知,Θi,j(x)=Li(x)也將隨之增加,所以xopt≤x。同理可證,若Li(x)≤Rj(x),則xopt≥x。

引理4表明通過反復比較Li(vl)和Rj(vl)可將xopt的取值范圍限定在某個區間[vk,vk+1]內。而在區間(vk,vk+1)內,避難點左右兩邊的最大擁堵點的位置將保持不變,這時的情形同文獻[4],可以在O(1)時間內直接求得xopt。

算法A:求TL(i,j),見表1。

表1 TL(i,j)的算法流程

表2 區間[vi,vj]的左動態表TL(i,j)

引理5 由動態表TL(i,j)可以在O(logn)時間內求得任意L(i,k)(i≤k≤j)。

綜上,即證由動態表TL(i,j)可在O(logn)時間內求得任意L(i,k)(i≤k≤j)。

定理1:問題P1可在O(nlogn)時間內求解。

運行一次算法A可得TL(i,j),同理可在O(nlogn)時間內求得TR(i,j)。由引理5可知,根據動態表TL(i,j)和TR(i,j)可在O(logn)時間內求得任意L(i,k)和R(k,j)。引理4表明采用二分法求解問題P1時,需要計算L(i,k)和R(k,j)的次數為O(logn)。所以在已求得動態表TL(i,j)和TR(i,j)后,可在O(log2n)時間內求解問題P1。

綜上,即證問題P1可在O(nlogn)時間內求解。

4 k-避難點選址問題

問題Pk滿足最優性原理,本節采用動態規劃算法進行求解。設區間[vi,vj]內建立的p個最優避難點的位置坐標為p維向量x*(p,i,j),與之對應的最優劃分點的角標為p-1維向量d*(p,i,j),它們的最小疏散完成時間為opt(p,i,j)。根據以上定義,求解問題Pk即求解opt(k,1,n)。根據第2節的討論有如下遞歸表達式:

(5)

用d(p,i,j)表示表示d*(p,i,j)的第p-1個值,即第p-1個和第p個避難點間的劃分點,同樣x(p,i,j)表示x*(p,i,j)的第p個值,x(1,d+1,j)表示區間[vd+1,vj]中的單個最優選址位置,則x*(p+1,i,j)和d*(p+1,i,j)可表示為:

x*(p+1,i,j)=(x*(p,i,d),x(1,d+1,j))

d*(p+1,i,j)=(d*(p,i,j),d(p+1,i,j))

(6)

由引理2可以推出d(p,i,j)和x(p,i,j)隨j的減小而減小,證明同Higashikawa等[17]。

性質1 對2≤p≤k,d(p,i,j-1)≤d(p,i,j)。

性質2 對1≤p≤k,x(p,i,j-1)≤x(p,i,j)。

由第2節的討論可知,求解單個避難點選址問題時,本文所采用的數據結構與Higashikawa等[17]不同。同樣這里在求解k-避難點選址問題時,采用的遞歸順序也與Higashikawa等[17]不同。根據建立的避難點的個數不同將它們分為k層,即1,2,…,k。在每一層采用與Higashikawa等[17]中所述的相反的順序遞歸,即按照opt(1,1,n),…,opt(1,1,1);opt(2,1,n),…,opt(2,1,1);…;opt(k,1,n)的順序求解。根據遞歸表達式(5),每一層都可以由上一層的所有值和opt(1,1,n),…,opt(1,n-1,n)求得。采用這樣的逆序遞歸時只需計算一次opt(1,2,n),…,opt(1,n-1,n),然后在每一層遞歸中都可以用到它們,如此減少了更新的次數,也減少了更新的方向。由于采用逆序,在每一層內更新時,只需由L(α,β)向2個方向更新,即L(α-1,β)和L(α,β-1)。具體如下圖5所示。

圖5 由opt(p,1,n)計算opt(p,1,n-1)示意圖

如圖5所示,用dn表示d(p,1,n),xn表示x(p,1,n),dn-1、xn-1同理。由第2節的討論可知,只要確定了xopt所屬的區間即可在O(1)時間內求解,而確定xopt所屬區間的方法是運用引理4,通過反復比較L(i,l)和R(l,j)得到。設xn屬于區間[vx[n],vx[n]+1),xn-1屬于區間[vx[n-1],vx[n-1]+1)。根據(5)式,由opt(p,1,n)計算opt(p,1,n-1)可以通過逐步比較opt(p-1,1,dn)和opt(1,dn+1,n-1),opt(p-1,1,dn+1)和opt(1,dn+2,n-1)等求得。根據遞歸順序可知,在計算第p層時,所有p-1層的opt(p-1,1,j)已知,只需進一步由opt(1,dn+1,n)計算opt(1,dn+1,n-1)…和opt(1,dn+2,n)…,根據第2節的分析,左邊涉及的是由L(dn+1,[xn])更新到L(dn-1+1,[xn-1])。由性質1和2可知,dn-1≤dn,xn-1≤xn。所以在整個計算過程中L(i,j)只有2個更新方向,即L(i-1,j)和L(i,j-1)。引理5已表明可由動態表TL(i,j)在O(logn)時間內求得任意L(i,k)。下面主要說明如何由動態表TL(i,j)在O(logn)時間內更新為TL(i-1,j)。

算法B1:由TL(i,j)求TL(i-1,j),見表3。

表3 TL(i-1,j)的算法流程

根據以上討論,設計如下算法B自下而上地求得opt(k,1,n)、d*(k,1,n)和x*(k,1,n),即是問題Pk的解。

算法B:求解問題Pk,見表4。

表4 Pk的算法流程

定理2采用遞歸算法B求解問題Pk的時間復雜度為O(knlogn)。

證明:第(1)步中,采用算法A建立動態表的時間復雜度為O(nlogn)。由定理1的證明可知根據動態表求解opt(1,1,n)的時間復雜度為O(log2n)。由引理2和4可知,x(1,1,n)≤x(1,2,n),所以按opt(1,2,n),…,opt(1,n-1,n)的順序求解時至多調用動態表O(n)次,由引理5可知,通過動態表計算L(i,j)和R(i,j)的時間復雜度為O(logn)。所以求解opt(1,1,n),…,opt(1,n-1,n)的時間復雜度為O(nlogn)。同理依次求解opt(1,1,n-1),…,opt(1,1,2)的時間復雜度也是O(nlogn)。即證第(1)步的計算時間復雜度為O(nlogn)。

第(2)步中,根據遞歸表達式(5),由前一次循環得到的opt(p-1,1,n),…,opt(p-1,1,p)和第1步得到的opt(1,1,n),…,opt(1,n-1,n)求解opt(p,1,n)的時間復雜度為O(logn)。在依次求解opt(p,1,n),opt(p,1,n-1),…,opt(p,1,p)時,由性質1和2可知:d(p,1,j-1)≤d(p,1,j),x(p,1,j-1)≤x(p,1,j)。所以d(p,1,j)至多變化n-p次,至多減少n-1,對于x(p,1,j)所屬的區間同理。根據建立的動態表TL(dn+1,[xn]),對于x(p,1,j)的變化不需更新表,可以直接在表中查找。對于d(p,1,j)的變化,d(p,1,j)每減少1則需要運行算法B1更新一次動態表。所以每一次循環中動態表TL(dn+1,[xn])需要更新O(n)次,即每一次循環中更新動態表的時間復雜度為O(nlogn)。由引理4可知,d(p,1,j)每減少1需要對L(dj+1,x[j])和R(x[j],j)進行O(1)次比較,所以每一次循環中需要通過動態表計算L(i,j)和R(i,j)的次數為O(n),即每一次循環中計算L(i,j)和R(i,j)的時間復雜度為O(nlogn)。第(2)步中有k-1次循環,即證算法B的第(2)步的時間復雜度為O(knlogn)。

綜上,采用遞歸算法B求解問題Pk的時間復雜度為O(knlogn),即證。

本文從理論上擴展了動態網絡中的k-避難點選址問題,在已有的權重確定、通行能力為1的直線型動態網絡(動態路徑網絡)中的k-避難點選址問題的基礎上,研究了權重確定、各路段通行能力為任意常數的動態路徑網絡中的k-避難點選址問題。目前針對權重確定、通行能力為1的直線型動態網絡中的k-避難點選址問題的研究中設計的求解算法復雜度分別為O(nlogkn)[15]、O(knlogn)[16]和O(kn)[17],本文主要在Higashikawa等[17]的基礎上進行了拓展,根據新問題的性質,放棄了原有的平衡二叉樹結構,采用一種動態鏈表結構儲存,求解思路仍是儲存擁堵點信息,只是本文因為各路段容量不同,而路段容量對最大擁堵點的計算有影響,所以儲存了決定對應疏散完成時間的容量值。如此,在采用動態規劃算法求解k-避難點問題時,雖然每次迭代更新時不能同文獻[17]一樣通過一次比較獲得,但根據動態表建立的序列可以在O(logn)時間內插入新的值。最后本文設計的求解k-避難點選址問題的遞歸算法的時間復雜度為O(knlogn)。

5 結語

現有研究避難點選址問題的文獻多沒有考慮道路通行能力及由此帶來的擁堵問題,使得避難點的選址位置在實際中并不適用。而在考慮道路通行能力的文獻中都假設各段道路通行能力相同,與實際的道路情況不符,且缺乏變通。本文在后者的基礎上,放松了各路段通行能力一致的假設,通過建立關于最大擁堵點的動態表結構,結合動態規劃的方法,解決了在確定型路徑動態網絡中建立k個避難點使得路徑上所有權重到達各自避難點的最大時間最小的問題,為后續研究提供參考。

[1] Toregas C,Swain R,Revelle C,et al.The location of emergency service facilities [J].Operations Research,1971,19(6):1363-1373.

[2] Adel A A,A.White J A.Probabilistic formation of the emergency service location problem [J].Journal of Operational Research Society,1978,29(12):1167-1179.

[3] Shier D,Dearing P.Optimal locations for a class of nonlinear single-facility location problems on a network [J].Operations Research,1983,31(2):292-303.

[4] Cheng S,Higashikawa Y,Katoh N,et al.Minimax regret 1-sink location problems in dynamic path networks [C]//Proceedings of the 10th International Comference on Theory and Applincations of models of Computation Hongkong,China,May 20-22,2013.

[5] Higashikawa Y,Augustine J,Cheng S,et al.Minimax regret 1-sink location problem in dynamic path networks[J].Theoretical Computer Science,2015,24-36.

[6] Wang Haitao.Minmax regret 1-facity location on uncertain path networks [J].European Journal of Operational Research,2014,239(3): 636-643.

[7] Li Hongmei,Xu Yinfeng,Ni Guanqun.Minimax regret vertex 2-sink location problem in dynamic path networks [J].Journal of Combinatorial Optimization,2014.

[8] Ni Guanqun,Xu Yinfeng,Dong Yucheng.Minimax regret k-sink location problem in dynamic path networks [C]//Proceedinys of the 10th International Conference on Algorithmic Aspects in Information and Management,Vancouver,Canada,July 8-11,2014.

[9] Arumugam G,Augustine J,Golin,et al.A polynomial time algorithm for minimax-regret evacuation on a dynamic path [J].Computerscience,2014,588(c):1404-5448.

[10] Bhattacharya B,Kameda T.Improved algorithms for computing minmax regret sinks on dynamic path and tree networks [J].Theoretical Computer Science,2014,607:411-425.

[11] Higashikawa Y,Golin M,Katoh N.Minimax regret sink location problem in dynamic tree networks with uniform capacity [M].2014,8344:125-137.

[12] Xu Xinfeng,Li Hongmei.Minimax regret 1-sink location problem in dynamic cycle networks [J].Information Processing Letters,2015,115(2):163-169.

[13] Naoyuki K,Katoh N.The universally quickest transshipment problem in a certain class of dynamic networks with uniform path-lengths [J].Discrete Applied Mathematics,2014,178:89-100.

[14] Li Hongmei,Xu Yinfeng.Minimax regret 1-sink location problem with accessibility in dynamic general networks [J].European Journal of Operational Research,2015,250(2):360-366.

[15] 倪冠群,徐寅峰,徐久平.考慮道路通行能力的應急避難點選址模型及算法 [J].中國管理科學,2015,23(1):82-88.

[16] Higashikawa Y,Golin M,Katoh N.Multiple sink location problems in dynamic path networks [J].Theoretical Computer Science,2015,607:2-15.

[17] Higashikawa Y,Golin M,Katoh N.Improved algorithms for multiple sink location problems in dynamic path networks [J].Computer Science,2014:1405-2014.

[18] Mamada S,Uno T,Makino K,et al.An algorithm for the optimal sink location problem in dynamic tree networks [J].Discrete Applied Mathematics,2006,154(16):2387-2401.

Abstract: From the viewpoint of disaster prevention from city planning and evacuation planning,it is important to establish effective evacuation planning systems against large scale disasters.Considering the different roads have different capacity,the k-sink location problem in dynamic network with different capacity is proposed.

In our model,each vertex supplies with a certain nonnegative value and each edge has a capacity representing the least upper bound for the units flowing into the edge per unit time.It is found that the time for a vertex weightωto go through the edgeewhich have a capacityceand a lengthleisleτ+ω/ce-ω/c,whereτis a constant representing the time required for traversing the unit distance of per unit weight and c is the flow ofω.Our goal is to findksinks andk-1 divides which minimize the maximum time for all units flowing into the corresponding sink that the divides have provided.First,the 1-sink location problem is analyzed and it is found the monotonicity and unimodality of the evacuation completion time.Then based on some new properties,the linked list data structure is used to store the completion time and the minimum road capacity on their way to the sink of the maximum congestion points,which make the solution process easier.On this basis,anO(nlogn) time algorithm is developed to solve the 1-sink location problem.Finally,anO(knlogn) time recursive algorithm is developed to solve thek-sink location problem based on dynamic programming,wherenis the number of vertices in the given network.

Since we are the first to analyze the sink location problem in dynamic network with different capacity,our research may be useful to the further research such as the sink location problem in dynamic tree network with different capacity and the sink location problem in dynamic network with interval weight and different capacity.

Keywords: road capacity;dynamic network;dynamic programming;emergency management

Min-max Multiple Sink Location Problem in Dynamic Path Networks with Different Traffic Capacity Constraint

ZHAORong,LIUKe-yan,RENPei-yu

(Business School,Sichuan University,Chengdu 610065,China)

C931;O221

A

1003-207(2017)09-0133-08

10.16381/j.cnki.issn1003-207x.2017.09.015

2016-02-25;

2016-06-08

國家自然科學基金資助項目(71371130,71501019);四川旅游發展研究中心項目(LYC16-16);賽爾網絡下一代互聯網技術創新項目

任佩瑜(1952-),男(漢族),重慶人,四川大學商學院教授,博士生導師,研究方向:管理科學與工程,E-mail:renpeiyu@scu.edu.cn.

主站蜘蛛池模板: 992Tv视频国产精品| 毛片基地美国正在播放亚洲 | 国产好痛疼轻点好爽的视频| 亚洲熟妇AV日韩熟妇在线| 毛片在线看网站| 日本黄网在线观看| 欧美精品v欧洲精品| 国产欧美综合在线观看第七页| 中文字幕无码制服中字| 无码精品福利一区二区三区| 国产主播一区二区三区| 99久久精品免费看国产电影| 国产精品99在线观看| 国产第八页| 亚洲伊人天堂| 在线精品亚洲国产| 欧美高清国产| 精品无码一区二区三区电影| 国产靠逼视频| 最新亚洲人成无码网站欣赏网 | 亚洲三级影院| 激情综合网址| 囯产av无码片毛片一级| 国产在线专区| 欧美一级一级做性视频| 九九九国产| 中国国语毛片免费观看视频| 午夜国产精品视频黄| 日韩资源站| 天天色天天综合| 日韩高清中文字幕| 91在线免费公开视频| 99视频在线观看免费| 国产精品一区不卡| 波多野结衣一区二区三区四区视频| 免费av一区二区三区在线| 午夜视频在线观看免费网站| 992tv国产人成在线观看| 亚洲三级视频在线观看| 亚洲成A人V欧美综合| 2019年国产精品自拍不卡| 九九热视频精品在线| 久久久久亚洲精品成人网 | 永久在线精品免费视频观看| 婷婷久久综合九色综合88| 亚洲欧美综合精品久久成人网| 91精品网站| 国产系列在线| 激情综合图区| 色老二精品视频在线观看| 日韩精品无码一级毛片免费| 91久草视频| 亚洲精品欧美日韩在线| 国产精品九九视频| 亚洲精品国产乱码不卡| 欧美色99| 国内精自线i品一区202| 澳门av无码| 国产成人无码AV在线播放动漫| 在线国产欧美| 国产三级毛片| 五月丁香伊人啪啪手机免费观看| 亚洲美女一级毛片| 亚洲二区视频| 欧美成人精品在线| 精品福利国产| 国产丝袜无码精品| 波多野一区| 亚洲国产亚洲综合在线尤物| 亚洲天堂首页| 国产成人一二三| 久久免费观看视频| 第一页亚洲| 四虎AV麻豆| 欧美成人在线免费| 免费又黄又爽又猛大片午夜| 日韩在线影院| 99久久精品国产自免费| 国产特级毛片aaaaaaa高清| 婷婷综合亚洲| 亚洲精品无码不卡在线播放| 91欧美亚洲国产五月天|