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

土地利用變化對雅礱江流域水文極值的影響研究

2020-06-11 12:43:56董立俊董曉華薄會娟
中國農村水利水電 2020年4期
關鍵詞:城鎮化模型

董立俊,董曉華,魏 沖,喻 丹,薄會娟,郭 靖

(1.三峽大學 水利與環境學院,湖北 宜昌 443002;2.水資源安全保障湖北省協同創新中心,武漢 430070;3.中國電建集團華東勘測設計研究院有限公司,杭州 310014)

全球環境變化是當前和未來人類社會發展面臨的共同問題, 全球環境變化包括氣候變化和人類活動引起的環境變化[1]。土地利用變化直接體現和反映了人類活動的影響程度。土地利用變化可以通過改變植物截留、蒸散發、下滲、地下水補給、基流和地表徑流等水文要素的方式直接影響流域水循環,它對水循環過程的影響結果將直接導致水資源供需關系發生變化,進而對流域水資源開發利用、生態環境、供水保障、防洪安全等方面產生重大影響[2]。

針對土地利用變化的水文響應,國內外學者進行了大量類似研究。歐春平等[1]使用 SWAT模型,定量分析了土地利用/覆被變化對蒸發、徑流等水循環要素的影響。張蘭影等[3]利用率定好的SWAT模型,結合情景分析方法,定量分析了古浪河流域氣候變化和土地利用/覆被變化對蒸散發、地表徑流等水循環要素的影響。劉芳等[4]基于GIS和Fragstats等軟件,分析了涇河流域土地利用類型和景觀格局變化及其對徑流的影響。Li等[5]、張凌[6]、竇小東[7]等分別研究了土地利用變化及氣候變化對藉河流域、黑河流域中上游、龍川江流域徑流的影響。羅開盛等[2]借助SWAT模型,并結合統計方法定量評估了黑河徑流對土地利用變化和氣候變化的敏感性。楊宏偉等[8]利用 SWAT 模型研究了東江流域典型子流域土地利用/覆被變化對地表徑流的影響。Zhou等[9]通過耦合CLUE-S 和SWAT模型,定量分析了土地利用/覆被變化對西苕溪流域徑流的影響。高玉琴等[10]通過建立研究區HEC-HMS水文模型,并設置3種土地利用情景(自然發展情景、林地限制情景、水田限制情景),研究了秦淮河流域高速城市化背景下土地利用變化的暴雨洪水響應機制。

土地利用/土地覆蓋變化(LUCC)受自然資源條件和社會經濟發展程度等因素約束,同時也直接影響著區域內各生態過程的變化[11]。在氣候變化、城鎮化進程加快以及流域生態環境問題日益凸顯等多重背景下,土地利用/覆被變化(LUCC)研究已成為全球環境變化研究領域的熱點[12]。目前土地利用/土地覆蓋變化模擬常用的模型主要包括回歸分析模型、Markov模型、CA(cellular automaton)模型、CLUE-S (the conversion of land use and its effects at small regional extent)模型[13,14]。CLUE-S 模型是荷蘭瓦赫寧根大學 “土地利用變化和影響”研究小組在CLUE 模型的基礎上開發的,CLUE-S模型通過土地利用類型與其自然條件及社會經濟等驅動因子之間的定量關系,來模擬土地利用變化的空間分布,進而對未來土地利用變化進行預測[15-17]。

目前雅礱江流域已開展的類似研究主要是針對歷史期徑流模擬及變化分析,陳昕等[18]運用北京師范大學水科學研究院自主開發的基于人工神經網絡方法的多模型輸出集合系統,對SWAT模型、BTOPMC模型、VIC模型和DTVG模型等4個水文模型模擬的雅礱江流域徑流結果進行了集合計算。黃艷艷等[19]采用相關分析和逐步多元回歸方法分析了2000-2014年雅礱江上游徑流變化及其影響因素,結果表明2000-2014年雅礱江上游徑流整體呈上升趨勢,且冬季積雪面積和汛期降水為主要影響因素。雅礱江流域為我國重要的水電基地,其規劃總裝機容量達3 000 萬kW,變化環境導致的徑流變化將直接影響雅礱江梯級水庫群運行安全和發電調度。量化分析土地利用/覆被變化對水文極值的影響對于雅礱江流域梯級水庫群防洪調度和確保流域防洪安全具有重要的理論價值和現實意義。

