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

小清河流域汛期多年降水變化趨勢

2018-02-01 15:15:10宋蘇林高曉曦左德鵬徐宗學韓先明蔡思揚
南水北調與水利科技 2018年6期

宋蘇林 高曉曦 左德鵬 徐宗學 韓先明 蔡思揚

摘要:以黃臺橋水文站為出口斷面的小清河流域作為研究區,根據流域內5個雨量站1977-2014年日降水資料,首先采用非參數Mann-Kendall法對流域內各站多年汛期降水變化進行趨勢分析,并用Morlet小波分析流域汛期降水的周期變化;其次使用Mann-Kendall檢驗法并結合滑動t檢驗、有序聚類法及Yamamoto法進行突變檢驗;最后使用Hurst指數法對流域各站點降水未來趨勢進行預測。研究結果表明:流域內各站點汛期多年降水變化呈現增加趨勢,但變化趨勢并不顯著;流域汛期降水變化存在22 a左右的主周期;突變分析表明汛期各站點降水的突變年份并不完全相同,而預計汛期降水量的未來變化將呈現出微弱的上升趨勢。

關鍵詞:小清河流域;汛期降水;趨勢變化;非參數檢驗;Hurst指數;Morlet小波

中圖分類號:P426文獻標志碼:A開放科學(資源服務)標識碼(OSID):宋蘇林Trend of precipitation variation in flood season in Xiaoqing River basin

SONG Sulin.1,GAO Xiaoxi.2,ZUO Depeng.2,XU Zongxue.2,HAN Xianming.2,CAI Siyang.2

(1.Ji′nan Hydrology Bureau,Ji′nan 250014,China;2.College of Water Sciences,Beijing Normal University,

Beijing Key Laboratory of Urban Hydrological Cycle and Sponge City Technology,Beijing 100875,China)

Abstract:We took the Xiaoqing River basin as the study area.Based on the daily precipitation data from 1977 to 2014 at five rainfall stations,we used the Mann-Kendall trend test to analyze the variation trend of flood-season precipitation at each station;then we used Morlet wavelet to analyze the periodical variation of flood-season precipitation.We detected the abrupt change points of the precipitation series using various methods,including the Mann-Kendall method,sliding t test,sequential clustering analysis,and Yamamoto method.Finally,we used the Hurst exponent method to predict the future trend of the flood-season precipitation at each station.The results showed that the flood-season precipitation tended to increase,but the abrupt change points were insignificant.The period of the flood-season precipitation was about 22 years.The abrupt change analysis showed that the abrupt change points of the flood-season precipitation detected by different methods were not consistent,and the future variation of the estimated flood-season precipitation showed a weak upward trend.

Key words:Xiaoqing River basin;flood-season precipitation;change trend;non-parameter test;Hurst index;Morlet wavelet

小清河流域是濟南市重要經濟區,在全市的經濟發展中占有舉足輕重的地位,同時流域內的小清河是濟南及沿岸地區唯一的排洪河道,擔負著濟南市、章丘、歷城等地的泄洪任務。近幾十年來,隨著全球氣候的變暖以及人類城市化進程的不斷加快,流域水循環發生了重要變化。同時,由于該流域降水年內分配極不均勻,且主要集中在汛期,汛期降水約占到全年降水的75.5%。因此分析該流域汛期降水的時

空變化規律及未來變化趨勢,對于探討城市化對水文氣象變化的影響具有重要意義。

近些年來,國內諸多學者針對小清河流域水文氣象變化進行了一系列研究,例如,竇實等[1]分析了小清河流域降水、徑流的時空變化,重點對濟南“070718”次特大暴雨洪水資料進行了分析,結果表明該暴雨歷時短、強度高,但空間分布相對均勻。黃國如等[2]根據流域內5個雨量站1977-2013年逐日降水資料分析了城市化對汛期降水的影響,結果表明汛期降水量和頻次均呈增加趨勢,且城區較郊區更為明顯。于翠松等[3]分析了小清河流域防洪現狀及問題,提出了加強河道、病險水庫治理,建立洪水預報系統等建議。李祥松等[4]根據小清河流域內7個雨量站1977-2013的逐日降水資料,主要分析多年汛期降水量趨勢變化和突變成分并對未來變化趨勢進行了預測,結果表明汛期降水呈增長趨勢,且未來將會繼續增加。

本文在對小清河流域的汛期時空演變分析的同時,結合多種突變方法對小清河汛期降水突變點進行識別,并采用Morlet小波分析不同尺度下多年汛期降水的周期分量,旨在為當地的水資源管理及應對汛期洪水危機提供指導。

