999精品在线视频,手机成人午夜在线视频,久久不卡国产精品无码,中日无码在线观看,成人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)氣候,多臺(tái)風(fēng)、暴雨等災(zāi)害天氣,年均降雨量1 879.8 mm、降雨日數(shù)145 d、暴雨11次左右,全年85%的雨量、90%的暴雨日數(shù)出現(xiàn)在4—9月。復(fù)雜的自然環(huán)境、氣候條件和人類活動(dòng),使深圳洪澇、崩塌、滑坡等災(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ǔ)平臺(tái)開發(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

猜你喜歡
方法
中醫(yī)特有的急救方法
中老年保健(2021年9期)2021-08-24 03:52:04
高中數(shù)學(xué)教學(xué)改革的方法
化學(xué)反應(yīng)多變幻 “虛擬”方法幫大忙
變快的方法
兒童繪本(2020年5期)2020-04-07 17:46:30
學(xué)習(xí)方法
用對方法才能瘦
Coco薇(2016年2期)2016-03-22 02:42:52
最有效的簡單方法
山東青年(2016年1期)2016-02-28 14:25:23
四大方法 教你不再“坐以待病”!
Coco薇(2015年1期)2015-08-13 02:47:34
賺錢方法
捕魚
主站蜘蛛池模板: 色欲色欲久久综合网| 日韩免费毛片视频| 日韩精品欧美国产在线| 日本精品视频| 中文字幕免费视频| 国产H片无码不卡在线视频| 国产日韩久久久久无码精品| 精品久久香蕉国产线看观看gif| 国产欧美日韩一区二区视频在线| 欧洲av毛片| 亚洲精品成人福利在线电影| 国产视频你懂得| 国产激情在线视频| 色爽网免费视频| 精品视频第一页| 国产美女免费| 久久特级毛片| 久久无码免费束人妻| 免费一极毛片| 日韩久草视频| 亚洲香蕉久久| 欧美在线国产| 免费av一区二区三区在线| 精品无码一区二区在线观看| 欧美亚洲激情| 亚洲国产精品成人久久综合影院| 中文字幕1区2区| 日韩精品一区二区三区免费在线观看| 青青操视频在线| 久久天天躁狠狠躁夜夜2020一| 黄色一级视频欧美| 18禁色诱爆乳网站| 国产欧美精品一区aⅴ影院| 日韩精品无码免费一区二区三区 | 国产美女91呻吟求| 九九热在线视频| 欧美另类图片视频无弹跳第一页 | 国产SUV精品一区二区| 丁香六月激情婷婷| 深夜福利视频一区二区| 99久久精品免费观看国产| 久久人体视频| 四虎国产精品永久一区| 免费A级毛片无码无遮挡| 国产亚洲视频免费播放| 国产网站一区二区三区| AⅤ色综合久久天堂AV色综合 | 国产成在线观看免费视频| 亚洲AⅤ永久无码精品毛片| 免费一级毛片不卡在线播放| 在线观看免费国产| 小说 亚洲 无码 精品| 2021国产乱人伦在线播放| 怡红院美国分院一区二区| 亚洲日本中文字幕乱码中文| 992tv国产人成在线观看| 国产一区二区网站| 中文字幕人妻无码系列第三区| 国产香蕉在线视频| 1024你懂的国产精品| 日本黄色不卡视频| 在线观看免费黄色网址| 99久久免费精品特色大片| 呦女亚洲一区精品| 成人小视频网| 人妻精品久久无码区| 91免费国产高清观看| a级毛片网| 国产在线视频福利资源站| 国产综合网站| 亚洲一区二区三区中文字幕5566| 亚洲成人77777| 人妻丰满熟妇啪啪| 欧美成人一级| 色香蕉影院| 全裸无码专区| 国产精品美女网站| 综合久久五月天| 青青青国产在线播放| 国产精品性| 国产亚洲成AⅤ人片在线观看| 亚洲第一在线播放|