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

天氣發生器MulGETS和k-NN對區域歷史氣象場特征重現能力的比較*

2019-06-17 09:33:16周凌峰孟耀斌伍甘霖張東妮宋昊政
中國農業氣象 2019年6期
關鍵詞:模型

周凌峰,孟耀斌**,逯 超,伍甘霖,張東妮,宋昊政,吳 丹

?

天氣發生器MulGETS和k-NN對區域歷史氣象場特征重現能力的比較*

周凌峰1,孟耀斌1**,逯 超1,伍甘霖1,張東妮1,宋昊政1,吳 丹2

(1.北京師范大學環境演變與自然災害教育部重點實驗室/北京師范大學應急管理部/教育部減災與應急管理研究院/北京師范大學地理科學學部,北京 100875;2.湖南瀏陽市水文水資源局,瀏陽 410300)

傳統單站點天氣發生器未考慮不同站點氣象變量間的空間相關性,導致其在區域影響評價中的應用受到限制,而多站點天氣發生器可以克服單站點天氣發生器的缺點,近年來得到迅速發展。評估和驗證多站點天氣發生器對區域歷史氣象場特征的重現能力是開展影響評價的前提和基礎。為此,本研究選取MulGETS(參數型)和k-NN(非參數型)發生器為代表模型,利用湘江流域12個氣象站點1981?2010年日序列降水量、最高氣溫、最低氣溫資料,通過均值、標準差、偏度、極值、空間相關系數、空間連接度和自相關系數等指標的對比,評估了MulGETS和k-NN模型的優缺點及適用性。結果表明:MulGETS和k-NN模型均較好地再現了原氣象場的均值、標準差和偏度,k-NN表現稍好于MulGETS。同時k-NN相比MulGETS在保持氣象要素空間相關性上具有優勢,特別是降水量的空間間歇性。由于算法本身的限制,k-NN無法模擬出超出歷史數據范圍的極值,而MulGETS具備一定的極值模擬能力。此外,MulGETS和k-NN在重現原始日尺度降水量的自相關性上均存在不足??傮w來看,兩個模型各具優勢和不足,MulGETS更適于極端氣象事件模擬,而k-NN可以更好地體現原始氣象場的空間差異,實際使用時應根據不同的研究目的選擇合適的模型。

多站點;天氣發生器;模型比較;湘江流域

天氣發生器(Weather Generators),又稱隨機天氣發生器(Stochastic Weather Generators),可用于研究某個地區天氣或氣候的一般特征,并根據這些統計特征隨機模擬該地區可能天氣序列。天氣發生器可用于隨機生成降水、最高氣溫、最低氣溫和輻射等多種氣象要素,已廣泛應用于水文循環、土壤侵蝕和作物生長等過程對氣候或氣象情景的響應評價[1]。

天氣發生器研究歷史悠久,幾十年來得到持續不斷的發展。以往使用比較廣泛的模型主要有WGEN[2]、LARS-WG[3]和CLIGEN[4]等,在國內也得到了廣泛應用[5?7]。但是,這些模型都屬于單站點模型,未考慮氣象要素的空間相關性,很難運用在區域影響評價中。將天氣發生器從單站點擴展到多站點是近年研究熱點和難點[8?9],Chen等[10]利用Wilks[11]提出的空間相關時間獨立隨機數方法保持降水量的空間相關性,開發了多站點天氣發生器MulGETS,擴展了傳統的WGEN型天氣發生器的應用范圍。Lall等[12?13]提出k近鄰(k-nearest neighbors,k-NN)重采樣方法,根據當前天氣狀況與歷史數據的相似性,從歷史相似數據中進行抽樣,實現了多站點隨機模擬,并保持了氣象要素的空間相關性。MulGETS和k-NN模型同屬于多站點天氣發生器,但建模原理差別較大,可分類為參數型方法和非參數型方法。MulGETS屬于傳統的參數型天氣發生器,降水是否發生由一階或二階馬爾可夫鏈模型來模擬,而降水量需要由假設的某種概率分布生成。k-NN屬于非參數型天氣發生器,不需假設分布,根據歷史數據與當前數據的相似性從歷史數據中重抽樣,以刻畫變量間的非線性特征見長。