本文通過建立雅礱江流域日尺度SWAT模型,并同時進行多出口率定,然后采用CLUE-S模型模擬了城鎮化發展情景下的土地利用空間分布,并在現狀土地利用基礎上進行退耕還林操作得到了退耕還林情景下的土地利用空間分布,進而模擬了3種土地利用情景(基準期、城鎮化情景及退耕還林情景)下的日徑流序列,并根據模擬的日徑流序列對水文極值變化情況進行了統計分析。

1 研究區概況及數據

1.1 研究區概況

雅礱江是長江上游金沙江的最大支流,其徑流量占金沙江總徑流量的41%。雅礱江流域水量豐沛、落差大,干支流蘊藏了豐富的水力資源,是我國十三大水電基地之一。雅礱江干流共規劃了21級水電站,規劃總裝機容量約3 000 萬kW,年發電量約1 500 億kWh。雅礱江流域位于青藏高原的東部,地理位置介于東經96°52′~102°48′、北緯26°32′~33°58′之間,流域大致在南北向呈條帶狀,流域總面積12.8 萬km2。雅礱江流域年降水量在河源區為500~600 mm,中下游區漸增為900~1 300 mm,中游地區最高可達1 500~1 800 mm。大致自北向南遞增,且東側多于西側。

1.2 數據來源及處理

數字高程模型DEM(Digital Elevation Model)是SWAT模型對研究區進行河流水系及流域邊界提取、坡度坡向計算、子流域劃分的依據。本文使用USGS(United State Geological Survey)官網下載的SRTM(Shuttle Radar Topography Mission)90 m×90 m數據作為建模的DEM數據,雅礱江流域DEM、水系圖及氣象水文站點位置見圖1。

本文采用的土壤基礎數據來源于聯合國糧農組織(FAO)和維也納國際應用系統研究所(IIASA)構建的世界協調土壤數據庫HWSD(Harmonized World Soil Database)。首先采用劃分的流域邊界對HWSD空間柵格數據進行裁剪,得到研究區土壤空間分布圖,并在此基礎上根據HWSD土壤屬性數據制作SWAT模型運行所需的土壤數據庫。

圖1 雅礱江DEM、水系圖及氣象站點位置Fig.1 Map of DEM, river system and hydrological stations in Yalong River Basin

土地利用數據來源于中國科學院資源環境科學數據中心下載的1∶10萬的1995年土地利用數據集,采用劃分的流域邊界進行裁剪,并根據SWAT模型中土地利用/植被覆蓋分類標準進行重分類,得到雅礱江流域土地利用空間分布圖。

SWAT模型所需的氣候變量包括天氣發生器、日降水量、日最高氣溫、日最低氣溫、太陽輻射、風速和相對濕度。其中,日降水量、日最高氣溫、日最低氣溫為實測數據,采用研究區內清水河、石渠、甘孜、色達、道孚、新龍、康定、木里、九龍、越西、昭覺、鹽源、西昌等13個國家氣象站1963-2004年的逐日數據。研究區氣象站點分布見圖 1。其余無觀測資料的氣象數據使用天氣發生器產生。天氣發生器參數庫根據美國國家環境預報中心(NCEP)基于全球氣候再分析系統(CFRS)提供的全球高精度氣候變量數據計算。

2 流域SWAT模型構建及率定

2.1 SWAT模型基本原理

SWAT(Soil and Water Assessment Tool)模型是美國農業部(USDA)農業研究中心開發的基于物理機制的、能夠連續模擬的流域尺度模型。該模型是基于GIS的半分布式水文模型,可模擬流域循環中多種水文物理化學過程,并可模擬和預測土地利用及氣候因素變化對流域水文循環的影響[20-22]。

2.2 SWAT模型率定和驗證結果

本文采用SWAT-CUP(SWAT Calibration and Uncertainty Programs)程序進行參數敏感性分析和模型率定。SWAT-CUP是瑞士聯邦水科學與技術研究所(Swiss Federal Institute of Aquatic Science and Technology)開發的,以SWAT模型模擬結果為基礎對參數進行敏感性分析、自動校準、不確定分析和驗證的分析工具[23-25]。評價模擬精度使用3個指標:Nash-Sutcliffe系數NS,確定性系數R2,相對誤差PBIAS。

