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

DEM和土地利用分辨率對SWAT水沙模擬效果的影響
——以涇河流域為例

2023-11-13 01:18:42李睿智杜佰林郭宗俊許鎏佳
水資源與水工程學報 2023年5期
關鍵詞:影響模型

李睿智, 吳 磊 , 杜佰林, 郭宗俊, 王 妍, 許鎏佳

(1.西北農林科技大學 旱區農業水土工程教育部重點實驗室, 陜西 楊凌 712100; 2.西北農林科技大學 黃土高原土壤侵蝕與旱地農業國家重點實驗室, 陜西 楊凌 712100; 3.西北農林科技大學 水利與建筑工程學院, 陜西 楊凌 712100)

1 研究背景

21世紀以來,隨著地理信息系統(Geographic Information System, GIS)技術與遙感(remote sensing,RS)技術的廣泛應用,各種水文模型被廣泛應用于水文和產沙機理過程模擬[1]以及水土保持措施的規劃與實施[2]研究中,如RUSLE(revised universal soil loss equation)模型[3]、WEPP(Water Erosion Prediction Project)模型[4]、DYRIM(digital Yellow River integrated model)模型[5]以及SWAT(soil and water assessment tool)模型[6]等。水文模型的模擬精度很大程度上取決于輸入數據的空間精確度,不同空間分辨率的數字高程數據(digital elevation model, DEM)和土地利用數據會對模型模擬結果造成一定的影響[7-8]。因此,有必要探究輸入數據的空間精度與模擬結果之間的關系,以提高模型的模擬精度。

基于此,國內外學者開展了與DEM分辨率相關的研究,發現DEM的分辨率會對流域地形細節的提取造成影響,一般而言,隨著分辨率的增加,地形細節會更加精確[8],然而高精度的數據并不總能提高模型性能。比如,郭偉玲等[9]通過研究發現高精度的DEM數據并不能提取出更準確的坡長;Gautam等[10]研究表明在尼泊爾的兩個流域使用分辨率為30 m或更精細的DEM數據能獲得最佳的徑流模擬性能;Tan等[11]在Kelantan River流域研究發現不同分辨率的DEM數據以及流域閾值均會對模型性能產生影響,DEM分辨率為20~60 m時可以取得更好的月徑流模擬效果;蔡朵朵等[12]發現過低的DEM分辨率會降低北洛河流域的徑流模擬效果。陳海濤等[13]發現DEM分辨率對潮河流域SWAT徑流模擬效果影響并不顯著,但會對總氮的模擬結果產生較大影響。另外,有研究發現土地利用數據質量對流域水文建模較為重要[7],Jin等[14]研究發現高分辨率的土地利用數據會提高1.1%~6.9%的徑流模擬性能,而Fisher等[15]研究表明過于精細的土地利用數據會大幅度增加經濟和時間成本。但也有部分研究表明,土地利用精度與模型模擬結果的相關性不強[16]。綜上,不同分辨率的輸入數據會對流域邊界、河網和坡度等基本參數產生影響[17],導致模型水沙的模擬計算過程產生偏差。現關于輸入數據對模型模擬效果方面的研究多集中在DEM或土地利用等單因素對流域SWAT模型模擬的影響,且輸入數據分辨率選取的數量不足,鮮有對嵌套流域多分辨率DEM和土地利用數據交互影響的多模型模擬效果的評價研究。

涇河是黃河流域十大水系之一,針對涇河流域開展水沙研究對于黃河流域的生態保護和高質量發展具有重要意義[18]。因此,本研究通過構建涇河流域的SWAT模型,探究DEM與土地利用分辨率對模型性能的交互影響,并利用TOPSIS(technique for order preference by similarity to an ideal solution)優選理論模型針對水沙模擬結果的多個指標進行綜合評價,評選出涇河流域及其嵌套流域(蒲河流域、涇河南支干流流域)DEM與土地利用分辨率最佳組合的SWAT水沙模擬模型,該結果可為后期涇河流域的水沙研究提供技術參考。