目前國內的天氣發生器研究主要集中在單站點天氣發生器[5?7],多站點天氣發生器的應用、效能評估和模型比較仍是空白,因此有必要對不同多站點天氣發生器的特點、優勢和劣勢進行比較,并對其在中國的適用性進行驗證?;诖耍狙芯恳院舷娼饔驗檠芯繀^,以均值、標準差、偏度、極值、空間相關系數、空間連接度和自相關系數為指標,對比兩類天氣發生器MulGETS模型(參數型)與k-NN模型(非參數型)的模擬結果,評估兩個模型的優缺點及適用性,以期為合理使用現有多站點模型提供參考。

1 資料與方法

1.1 研究區域和數據

湘江是長江的重要支流,同時也是匯入洞庭湖流域的最大河流。湘江流域位于中國中南部,流域面積約9.46萬km2(圖1)。屬亞熱帶季風性氣候,雨量充沛,年平均降水量1200~1700mm,雨季在4?6月,夏季炎熱潮濕,冬季溫暖濕潤,多年平均氣溫在16~18℃。從中國國家氣象局(http://data.cma. cn/)收集、整理湘江流域1981?2010年共30a的逐日降水量、最高氣溫和最低氣溫資料,包括馬坡嶺、株洲、雙峰、南岳、衡陽、攸縣、永州、常寧、道縣、嘉禾、郴州和桂東共12個站點(表1)。各站30a數據累積缺測時間最多不超過11d,少數缺測1~2a的站點數據由附近站點的平均值補充。

圖1 研究區域(湖南湘江流域)12個氣象站點的分布

表1 研究區各氣象站點基本信息

Table 1 Basic information of the meteorological stations in study area

1.2 模型原理介紹

1.2.1 MulGETS模型

多站點多變量天氣發生器MulGETS模型首先模擬降水量,采用一階二狀態馬爾科夫鏈確定某一天是干日還是濕日,假設某日是否降水僅取決于前一天的狀態。如果為濕日,則還要進行日降水量模擬。與單站點參數型降水量模擬方法不同的是,MulGETS采用了Wilks[11]提出的時間獨立空間相關的隨機數構造方法,該方法可以較好地保持不同站點降水量的空間相關性。同時,MulGETS還將降水發生指數KS與Gamma分布的均值m和標準差s建立關系。其中,某天第S站的降水量發生指標KS定義為

式中,O為S站點的周邊11個站點某天是否發生降水的行向量,如[1,0,1,1,1,0,1,1,0,0,0],1代表有降水,0代表無降水。U為全1行向量,即[1,1,1,1,1,1,1,1,1,1,1,1],CT為S站點實測的0?1序列(是否降水)與周邊11個站點對應序列的相關系數列向量,如本研究中馬嶺坡站與周邊站點的相關系數為[0.78,0.71,0.64,0.63,0.62,0.62,0.56,0.57,0.46,0.42,0.45,0.41]T。

可以看到,計算得到的KS數值在[0,1]之間,且包含了周邊站點是否降水的信息。MulGETS按照KS的數值大小分成不同的區間,每個區間內的降水量都服從Gamma分布,即整體服從多元Gamma分布。為了保證每個區間有足夠的數據參與擬合,MulGETS假設季節內KS與Gamma分布均值m和標準差s的關系保持不變,所以最終可以得到不同站點不同季節的KS與m和s的關系。根據模擬的降水量發生指數KS來生成降水量可以一定程度上解決降水量的空間間歇性問題。

最高/最低氣溫采用傳統的多變量一階自回歸模型模擬,為了考慮氣溫的空間相關性,同樣采用上述時間獨立空間相關隨機數方法,具體方法見文獻[10]。

1.2.2 k-NN模型

Lall等[12]最早引入k-NN方法并將其運用在月徑流時間序列的隨機模擬中,Rajagopalan等[13]隨后利用k-NN方法模擬包含降水量、氣溫和輻射等的多變量時間序列,Buishand等[14?15]又將其擴充為多站點模擬方法。多站點多變量天氣發生器k-NN模型的具體算法分為以下步驟:

(2)從歷史數據N年的1月1日中隨機選擇1日作為模擬第1天的天氣值。

(3)選擇窗口期為w天,則可以確定共L天的相鄰天,其中

假設當前模擬到第t天,以t =1為例,即1月1日。w設為14d,所以30a中每年的12月25日至第二年的1月8日共449d(不包括當天)為相鄰天。

(4)計算當天天氣與備選L天的天氣狀態的相似程度,以馬氏距離(Mahalanobis distances)di為指標,St為協方差矩陣。

(6)建立累計概率密度函數cpj,用于從k個最近相鄰天中抽樣,距離di越近越容易被選中。

(7)產出0~1均勻分布隨機數,從備選的k天潛在最近相鄰天中隨機抽一天作為當前天的最近相鄰天,將該最近相鄰天的后一天作為第t+1天的模擬天氣值。

(8)重復步驟4?步驟7,直至達到模擬天數要求。

1.3 模擬效果評估

基于12個站點1981?2010年3個變量(降水量、最高氣溫和最低氣溫)組成的歷史氣象場,運行MulGETS和k-NN模型各生成20個30a的模擬氣象場,針對每個情景分別計算各氣象要素的均值、標準差、偏度、極值、自相關系數、空間相關系數和空間連接度指標,并與歷史氣象場對應指標進行比較,評估兩個模型重現歷史氣象場統計特征的能力?;窘y計量(均值、標準差和偏度)和極值分別考察模型重現氣象場平均態和極端狀態的能力??臻g相關系數和空間連接度用于評估模型重現氣象場空間結構特征能力。氣象場時間結構特征則由自相關性系數來表征。

由于MulGETS對降水的模擬是在季節尺度上建立的,為了達到綜合和公平評價的目標,按4個季節分別統計日均值、日標準差和日偏度,并采用平均相對誤差絕對值(mean absolute relative error,MARE)作為指標,MARE定義為

式中,Oi,j代表第i站(共12站)第j個季節(共4季)的觀測數據的統計指標(均值、標準差和偏度),Si,j,k代表對應的第i站(共12站)第j個季節(共4季)第k次(共20次)的模擬數據的統計指標(均值、標準差和偏度)。

氣象變量的空間相關系數(Spatial correlation,簡稱SC)指某氣象變量兩兩站點間的線性相關系數,可表示為

式中,x表示某氣象變量(降水量、最低最高氣溫),i和j代表不同的站點。

由于降水在空間和時間上存在高度的不連續性[11],是各氣象變量模擬中最難也是最重要的部分。為了衡量天氣發生器對降水場空間間歇性的模擬好壞,Wilks[11]定義了空間連接度指標(Continuity Ratio,簡稱CR)。

式中,x指降水量,i和j代表不同站點。空間連接度指標是一個比值,式(8)中分子表示j站點不下雨且i站點下雨時i站點雨量的均值,分母為j站點和i站點均下雨時i站點雨量的均值,比值大表示站點間是/否降水的相關性小,比值小表示站點間是/否降水具有一致性。

自相關性是降水時間序列的基本屬性,但在以往天氣發生器評價研究中考察較少。為了衡量天氣發生器對降水序列時間結構特征的重現能力,定義降水量序列的樣本自相關系數(Autocorrelation,簡稱AC)為

式中,xt表示某氣象變量的時間序列,xt?l表示該氣象變量滯后l個步長的時間序列,l表示延遲的步長數。

2 結果與分析

2.1 兩模型對歷史氣象場平均態重現能力的評估

圖2和圖3展現了MulGETS模型和k-NN模型對各站點日降水量、最高和最低氣溫均值、標準差和偏度的模擬值與實測值的對比,其中最高/最低氣溫的偏度很小所以不作對比??傮w來看(表2),k-NN相比于MulGETS能較好地再現日值降水量的統計特征。k-NN最大偏差出現在降水的偏度上(MARE=8.4%),MulGETS對應的降水偏度MARE為12.6%且一定程度上低估了降水偏度較大的站點(圖2c)。從模擬的氣象變量來看,氣溫場的MARE在0.6%~13.3%,而降水量場的MARE在4.9%~12.6%,可見兩個模型對氣溫場的模擬效果均好于降水量場。

圖2 MulGETS和k-NN對日降水量均值、標準差和偏度的模擬值與實測值的對比

圖3 MulGETS和k-NN對日最高/最低氣溫均值和標準差的模擬值與實測值的對比

表2 MulGETS和k-NN模型對氣象變量均值、標準差和偏度的平均相對誤差絕對值(MARE)(%)

Table 2 The MARE value of observed and model (MulGETS and k-NN) simulated mean, standard deviation(Std) and skewness of precipitation, maximum and minimum temperature(%)

2.2 兩模型對歷史氣象場極端狀態重現能力的評估

極端天氣狀態模擬是考核天氣發生器的重要指標。兩模型對各站點30a降水量、最高氣溫和最低氣溫的極值模擬與實測對比表明(表3),k-NN模型整體低估了各氣象要素的極值,這是k-NN模型重采樣方法本身的問題。若不討論模擬極值范圍的合理性,MulGETS模型具備一定的極值模擬能力。對于溫度,各站點MulGETS模擬得到的結果均超過歷史觀測數據范圍。但對于降水量,MulGETS模擬的各站點極端降水量可以超過低值站點(如S2、S8),無法超過高值站點(S1、S4),總體上有平均化的趨勢。

表3 MulGETS和k-NN模型對氣溫和降水量極值的模擬值與實測值的對比

Table 3 Comparison of observed value and simulated value by MulGETS and k-NN of maximum precipitation(Pcp), maximum and minimum temperature

注:O代表觀測值,M代表MulGETS模擬值,K代表k-NN模擬值,S1?S12為站點。

Note: O is the observation value, M and K refer to the simulation value of MulGETS and k-NN model, respectively. S1?S12 are stations.

2.3 兩模型對歷史氣象場空間結構特征重現能力的評估

2.3.1 氣象要素的空間相關性

空間相關性是評價多站點天氣發生器的主要指標,氣象要素不同站點相關性模擬結果的好壞直接決定了后續影響評價模型的可信度。從MulGETS模型和k-NN模型空間相關系數的模擬值與實測值的對比可以看到(圖4),雖然站點間降水量的空間相關性系數差異很大,從0.21到0.74不等,但MulGETS和k-NN的模擬結果均較好地保持了原始降水量觀測場的空間相關性系數。同樣地,兩個模型對于溫度場的空間相關性也模擬較好,但k-NN要優于MulGETS。

2.3.2 降水量的空間間歇性

降水量與溫度相比具有較強的空間間歇性,保持原始降水量場的空間間歇性是多站點降水量模擬的重要需求。從降水量場空間連接度指標CR的模擬與實測可見(圖5)。CR模擬值與實測值的相關系數由MulGETS的0.27顯著提高到了k-NN的0.97,可以看出,k-NN相比于MulGETS更好地把握住了降水量場的空間間歇性。盡管在空間相關性系數指標上表現較好,但MulGETS模擬的空間連接度指標在低值區有所高估,在高值區又有所低估,總體上有平均化空間間歇性的趨勢??赡茉蚴窃跀M合降水發生指數KS與均值m和標準差s的關系式時,仍存在一定的誤差。圖6以南岳站為例,其夏季的標準差和秋季均值仍存在一些偏離較大的散點,這些偏差會導致模擬降水量場的空間連接度指標偏離實測值。

圖4 MulGETS和k-NN模型對氣象要素的空間相關性系數模擬值與實測值的對比

圖5 MulGETS和k-NN模型對日降水量的空間連接度指標模擬值與實測值的對比

2.4 兩模型對歷史氣象場時間結構特征重現能力的評估

自相關性是降水時間序列的基本屬性,但在以往天氣發生器評價研究中考察較少。為了衡量天氣發生器對降水場時間結構特征的重現能力,以永州站(S8)第10次(共20次)降水量模擬為例,評估MulGETS模型和k-NN模型重現不同尺度(日、月、季)降水量自相關性的能力。如圖7所示,在0.05的顯著性水平下,觀測的日尺度、月尺度和季尺度降水量,均存在不同的落后時間自相關系數超出的臨界值情況,日值的最大自相關系數(AC)為0.20,出現在滯后1d時(圖7a);月值的最大自相關系數(AC)為0.33,出現在滯后60個月(5a)時(圖7d)。季值的最大自相關系數(AC)為0.47,出現在滯后20季(5a)時(圖7g),由此可見,永州站降水量序列存在一個顯著的5a周期。以永州站自相關系數模擬結果為例,MulGETS在日、月和季尺度的模擬值與實測值的相關系數分別為0.30、0.82和0.88,k-NN在日、月和季尺度分別為0.34、0.81和0.83。總體來看,MulGETS和k-NN在月尺度和季尺度的模擬效果較好,在日尺度自相關模擬上還有不足,特別是低估了降水量在滯后1d時的最強自相關性。此外,從圖7b和圖7c可以發現,觀測中部分顯著的自相關系數被低估,在MulGETS和k-NN模擬中都接近于0。

圖6 MulGETS模型在南岳站對四季降水發生指數與降水量均值(1)和標準差(2)的擬合

圖7 MulGETS和k-NN對永州站日尺度(a)、月尺度(b)和季尺度(c)降水量自相關系數的模擬值與實測值的對比

3 結論與討論

MulGETS和k-NN兩個模型對氣溫場的模擬效果均好于降水量場,降水量場隨機模擬仍然是現有天氣發生器需要克服的難點。MulGETS和k-NN模型均較好地再現了原氣象場的均值、標準差和偏度,且k-NN表現稍好于MulGETS。MulGETS在重現降水量的空間間歇性方面存在不足,但具備一定的極值模擬的能力。k-NN在保持氣象要素空間相關性上也具有優勢,特別是降水量場的空間間歇性。但由于k-NN重采樣算法上的限制,k-NN無法模擬出歷史數據范圍外的極值。此外,MulGETS和k-NN模型在重現原始氣象場自相關性上均存在不足,特別是日尺度降水量,今后應加強這方面研究。

MulGETS和k-NN模型不同的算法決定了模擬得到的氣象序列的不同統計特性,兩個模型各有優勢和不足,MulGETS更適合用于極端氣象事件模擬,而K-NN可以更好地體現原始氣象場的空間差異,應根據不同的研究目的選擇合適的模型。

[1] Wilks D S,Wilby R L.The weather generation game:a review of stochastic weather models[J].Progress in Physical Geography, 1999,23(3):329-357.

[2] Richardson C W.Stochastic simulation of daily precipitation, temperature,and solar radiation[J].Water Resources Research, 1981,17(1):182-190.

[3] Semenov M,Brooks R,Barrow E,et al.Comparison of the WGEN and LARS-WG stochastic weather generators for diverse climates[J].Climate Research,1998,10(2):95-107.

[4] Nicks A D,Lane L J,Gander G A.In USDA-water erosion prediction project:hillslope profile and watershed model documentation[R].USDA-ARS National Soil Erosion Research Laboratory,West Lafayette,USA,1995:2.1-2.22.

[5] 王幼奇,樊軍,邵明安.LARS-WG天氣發生器在黃土高原的適應性研究[J].中國水土保持科學,2007,5(3):24-27.

Wang Y Q,Fan J,Shao M A.Adaptability of climate generator of LARS-WG on the Loess Plateau[J].Science of Soil & Water Conservation,2007,5(3):24-27.(in Chinese)

[6] 李明財,郭軍.CLIGEN在華北地區生成降水的效能評估[J].中國農業氣象,2010,31(2):183-187.

Li M C,Guo J.Efficiency evaluation of precipitation generated by CLIGEN in North China[J].Chinese Journal of Agrometeorology, 2010,31(2):183-187.(in Chinese)

[7] 廖要明,劉綠柳,陳德亮,等.中國天氣發生器模擬非降水變量的效果評估[J].氣象學報,2011,69(2):310-319.

Liao Y M,Liu L L,Chen D L,et al.An evaluation of the BCC/RCG-WG’s performance in simulating daily non- precipitation variables in China[J].Acta Meteorologica Sinica, 2011,69(2):310-319.(in Chinese)

[8] Vallam P,Qin X S.Multi-site rainfall simulation at tropical regions:a comparison of three types of generators[J]. Meteorological Applications,2016,23(3):425-437.

[9] Steinschneider S,Brown C.A semiparametric multivariate, multisite weather generator with low-frequency variability for use in climate risk assessments[J].Water Resources Research, 2013,49(11):7205-7220.

[10] Chen J,Brissette F P,Zhang X J.A multi-site stochastic weather generator for daily precipitation and temperature[J]. Transactions of the ASABE,2014,57(5):1375-1391.

[11] Wilks D S.Multisite generalization of a daily stochastic precipitation generation model[J].Journal of Hydrology,1998, 210(1):178-191.

[12] Lall U,Sharma A.A nearest neighbor bootstrap for resampling hydrologic time series[J].Water Resources Research,1996, 32(3):679-693.

[13] Rajagopalan B,Lall U,Tarboton D G,et al.Multivariate nonparametric resampling scheme for generation of daily weather variables[J].Stochastic Hydrology and Hydraulics, 1997,11(1):65-93.

[14] Buishand T A,Brandsma T.Multisite simulation of daily precipitation and temperature in the Rhine basin by neares-neighbor resampling[J].Water Resources Research, 2001,37(11):2761-2776.

[15] Yates D,Gangopadhyay S,Rajagopalan B,et al.A technique for generating regional climate scenarios using a nearest- neighbor algorithm[J].Water Resources Research,2003,39(7): 1199-1213.

Reproducibility Evaluation of Multi-Site Stochastic Weather Generators: a Comparison between a Typical Parametric Model and a Non-Parametric Model

ZHOU Ling-feng1, MENG Yao-bin1, LU Chao1, WU Gan-lin1, ZHANG Dong-ni1, SONG Hao-zheng1, WU Dan2

(1. Key Laboratory of Environmental Change and Natural Disaster of Ministry of Education/Academy of Disaster Reduction and Emergency Management/Ministry of Emergency Management & Ministry of Education/Faculty of Geographical Science, Beijing Normal University, Beijing 100875, China; 2. Hunan Liuyang Hydrology and Water Resources Bureau, Liuyang 410300)

Many impact models (e.g., hydrological and agricultural models) require simulations of weather variables reflecting the spatial and temporal dependence of observed meteorological fields. New techniques are recently available to generate weather variables simultaneously at multiple locations. This paper presents a comparison of two types of multi-site stochastic weather generators (MulGETS model and k-NN model) for simulation of precipitation and temperature at a network of 12 stations in Xiang River Basin, China. These two models were evaluated for their ability to reproduce the statistical features of the historical meteorological field. The results showed that both MulGETS and k-NN model were successful in reproducing the mean, standard deviation, and skewness of the weather variables, while the performance of k-NN was generally superior to that of MulGETS. The k-NN model was found to perform satisfactorily in preserving the spatial structure of the weather variable, especially the spatial intermittence. Only MulGETS model could generate extreme values out of the historical range. New technology is needed because both MulGETS and k-NN model have the limitation in representing temporal dependence of weather sequence, especially the autocorrelation of daily precipitation.

Multi-site;Weather generator;Model comparison;Xiang River Basin

10.3969/j.issn.1000-6362.2019.06.001

收稿日期:2018?12?11

通訊作者。E-mail:yaobin-meng@bnu.edu.cn

北京師范大學地表過程與資源生態國家重點實驗室方向性項目(2017-FX-07)

周凌峰(1992?),博士生,主要從事氣象水文學研究。E-mail:zhoulf@mail.bnu.edu.cn

周凌峰,孟耀斌,逯超,等.天氣發生器MulGETS和k-NN對區域歷史氣象場特征重現能力的比較[J].中國農業氣象,2019,40(6):341-349

猜你喜歡
模型
一半模型
一種去中心化的域名服務本地化模型
適用于BDS-3 PPP的隨機模型
提煉模型 突破難點
函數模型及應用
p150Glued在帕金森病模型中的表達及分布
函數模型及應用
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
3D打印中的模型分割與打包
主站蜘蛛池模板: 99r在线精品视频在线播放| 国产成人91精品免费网址在线 | 夜夜操国产| 一级爱做片免费观看久久| 国产精品自在线拍国产电影| 免费jjzz在在线播放国产| 四虎国产在线观看| 国产成人综合久久| 亚洲欧洲自拍拍偷午夜色无码| 亚洲码在线中文在线观看| 欧美在线网| 国内熟女少妇一线天| 欧美亚洲激情| 天天色天天综合网| 国产精品偷伦视频免费观看国产 | 国产中文在线亚洲精品官网| 色国产视频| 国产三级视频网站| 日韩a级毛片| 在线观看亚洲精品福利片| 狠狠干综合| 91色爱欧美精品www| 亚洲欧美一区二区三区麻豆| vvvv98国产成人综合青青| 国产美女视频黄a视频全免费网站| 日韩天堂网| 日本精品一在线观看视频| 成人在线综合| 最近最新中文字幕在线第一页| 一级爆乳无码av| 欧美精品伊人久久| 亚洲妓女综合网995久久| 乱码国产乱码精品精在线播放| 国产成人精品男人的天堂| 中文字幕av无码不卡免费| 日本国产在线| 欧美精品一二三区| 国产精品私拍99pans大尺度| 国产精品亚欧美一区二区| 欧美区一区| 99久久国产精品无码| 亚洲欧美另类中文字幕| 2021最新国产精品网站| 波多野结衣久久高清免费| 久久精品丝袜| 91色国产在线| 日韩美毛片| 国产精品美女免费视频大全| 日韩精品一区二区三区大桥未久| 综合色区亚洲熟妇在线| 欧美人人干| 國產尤物AV尤物在線觀看| 亚洲va欧美ⅴa国产va影院| 草草线在成年免费视频2| 午夜欧美理论2019理论| 欧美日韩精品在线播放| 久久久久久久久久国产精品| 综合色天天| 18黑白丝水手服自慰喷水网站| 国产又色又刺激高潮免费看| 18禁色诱爆乳网站| 欧美在线视频a| 久久人妻系列无码一区| 亚洲国产天堂久久综合| 欧美综合区自拍亚洲综合天堂| 国产99精品久久| 国产精品福利尤物youwu | 国产91高跟丝袜| 亚洲永久精品ww47国产| 久久96热在精品国产高清| 婷婷色狠狠干| a级高清毛片| 日本欧美成人免费| 99re热精品视频中文字幕不卡| 国产精品网拍在线| 欧美日韩高清| 国产白浆一区二区三区视频在线| 日韩在线视频网| 中国国产一级毛片| 久久免费看片| 色综合狠狠操| 亚洲啪啪网|