參數敏感性分析基于SWAT-CUP的LH-OAT全局方法進行,參數率定通過SUFI-2優化算法實現。SWAT模型中與徑流模擬有關的參數共有28個,根據LH-OAT方法得出的P-value值和t-Stat值對參數敏感性進行排序,P-value值越小,t-Stat值的絕對值越大,表明參數越敏感;反之,參數越不敏感[26]。根據參數敏感性分析結果,結合雅礱江流域特點(流域北部氣溫較低,常年有積雪,因此加入與融雪模塊相關的參數),以及本文研究目的(重點考慮與水文極值模擬有關的參數),選取敏感性高且對流域水文過程影響較大的20個參數進行率定,參數率定結果見表1。

為驗證SWAT模型在雅礱江流域的空間適用性,進而分析土地利用變化對流域不同區域的水文極值的影響程度,結合徑流資料收集情況,本文選取雅江站、洼里站以及小得石站3個出口率定雅礱江流域SWAT模型,同時模擬雅江站、洼里站以及小得石站1981-1997年的日徑流。其中以1981-1982年為預熱期,1983-1990年為率定期,1991-1997為驗證期。本研究區構建的SWAT模型子流域及水文響應單元劃分參數設置:最小集水面積值閾設置為170 000 ha,設置四級坡度,分別為0%~15%、15%~25%、25%~40%、≥40%,拆分合并到其他類型的閾值分別設置為:坡度閾值10%,土地利用閾值10%,土壤類型閾值15%。研究區域劃分為41個子流域,641個水文相應單元。模擬結果見圖2~圖4,模型精度統計指標見表2。

根據上述SWAT模型模擬結果,模擬日徑流與實測日徑流吻合較好,特別是對于年內最大日流量,基本上3個站點都有較好響應,基本能夠重現最大日流量的出現時間和洪水量級等主要特征,這對水文極值變化的模擬和評價是十分有利的。由表 2可知,模型率定期和驗證期的Nash-Sutcliffe系數NS、確定性系數R2均在0.70以上,除洼里站的驗證期外,PBIAS的絕對值均小于10%,模擬精度較高。

表1 參數率定結果Tab.1 Results of parameter calibration

注:率定值中…/…表示 1~21號子流域率定值/22~41號子流域率定值。

圖2 雅江站日徑流模擬結果Fig.2 Daily runoff simulation at Yajiang station

圖3 洼里站日徑流模擬結果Fig.3 Daily runoff simulation at Wali station

圖4 小得石站日徑流模擬結果Fig.4 Daily runoff simulation at Xiaodeshi station

3 土地利用情景設置

3.1 CLUE-S模型

CLUE-S(Conversion of Land Use and its Effects at Small Region Extent)模型通過量化土地利用與自然條件與社會經濟等驅動因子之間的關系來模擬土地利用空間分布,從而預測未來土地利用。CLUE-S模型主要由非空間需求模塊和空間分配模塊組成,其中非空間模塊是獨立于模型而運行的,通過對人口增長、社會經濟以及區域綜合發展規劃等因素的分析,計算研究區不同土地利用類型每年的數量需求。空間模塊則是以柵格化空間數據為基礎,根據每個柵格單元出現某種土地類型的概率、土地利用轉換規則等,從而對模擬年份的土地利用需求進行空間位置分配[27]。

表2 模型模擬精度Tab.2 Accuracy of model simulation

3.2 驅動因子選擇

本文在考慮資料的可獲取性,因子是否能定量化,在研究區內部存在是否空間差異,與研究區土地利用變化的相關性,自然條件與社會經濟因子綜合選取等原則的基礎上,從地形條件、土壤條件、氣象條件以及社會條件4個方面選擇了9個驅動因子[28](見圖5)。

圖5 雅礱流域CLUE-S模型驅動因子選擇Fig.5 Selection of driving factors for CLUE-S model in Yalong River Basin

CLUE-S模型使用Logistic回歸分析對不同土地利用類型與各種驅動因素之間的關系進行定量分析[29]。本文使用二元Logistic回歸對每一柵格單元中可能出現某種土地利用類型的概率進行分析,二元Logistic回歸方程的基本形式見。雅礱江流域各種土地利用類型的Logistic二元回歸方程系數及常量見表3。

(1)