1研究區概況

濟南市小清河流域,南依泰沂山脈,北界黃河,位于魯中山區與華北平原的過渡地帶,屬暖溫帶半濕潤大陸性季風氣候區。小清河濟南段全長70.5 km,自西向東流經槐蔭、天橋、歷城、章丘四區(市),流域面積為2 792 km.2。流域內地形地勢復雜,地形南高北低,自西向東北傾斜,由南到北依次為山區丘陵、平原、洼地。流域內大小支流有20條左右,大多分布在河的右岸河流,且都是山區性河流坡陡流急。流域多年平均降水量617.2 mm,降水年際變化大,時空分配不均,汛期降水量要占到全年降水量的絕大多數。濟南小清河流域位置見圖1。

2研究方法和數據來源

2.1研究方法

2.1.1趨勢分析

Mann-Kendall法是一種基于秩的非參數統計[CM(22]檢驗法,常用來檢測水質、氣溫、降水等水文氣象序[CM)]

列資料的趨勢變化[5]。其優點是樣本可以不遵循某種特定的分布,且不受異常值的干

擾,計算相對簡便[6]。Mann-Kendal法中有兩個參數比較特殊,一個是代表序列趨勢顯著性水平的Z值,另一個是表示序列變化趨勢傾斜度的β值。其中利用Z判斷序列的顯著性水平,當Z>0時,表明序列呈上升趨勢;當Z<0時,序列呈遞減趨勢;而變化趨勢的大小可根據β值確定,當β>0時,表示序列變化呈遞增趨勢;β<0時則相反,且β絕對值越大,表示變化程度越大[7]。

2.1.2周期分析

小波分析具有多分辨率的功能,可以揭示水文序列的多時間尺度特性。小波分析可把時間序列同時在時域和頻域上展開,對時頻進行更精確的分析,其在時域和頻域上同時具有的局部化功能,能更清楚的看出時間序列的周期變化情況[8]。本文采用Morlet小波分析小清河流域降水時間序列的多時間尺度特征及其周期。

2.1.3突變檢驗

由于突變檢測的方法還不是很成熟,各種檢測方法有其固有的優缺點,很難只根據一種方法就能對一時間序列做出準確的突變分析。若方法不當,很可能會得出錯誤的結論。因此,在確定某氣候系統或過程發生突變的現象時,最好使用多種方法進行比較分析[9]。本研究使用Mann-Kendall法,并結合滑動t檢驗、有序聚類法法及Yamamoto法進行綜合對比和驗證。

(1)Mann-Kendall檢驗法。 Mann-Kendall突變檢驗為非參數檢驗法,能夠客觀地反映樣本時間序列的變化趨勢。進行突變分析時,給定一定的顯著性水平并確定臨界值,計算并繪制UF和UB曲線。若UF或UB的值大于0,表明序列呈上升趨勢,反之,則相反。如果UF和UB在置信區間內有交點,可認為該點為突變點[10]。

(2)滑動t檢驗。 對于時間序列x1,x2,x3,…,xn,設定某一時刻為基準點,計算基準點前后兩段子序列的平均值和方差,并計算得到統計量t。對于給定的顯著性水平α,根據t分布表得到tα。當|ti|≥tα時,說明存在顯著性差異,否則認為兩子序列均值無顯著差異。由于子序列長度選擇帶有人為性,可能會造成突變點的飄移,因此在具體使用時,應該反復變動子序列長度進行試驗比較,提高計算結果的可靠性[9]。

(3)有序聚類法。 對存在變化趨勢的水文序列利用有序聚類法進行跳躍性檢驗,判斷序列中的可能干擾點,該點使同類之間的離差平方和最小,類與類之間的離差平方和最大。對存在多個階段變化過程的時間序列來說,會形成多個谷底,可根據谷底發生時間確定突變點[11]。

(4)Yamamoto法。 Yamamoto法與滑動t檢驗原理類似,也是通過檢驗兩子序列均值差異來判別突變的,但是比滑動t檢驗更簡潔[12]。通過人為設置基準點,計算前后兩段序列的均值和標準差進而確定信噪比,在給定的α下確定有無突變點[9]。

2.1.4未來變化趨勢分析

