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

        ?

        五種常用降水量插值方法誤差時(shí)空分布特征研究——以深圳市為例

        2010-12-28 03:19:00吳小娟肖晨超
        地理與地理信息科學(xué) 2010年3期
        關(guān)鍵詞:插值法雨量插值

        鄔 倫,吳小娟,肖晨超,田 原,3*

        (1.北京大學(xué)遙感與地理信息系統(tǒng)研究所,北京 100871;2.武漢大學(xué)遙感信息工程學(xué)院,湖北武漢430079; 3.香港理工大學(xué)土地測量及地理資訊學(xué)系,香港)

        五種常用降水量插值方法誤差時(shí)空分布特征研究
        ——以深圳市為例

        鄔 倫1,吳小娟2,肖晨超1,田 原1,3*

        (1.北京大學(xué)遙感與地理信息系統(tǒng)研究所,北京 100871;2.武漢大學(xué)遙感信息工程學(xué)院,湖北武漢430079; 3.香港理工大學(xué)土地測量及地理資訊學(xué)系,香港)

        在洪水、滑坡等地質(zhì)災(zāi)害預(yù)警預(yù)報(bào)中,通常需要對多個(gè)時(shí)點(diǎn)、多個(gè)站點(diǎn)的降水量觀測數(shù)據(jù)進(jìn)行高精度插值,雨量插值精度對災(zāi)害預(yù)警預(yù)報(bào)的可靠度具有很大的影響,因此研究降水量插值方法誤差的時(shí)空分布特征具有重要的科研和實(shí)用價(jià)值。該文以深圳市2008年6月12日至14日百年一遇的強(qiáng)降水過程為例,采用交叉驗(yàn)證方法對反距離權(quán)重法、普通克里金法、全局多項(xiàng)式法、局部多項(xiàng)式法和徑向基函數(shù)法五種常用空間插值方法誤差的時(shí)空分布特征進(jìn)行分析,研究成果可為根據(jù)雨量時(shí)空分布特點(diǎn)選取適用雨量插值模型提供相關(guān)依據(jù),并為相關(guān)研究提供借鑒。

        降水量;插值誤差;誤差時(shí)空分布特征;地質(zhì)災(zāi)害;預(yù)警預(yù)報(bào)

        0 引言

        高質(zhì)量的降水?dāng)?shù)據(jù)在水資源管理、水文建模、氣候變化、地質(zhì)災(zāi)害預(yù)警預(yù)報(bào)等方面具有重要應(yīng)用[1]。實(shí)踐中,雨量觀測站點(diǎn)的布設(shè)受到地形、人力、財(cái)力等因素的限制,直接觀測數(shù)據(jù)無法精確表達(dá)降水量的時(shí)空分布,通常需要采用空間插值方法對站點(diǎn)觀測數(shù)據(jù)進(jìn)行插值,以獲取整個(gè)區(qū)域內(nèi)降水量的總體分布數(shù)據(jù)。目前,空間插值方法在降水量插值方面得到了廣泛應(yīng)用。徐超在山東省境內(nèi)分別采用反距離權(quán)重法、徑向基函數(shù)法和普通克里金法對多年氣象要素進(jìn)行了空間插值分析,發(fā)現(xiàn)普通克里金法的插值效果更理想[2];朱芮芮等對日降雨量的時(shí)空變異特征進(jìn)行分析,得出普通克里金法和反距離權(quán)重法總體效果較好[3];Bussires等在日累計(jì)降水量的插值研究中發(fā)現(xiàn)地統(tǒng)計(jì)學(xué)克里金法優(yōu)于簡單的泰森多邊形法和反距離權(quán)重法[4];Dirks等比較了克里金法、反距離權(quán)重法、泰森多邊形法在年、月、日、時(shí)四種時(shí)間分辨率情況下的插值結(jié)果,發(fā)現(xiàn)克里金法插值效果最好[5]。

        前人研究表明,由于時(shí)間及空間尺度不同,適用的空間插值方法和降水模型也不相同[2-5]。目前降水量插值方法的空間和時(shí)間誤差分布研究主要集中在空間插值方法對區(qū)域總體年、月、日累計(jì)降水量的插值效果方面,對于洪水、滑坡等地質(zhì)災(zāi)害預(yù)警預(yù)報(bào)工作需要的更高時(shí)間和空間分辨率下插值方法的誤差時(shí)空分布很少涉及。實(shí)踐中,一次特大洪水的形成不但與前期累計(jì)雨量有關(guān),也與致洪暴雨的時(shí)空分布特別是與暴雨量持續(xù)時(shí)間的長短有關(guān)[6];降水誘發(fā)的滑坡、崩塌等地質(zhì)災(zāi)害發(fā)生的時(shí)間和區(qū)域與強(qiáng)降雨的時(shí)空分布也具有耦合性[7-9]。在這些應(yīng)用中,不僅需要獲得少數(shù)重點(diǎn)時(shí)點(diǎn)的累計(jì)雨量和分布情況,更需要獲得高空間和時(shí)間分辨率的全過程雨量分布數(shù)據(jù)。因此,有必要對常見雨量插值方法在高空間和時(shí)間分辨率下的誤差分布進(jìn)行研究。

        1 插值方法及其精度評價(jià)

        1.1 插值方法

        空間插值是利用已知的部分空間樣本信息對未知的地理空間特征進(jìn)行估計(jì),也是 GIS的重要功能模塊之一[10]。Dirks將空間插值方法歸結(jié)為整體插值法和局部插值法[5],Vicente-Serrano等將插值方法細(xì)分為整體插值法、局部插值法、地學(xué)統(tǒng)計(jì)法和混合插值法[11]。雨量觀測數(shù)據(jù)均為點(diǎn)數(shù)據(jù),在對點(diǎn)數(shù)據(jù)進(jìn)行插值時(shí),常用的插值方法有五種:1)反距離權(quán)重插值法(Inverse Distance Weighting,IDW),以插值點(diǎn)與樣本點(diǎn)間的距離為權(quán)重進(jìn)行加權(quán)平均,離插值點(diǎn)越近的樣本點(diǎn)賦予的權(quán)重越大(本文采用反距離平方加權(quán)法[12]);2)普通克里金插值法(Ordinary Kriging,OK),是區(qū)域化變量的線性估計(jì),但其考慮了空間相關(guān)性,更加符合空間數(shù)據(jù)的特點(diǎn)(本文用球形模型[13,14]);3)全局多項(xiàng)式插值法(Global Polyno-m ial,GP),以整個(gè)研究區(qū)的樣點(diǎn)數(shù)據(jù)集為基礎(chǔ),用一個(gè)多項(xiàng)式計(jì)算預(yù)測值(本文采用二次多項(xiàng)式進(jìn)行擬合);4)局部多項(xiàng)式插值法(Local Polynomial,LP),采用處在特定重疊鄰近區(qū)域內(nèi)的多個(gè)多項(xiàng)式進(jìn)行插值,其產(chǎn)生的曲面更依賴于局部數(shù)據(jù)的變異(本文采用二次多項(xiàng)式進(jìn)行擬合);5)徑向基函數(shù)插值法(Radial Basis Function,RBF),是多個(gè)數(shù)據(jù)插值方法的組合,即經(jīng)過各個(gè)已知樣點(diǎn)生成一個(gè)圓滑曲面,并使表面的總曲率最小[15](本文采用張力樣條函數(shù)作為基函數(shù),并設(shè)置表面光滑度參數(shù)為2)。

        1.2 插值精度評價(jià)

        本文采用交叉驗(yàn)證法[16]評價(jià)各種方法的插值效果。研究中將50個(gè)站點(diǎn)隨機(jī)排列并分為5組,每組10個(gè)站點(diǎn),依次假設(shè)其中一組站點(diǎn)的降水量未知,用其余40個(gè)站點(diǎn)的降水量進(jìn)行插值,得到所有站點(diǎn)所有時(shí)段的降水估計(jì)值,最后計(jì)算每個(gè)站點(diǎn)每時(shí)段實(shí)際觀測值與估算值的誤差。

        誤差評價(jià)指標(biāo)采用平均絕對誤差(MAE)、平均相對誤差(M RE)和均方根誤差(RM SE)。其中MAE反映估計(jì)值的誤差范圍,M RE反映估計(jì)值對于觀測值的準(zhǔn)確度,RM SE則反映估計(jì)值的靈敏度和極值情況[1,17]。表達(dá)式為:

        其中:Za,i、Ze,i分別為第i個(gè)站點(diǎn)的實(shí)際觀測值和插值預(yù)測值,n為驗(yàn)證站點(diǎn)數(shù)。

        分別從時(shí)間和空間兩個(gè)角度對上述插值方法的誤差進(jìn)行統(tǒng)計(jì)分析。首先計(jì)算全區(qū)各時(shí)點(diǎn)雨量觀測值與插值結(jié)果之間的誤差,以分析各插值方法誤差在時(shí)間序列上的分布特征;之后計(jì)算各站點(diǎn)在整個(gè)降水過程中各時(shí)點(diǎn)觀測值與估計(jì)值間的誤差,以分析不同插值方法誤差的空間分布特征。研究中繪制了MAE、M RE和RM SE在時(shí)間和空間序列上的變化曲線,并計(jì)算MAE、MRE和RM SE的均值、標(biāo)準(zhǔn)差、最大值和最小值,用以表達(dá)各種插值方法誤差的時(shí)空分布特征,以及各插值方法誤差與雨量極值出現(xiàn)的時(shí)間與空間位置間的關(guān)系。

        2 實(shí)例研究

        2.1 研究區(qū)概況及研究數(shù)據(jù)

        深圳地處我國東南沿海,屬亞熱帶季風(fēng)氣候,多臺風(fēng)、暴雨等災(zāi)害天氣,年均降雨量1 879.8 mm、降雨日數(shù)145 d、暴雨11次左右,全年85%的雨量、90%的暴雨日數(shù)出現(xiàn)在4—9月。復(fù)雜的自然環(huán)境、氣候條件和人類活動,使深圳洪澇、崩塌、滑坡等災(zāi)害頻發(fā)[18]。

        本研究以深圳市2008年6月12日至14日百年一遇強(qiáng)降水過程為例,選取降水較為集中的深圳市西部地區(qū)作為研究區(qū),以研究區(qū)內(nèi)工作狀態(tài)良好、空間分布相對均勻的50個(gè)雨量站的觀測數(shù)據(jù)作為雨量觀測數(shù)據(jù)(圖1)。降水持續(xù)36 h,以1 h為間隔共劃分為36個(gè)時(shí)段。降水量最大值出現(xiàn)在第24個(gè)時(shí)段,高達(dá)1 602.2 mm;最小值出現(xiàn)在第34個(gè)時(shí)段,為13 mm。為了分析誤差的空間分布特征,對50個(gè)站點(diǎn)按自西向東的順序依次編號,其中13號站點(diǎn)累計(jì)降水量最大(526.9 mm),31號站點(diǎn)累計(jì)降水量最小(105.3 mm)

        圖1 研究區(qū)位置與雨量計(jì)分布Fig.1 The distribution of rainfall gauges in the study area

        2.2 插值過程及結(jié)果分析

        筆者以A rcGIS軟件為基礎(chǔ)平臺開發(fā)了交叉檢驗(yàn)程序。使用的插值模塊為A rcGIS的 Geostatistical Analyst,按照交叉檢驗(yàn)方案對每個(gè)時(shí)點(diǎn)的觀測數(shù)據(jù)進(jìn)行了插值。

        2.2.1 插值誤差時(shí)間序列分布特征 利用交叉檢驗(yàn)方法對各個(gè)時(shí)點(diǎn)的插值結(jié)果進(jìn)行誤差統(tǒng)計(jì),計(jì)算出五種方法插值結(jié)果的M AE、M RE和RM SE(圖2)。表 1列出了五種插值方法在時(shí)間序列上的MAE、M RE和RM SE與降水量的相關(guān)系數(shù)及各自的平均值、標(biāo)準(zhǔn)差、最大值和最小值。

        對MAE的分析可以得出:OK方法的平均值、標(biāo)準(zhǔn)差和最大值均最小,說明其插值誤差范圍最小; RBF、IDW和LP方法次之;GP方法誤差范圍最大。對MRE的分析可以得出:LP方法的估值誤差與觀測值的比值最小,OK方法的比值最大。RMSE的平均值、標(biāo)準(zhǔn)差和最大值的比較結(jié)果具有一致性,即GP>LP>IDW>RBF>OK。表1中相關(guān)系數(shù)表明:所有方法的MA E與RM SE和降水量呈正相關(guān),表明降雨量大的時(shí)刻誤差范圍和誤差極值都較大;而MRE與降水量呈負(fù)相關(guān),表明降雨量大的時(shí)刻雨量估計(jì)值相對于觀測值的準(zhǔn)確度較高。

        圖2 五種插值方法插值誤差的時(shí)間序列分布Fig.2 Temporal distributions of errors of five interpolation models

        表1 五種插值方法時(shí)間序列誤差指標(biāo)分析Table 1Analysis on temporal distributions of errors of five interpolation models

        進(jìn)一步比較各種插值方法在降水量極值時(shí)點(diǎn)的插值效果,表2列出各種方法的 M AE、M RE和RMSE。在降水量最大時(shí)段,OK方法誤差最小,GP方法誤差最大;在降水量最小時(shí)段,GP方法誤差最小,而RBF方法誤差最大。

        綜合上述分析,從時(shí)間序列看五種插值方法中OK方法的誤差總體最小,這與文獻(xiàn)[3-5]的研究結(jié)果一致,且其在降水量最大的時(shí)段誤差遠(yuǎn)小于其他4種方法;RBF方法和IDW方法的誤差居中;LP方法和GP方法的誤差整體偏大,但較適合于降水量小的時(shí)段。

        2.2.2 插值方法誤差空間分布特征 依次對每個(gè)站點(diǎn)各時(shí)點(diǎn)插值結(jié)果進(jìn)行誤差統(tǒng)計(jì),計(jì)算出每個(gè)站點(diǎn)在整個(gè)降水過程中的MAE、MRE和RM SE(圖3)。表3列出了各種插值方法插值結(jié)果在各個(gè)站點(diǎn)MAE、M RE和RM SE與降水量的相關(guān)系數(shù)及各自的平均值、標(biāo)準(zhǔn)差、最大值和最小值。各插值方法在各站點(diǎn)MA E、M RE和RM SE的誤差分布表現(xiàn)與其在時(shí)間序列上的誤差表現(xiàn)基本一致:OK方法的誤差最小;RBF、IDW和LP方法誤差居中;GP方法插值誤差最大,且極值情況最為嚴(yán)重。表3中相關(guān)系數(shù)總體規(guī)律與表1類似:所有插值方法在空間分布上的MAE和RM SE均與對應(yīng)站點(diǎn)的降水量呈正相關(guān),表明在降水量較大的站點(diǎn)誤差范圍和誤差極值都較大;而MRE與對應(yīng)站點(diǎn)的降水量呈負(fù)相關(guān),表明在降水量較大的站點(diǎn)雨量估計(jì)值相對于觀測值的準(zhǔn)確度較高。

        為了進(jìn)一步分析各插值方法在雨量極值站點(diǎn)的插值表現(xiàn),表4列出了各種插值方法在降水量最大的站點(diǎn)和降水量最小的站點(diǎn)上的MAE、MRE和 RM SE。在降水量最大的站點(diǎn),OK方法總體誤差最小,GP方法誤差最大;在降水量最小的站點(diǎn),IDW方法誤差最小,GP方法誤差最大。

        表3 五 種插值方法空間分布誤差指標(biāo)分析Table 3Analysis on spatial distributions of errors of five interpolation models

        表4 五種插值方法在降水量最大、最小站點(diǎn)誤差值Table 4 Errors of five interpolation modelsat themaximum and m inimum precipitation locations

        3 結(jié)論與展望

        對五種常見雨量插值方法誤差的時(shí)空分布特征進(jìn)行對比分析,可得出以下結(jié)論:1)在時(shí)間序列和空間點(diǎn)位分布上,OK方法總體誤差最小,RBF和IDW方法總體誤差次之,LP和 GP方法總體誤差最大。2)在時(shí)間序列和空間點(diǎn)位分布上,各種插值方法的MAE和RM SE均與對應(yīng)時(shí)刻或站點(diǎn)的降水量呈正相關(guān),M RE則與之呈負(fù)相關(guān)。3)在時(shí)間序列上,在降水量最大的時(shí)段,OK方法誤差最小,GP方法誤差最大;在降水量最小的時(shí)段,GP方法誤差最小,而RBF方法誤差最大。4)在空間分布上,在降水量最大的站點(diǎn),OK方法誤差最小,GP方法誤差最大,這與各插值方法的雨量最大值時(shí)刻的誤差對比結(jié)果一致;在降水量最小的站點(diǎn),IDW方法誤差最小,GP方法誤差最大。

        五種常見雨量插值模型的插值誤差有著不同的時(shí)空分布特征,總體而言,O K方法適用于降水量較大的時(shí)段和區(qū)域,IDW方法則適用于降水量較小的時(shí)段和區(qū)域,在選用時(shí)可以綜合考慮。為取得更好的插值精度,還可以進(jìn)一步考慮在一個(gè)降水過程中針對不同的降雨時(shí)段或區(qū)域采用不同的插值方法。

        [1] 宋麗瓊.日降水量的空間插值方法與應(yīng)用對比分析——以深圳市為例[J].地球信息科學(xué),2008,10(5):566-572.

        [2] 徐超,吳大千,張治國.山東省多年氣象要素空間插值方法比較研究[J].山東大學(xué)學(xué)報(bào)(理學(xué)報(bào)),2008,43(3):1-5.

        [3] 朱芮芮,李蘭,王浩,等.降水量的空間變異性和空間插值方法的比較研究[J].中國農(nóng)村水利水電,2004(7):25-28.

        [4] BUSSIRESN,HOGGW.The objective analysisof daily rainfall by distance weighting schemes on a mesoscale grid[J].A tmosphere-Ocean,1989,27(3):521-541.

        [5] DIRKS K N,HAY J E,STOW CD,et al.High-resolution studies of rainfall on Norfolk Island.Part II:Interpolation of rainfall data[J].Hydrology,1998,208:187-193.

        [6] 趙玉春,王仁喬.一次致洪暴雨的中尺度分析[J].氣象科技, 2005,33(3):245-249.

        [7] BRAND E W,PREMCH ITT J,PH ILL IPSON H B.Relationship between rainfall and landslides in Hong Kong[A].Proceedings of the Fourth International Symposium on Landslides, To ronto,1984,1:377-384.

        [8] JIA G Y,TIAN Y,L IU Y,et al.A static and dynamic factorscoupled forecasting model of regional rainfall-induced landslides:A case study of Shenzhen[J].Science in China Series ETechnological Sciences,2008,51:164-175.

        [9] TIAN Y,XIAO C C,L IU Y,et al.Effects of raster resolution on landslide susceptibility mapping:A case study of Shenzhen [J].Science in China Series E-Technological Sciences,2008,51: 188-198.

        [10] LAM N.Spatial interpolation methods:A review[J].The American Cartographer,1983,10(2):129-150.

        [11] V ICENTE-SERRANO S M,SAZ-SáNCHEZ M A,CUADRAT J M.Comparative analysisof interpolationmethods in themiddle Ebro Valley(Spain):Application to annual p recipitation and temperature[J].Climate Research,2003,24(2):161-180.

        [12] PATRICK M B,KELLER C P.Multivariate interpolation to incorporate thematic surface date using inverse distance weighting (IDW)[J].Computer&Geosciences,1996,22(7):795-799.

        [13] PARDO-IGUZQU IZA E,DOWD P A.The second-order stationary universal Kriging model revisited[J].Mathematical Geology,1998,30(4):347-378.

        [14] KARN IEL IA.Application of Kriging technique to areal p recipitation mapping in A rizona[J].GeoJournal,1990,22(4): 391-398.

        [15] 王樹良,史文中,李德毅,等.基于張力樣條插值函數(shù)的土地?cái)?shù)據(jù)挖掘[J].計(jì)算機(jī)工程與應(yīng)用,2003,39(25):5-7.

        [16] SEAMAN R S.Objective analysis accuracies of statistical in-terpolation and successive correction schemes[J].Australian Meteorological Magazine,1983,31:225-240.

        [17] 高華喜,殷坤龍.降雨與滑坡災(zāi)害相關(guān)性分析及預(yù)警預(yù)報(bào)閾值之探討[J].巖土力學(xué),2007,28(5):1055-1060.

        [18] 李朝奎,陳良,王勇.降水量分布的空間插值方法研究——以美國愛達(dá)荷州為例[J].礦產(chǎn)與地質(zhì),2007,21(6):684-687.

        On Temporal and Spatial Error Distributions of Five Precipitation Interpolation M odels:A Case of Shenzhen

        WU Lun1,WU Xiao-juan2,XIAO Chen-chao1,TIAN Yuan1,3
        (1.Institute of RS and GIS,Peking University,Beijing 100871;
        2.Institute of Remote Sensing and Inform ation Engineering,W uhan University,W uhan 430079;
        3.Department of Land Surveying and Geo-Informatics,Hong Kong Polytechnic University,Kow loon,Hong Kong,China)

        Precipitation isan impo rtantmeteorological elementw idely used in the fo recasting of many kindsof natural disasters, such as floods and landslides.Typically,interpolation models are app lied to calculate the p recipitation distribution over the wo rking area based on limited observation data.It is important to study the tempo ral and spatial error distribution of p recipitation interpolation models to get p recise p recipitation data and to make the disaster fo recasting more credible.Five typical interpolationmodels,Inverse Distance Weighted,Ordinary Kriging,Global Polynomial,Local Polynomial,and Radial Basis Function,are chosen to carry a case study of Shenzhen,based on observation data of a rainfall p rocess between June 12 and June 14,2008.It can be concluded that the erro rsof all the fivemethods are highly relevant to p recipitation.Considering both the time and space series distribution,the erro rsof Ordinary Krigingmethod areminimum of all and the erro rsof Local Polynomial and Global Polynomial are themaximum.A t the moment of maximum p recipitation and at the locations of maximum p recipitation,Ordinary Kriging method has theminimum errors w hile Global Polynomialmethod has themaximum errors.But at themoment of minimum p recipitation,the errors of Global Polynomial,Local Polynomial and Inverse Distance Weighted are smaller than that of Ordinary Kriging,w hile erro rsof Radial Basis Function is themaximum.On the location of minimum p recipitation,the errorsof Inverse Distance Weighted and Global Polynomial are maximum and minimum respectively.The conclusions of this study may p rovide useful guidance on choosing suitable interpolation model.

        p recipitation;interpolation erro r;tempo ral and spatial error distribution;geological disaster;forecasting and warning

        P208;P426.6

        A

        1672-0504(2010)03-0019-06

        2010-01-22;

        2010-03-31

        國家科技支撐計(jì)劃課題(2008BAJ11B04、2006BAJ14B04);國家自然科學(xué)基金項(xiàng)目(40701134、40928001);國家高科技研究發(fā)展計(jì)劃項(xiàng)目(2007AA 120502);香港理工大學(xué)研究基金(Project No.G-U632)

        鄔倫(1964-),男,教授,博士生導(dǎo)師,研究方向?yàn)榈乩硇畔⒖茖W(xué)。*通訊作者E-mail:tianyuanpku@pku.edu.cn

        猜你喜歡
        插值法雨量插值
        寧夏紅柳溝流域水沙變化及產(chǎn)沙分析
        《計(jì)算方法》關(guān)于插值法的教學(xué)方法研討
        基于小波去噪的稱重雨量數(shù)據(jù)分析
        基于Sinc插值與相關(guān)譜的縱橫波速度比掃描方法
        SL—1 型雨量傳感器故障分析排除和維護(hù)
        西藏科技(2016年5期)2016-09-26 12:16:40
        一種改進(jìn)FFT多譜線插值諧波分析方法
        基于四項(xiàng)最低旁瓣Nuttall窗的插值FFT諧波分析
        基于二次插值法的布谷鳥搜索算法研究
        Newton插值法在光伏發(fā)電最大功率跟蹤中的應(yīng)用
        Blackman-Harris窗的插值FFT諧波分析與應(yīng)用
        蜜桃视频免费在线视频| av一区二区三区人妻少妇| 18级成人毛片免费观看| 国产永久免费高清在线观看视频| 免费观看一区二区三区视频| 精品国产午夜肉伦伦影院| 一性一交一口添一摸视频| 亚洲另类欧美综合久久图片区 | 三级日本午夜在线观看| 亚洲av久播在线一区二区| 精品9e精品视频在线观看| 中文字幕不卡在线播放| 日韩女同一区在线观看| 精品人妻一区三区蜜桃| 久久精品国产亚洲av四虎| 久久99精品久久久久久齐齐百度 | 中文字幕人妻互换av| 国产综合久久久久久鬼色| 中文字幕无码无码专区| 久久久9色精品国产一区二区三区 国产三级黄色片子看曰逼大片 | 国产精品入口牛牛影视| 国产91九色视频在线播放| 久久久精品国产性黑人| 国产国拍亚洲精品mv在线观看| 久久精品中文字幕极品| 中文字幕乱码在线婷婷| 国产爆乳美女娇喘呻吟| 亚洲欧美日韩综合久久| 国产成人综合亚洲av| 激情久久黄色免费网站| 久久精品无码一区二区三区免费| 亚洲视频天堂| 一本之道加勒比在线观看| 国产自国产自愉自愉免费24区| 国模少妇一区二区三区| 尤物蜜芽福利国产污在线观看| 日本九州不卡久久精品一区 | 制服丝袜人妻中文字幕在线| 丰满人妻一区二区乱码中文电影网 | 亚洲国产精品亚洲一区二区三区 | 婷婷色国产精品视频二区|