式中:Pi為每個柵格可能出現某種土地利用類型i的概率;Xn,i為驅動因子;β為回歸系數,其中β0為常量。

表3 Logistic二元回歸方程系數及常量Tab.3 Accuracy of model simulation

由表 3可知,影響雅礱江流域土地利用類型的主要驅動因子有坡度、有機質含量、排水條件、土壤pH、總用水量及人口分布。采用ROC(Relative Operating Characteristics)方法對二元Logistic回歸結果進行檢驗,雅礱江流域主要的6種土地利用類型的ROC值為0.75~0.86,均大于0.7,表明表3中的Logistic回歸方程對雅礱江流域的土地利用具有較強的解釋能力,可以用來模擬雅礱江流域的土地利用。

3.3 CLUE-S模擬及驗證

以雅礱江流域1995年土地利用為基準期,2005年土地利用為目標期,率定CLUE-S模型參數,用Kappa指數驗證模擬效果[30]。Kappa指數主要是對圖像的精確性及一致性進行評價,其計算公式如下。

Kappa=(Po-Pc)/(Pp-Pc)

(2)

式中:Po為正確模擬的比例;Pc為隨機情況下期望的正確模擬比例,研究區內總共土地利用類型有6種,即Pc=1/6;Pp為理想分類情況下正確模擬的比例,取值為1。

Kappa指數的值位于0~1之間,值越接近1,說明模擬結果越好。當Kappa指數值大于0.7時,表明結果較好。

經計算,整個雅礱江流域共有127 885個柵格,正確模擬的柵格為126 788個,即P0=126 788/127 885=0.99,Kappa=0.988,可知該模型模擬的精度很好,可以用于雅礱江流域土地利用狀況的模擬。

3.4 土地利用情景設置及模擬

雅礱江流域整個流域林地和草地合計占總面積的比例超過85%,且流域內生態環境脆弱,一方面開墾適宜耕作的后備土地資源的難度較大,另一方面具有生態保護功能的林地和草地屬于保護用地,面積應保持相對穩定。根據流域內地方政府(主要是四川省甘孜藏族州和涼山彝族州)土地利用規劃,結合城鎮化發展和生態保護的要求,主要設置兩種土地利用情景,一種是城鎮化發展情景,即按照1990-2015年城鎮用地的年平均增長速率線性外延至2050年,未利用地和耕地面積減小,其他土地類型面積不變;另一種是退耕還林情景,即以2015年土地利用現狀為基礎,大于25°的耕地退耕還林,15°~25°的耕地退耕還草。根據以上設置的土地利用發展情景,各土地利用類型的面積和比例見表4。

表4 雅礱江流域土地利用情景Tab.4 Land use scenario in Yalong River Basin

根據表4可知,與基準期相比,城鎮化情景和退耕還林情景土地利用變化總體較為平緩,其中城鎮化情景下城鎮面積占比由0.12%增加到0.37%,耕地和未利用地面積占比略有減少;退耕還林情景下耕地面積占比由5.43%降至2.85%,林地和草地面積占比相應增加。

根據各土地利用需求量的設置,采用已率定的CLUE-S模型模擬在城鎮化發展情景下的雅礱江土地利用空間分布,并對現狀土地利用空間分布進行退耕還林操作,得到退耕還林情景下的雅礱江土地利用空間分布。城鎮化情景及退耕還林情景下的雅礱江流域土地利用空間分布見圖6。

4 徑流模擬及水文極值分析

基于率定好的SWAT模型,分別輸入城鎮化情景、退耕還林情景下的土地利用空間分布圖(圖6),得到城鎮化情景、退耕還林情景下的雅礱江流域SWAT模型,然后保持其他參數不變分別模擬3種情景(基準期土地利用情景、城鎮化土地利用情景以及退耕還林土地利用情景)下雅礱江流域1963-2004年日徑流,其中1963-1964年為預熱期,同時對雅江站、麥地龍站和洼里站進行日徑流模擬。這是因為,雅江站為在建兩河口電站的設計代表站,且位于雅礱江流域上中游交界處;麥地龍為在建楊房溝電站的設計代表站;洼里站為已建成的錦屏一級水庫的入流代表站,且位于雅礱江流域中下游交界處。根據模擬日徑流序列,分別提取雅江站、麥地龍站以及洼里站每年的最大日流量,3個站點的各年度最大日流量見圖7,最大日流量序列統計參數見表5。