Hurst指數法是由英國水文專家H.E.Hurst在研究尼羅河徑流資料時發現的,經國內外學者的不斷發展,Hurst指數已經成為對時間序列的未來變化趨勢有很強預測能力的方法。它不需要假定R/S時間序列的分布特征[13]。根據R/S分析計算的Hurst指數大小,可以判斷序列中是否存在趨勢性成分及其趨勢強度。當H=0.5,表明序列是隨機分布的,未來變化趨勢和過去沒有關系;當0≤H<0.5時,表明時間序列具有反持續性,未來變化趨勢將與過去相反;當0.5H≤1時,表明時間序列具有長期持續性,未來變化趨勢將與過去相同[14],且H越大,持續性越強。

2.2數據來源

本文選擇的研究區是以黃臺橋為出口斷面的小清河流域。研究采用的數據是小清河流域及其周邊區域6個雨量站汛期日降水序列資料,由濟南市水文局提供。由于各站點降水序列起始年份不一,序列長度不同,為保證各站序列的一致性,本研究采用1977-2014年作為研究時段。

3結果分析

3.1汛期降水空間分布

為研究小清河流域汛期降水的空間分布,根據流域內6個雨量站1977-2014年的汛期平均降水量,采用克里金法插值得到流域汛期降水空間變化,見圖2。從圖中可以看出,汛期降水從西北部向東部及東南部逐漸遞增。此結果與遲竹萍[15]對山東省夏季降水空間變化研究得出的結果相吻合。其中流域汛期平均降水量最低的是劉家莊站,最高的是興隆站,其降水量分別為417.3 mm,479.0 mm,兩者之間的插值為61.7 mm,空間變化較大。

3.2汛期降水趨勢分析

使用Mann-Kendall法對小清河流域雨量站多年汛期降水變化進行趨勢檢驗,同時對汛期降水的多年平均變化率進行計算,最終統計結果見表1。由表1可知,流域內各站點汛期降水多年變化呈上升趨勢,除黃臺橋站通過了置信度90%的檢驗外,其余站點變化均不顯著。各站點中增長速度最快的是劉家莊站1.60 mm/a,其次是黃臺橋站1.39 mm/a和東紅廟站1.32 mm/a。許多學者[16-19]研究表明城市化會使城市下風區降水增加,而李祥松等[4]對小清河流域在1984年、1992年、2002年和2013年四期土地利用的統計表明,研究區的城市建設用地近些年來急劇增加,[HJ2.2mm]其它用地面積迅速減少。而研究時段正好是在汛期(6月-8月),偏南風頻次較多,而劉家莊與黃臺橋正處于城市的下風口,可能受城市化影響,導致降水量增長率相較其他站點相對要快。

3.3汛期降水突變分析

本文采用Mann-Kendall突變檢驗、滑動t檢驗、有序聚類法和Yamamoto法綜合對比識別小清河流域各雨量站多年汛期降水量的突變點,這里只以東紅廟雨量站為例進行分析,在給定顯著性水平α=0.05下,對東紅廟多年汛期降水采用Mann-Kendall法進行突變點識別,見圖3。由圖3中UF曲線可以看出,汛期降水在1977-1980年和1988-1989年左右呈上升趨勢;在1980-1986年和1989年左右呈下降趨勢。總的來說,在1977-1990年汛期降水變化有增有減,呈波動狀態,而在1990年之后則一直呈上升趨勢,但變化范圍處于±1.96的信度線之間,說明多年變化并不顯著。從UF和UB曲線位于置信區間內的交點來看,突變發生的時刻有1978年、1987年左右,但并不能直接下結論說這幾個點就一定是突變點,也可能存在虛假的突變點,應結合其它檢驗方法綜合分析,盡量避免錯誤產生。

結合滑動t檢驗法繪制t統計量序列圖,見圖4。其中滑動長度取n1=n2=7,給定顯著性水平0.05。從圖4中可以看出,一共有兩個點超過置信水平,其中一處是正值(1996年),一處是負值(2003年),說明在這個階段內,降水出現了兩次明顯的突變。在1996年左右,汛期降水出現由多轉少的突變;在2003年左右汛期降水由少轉多的突變。由于人為選擇滑動長度不同所導致的突變點漂移問題,除了要反復變動滑動長度進行對比,還可采用多種方法對比可提高最終檢測結果的可靠性。故本文同時選擇了有序聚類法對汛期降水進行突變檢驗。

汛期降水序列離差平方和S(τ)曲線見圖5。從序列圖中可以看出,整個序列的最低點出現在1986年,局部的極小點有1978年、1989年、1993年、1998年、2002年、2010年。可以看出,有序聚類檢驗在東紅廟站并未檢驗出有效的突變點。

