亚洲免费av电影一区二区三区,日韩爱爱视频,51精品视频一区二区三区,91视频爱爱,日韩欧美在线播放视频,中文字幕少妇AV,亚洲电影中文字幕,久久久久亚洲av成人网址,久久综合视频网站,国产在线不卡免费播放

        ?

        中國和巴西港口擁塞指數(shù)多重分形研究

        2016-05-25 00:37:02林國龍馬理博張辰彥何紅弟
        關鍵詞:長程岬型分形

        林國龍,馬理博,張辰彥,何紅弟

        (上海海事大學 物流研究中心,上海 201306)

        中國和巴西港口擁塞指數(shù)多重分形研究

        林國龍,馬理博,張辰彥,何紅弟

        (上海海事大學 物流研究中心,上海 201306)

        運用多重分形消除趨勢波動分析法(MF-DFA)與多重分形消除趨勢波動交叉相關分析法(MF-DCCA)對全球港口擁塞指數(shù)的中國港口ETA(Estimated Time of Arrival)海岬型船與巴西港口ETA海岬型船港口擁塞指數(shù)進行研究。研究表明:中國港口ETA海岬型船與巴西港口ETA海岬型船時間序列廣義Hurst指數(shù)及MF-DCCA廣義Hurst指數(shù),均單調遞減且值域在(0.5,1)之間,標度指數(shù)均表現(xiàn)凹向橫軸的函數(shù),分形譜密集鐘形分布;兩組時間序列均具有明顯的多重分形特征,同時具有長期記憶性與多標度性,且二者多重分形交叉相關,表現(xiàn)長期相互影響。

        交通運輸工程;港口擁塞指數(shù);多重分形;消除趨勢波動分析;交叉相關分析

        0 引 言

        由全球港口〔G-ports (UK) Ltd〕發(fā)布的全球港口擁塞指數(shù)報告(The Global Ports Congestion Index,GPCI)是衡量全球港口擁堵狀況的有效參考。GPCI以周報的形式發(fā)布,詳細及時地統(tǒng)計反映全球滯留主要港口的煤炭、礦石等干散貨船數(shù)量,分析由于船舶壓港、港口擁塞對干散貨市場供求關系以及散貨船租金所帶來的影響,及時追蹤全球主要港口煤炭、礦石等干散貨船滯期情況,覆蓋了全球不同國家80個煤炭、礦石碼頭的擁塞情況,包括澳大利亞、巴西、中國、印度以及南非等國。GPCI能夠反映全球港口擁堵現(xiàn)象及程度。報告中港口擁塞指數(shù)是非平穩(wěn)的時間序列。通過研究港口擁塞指數(shù)時間序列的變化情況,探索港口擁塞指數(shù)間的內在相關性,發(fā)現(xiàn)港口擁塞指數(shù)的發(fā)展特征,從而能夠預測港口擁堵現(xiàn)象的發(fā)展方向。

        港口擁堵問題備受全球關注,源于2002年美西碼頭工人罷工導致洛杉磯、長灘碼頭擁堵。近年來進出口貿易往來頻繁對國際海運航運行業(yè)刺激加劇,港口擁堵問題也隨之加劇且成為全球性問題。尤其是集裝箱貿易運輸,運價運力皆不斷攀升。張榮忠[1]在全球性港口擁堵現(xiàn)象透析文中就港口擁堵現(xiàn)狀的原因及解決辦法給出過介紹。尹銘鈺[2-3]依據(jù)全球港口擁塞指數(shù)報告,詳細分析全球干散貨港口擁堵現(xiàn)狀,定性指出滯港船舶情況呈現(xiàn)區(qū)域性不均衡、巴西、澳大利亞港口擁堵嚴重、中國港口擁堵集中于鐵礦石運輸港的港口擁堵問題;隨后針對巴西港口擁堵情況分析得出,2014年港口擁堵情況略有好轉,但整體而言,巴西港口的擁堵現(xiàn)狀明顯要嚴峻于世界平均水平。徐劍華[4]指出美西港口擁堵的深層原因包括船舶大型化、班輪公司超級聯(lián)盟和承運人為節(jié)約成本出售自有底盤車隊。 P.T.LEACH[5-6]就美西港口擁堵問題仍對航運市場造成的影響,提出重新定義航運業(yè)務旺季。

        分形(fractal)是20世紀70年代Mandelbort表征復雜圖形和過程時首先引入于科學領域,表示不規(guī)則的物體。C.K.PENG等[7]在研究DNA序列時首次提出消除趨勢波動分析方法(detrended fluctuation analysis,DFA),該方法能有效計算一個非平穩(wěn)時間序列地長程相關性,即非平穩(wěn)時間序列的內在自相似性。J.W.KANTELHARDT等[8]提出檢驗一個非平穩(wěn)時間序列是否具有多重分形特征的方法,多重分形消除趨勢波動分析法(multifractal detrended fluctuation analysis,MF-DFA)是在DFA基礎上更加準確有效的判斷一個時間序列長程相關性的定量方法。B.PODOBNIK等[9]首次提出定量分析兩列非平穩(wěn)時間序列相關性的科學分析方法消除趨勢波動交叉相關分析法(Detrended Cross-Correlation Analysis,DCCA)。L.Y.HE等[10]在檢驗中美兩國農產品期貨市場的量價關系時,首次采用多分形消除趨勢波動交叉相關分析法(Multifractal Detrended Cross-Correlation Analysis,MF-DCCA),并證實其多重分形特征及交叉相關性。近年來,DFA,MF-DFA,DCCA和MF-DCCA被應用于多個領域時間序列研究,包括生命自然科學[11]、證券黃金期貨市場[12-14]、原油集裝箱市場[15-16]、水文地理地質[17-18]、設備故障診斷[19]、能源[20]、交通[21]等多方領域。如此廣泛地應用足以說明消除趨勢波動分析法與消除趨勢波動交叉相關分析法在研究復雜非平穩(wěn)時間序列內在自相關和互相關的價值與作用。

        采用MF-DFA與MF-DCCA對全球港口擁塞指數(shù)的中國港口ETA(Estimated Time of Arrival)海岬型船(100000mts+)與巴西港口ETA海岬型船港口擁塞指數(shù)進行研究,MF-DFA結果發(fā)現(xiàn)中國港口ETA海岬型船與巴西港口ETA海岬型船港口擁塞指數(shù)時間序列具有明顯的多重分形特征,廣義Hurst指數(shù)均大于0.5小于1,尺度函數(shù)均為凹向橫軸的函數(shù),多重分形譜呈密集鐘形分布。MF-DCCA結果發(fā)現(xiàn)中國港口ETA海岬型船與巴西港口ETA海岬型船港口擁塞指數(shù)時間序列具有明顯的多重分形交叉相關特征, 廣義Hurst指數(shù)大于0.5小于1,尺度函數(shù)為凹向橫軸的函數(shù),多重分形譜呈密集鐘形分布。

        1 MF-DFA與MF-DCCA基本原理

        多重分形消除趨勢波動交叉相關分析法(MF-DCCA)是計算兩個非平穩(wěn)時間序列內在交叉相關性的一種定量方法,理論基礎是多重分形消除趨勢波動分析法(MF-DFA),MF-DFA是在DFA[22]的基礎上檢驗一個非平穩(wěn)時間序列多重分型特性的方法。對長度為N的時間序列{xi},i=1,2,…,N和等長度時間序列{yi},i=1,2,…,N,MF-DCCA的具體步驟:

        1)計算兩個時間序列的累積離差:

        (1)

        (2)

        2)將累積離差X(t),Y(t)分別等分成Ns個互不相關時間長度為s的等時距區(qū)間。由于時間序列總長N不一定是長度s的整數(shù)倍,為保證時間序列信息完整,將時間序列反向再計算一次累積離差。因此求得兩組2Ns個長度相同子區(qū)間,記每個子區(qū)間為v=(1,2,…,2Ns)。

        3)對兩組時間序列每個子區(qū)間的時間數(shù)據(jù)采用最小二乘法分別進行多項式回歸擬合,得擬合多項式xv(t)和yv(t),求局部協(xié)方差函數(shù)。

        當v=1,2,…,2Ns時:

        {Y[(v-1)s+i]-yv(i)}

        (3)

        當v=Ns+1,…,2Ns時:

        (4)

        4)把所有子區(qū)間的局部協(xié)方差取均值,計算所有區(qū)間的平均值,得到時間序列的q階波動函數(shù)

        (5)

        當q=0時,時間序列的波動函數(shù)為

        (6)

        當q=2時,對一列時間序列的MF-DFA退化成標準DFA。

        5)若時間序列長程冪律交叉相關,F(xiàn)q(s)與s應滿足冪律關系式:

        Fq(s)~sHxy(q)

        (7)

        將式(7)兩邊同時取對數(shù)得

        ln[Fq(s)]~lns

        (8)

        將時間序列在對數(shù)坐標中的散點圖用最小二乘法進行擬合,擬合函數(shù)斜率即為q階廣義Hurst指數(shù)Hxy(q)。

        6)經(jīng)過MF-DCCA的計算得到對應不同q下的廣義Hurst指數(shù)Hxy(q),我們也稱其為標度指數(shù),這個指數(shù)是衡量時間序列是否長程冪律交叉相關的參數(shù),標度指數(shù)的變化特征描述為:①Hxy(q)是常數(shù)時,時間序列具有單分形特征;②是q的單減函數(shù),時間序列具有多重分型特征;③Hxy(q)值域為(0,1)時,可判斷時間序列是否長程冪律交叉相關;④Hxy(q)值域為(0,0.5)時,時間序列具有反長程冪律交叉相關性;⑤Hxy(q)為常數(shù)0.5時,時間序列隨機游走;⑥Hxy(q)值域為(0.5,1)時,時間序列具有正長程冪律交叉相關性。

        7)廣義Hurst指數(shù)Hxy(q)與質量指數(shù)τ(q)的關系式為

        τ(q)=qHxy(q)-1

        (9)

        若實驗時間序列表現(xiàn)多重分形特征,則質量指數(shù)為凹向橫軸的函數(shù),質量指數(shù)可進一步檢驗時間序列是否具有長程冪律交叉相關的特性。

        8)經(jīng)過Legendere變換,可得到多標度時間序列的多重分形譜f(α):

        (10)

        f(α)=q[α-Hxy(q)]+1

        (11)

        式中:α為時間序列區(qū)間上的奇異指數(shù),α值越小表示時間區(qū)間奇異性越大;f(α)為多重分形譜,多重分形譜能說明實驗時間序列發(fā)展趨勢,分形譜越寬則實驗時間序列的波動越大,分形譜曲線極大值為1,若左端高于右端說明實驗時間序列出現(xiàn)高位置概率更大,反之,若右端高于左端則說明實驗時間序列出現(xiàn)低位置概率更大,而且圖像左右高度距離差值越大,時間序列出現(xiàn)高低位置概率相差也就越大。

        2 港口擁塞指數(shù)多重分形研究

        筆者用全球港口擁塞指數(shù)報告中時間序列數(shù)據(jù)進行MF-DFA,MF-DCCA研究,發(fā)現(xiàn)中國港口ETA海岬型船與巴西港口ETA海岬型船港口擁塞指數(shù)時間序列具有較為明顯的多重分形特征,且二者明顯多重分形交叉相關。

        2.1 數(shù)據(jù)來源

        選擇地理位置上無明顯相關關系的中國和巴西港口,選用由全球港口公司發(fā)布的全球港口擁塞指數(shù)報告的中國港口ETA海岬型船和巴西港口ETA海岬型船港口擁塞指數(shù),2008年第2周至2015年第17周共380個周數(shù)據(jù)進行研究,其原始數(shù)據(jù)分布如圖1。

        圖1 中國和巴西港口ETA海岬型船港口擁塞指數(shù)波動Fig.1 Congestion index fluctuations of Chinese and Brazil ETA Capesize port

        2.2 研究分析

        運用MF-DCCA對港口擁塞指數(shù)時間序列進行實驗,將實驗數(shù)據(jù)導入已編碼的MF-DCCA MATLAB R2013a程序并運行,提取運行數(shù)據(jù)用OriginPro8.0作圖,見圖2。

        圖2 Fq(s)與s的對數(shù)關系圖Fig.2 The logarithmic relation diagram of Fq(s) and s

        由圖2可見,對應式(1)~式(6),表示q變化時波動函數(shù)Fq(s)與s呈雙對數(shù)關系,圖中從上到下3條線表示q=5,0,-5時時間序列的擬合狀況,隨著q的增大擬合函數(shù)圖像單調遞增且斜率減小,因為擬合函數(shù)的斜率即是為廣義Hurst指數(shù)h(q),因此實驗結果能夠初步說明廣義Hurst指數(shù)h(q)隨q的增大而減小。

        圖3為對應式(7)和式(8)給出MF-DFA與MF-DCCA廣義Hurst指數(shù)。

        圖3 港口擁塞指數(shù)時間序列廣義Hurst指數(shù)Fig.3 Generalized Hurst index of port congestion index in time series

        由圖3可見,q增大h(q)單調遞減,h(q)值域在(0.5,1)內,充分說明巴西港口ETA海岬型船與中國港口ETA海岬型船港口擁塞指數(shù)時間序列均具有多重分形特征,均表現(xiàn)正的長程冪律相關。巴西港口ETA海岬型船港口擁塞指數(shù)時間序列廣義Hurst指數(shù)全部位于中國港口ETA海岬型船港口擁塞指數(shù)時間序列廣義Hurst指數(shù)下方,且值域更大,說明巴西港口ETA海岬型船港口擁塞指數(shù)序多重分形特征更為明顯,其函數(shù)圖像的斜率下降更快。MF-DCCA分析結果表明,廣義Hurst指數(shù)Hxy(q)位于MF-DFA Hurst指數(shù)h(q)結果之間且單調遞減。值域在(0.5,1)內,說明巴西港口ETA海岬型船與中國港口ETA海岬型船港口擁塞指數(shù)時間序列之間具有多重分形交叉相關特性,且呈正的冪律相關關系。

        為了進一步驗證時間序列多重分形特征,同時給出港口擁塞指數(shù)時間序列尺度函數(shù)如圖4。

        圖4 港口擁塞指數(shù)時間序列尺度函數(shù)Fig.4 Port congestion index time series scale function chart

        由圖4可見,巴西港口ETA海岬型船與中國港口ETA海岬型船港口擁塞指數(shù)尺度函數(shù)τ(q)均為單調遞增的凹向橫軸的函數(shù),證明實驗時間序列是多重分形的。MF-DCCA對比結果顯示,三者為基本重合且凹向橫軸的單調遞增曲線,更進一步說明二者之間的多重分形交叉相關特性。

        此外,根據(jù)式(10)和式(11)計算給出港口擁塞指數(shù)時間序列分形譜,如圖5。

        圖5 港口擁塞指數(shù)時間序列分形譜Fig.5 Fractal spectrum of port congestion index in time series

        由圖5可知,巴西港口ETA海岬型船港口擁塞指數(shù)時間序列比中國港口ETA海岬型船港口擁塞指數(shù)時間序列分形譜更寬,且左端明顯高于右端,幅度更大。說明巴西港口ETA海岬型船港口擁塞指數(shù)時間序列波動比中國港口ETA海岬型船港口擁塞指數(shù)時間序列波動較大,兩者港口擁塞指數(shù)均表現(xiàn)為出現(xiàn)高位置時概率大,而巴西港口此特征更為明顯。MF-DCCA分析結果表示,兩組實驗時間序列多重分形交叉相關,因此港口擁塞指數(shù)時間序列受自身歷史數(shù)據(jù)影響,同時受另外一組港口擁塞指數(shù)實驗時間序列歷史數(shù)據(jù)影響。

        3 結 語

        通過對中國港口ETA海岬型船與巴西港口ETA海岬型船港口擁塞指數(shù)進行多重分形分析,結果發(fā)現(xiàn)中國港口和巴西港口擁塞指數(shù)的廣義Hurst指數(shù)均單調遞減且值域在(0.5,1)之間,并且巴西港口的擁塞指數(shù)的分形譜明顯大于中國港口的多重分形譜。這說明兩個國家的港口擁塞指數(shù)都具有多重分形的特征和長程相關的特征,并且巴西港口的多重分析特征相對更明顯。此外,還對兩個國家港口擁塞指數(shù)的互相關行為進行研究,結果發(fā)現(xiàn)兩者之間具有明顯的長程互相關特征。這些結果對研究不同港口擁塞指數(shù)的相互行為具有很好的參考價值,從而為研究港口擁堵問題提供依據(jù)。

        [1] 張榮忠. 全球性港口擁堵現(xiàn)象透析[J]. 交通企業(yè)管理,2005,20(5):31-32. ZHANG Rongzhong. Analysis of the global port congestion phenomenon[J].TransportationEnterpriseManagement,2005,20(5):31-32.

        [2] 尹銘鈺. 全球干散貨港口擁堵現(xiàn)狀分析[J]. 中國港口,2013(11):4-7. YIN Mingyu. Analysis of the current situation of global dry bulk port congestion[J].ChinaPorts,2013(11):4-7.

        [3] 尹銘鈺. 巴西港口在探索變革中成長[J]. 中國港口,2014(8):61-62. YIN Mingyu. The development of Brazil port in the exploration of the changes[J].ChinaPorts,2014(8):61-62.

        [4] 徐劍華. 美西港口擁堵迷局[J]. 中國船檢,2015(3):18-22. XU Jianhua. West Port congestion puzzles[J].ChineseShipSurvey,2015(3):18-22.

        [5] LEACH P T. 美西港口擁堵余波未平[J]. 中國遠洋航務,2015(4):60-61. LEACH P T. West Port congestion are still repercussions[J].MaritimeChina,2015 (4):60-61.

        [6] LEACH P T. 重新定義航運旺季[J]. 中國遠洋航務,2015(8):48-49. LEACH P T. Redefining shipping season[J].MaritimeChina,2015(8):48-49.

        [7] PENG C K, BULDYREV S V, HAVLIN S,et al. Mosaic organization of DNA nucleotides[J].PhysicalReviewE,1994,49(2):1685-1689.

        [8] KANTELHARDT J W, ZSCHIEGNER S A, KOSCIELNY B E, et al. Multifractal detrended fluctuation analysis of nonstationary time series[J].PhysicaAStatisticalMechanics&ItsApplications,2002,316(1/2/3/4):87-114.

        [9] PODOBNIK B, STANLEY H E. Detrended cross-correlation analysis: a new method for analyzing two nonstationary time series[J].PhysicalReviewLetters,2007,100(8):38-71.

        [10] HE L Y, CHEN S P. Nonlinear bivariate dependency of price-volume relationships in agricultural commodity futures markets: A perspective from Multifractal Detrended Cross-Correlation Analysis[J].PhysicaAStatisticalMechanics&ItsApplications,2011,390(2):297-308.

        [11] WANG Jun, ZHAO Daqing. Detrended cross-correlation analysis of electroencephalogram[J].ChinesePhysicsB,2012,21(2):577-580.

        [12] 曾志堅, 張倩倩. 基于MF-DCCA方法的證券市場間交叉相關性研究[J]. 財經(jīng)理論與實踐,2013,34(6):45-49. ZENG Zhijian, ZHANG Qianqian. Study on the cross correlation of securities market based on MF-DCCA method[J].TheTheoryandPracticeofFinanceandEconomics,2013,34(6):45-49.

        [13] 曹建軍,顧榮寶. 基于MF-DFA方法的黃金市場多重分形分析[J]. 黃金,2011,32(3):5-7. CAO Jianjun, GU Rongbao. Multi fractal analysis of gold market based on MF-DFA method[J].Gold,2011,32 (3):5-7.

        [14] 郭麗娜,王宏勇. 大宗商品期貨市場收益率的多重分形分析 [J]. 南京財經(jīng)大學學報,2013,20(1):43-50. GUO Li’na, WANG Hongyong. Return rate of commodity futures market multifractal analysis[J].JournalofNanjingUniversityofFinanceandEconomics,2013,20(1):43-50.

        [15] 張德園,王雪飛. 基于多重分形理論的原油市場耦合關系研究 [J]. 成都理工大學學報(社會科學版),2014,22(5):58-64. ZHANG Deyuan, WANG Xuefei. Research on the coupling relationship of crude oil markets based on multifractal theory[J].JournalofChengduUniversityofTechnology(SocialSciences) ,2014,22(5):58-64.

        [16] 劉俊超,陳秀乾. 上海集裝箱運價衍生品市場多重分形特征研究[J]. 航海,2014(1):47-51. LIU Junchao, CHEN Xiuqian. Research on multi fractal characteristics of container transport market in Shanghai[J].Navigation,2014(1):47-51.

        [17] SEO S K, KIM K, CHANG K H, et al. Determination of the dynamical behavior of rainfalls by using a multifractal detrended fluctuation analysis[J].JournaloftheKoreanPhysicalSociety,2012,61(4):658-661.

        [18] LIM G, KIM K, PARK J K, et al. Dynamical analyses of the time series for the temperature and the humidity[J].JournaloftheKoreanPhysicalSociety,2013,62(1):193-196.

        [19] 胡江,蘇懷智,馬福恒,等. MF-DFA在大壩安全監(jiān)測序列分析和整體性態(tài)識別中的應用[J]. 水利水電科技進展,2014,34(3):50-55. HU Jiang, SU Huaizhi, MA Fuheng, et al. The application of MF-DFA in time series analysis and global state recognition of dam safety monitoring[J].AdvancesinScienceandTechnologyofWaterResources,2014,34(3):50-55.

        [20] WANG F, LIAO G P, ZHOU X Y, et al. Multifractal detrended cross-correlation analysis for power markets[J].NonlinearDynamics,2013,72(1/2):353-363.

        [21] XU N, SHANG P, KAMAE S. Modeling traffic flow correlation using DFA and DCCA[J].NonlinearDynamics,2010,61(1/2):207-216.

        [22] 魏宇,黃登仕. 中國股票市場波動持久性特征的DFA分析[J]. 中國管理科學,2004,12(4):13-20. WEI Yu, HUANG Dengshi. A DFA study on the persistence of fluctuations in china stock market[J].ChineseJournalofManagementScience,2004,12(4):13-20.

        Multifractal Research on Congestion Index of Chinese and Brazil Port

        LIN Guolong, MA Libo, ZHANG Chenyan, HE Hongdi

        (Logistics Research Center, Shanghai Maritime University, Shanghai 201306, P.R.China)

        The multifractal detrended fluctuation analysis method (MF-DFA) and multifractal detrended fluctuation cross correlation analysis method (MF-DCCA) were used to study the congestion index of ETA (Etimated Ttime of Arrival) capesize in Chinese ports and in Brazil ports, which were components of GPCI (the global port congestion index). It is found that generalized Hurst index of ETA capesize in Chinese ports and in Brazil ports in time series and MF-DCCA generalized Hurst index are both monotonically decreasing and the range of their values is between 0.5 and 1; scaling index shows a concave function and fractal spectrum shows a intensive bell shaped distribution. Research results show that the two groups of experimental time series both have obvious multifractal characteristics and also have a long memory and multi-marked degree; moreover, the two multifractals are cross-correlated, and their performances are mutually influenced in a long-term.

        traffic and transportation engineering; port congestion index; multifractal; detrended fluctuation analysis; cross correlation analysis

        2015-08-21;

        2016-03-15

        林國龍(1951—),男,浙江象山人,教授,博士生導師,主要從事國際航運、自由貿易區(qū)、綜合物流環(huán)境與供應鏈一體化管理方面的研究。E-mail:LingLzm@163.com。

        10.3969/j.issn.1674-0696.2016.06.30

        U691+.32

        A

        1674-0696(2016)06-148-05

        猜你喜歡
        長程岬型分形
        長程動態(tài)心電圖對心律失常的檢出率分析
        感受分形
        分形之美
        分形空間上廣義凸函數(shù)的新Simpson型不等式及應用
        長程電子關聯(lián)對聚合物中激子極化率的影響
        基于分形理論的一種新的機器學習方法:分形學習
        亚洲婷婷丁香激情| 中国黄色一区二区三区四区| 中文字幕一区二区综合| 亚洲深深色噜噜狠狠网站| 97人伦色伦成人免费视频| 欧美亚洲国产片在线播放| 麻豆五月婷婷| 抖射在线免费观看视频网站 | 在线观看免费视频发布白白色| 精品人妻一区二区三区视频| 护士的小嫩嫩好紧好爽| 亚洲av无码成人专区片在线观看| 国产精品99久久精品爆乳| 久久精品综合国产二区| av国产免费在线播放| 精品一区二区三区芒果| 日韩乱码人妻无码中文字幕久久 | 午夜精品久久久久久| 中文字幕一区二区人妻出轨 | 最新国产乱人伦偷精品免费网站| 久久精品人人做人人爽电影蜜月| 日韩精品欧美激情亚洲综合| 免费女同毛片在线不卡| 精品人妻一区二区三区视频 | 中文字幕+乱码+中文字幕一区| 久久亚洲国产中v天仙www| 妺妺窝人体色www在线直播| 青青草伊人视频在线观看| 久久久精品人妻一区二区三区游戏| 国产精品白浆在线观看免费| 久热综合在线亚洲精品| 亚洲www视频| 免费国产不卡在线观看| 国模gogo无码人体啪啪| 在线综合亚洲欧洲综合网站 | 无码啪啪熟妇人妻区| 亚洲国产国语对白在线观看| 国产亚洲精品美女久久久m| 日本丰满人妻xxxxxhd| 九月色婷婷免费| 日韩五码一区二区三区地址|