2 數據來源與研究方法

2.1 研究區概況

涇河流域(106°14′~108°42′E, 34°46′~37°19′N)位于黃土高原中部,流經陜、甘、寧部分地區,流域面積約為45 000 km2,河網長度為455 km,地勢西北高、東南低。涇河發源于寧夏回族自治區涇源縣六盤山東麓,在陜西省高陵縣注入渭河,是渭河的一級支流、黃河的二級支流。流域內地形復雜,最小高程為219 m,最大高程為2 890 m,是黃土高原水土流失最嚴重的區域之一。流域處于半濕潤半干旱氣候的過渡地帶,為典型的溫帶大陸性氣候,多年平均氣溫為8 ℃,年降雨量在350~600 mm之間,降雨大多集中在7、8月份,徑流年際變化明顯。流域土地利用類型以草地和耕地為主,土壤以石灰性雛形土和鈣積潛育土為主。涇河流域的地理位置及水系、高程分布如圖1所示。

圖1 涇河流域高程、水系及水文氣象站點分布

2.2 研究數據及來源

2.2.1 DEM數據 在本研究中,選用了ASTER GDEM V2 DEM數據,數據來源于地理空間數據云(https://www.gscloud.cn),空間分辨率為30 m×30 m。為獲得研究所需的不同分辨率的數據,基于最鄰近分配法將 DEM 數據分別重采樣至60、90、150、300、500、750、1 000、2 000、3 000 m共9組不同空間分辨率的 DEM 數據。其中,3種不同量級(30、300、3 000 m)分辨率的DEM數據如圖2所示,圖中還顯示了同一局部位置的放大圖。

圖2 涇河流域3種不同量級分辨率的DEM數據圖

2.2.2 土地利用數據 土地利用數據來源于地理國情監測云平臺(http://www.geodata.cn/),空間分辨率為30 m×30 m。自20世紀80年代起,涇河流域開始實施大規模的水土保持措施,為降低淤地壩、水庫等客觀因素對模擬結果的影響,本研究選用1980年的土地利用數據。與DEM數據處理方法一致,采用最鄰近分配法將土地利用數據分別重采樣至60、90、150、300、500、750、1 000、2 000、3 000 m共9組不同空間分辨率的土地利用數據,其中,不同量級(30、300、3 000 m)分辨率的土地利用數據如圖3所示,圖中也顯示了同一局部位置的放大圖。

圖3 涇河流域3種不同量級分辨率的土地利用數據圖

2.2.3 土壤數據 本研究所需的土壤數據來源于HWSD(Harmonized World Soil Database)數據集(http://westdc.westgis.ac.cn/),空間分辨率為1 km×1 km。

2.2.4 氣象與水文數據 涇河流域內固原、長武、西峰、環縣4個氣象站的逐日氣象數據來源于中國氣象數據網,毛家河、楊家坪及張家山水文站逐日徑流泥沙數據來源于水文年鑒及黃土高原科學數據中心,該3個水文站分別為蒲河流域、涇河南支干流以及整個涇河流域的控制站。為降低人類活動等客觀因素造成的影響,并與土地利用數據的年份相匹配,采用1972—1986年的水文數據校準SWAT模型。

2.3 SWAT模型的構建與校準

SWAT模型是一個基于物理的半分布式流域水文模型,該模型可從DEM數據中提取坡度、河流的長度與寬度等水文參數。當流域內不同柵格的土地利用、土壤以及坡度類型相匹配時,SWAT會將這些柵格劃分為一個水文響應單元(hydrological response unit,HRU)進行處理。并計算HRU水平上的水文參數、徑流量、產沙量等水文水質信息,HRU的數量以及分布與輸入數據的空間分辨率相關。該模型被廣泛應用于評估土地管理和氣候變化對流域水質水量的影響[19]。

本研究基于土壤、氣象數據以及不同分辨率的DEM和土地利用數據構建了SWAT模型,并根據張家山站實測的水沙數據來確定模型參數。在分析DEM與土地利用分辨率對模型的影響時,控制土壤、氣象數據以及模型的參數類型保持不變。在建模的過程中,面積閾值統一設定為80 000 hm2,坡度統一劃分為3個等級(0°~7°,7°~25°,大于25°),HRU定義閾值統一設定為土地利用5%、土壤10%、坡度10%。使用SUFI-2算法程序對模型的參數進行校準,并選取該程序中的全局敏感性分析方法對模型的參數進行敏感性分析,其中t-Stat代表參數敏感性程度,P-Value代表參數敏感的顯著性程度,t-Stat越大、P-Value越接近于0,則參數的敏感性越強。

選取毛家河、楊家坪和張家山站實測的徑流、泥沙數據對模型逐一進行水沙同步校準。選取納什效率系數NSE和決定系數R2作為3個水文站點的模擬結果評價指標。其中,NSE為小于1的數值,可為負值;R2為小于1的正數,該兩個值的上限數值均為1,兩者數值大小越接近上限值1則模型的模擬結果越好,模擬值與實測值越接近。所構建的100個校準和不確定分析程序(calibration and uncertainty programs, CUP)工程均包含3個水文站的率定,每個站點均迭代3 000次,并對校準后的結果進行討論分析。

2.4 TOPSIS模型優選評價

TOPSIS是多目標決策方法的一種,其基本原理是計算評價對象到正理想解和負理想解的距離,然后根據正、負理想解的距離計算出相對貼近度并進行排序[20]。該方法的計算步驟如下:

首先采用平方和歸一化對原始數據進行無量綱處理并構建評價矩陣Z=(zij)m×n(i=1,2,…,m;j=1,2,…,n);zij為第i個評估對象的第j個評價指標。平方和歸一化的公式如下:

(1)

式中:xij為第i個評估對象第j個評價指標的原始數據。

然后確定模型的正理想解Z+和負理想解Z-,其分別代表理想中的最優方案和最差方案,計算公式如下:

Z+=(maxzi1,maxzi2,…,maxzin)

(2)

Z-=(minzi1,minzi2,…,minzin)

(3)

式中:maxzij、minzij分別為第j個評價指標中的最大值、最小值。

(4)

(5)

計算每個待選方案的相對貼近度Ci,計算公式如下:

(6)

最終根據相對貼近度Ci值的大小對待選方案進行排序,評選出最優的模型。

3 結果與分析

3.1 水文邊界與河網

受空間分辨率的影響,基于較低空間分辨率的DEM數據生成的河網和水文邊界可能與實際有所偏差。為解決建模中DEM空間分辨率低于1 000 m時河網提取與實際差別較大、不連續等問題,本研究采用SWAT模型中自帶的Burn in工具,對1 000 m分辨率DEM數據產生的河網文件進行修正,以獲得連續河網。不同分辨率DEM數據提取的河網與流域邊界如圖4所示。由圖4可以看出,當DEM分辨率在1 000 m及以下時,涇河流域水文邊界的劃分并未產生較大的差異,提取的河網接近實際河網。而當DEM的分辨率在2 000和3 000 m時,北支涇河流域河網的提取在部分支流中出現了較大的偏差,產生了一些偽河道(圖4 (a)、4 (b)),流域西南部部分邊界發生明顯的凹陷(圖4 (c))。DEM分辨率的降低會減弱水文邊界以河網的提取能力,在涇河流域,30~1 000 m分辨率的DEM數據提取河網的能力能基本滿足SWAT建模的需要。

圖4 涇河流域不同分辨率DEM數據提取的河網與流域邊界

3.2 水文特征與HRU數量分析

為進一步表明涇河流域在不同DEM分辨率下的水文特征情況,在土地利用數據分辨率為30 m的條件下,針對建模過程中所提取的子流域數量、流域面積、坡度及高程等展開分析,如表1所示。

表1 涇河流域不同DEM分辨率數據提取的流域水文特征

由表1可知,DEM分辨率的變化對流域面積有著顯著的影響,模型中涇河流域的計算面積與DEM分辨率呈負相關關系。當DEM分辨率由30 m降低至3 000 m時,流域面積由44 420 km2減小至38 970 km2,面積減少了12.27%。當DEM的分辨率為30~150 m時,子流域數量穩定不變,當DEM的分辨率低于150 m時,子流域數量呈現減少趨勢。更粗糙的DEM分辨率會導致高程信息的誤差增加,與30 m分辨率的DEM數據相比,3 000 m分辨率的DEM最大高程低估了312 m,最小高程高估了139 m,平均高程降低了14.54 m,標準差降低了29.02 m。DEM分辨率的降低會導致高程范圍以及標準差的減小,這是由地形特征信息的丟失所導致的。DEM分辨率對于坡度的影響更為顯著,在30 m分辨率DEM的模型中,46.63%的地區的坡度在25°以上,僅有11.31%的地區的地形比較平坦,坡度在7°以下。而當分辨率低于500 m時,25°以上的坡度信息幾乎完全丟失,當分辨率低于2 000 m時,所有的坡度信息幾乎完全丟失,99%以上地區的坡度均呈現為7°以下。

此外,DEM和土地利用的分辨率還會對HRU數量產生交互影響,圖5給出了兩者分辨率對HRU數量的交互影響結果。圖5表明,在相同的土地利用分辨率下,隨著DEM分辨率的降低,HRU數量呈明顯的減少趨勢,但在相同的DEM分辨率下,土地利用分辨率的變化對HRU數量的影響有限。

圖5 DEM和土地利用數據的分辨率對HRU數量的交互影響

3.3 參數敏感性分析

采用SWAT-CUP的SUFI-2算法,對不同分辨率的DEM和土地利用數據構建的SWAT模型進行參數敏感性分析,探究DEM和土地利用數據對參數敏感性的影響。本研究以5個模型為例,分別對每個模型最敏感的前8個參數進行分析,各參數及其物理意義見表2,不同模型的敏感性參數匯總如表3所示。

表2 模型各參數及其物理意義

表3 不同分辨率模型的參數敏感性分析結果匯總

由表3可知,在所有5個模型中,USLE_C、EPCO、ALPHA_BNK以及CN2均為最敏感的參數。對比模型30D30L、90D30L、3000D30L來探究DEM分辨率對參數敏感性的影響,DEM分辨率變化較小時,參數的敏感性幾乎不發生變化,當DEM分辨率下降到3 000 m時,3個產匯流參數OV_N、SMTMP、SURLAG以及一個產沙參數USLE_P不再是最敏感的參數,此時新增的敏感性參數包括3個產沙參數BIOMIX、 LAT_SED、USLE_K 和一個匯流參數CH_N1。對比模型30D30L、30D90L、30D3000L可發現土地利用分辨率的變化對參數敏感性的影響不大,僅個別參數的敏感性排序略有改變。

3.4 DEM和土地利用分辨率對水沙模擬效果的交互影響

利用不同分辨率的DEM和土地利用數據共構建了100個SWAT模型,經過參數校準后不同水文站點的徑流或泥沙較優模擬效果(NSE系數較高)的分辨率組合如表4所示。在所有的模型中毛家河站徑流和泥沙的NSE系數均大于0.63,均能滿足模擬精度要求。由表4可見,毛家河站徑流模擬的NSE系數最佳為0.80,此時DEM和土地利用分辨率分別為2 000和150 m;泥沙模擬的NSE系數最佳為0.81,此時DEM和土地利用分辨率分別為500 和2 000 m。楊家坪站徑流模擬的NSE系數最佳為0.76,此時DEM和土地利用分辨率分別為30和90 m;泥沙模擬的最佳NSE系數為0.79,此時DEM和土地利用分辨率分別為30 和2 000 m,有3個SWAT模型徑流或泥沙模擬的NSE系數小于0.50,不滿足模擬要求。張家山站徑流模擬的最佳NSE系數為0.73,此時DEM分辨率為90 m、土地利用分辨率分別為30、60、90 m;泥沙模擬的NSE系數最佳為0.66,此時DEM和土地利用的分辨率均為2 000 m,該站共有22個SWAT模型泥沙模擬的NSE系數小于0.50。與其他兩個水文站點相比,張家山站的模擬效果較差,這是由于建模所用到的4個氣象站多分布在西南部,無法有效覆蓋北部支流,因而降低了北部支流的模擬精度,影響了張家山站的模擬結果。

表4 涇河流域不同水文站點水沙模擬效果較優的分辨率組合

以張家山站為例,徑流和泥沙模擬的NSE系數與DEM及土地利用分辨率的交互影響關系如圖6所示。通常模擬結果隨著DEM和土地利用分辨率的降低而變差,然而圖6顯示出徑流和泥沙的模擬結果隨著土地利用分辨率的變化并沒有發生明顯的改變,而隨著DEM的分辨率的降低呈現了“增大→降低→增大→降低”的不規則變化。

圖6 徑流和泥沙模擬的NSE系數與DEM及土地利用分辨率的交互影響關系(以張家山站為例)

3.5 模型綜合優選評價

采用TOPSIS對涇河流域3個水文站模型模擬效果進行綜合評價與方案優選。取徑流和泥沙模擬的NSE和R2共4個參數作為評價指標。對評價指標進行平方和歸一化處理,計算出正理想解距離、負理想解距離和相對貼近度C,據此對模型進行排序,相對貼近度越大,表明該模型的模擬效果越好,3個水文站點的相對貼近度計算結果如圖7所示。

圖7 基于TOPSIS理論的模型綜合優選評價

由圖7可見,張家山站模擬效果最好的DEM和土地利用分辨率分別為150和3 000 m,相對貼近度為0.85(徑流模擬:R2=0.71、NSE=0.69;泥沙模擬:R2=0.64、NSE=0.62),模擬效果最差的DEM和土地利用分辨率分別為300 和30 m。采用DEM和土地利用均為30 m分辨率的數據構建的模型排序僅為49,這也進一步說明并非DEM和土地利用的分辨率越高模型的模擬效果越好。毛家河站和楊家坪站模擬效果最優的DEM和土地利用分辨率組合分別為750和300 m以及30 和90 m,相對貼近度分別為0.86(徑流模擬:R2=0.81、NSE=0.78;泥沙模擬:R2=0.79、NSE=0.78)和0.95(徑流模擬:R2=0.76、NSE=0.76;泥沙模擬:R2=0.79、NSE=0.77)。

4 討 論

本研究表明,隨著DEM數據分辨率的降低,水文邊界的提取會產生偏差。這是由于DEM數據的分辨率會對地形參數的垂直精度產生影響[21],低分辨率的DEM數據會對高程信息產生概化作用[13],降低了提取數據的準確性,導致流域邊界的偏移及流域面積的缺失[22]。另外,SWAT模型中河道的提取基于D8算法[23],高程信息的偏差會對水流方向的計算產生較大影響。而且,DEM數據分辨率的降低會造成黃土丘陵區及殘塬區為主的涇河流域中部和北部地區[24]地形細節損失,更容易產生偽河道。

在相同的土地利用、土壤類型及坡度劃分閾值下,HRU的數量與流域柵格數、子流域數量和子流域中土地利用、土壤類型及坡度的分類數有關[25]。本研究中,DEM分辨率對HRU數量的影響占主導作用,這是由三方面原因導致的。第一,整個流域的柵格數量是由DEM分辨率決定的,DEM的分辨率越高則柵格的數量越多,HRU劃分的過程就會越精細。第二,DEM分辨率的降低會導致子流域數量減少,從而引起HRU數量的降低。第三,土地利用數據和土壤數據會被投影在DEM圖形柵格上,DEM分辨率的降低還會導致土地利用、土壤以及坡度的分類丟失[26],這也導致HRU數量受DEM分辨率的影響更大。此外,由于涇河流域的地形較為復雜,相較于地形特征的變化,相鄰柵格數據之間的土地利用變化較為平滑,從而導致了土地利用分辨率的變化對HRU數量的影響遠小于DEM。當DEM分辨率為30 m時,HRU的數量多達612個;當DEM分辨率為3 000 m時,HRU的數量減少至154個,同比減少了74.8%。

DEM數據分辨率的降低會造成坡度信息的丟失,地形趨于平緩,匯流時間增加,從而導致了部分參數的敏感性發生變化。曼寧系數OV_N和水土保持措施因子USLE_P的敏感性均與坡度相關,坡度越大則敏感性越強,而地表徑流滯后系數SURLAG會隨匯流時間的增大而更加敏感[27]。在DEM分辨率變化較小時,坡度和匯流時間變化幅度也較小,對參數的敏感性影響有限;當DEM分辨率變化較大時,OV_N、USLE_P以及SURLAG等參數的敏感性降低。

土地利用分辨率對模型水沙模擬效果和模型參數的敏感性的影響要小于DEM分辨率。具體來說,土地利用分辨率的變化會影響土地利用的類型及空間分布,但對土地利用類型的面積變化影響有限,在由30 m分辨率DEM數據構建的所有模型中,土地利用數據分辨率由30 m降低至3 000 m時,僅造成土地利用類型“牧草”的面積占比減少了1.57%以及土地利用類型“耕地”的面積占比增加了1.26%,其余土地類型的變化幅度均不超過0.2%。因此對整個流域的水沙模擬結果及參數的敏感性影響并不顯著。然而,DEM的分辨率會影響模型水文參數的提取,導致水文循環及泥沙侵蝕的模擬過程產生變化,從而對水沙模擬效果產生較大的影響。另外,30 m分辨率的DEM數據并未構建出精度最高的模型。這是因為過于精細的DEM數據會生成更多數量的HRU,導致SWAT-CUP的運算量大幅度增加,不僅延長了模型運行時間,還有可能導致優化算法陷入局部最優解,無法計算出最佳參數[28-29]。低分辨率數據構建的模型也可能具有良好的適用性,比如,在模型綜合評價中,毛家河站與張家山站最優模型的DEM和土地利用分辨率組合為750 和300 m以及150 和3 000 m,在水沙模擬的過程中均具有較好的適用性,這也與韓振宇等[30]的研究結果一致。

本研究基于TOPSIS理論優選出涇河流域SWAT水沙模擬模型的DEM和土地利用分辨率組合,可為涇河流域的生態保護和水沙治理提供技術參考。然而,由于各流域的下墊面情況不盡相同,該流域的建議輸入數據不一定適用于其他流域,各流域可根據自身條件進行模型模擬評價,優選出不同分辨率輸入數據的最佳組合。此外,后續研究可以考慮數據源和土壤數據的分辨率對模擬效果的影響。

5 結 論

(1) DEM分辨率的變化會影響水文邊界與河網的提取,30~1 000 m分辨率的DEM均能滿足建模需求。當DEM分辨率低于1 000 m時,會產生偽河道和不連續的河網,河網邊界也會發生凹陷。隨著DEM的分辨率的降低,流域的計算面積累積減少了12.27%,最大高程低估了312 m,最小高程高估了139 m。坡度受分辨率變化的影響更大,當分辨率低于2 000 m時,坡度信息幾乎全部丟失。HRU的數量會隨著DEM分辨率的降低而明顯減少,土地利用對HRU數量的影響較小。

(2) 涇河流域毛家河和楊家坪水文站的模擬效果優于張家山站,這與氣象站的坐標以及流域面積有關,參數自身的不確定性會影響模型的模擬效果。因此,低分辨率的模型也可能具有較好的適用性。在所有的模型中,USLE_C,EPCO、ALPHA_BNK以及CN2均為最敏感的參數,DEM對參數敏感性的影響大于土地利用。

(3) 由于模型和參數的不確定性,DEM和土地利用數據的分辨率對模型水沙模擬效果的影響并不明顯。不同流域的最佳模型的DEM和土地利用分辨率組合不同,蒲河流域、涇河南支干流流域以及整個涇河流域的最佳模型DEM和土地利用分辨率組合分別是750和300 m、30和90 m以及150和3 000 m。

猜你喜歡
影響模型
一半模型
是什么影響了滑動摩擦力的大小
哪些顧慮影響擔當?
當代陜西(2021年2期)2021-03-29 07:41:24
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
沒錯,痛經有時也會影響懷孕
媽媽寶寶(2017年3期)2017-02-21 01:22:28
3D打印中的模型分割與打包
擴鏈劑聯用對PETG擴鏈反應與流變性能的影響
中國塑料(2016年3期)2016-06-15 20:30:00
基于Simulink的跟蹤干擾對跳頻通信的影響
FLUKA幾何模型到CAD幾何模型轉換方法初步研究
主站蜘蛛池模板: 久久精品国产国语对白| 国产成人AV综合久久| 天天操天天噜| 中文成人无码国产亚洲| 国产毛片片精品天天看视频| 思思99思思久久最新精品| 四虎永久免费在线| 亚洲欧美综合在线观看| 尤物亚洲最大AV无码网站| av在线5g无码天天| 日韩精品资源| 波多野结衣国产精品| 91国语视频| 欧美亚洲一区二区三区导航| 国产在线观看91精品| 国产人人乐人人爱| 亚洲日本中文字幕乱码中文| 国产v精品成人免费视频71pao| 日韩视频福利| 亚洲中文在线看视频一区| 毛片网站在线播放| 亚洲人人视频| 日本www色视频| 91精品久久久无码中文字幕vr| 欧美精品成人一区二区视频一| 亚洲欧美日韩高清综合678| 日本不卡在线| 国产成人精品亚洲日本对白优播| 日韩av电影一区二区三区四区| 国产精品人莉莉成在线播放| 亚洲国产系列| 国产精品视频导航| 国产无人区一区二区三区| 一级毛片高清| 91青青草视频在线观看的| 日韩在线1| www.国产福利| 亚洲六月丁香六月婷婷蜜芽| 亚洲成人精品| 伊人国产无码高清视频| 91无码网站| 福利在线一区| 一本一道波多野结衣av黑人在线| 国产不卡国语在线| 极品av一区二区| 四虎亚洲国产成人久久精品| 久久伊人色| 园内精品自拍视频在线播放| 亚洲国产看片基地久久1024| 人妻一本久道久久综合久久鬼色| A级毛片无码久久精品免费| 日本午夜精品一本在线观看| 日韩精品视频久久| 久久精品视频亚洲| 丰满人妻中出白浆| 在线播放国产99re| 欧美日韩国产综合视频在线观看| 久久人人爽人人爽人人片aV东京热| 红杏AV在线无码| 日韩视频福利| 伊人天堂网| 91青青草视频| 欧美性精品不卡在线观看| 视频在线观看一区二区| 亚洲综合久久一本伊一区| 女人18毛片一级毛片在线 | 国产精品30p| 超碰aⅴ人人做人人爽欧美 | 亚洲精品国产精品乱码不卞| 精品三级网站| 国产在线视频导航| 成人免费黄色小视频| 四虎永久免费在线| 久久6免费视频| 国产拍在线| 亚洲天堂首页| 成人国产精品网站在线看| 99这里只有精品免费视频| 成人在线第一页| 亚洲中文字幕97久久精品少妇| 国产av剧情无码精品色午夜| 欧美成人精品在线|