最后選擇Yamamoto法,取n1=n2=10,滑動計算信噪比SNR的值,但繪制出的曲線變化范圍始終小于1,并沒有檢測出有突變點,見圖6。

對比以上四種突變方法的檢驗結果,說明各方法在檢測突變點時各有其局限性和不足。對各種方法檢測出的突變點進行對比分析,最后確定東紅廟突變年份為1986年、1989年。最終確定的各站點突變年份見表2。

3.4汛期降水周期分析

通過泰森多邊形法計算流域的面雨量,采用Morlet小波并結合Matlab對小清河流域汛期降水量進行多尺度時間分析。根據Morlet小波變換得到的小波變換實部等值線圖可以清晰的顯示汛期降水的時間尺度變化及位相結構。汛期降水小波變換實部等值線圖見圖7。從中可以看出,在1977-2014年時間序列中,存在兩個比較明顯的震蕩周期。其中3~8 a的尺度表現的較強,其中心尺度在5 a左右,且呈現正負位相交替出現的現象。而11~27 a的尺度表現也十分明顯,其中心尺度分別是15 a和22 a左右。從時間尺度為11~27 a,中心尺度為15 a這一時間來看,1977-1979年,1986-1989年、1995-1998年、2003-2006年各時段降水偏多,其余時間段則偏少。從2014年以后未閉合的等值線圖來看,未來幾年汛期降水將很可能繼續增加,同時可能存在更長時間尺度的時間周期,需要增加時間序列的長度以期進一步驗證。

小波方差能反映信號波動能量隨尺度的分布,因此可以確定信號中不同尺度擾動的相對強度和存在的主要時間尺度[20]。小清河流域汛期降水小波方差見圖8。從圖中可以看出存在三個明顯的波峰,對應的時間尺度分別是5 a、15 a和22 a,其中最高點出現在22 a左右,說明汛期降水在22 a左右的震蕩周期最大,可認為是小清河流域多年汛期降水的主周期,第二周期的峰值僅次于第一周期,為15 a,第三周期為5 a。從未閉合的小波變換實部等值線圖來看,以22 a為周期,未來幾年降水會減少;以15 a周期,未來幾年降水將繼續增加;以5 a為周期,未來將呈現波動變化,由于第一周期與第二周期峰值相差不大,并不能準確判斷未來的降水趨勢變化,需進一步分析。

3.5汛期降水未來變化趨勢預測

使用R/S法繪制的小清河流域各雨量站點多年汛期降水的Hurst指數見圖9,具體統計見表3。從中可以看出各站點Hurst指數均大于0.5,表明汛期降水的未來趨勢與過去相同,即將繼續呈現增長趨勢,這與以15 a為周期時,小波分析得出的結論一致。其中Hurst指數最大的是興隆,其次是劉家莊、邵而。從中發現未來增長趨勢最顯著的三個站點也都是分布在城鎮的周圍,分析其原因可能與城市化對降水的影響有關。城市化水平的不斷提高使得降水不斷增加,但降水增加的趨勢最終會受到城市化自身發展的影響,不可能造成大范圍降水的增加,所以各站點未來汛期降水增長趨勢要相對較弱[21]。

4結論

(1)小清河流域汛期降水量由西北部向東部及東南部逐漸遞減。其中降水量最高、最低的站點分別為興隆站和劉家莊站,變化差值為617 mm,空間變化差異較大。

(2)流域各站點多年汛期降水變化趨勢基本相同,各站點突變年份也基本一致,都在1989年左右檢測出了突變點。

(3)流域汛期降水時間序列存在多個周期,時間尺度不一,其中第一主周期為22 a,第二、第三周期分別為15 a和5 a。

(4)未閉合的小波等值線圖與各站點Hurst指數綜合分析表明流域未來汛期降水繼續呈增長的趨勢。

參考文獻(References):

[1]竇實,曹升樂.濟南市小清河流域降雨徑流及暴雨洪水特性分析[C].// 第五屆中國水論壇文集,2007:710-713.( DOU S,CAO S L.Characteristics of storm flood in Xiaoqing River Basin of Jinan City[C].// Water Forum of China,2007:710-713.(in Chinese))

[2]黃國如,何泓杰.城市化對濟南市汛期降雨特征的影響[J].自然災害學報,2011,20(3):7-12.(HUANG G R,HE H J.Impact of urbanization on features of rainfall during flood period in Jinan City[J].Journal of Natural Disasters,2011,20(3):7-12.(in Chinese))