圖6 雅礱江流域土地利用對比圖Fig.6 Land use contrast map of Yalong River Basin

根據圖7及表5可知,在城鎮化土地利用情景下,雅礱江流域多年日最大流量整體上增加,雅江站、麥地龍站及洼里站的日最大流量序列均值增幅依次為4.6%、3.4%及2.3%,增幅從上游至下游逐步減小。在退耕還林土地利用情景下,雅礱江流域的日最大流量略有減少,但減少幅度十分有限,洼里站減少幅度最大,也僅為-1.2%。

圖7 模擬最大日流量對比Fig.7 Comparison of simulatedmaximum daily flow

表5 雅礱江流域最大日流量序列統計參數
Tab.5 Land use scenario in Yalong River Basin

項 目雅江站基準期城鎮化情景退耕還林情景麥地龍站基準期城鎮化情景退耕還林情景洼里站基準期城鎮化情景退耕還林情景均值/(m3·s-1)2 7142 8402 6984 1454 2874 1255 0295 1444 966均值變幅/%4.6-0.63.4-0.52.3-1.2離差系數0.23 0.23 0.23 0.19 0.19 0.19 0.18 0.18 0.19

上述結果表明:城鎮化發展會對雅礱江流域的水文極值造成一定影響。退耕還林對雅礱江流域的水文極值的影響則十分微弱,其主要是發揮生態保護作用。在本文設置的城鎮化發展情景下,城鎮面積增加,即地表不透水面積增加,一方面減少了徑流下滲,降雨形成地表快速徑流的比例增加,另一方面縮短了地表徑流匯流時間,兩方面因素影響導致了日最大流量增加。

城鎮化土地利用情景會導致雅礱江流域日最大流量增加,對流域防洪安全不利,但在雅礱江流域,人類活動的另一方面(即水庫調蓄)同樣也對水文極值有著直接影響。雅礱江流域作為我國重要的水電基地,中下游已建成錦屏一級、二灘等大型調節水庫,在建兩河口調節水庫,通過梯級調度,充分發揮上述水庫的調蓄作用,可以削減洪峰,有效控制出庫流量,對最大出庫流量以及洪水過程線都有著較好控制作用,可使洪水過程線有陡峭變為平緩,這對防洪安全是十分有利的。

本文設置的城鎮化情景預計城鎮面積占比為0.37%,小于涼山彝族州規劃的城鎮面積占比1.01%,如后續城鎮化迅速發展,城鎮面積占比明顯增加,其對雅礱江流域水文極值的影響會更加顯著,再加之氣候變化的影響,二者疊加后對雅礱江流域水文極值以及雅礱江流域梯級水庫群防洪調度的影響不容忽視。因此在后續流域綜合規劃修編過程中,應按照集約用地、適當控制城鎮建設規模的原則,科學規劃城鎮建設布局,合理控制城鎮建設用地增長速度;同時應進一步優化調度,充分發揮流域梯級水庫群的調蓄作用,盡可能降低氣候變化及土地利用變化對流域防洪安全造成的不利影響。

5 結 論

本文建立了雅礱江流域日尺度SWAT模型,并對其進行多出口率定,然后采用CLUE-S模型模擬了城鎮化發展情景下的土地利用空間分布,并在現狀土地利用基礎上進行退耕還林操作得到了退耕還林情景下的土地利用空間分布,進而模擬了3種土地利用情景(基準期、城鎮化情景及退耕還林情景)下的日徑流序列,并根據模擬的日徑流序列對水文極值進行了統計分析,主要結論如下:

(1)建立的SWAT模型能較好模擬和重現雅礱江流域基準期水文極值,基本能夠重現雅江站、洼里站及小得石站水文極值的發生時間和量級等主要特征,模型率定期和驗證期的Nash-Sutcliffe系數NS、確定性系數R2均在0.70以上,除洼里站的驗證期外PBIAS的絕對值小于10,模擬精度較高。

(2)城鎮化土地利用情景下,城鎮面積占比由0.12%增加到0.37%,耕地和未利用地面積占比略有減少;雅江站、麥地龍站及洼里站的日最大流量整體上增加,日最大流量序列均值增幅依次為4.6%、3.4%及2.3%,增幅從上游至下游逐步減小。