[3]于翠松,王艷玲.小清河流域的防洪問題及對策[J].水利發展研究,2002,2(5):39-41.(YU C S,WANG Y L.Flood protection and countermeasures in Xiaoqing River Basin[J].Water Resources Development Research,2002,2(5):39-41.(in Chinese))

[4]李祥松,于翠松,曹升樂,等.濟南市小清河流域汛期降水時空演變規律分析[J].水資源與水工程學報,2016,27(6):72-78.(LI Y S,YU C S,CAO S L.The spatial-temporal evolution law of flood season precipitation in Xiaoqing River Basin of Jinan[J].Journal of Water Resources and Water Engineering,2016,27(6):72-78.(in Chinese)) DOI:10.11705/j.issn.1672-643X.2016.06.13.

[5]徐宗學.水文模型[M].北京:科學出版社,2009.(XU Z X.Hydrological Models[M].Beijing:Science Press,2009.(in Chinese))

[6]徐宗學,隋彩虹.黃河流域平均氣溫變化趨勢分析[J].氣象,2005,31(11):7-10.(XU Z X,SUI C H.Long-term trend of temperature in the Yellow River Basin[J].Meteorological Monthly,2005,31(11):7-11.(in Chinese)) DOI:10.3969/j.issn.1000-0526.2005.11.002.

[7]傅麗昕.近57年來和豐縣氣溫和降水量的趨勢性及突變特征[J].南水北調與水利科技,2014,12(4):38-41.(FU L X.Tendency and mutation analysis of annual temperature and precipitation of Hefeng County in recent 57 years[J].South-to-North Water Transfers and Water Science & Technology,12(4):38-41.(in Chinese)) DOI:10.13476/j.cnki.nsbdqk.2014.04.009.

[8]王文圣,丁晶,李躍清.北京:水文小波分析[M].化學工業出版社,2005.(Wang W S,Ding J,Li Y Q.Hydrology wavelet analysis.Beijing:Chemical Industry Press,2005.(in Chinese))

[9]魏鳳英.現代氣候統計診斷與預測技術[M].北京:氣象出版社,1999.(WEI F Y.Beijing:The technology of statistical diagnosis and prediction for modern climate[M].Beijing:China Meteorological Press,1999.(in Chinese))