(3)退耕還林土地利用情景下,耕地面積占比由5.43%降至2.85%,林地和草地面積占比相應增加;雅江站、麥地龍站及洼里站的日最大流量略有減少,但減少幅度十分有限,洼里站減少幅度最大,也僅為-1.2%。

(4)雖然城鎮化發展會對雅礱江流域水文極值的影響有限,但如果后續城鎮化迅速發展(本文設置的城鎮化情景城鎮面積占比為0.37%,小于涼山彝族州規劃的城鎮面積占比1.01%),再加之氣候變化的影響,二者疊加后對雅礱江流域水文極值的影響不容忽視。因此在后續流域綜合規劃修編過程中,應按照集約用地、適當控制城鎮建設規模的原則,科學規劃城鎮建設布局,合理控制城鎮建設用地增長速度;同時應進一步優化調度,充分發揮流域梯級水庫群的調蓄作用,盡可能降低氣候變化及土地利用變化對流域防洪安全造成的不利影響。

猜你喜歡
城鎮化模型
一半模型
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
3D打印中的模型分割與打包
讓老年農民挑起城鎮化的重擔?
FLUKA幾何模型到CAD幾何模型轉換方法初步研究
堅持“三為主” 推進城鎮化
學習月刊(2015年14期)2015-07-09 03:37:50
城鎮化面臨的突出問題和應對之道
全球化(2015年2期)2015-02-28 12:38:55
城鎮化
江蘇年鑒(2014年0期)2014-03-11 17:09:40
加快推進以人為本的新型城鎮化
主站蜘蛛池模板: 亚洲综合婷婷激情| 国产精品网曝门免费视频| 一边摸一边做爽的视频17国产| 亚洲男人的天堂视频| 日韩a级片视频| 亚洲天堂网站在线| 狠狠做深爱婷婷综合一区| 国产成人精品一区二区| 中文字幕欧美日韩高清| 91人妻日韩人妻无码专区精品| 亚洲AV无码乱码在线观看代蜜桃| 深爱婷婷激情网| 国产精品无码AV片在线观看播放| 欧美日韩午夜视频在线观看 | 免费啪啪网址| 欧美一级视频免费| 精品少妇人妻av无码久久| vvvv98国产成人综合青青| 91无码人妻精品一区二区蜜桃| a在线观看免费| 五月天久久婷婷| 国产福利一区视频| 久久超级碰| 91麻豆精品国产高清在线| 国产福利一区二区在线观看| 婷婷成人综合| 免费一极毛片| 国产自在线播放| 久久久久青草大香线综合精品| 久久久四虎成人永久免费网站| 国产高清在线精品一区二区三区| 欧美在线国产| 亚洲国产理论片在线播放| 毛片基地美国正在播放亚洲 | 999精品在线视频| 91无码国产视频| 国产亚洲精品在天天在线麻豆 | 亚洲AⅤ波多系列中文字幕| 日韩av无码精品专区| 亚洲成aⅴ人在线观看| 99性视频| 天天摸天天操免费播放小视频| 亚洲第一视频免费在线| 99九九成人免费视频精品| 波多野结衣第一页| 精品一区二区三区视频免费观看| 亚洲国产欧美国产综合久久 | 精品人妻无码中字系列| 国产熟女一级毛片| 国产大片喷水在线在线视频| 自拍欧美亚洲| 国产精品久线在线观看| 国产精品无码翘臀在线看纯欲| 日韩无码黄色| 久久婷婷国产综合尤物精品| 99re在线观看视频| 成人一级免费视频| 亚洲人视频在线观看| 国产亚洲欧美在线人成aaaa| 亚洲福利一区二区三区| 亚洲嫩模喷白浆| 国产原创自拍不卡第一页| 婷婷丁香色| 日韩在线永久免费播放| 国产欧美视频在线| 日本在线国产| 91精品国产一区自在线拍| 亚欧乱色视频网站大全| 日韩av高清无码一区二区三区| 精品欧美一区二区三区久久久| 国产精品视频观看裸模| 91精品国产综合久久香蕉922| 色综合中文| 不卡无码网| 亚洲va在线∨a天堂va欧美va| 无码乱人伦一区二区亚洲一| 999精品色在线观看| 91啪在线| 国产精品天干天干在线观看| 成人无码一区二区三区视频在线观看| 国产91丝袜在线观看| 99国产在线视频|