[10][ZK(#]蔡霞,蔡琳,李春華,等.晉北地區降水量時空變化及突變分析[J].干旱地區農業研究,2012,30(2):247-254.(CAI X,CAI L,LI C H,et al.Analysis of characteristics of precipitation variation and mutation in Northern Shanxi[J].Agricultural Research in the Arid Areas,2012,30(2):247-254.(in Chinese)) DOI:10.3969/j.issn.1000-7601.2012.02.042.

[11]張學真,劉燕.灞河出山徑流序列變化的小波分析[J].水資源保護,2006,22(3):12-15.( Zhang X Z,Liu Y.Wavelet analysis on runoff sequences from mountainous watershed of Bahe River[J].Water Resources Protection,2006,22(3):12-15.(in Chinese)) DOI:10.3969/j.issn.1004-6933.2006.03.004.

[12]李珍,姜逢清.1961-2004年新疆氣候突變分析[J].冰川凍土,2007,29(3):351-359.( ZHEN L I,JIANG F Q.A study of abrupt climate change in Xinjiang Region during 1961—2004[J].Journal of Glaciology & Geocryology,2007,29(3):351-359.(in Chinese)) DOI:10.3969/j.issn.1000-0240.2007.03.003.

[13]燕愛玲,黃強,劉招,等.R/S法的徑流時序復雜特性研究[J].應用科學學報,2007,25(2):214-217.(YAN A L,HUANG Q,LIU Z,et al.Complicated property of runoff time series studied with R/S method[J].Journal of Applied Sciences,2007,25(2):214-217.(in Chinese)) DOI:10.3969/j.issn.0255-8297.2007.02.021.

[14]馮新靈,羅隆誠,邱麗麗,等.青藏高原至中國東部年雨日變化趨勢的分形研究[J].地理研究,2007,26(4):835-843.(FENG X L,LUO L C,QIU L L,et al.Fractal research of rainy day changing trend from Tibetan Plateau to Eastern China[J].Geographical Research,2007,26(4):835-843.(in Chinese)) DOI:10.3321/j.issn:1000-0585.2007.04.021.

[15]遲竹萍.近45年山東夏季降水時空分布及變化趨勢分析[J].高原氣象,2009,28(1):220-226.(CHI Z P.Spatial and temporal distributions and climatic change of summer precipitation in Shandong[J].Plateau Meteorology,2009,28(1):220-226.(in Chinese))

[16]HUFF F A,CHANGNON S A J.Climatological assessment of urban effects on precipitation:Final report part I[J].Journal of Applied Meteorology,1972,11(5).

[17]CHANGNON S A.Rainfall Changes in Summer Caused by St.Louis[J].Science,1979,205(4404):402.

[18]SHEPHERD J M,PIERCE H,NEGRI A J.Rainfall modification by major urban areas:Observations from spaceborne rain radar on the TRMM satellite.[J].Journal of Applied Meteorology,2007,41(7):689-701.DOI:10.1175/1520-0450(2002)041<0689:RMBMUA>2.0.CO;2.

[19]北京市氣象局氣候資料室.北京城市氣候[M].氣象出版社,1992:50-57.(Beijing Meteorological Bureau climate data room.Beijing City Climate[M].China Meteorological Press,1992:50-57.(in Chinese))

[20]常浩娟,劉衛國,吳瓊.60年瑪納斯河紅山嘴徑流規律特征分析[J].水土保持研究,2016,23(6):128-134.(CHANG H,LIU W,QIONG W U,et al.Runoff characteristics of Hongshanzui Hydrologic Station of Manas River in the past 60 years[J].Research of Soil & Water Conservation,2016,23(6):128-134.(in Chinese))

[21]周翠寧,任樹梅,楊培嶺,等.城市化對降雨特征影響研究[J].水利水電技術,2007,38(10):62-65.(ZHOU C N,REN S M,YANG P L,et al.Study on impact from urbanization on rainfall characteristics[J].Water Resources & Hydropower Engineering,2007,38(10):62-65.(in Chinese)) DOI:10.3969/j.issn.1000-0860.2007.10.019.

主站蜘蛛池模板: 亚洲AV无码久久精品色欲| 色首页AV在线| 欧美色图第一页| 日韩欧美国产综合| 久热99这里只有精品视频6| 国产精品香蕉| 中文字幕在线免费看| 亚洲天堂首页| 首页亚洲国产丝袜长腿综合| 9999在线视频| 亚洲一区二区日韩欧美gif| 免费一级全黄少妇性色生活片| 国产xx在线观看| 国产精品嫩草影院视频| 亚洲男人在线| 爽爽影院十八禁在线观看| 久久国产精品麻豆系列| 欧美国产菊爆免费观看| 91青青草视频| 日本久久免费| 国产女人在线观看| 成人毛片在线播放| 国产精品亚洲精品爽爽| 欧美一级专区免费大片| 国产9191精品免费观看| a级毛片毛片免费观看久潮| 五月婷婷综合网| 国产精品开放后亚洲| 伦精品一区二区三区视频| 国产成人欧美| 最新国产网站| 欧美亚洲一二三区| 国产精品对白刺激| 五月婷婷欧美| 在线精品亚洲一区二区古装| 波多野结衣久久高清免费| 看你懂的巨臀中文字幕一区二区| 精品国产成人三级在线观看| 性欧美久久| 国产福利2021最新在线观看| 亚洲性影院| 国产国产人免费视频成18| 亚洲国产精品成人久久综合影院| 国产成年女人特黄特色大片免费| 欧美亚洲国产一区| 69免费在线视频| 欧美亚洲国产一区| 青青久久91| 99精品视频九九精品| 国产福利在线免费观看| 一本无码在线观看| 国产幂在线无码精品| 在线观看亚洲精品福利片| 91福利国产成人精品导航| 亚州AV秘 一区二区三区| 国产成人精品亚洲日本对白优播| 亚洲性日韩精品一区二区| 亚洲成人动漫在线观看| 国产亚洲第一页| 久久99蜜桃精品久久久久小说| 制服无码网站| 在线看片国产| 在线免费亚洲无码视频| 91午夜福利在线观看| 干中文字幕| 色综合a怡红院怡红院首页| 综合色在线| 中文字幕在线播放不卡| 人妻夜夜爽天天爽| 日本在线国产| 99在线视频精品| 五月天香蕉视频国产亚| 成人夜夜嗨| 国内熟女少妇一线天| 又爽又大又黄a级毛片在线视频 | 国产亚洲高清视频| 欧亚日韩Av| 国产91精选在线观看| 亚洲专区一区二区在线观看| 欧美精品亚洲精品日韩专区va| 欧美性久久久久| 欧美日本二区|