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

甘肅省高寒草甸植被覆蓋度反演及其時空變化研究

2021-04-09 03:36:06何國興柳小妮張德罡蒲小鵬劉志剛關文昊楊軍銀韓天虎潘冬榮
草地學報 2021年3期
關鍵詞:模型研究

何國興, 柳小妮*, 張德罡, 李 強, 蒲小鵬, 劉志剛,關文昊, 楊軍銀, 韓天虎, 孫 斌, 潘冬榮

(1. 甘肅農業大學草業學院,甘肅 蘭州 730070;2. 草業生態系統教育部重點實驗室(甘肅農業大學),甘肅 蘭州 730070;3. 甘肅省草原技術推廣總站,甘肅 蘭州 730000)

作為我國重要的國土資源和可再生性的自然資源,草地具有重要的生態價值和經濟價值[1]。FVC不僅可以直觀反映草原地上植被生長狀況,同時也是指示草地健康狀況及監測環境變化的常用評估指標之一[2]。地面實際測量和遙感反演是獲取植被覆蓋度數據的兩種主要方法[3-4]。地面測量通常采用針刺法和目測法[5],但地面測量具有較強的主觀性、效率較低、成本較高等缺點,導致其無法大面積推廣[6]。隨著遙感技術的不斷發展,利用遙感影像反演植被覆蓋度已成為較成熟的技術,其主要方法有像元二分法、經驗回歸模型法和混合像元分解法[7]。由于遙感反演具有效率高、成本低、客觀性較強等優點,已成為在大空間尺度和長時間序列范圍內獲取草地覆蓋度數據和監測評估草地健康狀況的重要手段[8]。因此,構建遙感反演模型提取區域植被覆蓋度的方法不僅是當前研究區域生態模型的基礎工作,也是今后開展生態工程評價的主要方向。

中分辨率成像光譜儀(moderate resolution imaging spectroradiometer,MODIS)是地球觀測系統(earth observation system,EOS)的核心傳感器,它的任務主要是每日獲取4次包括大氣、海洋、陸地等在內的地球系統相關要素變化數據[9]。游浩妍等[10]基于MODIS數據對研究區的植被覆蓋度利用像元二分模型進行估算,發現基于NDVI反演模型擬合效果最佳;周偉等[11]基于MODIS NDVI數據反演了中國草地覆蓋度的空間格局和變化趨勢,揭示了草地覆蓋度變化機理;宋清潔等[12]基于MODIS增強型植被指數(enhanced vegetation index,EVI)和NDVI建立了研究區草地植被覆蓋度的反演模型,結果表明基于EVI構建的對數模型為最優反演模型;Wiesmair等[13]使用隨機森林回歸建立了估算格魯吉亞高加索草地FVC變化的遙感模型。

甘肅省作為全國重點牧區之一,草原面積2.57×108hm2[14],其中高寒草甸是一重要類型,具有調節氣候、涵養水源、保持水土、凈化空氣等重要的生態功能[15]。高寒草甸作為全球公認的氣候變化敏感區和生態脆弱區,它的形成與青藏高原獨特的地形和氣候條件有著密不可分的關系[16]。同時,高寒草甸也是調節高寒區草地生態平衡和氣候變化的基質,不僅有助于畜牧業的發展,而且對維持生態系統穩定性具有十分重要的作用[17]。目前利用遙感技術獲取植被覆蓋度數據、監測草地生長狀況的研究中大多使用的是基于像元二分法或直接利用NDVI去反映植被覆蓋度,但對于某一特定草地類型覆蓋度反演的研究較少,不同草地類型對植被指數的敏感程度和模型的精度又不盡相同[2]。因此利用經驗回歸模型法將地表實測數據和VI相結合構建甘肅省高寒草甸FVC遙感反演模型,可以提高甘肅省高寒草甸FVC的遙感反演精度,對高寒草甸草地生態可持續發展和恢復重建具有十分重要的意義。

本研究基于地面實測樣點數據和MODIS產品數據,在研究高寒草甸不同植被覆蓋度遙感監測模型基礎上,模擬分析2000—2019年甘肅省高寒草甸生長季植被覆蓋度的時空變化特征,以期為高寒草甸草地資源的動態監測和健康評價提供理論依據和技術支撐。

1 材料與方法

1.1 研究區概況

以甘肅省高寒草甸草原為研究區(圖1),行政區域包括甘南藏族自治州、天祝藏族自治縣、肅南裕固族自治縣和山丹縣,位于97°23′~104°45′ E,33°06′~39°43′ N。境內地形十分復雜,山地、高原占總面積的95%以上,平均海拔3 000 m,從西北至東南地形逐漸降低。研究區屬于高原大陸性氣候,年降水量400~800 mm,年均溫1~3℃,年日照時數2 200~2 400 h,自西北至東南逐漸降低[12]。研究區高寒草甸面積40 510 km2,占總面積的36%。

高寒草甸是指以冷中生多年生草本植物為優勢形成的植物群落,該類草地主要以嵩草屬植物為主,草層低矮、叢生,結構簡單,層次分化不明顯,葉片縮小,被絨毛或呈墊狀;主要分布在甘南高原和祁連山山地并向東延伸到海拔3 000 m以上的西秦嶺山地。

1.2 數據來源及預處理

1.2.1草地類型數據 將來源于中國草業與生態網絡(http://www.ecograss.com.cn/)的《中國1∶100萬草地資源圖》經過校正后得到本研究中的草地類型數據[18];甘肅省1∶100萬行政區劃圖來源于地球科學數據共享平臺(http://www.geodata.cn)。

圖1 研究區及采樣點示意圖Fig.1 Study area and distribution of sampling points

1.2.2地面實測數據 2014年研究區草原調查數據,來源于甘肅省草原總站,共370個樣點(圖1),包括經緯度、海拔高度、草地類型、草地覆蓋度等。所有樣點的調查時間集中在2014年8—9月,該時期是草地植被生長旺盛期。將樣點數據在R 3.6.1中隨機抽樣70%和30%,分別作為建模數據組和檢驗數據組。

1.2.3遙感數據 本研究選用美國航空航天局(national aeronautics and space administration,NASA) 對地觀測系統數據共享平臺(https://lpdaacsvc.cr.usgs.gov/appeears/explore)提供的2000—2019年MOD13Q1產品,包括NDVI和EVI,空間分辨率250 m,時間分辨率16 d。紅光波段和近紅外波段數據為MOD09 GQ. 006產品,空間分辨率250 m,時間分辨率8 d。

1.3 研究方法

1.3.1植被指數選擇 差值植被指數(difference vegetation index,DVI)雖與植被覆蓋度相關性較高,但土壤背景較強的影響該指數。通常在較低覆蓋度時,RVI隨著植被的增加而增加,當植被覆蓋度大于80%后,RVI的敏感程度呈現下降趨勢[2]。NDVI具有植被探測范圍較寬、靈敏度較高的優點,是一種被廣泛使用的植被指數,該植被指數可以在不同程度上克服地勢、太陽高度角、大氣以及群落結構等因素的影響[19]。EVI能分離紅光和藍光通過的氣溶膠,補償氣體殘留單體對紅光的吸收,能夠反映植被的季節性[20]。土壤調整植被指數(soil adjusted vegetation index,SAVI)由于加入土壤調節系數L,可以減弱土壤背景的影響。修正土壤調整植被指數(modified soil adjusted vegetation index,MSAVI)也能夠有效減弱土壤背景的影響[21]。基于此,選擇以上植被指數(表1),進行研究區高寒草甸FVC反演。

表1 植被指數計算公式Table 1 Vegetation index formula

運用MRT (modis re-projection tool)讀取數據,同時轉換遙感影像的數據格式和地圖投影;利用ENVI 5.3中的Band math工具將MOD13Q1數據和MOD09 GQ.006數據的每個像元值縮小10 000倍,得到真實的NDVI、EVI值和紅光波段、近紅外波段值[22]。根據計算公式對DVI,RVI,NDVI,EVI,SAVI和MSAVI進行波段組合計算。將植被指數數據采用最大值合成法(maximum value composite,MVC)[23]生成2000—2019年8月份月植被指數,以減少云、積雪、大氣和太陽高度角等外界因素的影響。

1.3.2FVC反演模型構建

(1)經驗回歸模型法 利用ArcGIS 10.2中Spatial Analyst工具提取各經緯度點植被指數值,將FVC實測值與提取的各植被指數在R 3.6.1中進行Pearson相關性分析。選擇線性模型、對數模型、指數模型、二次多項式模型和冪函數模型等具有代表性的常用回歸模型[16],構建FVC與各植被指數間的模型,再根據精度結果篩選最佳模型。

式中:FVC為植被覆蓋度;NDVIsoil為裸地土壤;NDVIveg為有植被完全覆蓋區域的NDVI值。

(3) 模型精度檢驗采用回歸方程的決定系數R2、均方根誤差(root mean square error,RMSE)來評價模型精度。RMSE常用來衡量模型的估測能力、描述模型的不確定性以及量化模型精度,RMSE數值越低,表明回歸模型精確度越高。

以30%草原調查中的草地FVC為真實值,對所建立的回歸模型進行精度評價,分別計算各反演模型估測的草地FVC與真實值之間的RMSE(公式如下),篩選出研究區不同草地類型草地FVC監測的最佳遙感反演模型。

綜合模型的R2,RMSE、實測值與預測值回歸方程的R2和斜率等因素,確認最優模型,并應用其反演并分析2000—2019年高寒草甸FVC動態變化。

1.3.3FVC分級 依據全國《草場資源調查技術規程》及植被劃分標準[25-27],并結合甘肅省高寒草甸植被特點,將FVC分為5個等級:低FVC(0~20%)、較低FVC(20~40%)、中FVC(40~60%)、較高FVC(60~80%)、高FVC(80~100%)。

1.3.4FVC趨勢分析 采用一元線性回歸分析法和趨勢法(Slope)進行像元尺度的計算,分析20年間高寒草甸FVC的變化趨勢和變化速率。計算公式如下[28]:

式中:θSlope為草地FVC變化趨勢;FVCi為第i年的FVC;i為年變量。當θSlope<0時,FVC減少;θSlope>0時,FVC增加。

采用F檢驗法進一步分析草地FVC隨時間發生的變化趨勢。計算公式如下[29]:

式中:n為時間序列長度20;i為2000—2019年的年份序號1~20;i為年份序號的平均值;FVCi為第i年的FVC值;FVC為2000—2019年FVC的平均值。

參照已有的研究[29-31],根據θslope值的大小將研究區高寒草甸FVC變化劃分為7個等級:顯著增加、中度增加、輕度增加、基本不變、輕度減少、中度減少、顯著減少,并統計各區間的面積百分比。

1.3.5植被覆蓋度變異系數 變異系數(coefficient of variation,Cv)可以反映高寒草甸FVC多年的變化波動情況,Cv≤0.15說明高寒草甸FVC變化較小,草地生長狀況相對穩定;Cv>0.15說明FVC變化較大,草地生長狀況不夠穩定[32]。本研究通過計算像元尺度的變異系數,進而了解研究區高寒草甸FVC的波動情況。計算公式如下[33]:

參照已有的研究[34],根據Cv值的大小將高寒草甸FVC穩定性分為5個等級:十分穩定(Cv≤0.10)、穩定(0.100.30),并統計各區間的面積百分比。

2 結果與分析

2.1 FVC遙感反演模型的構建

2.1.1FVC與各植被指數間相關性 由6種植被指數DVI,RVI,NDVI,EVI,SAVI和MSAVI與FVC實測值的相關性分析(圖2)結果可知,各植被指數與FVC均顯著正相關(P<0.05)。

其中,NDVI與FVC間的相關系數最高,達到0.90;EVI和SAVI與FVC間的相關性高于DVI和RVI。6種植被指數中RVI與FVC間的相關系數最低(0.59),其余植被指數與FVC的相關性系數均高于0.65,因此選擇DVI,NDVI,EVI,SAVI和MSAVI與FVC進行經驗模型的構建。

2.1.2FVC與植被指數間反演模型的建立 基于FVC與DVI,NDVI,EVI,SAVI和MSAVI構建反演模型,選擇R2>0.60的模型[35]和像元二分模型計算RMSE(表2)。

5個植被指數模型中,NDVI模型優于其他模型。其中NDVI的指數模型、二項式模型和冪函數模型的R2均高于0.80,而NDVI的二項式模型、冪函數模型和線性模型的RMSE均小于其他模型。EVI的二項式模型(R2=0.64)優于冪函數模型(R2=0.63)。SAVI的對數模型(R2=0.64)和冪函數模型(R2=0.68)優于MSAVI的二項式模型(R2=0.61)。

2.1.3FVC與植被指數間反演模型精度檢驗 NDVI的線性模型、指數模型、對數模型、二項式模型和冪函數模型,以及EVI的二項式模型,SAVI的對數模型和冪函數模型、像元二分模型共9個模型的檢驗結果如圖3所示。

圖2 Pearson相關性矩陣Fig.2 Pearson correlation matrix注:相關性矩陣分為兩部分,上三角為相關系數與顯著性檢驗,下三角為散點圖與擬合曲線,**與***代表在0.05和0.01水平上差異顯著。各變量間相關關系相交于變量的行列處Note:The correlation matrix is divided into two parts,the upper triangle is the correlation coefficient and the significance test,the lower triangle is the scatter chart and the fitted curve,** and *** indicate significant difference at the 0.05 and 0.01 level,respectively. The correlations between variables can be found at the row and column intersection of the two variables

表2 植被覆蓋度與各植被指數間的反演模型Table 2 Inversion model of vegetation fractional cover and vegetation indices

圖3 FVC實測值與預測值比較Fig.3 Comparison of measured and predicted values of FVC注:(a)NDVI線性模型;(b)NDVI指數模型;(c) NDVI對數模型;(d) NDVI二項式模型;(e) NDVI冪函數模型;(f) EVI二項式模型;(g) SAVI對數模型;(h) SAVI冪函數模型;(i)像元二分法。圖中虛線為預測值等于真實值,實線為預測值與真實值間的回歸直線,實線離虛線越近說明預測結果越準確Note:(a) NDVI Linear model;(b) NDVI Exponential model;(c) NDVI Logarithmic model;(d) NDVI Quadratic polynomial model;(e) NDVI Power function model;(f) EVI Quadratic polynomial model;(g) SAVI Logarithmic model;(h) SAVI Power function model;(i) Pixel dichotomy. In the figure the dashed line is the predicted value equal to the measured value. The solid line is the regression line between the predicted value and the measured value. The closer the solid line is to the dotted line,more accurate the prediction

NDVI二項式模型的擬合R2=0.89,高于其他模型。像元二分模型線性擬合方程的斜率為0.90,高于其他模型。

綜合模型R2,RMSE,實測值與預測值回歸方程的R2和斜率等因素,確認NDVI的二項式模型:y=—0.65x2+1.97x—0.23為最優模型。

2.2 2000—2019年高寒草甸FVC的時空變化

2.2.1FVC的空間分布特征 應用NDVI二項式模型y=—0.65x2+1.97x—0.23反演2000—2019年高寒草甸FVC的動態變化。由20年高寒草甸均值空間分布(圖4a)可知,高寒草甸FVC整體較高,其中FVC>80%的面積最大(47.44%),60

從區域分布來看,研究區南北跨越較大,高寒草甸FVC呈現明顯變化趨勢。瑪曲縣、碌曲縣、夏河縣、合作市和天祝藏族自治縣等海拔較高且人煙稀少的區域,高寒草甸FVC較高(>80%);60

2.2.2FVC的時間變化規律 由高寒草甸FVC年際變化圖(圖4b)中可以看出,2000—2019年高寒草甸年均FVC呈現波動上升趨勢(P<0.05),分別在2001年、2007年和2012年FVC出現峰值,隨后FVC下降。線性擬合發現,高寒草甸FVC以0.15%·a-1(R2=0.58)的速率呈顯著上升趨勢(P<0.05)。從整體變化來看,2002—2007年、2008—2013年和2015—2019年這3個時間段內FVC均呈現上升趨勢。高寒草甸FVC大約每5年會出現明顯波動。

將高寒草甸2000—2019年每年FVC等級進行統計(圖4c)后發現,高寒草甸年均FVC大部分以中、高覆蓋度為主,說明高寒草甸整體上植被覆蓋度較高,生長狀況良好。2000年極高覆蓋度水平比例最小,極低和低覆蓋度水平比例達到研究年份的最大。其余年份內低和極低覆蓋度水平的比例基本趨于穩定,高和極高覆蓋度水平的比例呈現此消彼長的變化趨勢。

圖4 2000—2019年高寒草甸FVC空間分布及年際變化趨勢Fig.4 Alpine meadow FVC spatial distribution and annual change trends from 2000 to 2019

2.2.3FVC空間趨勢變化分析 趨勢分析結果表明(圖5):2000—2019年FVC呈上升趨勢的面積為21 078 km2,占高寒草甸總面積的52.76%;穩定趨勢的面積為11 017 km2,面積占27.58%;下降趨勢的面積為7 852 km2,面積占19.66%。其中輕度增加的面積比例最高,占33.14%;顯著減少所占的面積比例最低,占0.82%。

圖5 2000—2019年高寒草甸FVC變化趨勢分級Fig.5 Alpine meadow FVC change trends classification from 2000 to 2019

從空間分布來看,FVC顯著增加的區域主要分布在肅南裕固族自治縣;中度和輕度增加的區域主要分布在夏河縣、合作市、天祝藏族自治縣和肅南裕固族自治縣;輕度和中度減少的區域主要分布在碌曲縣的中部和瑪曲縣的東南部;顯著減少的區域零散分布于瑪曲縣的邊緣和碌曲縣的中部。

2.2.4FVC穩定性分析 研究區2000—2019年高寒草甸Cv≤0.15的面積占總面積的86.48%,高寒草甸FVC這20年內波動較小;Cv<0.10的面積占49.59%;0.100.15)區域主要分布在肅南裕固族自治縣、民樂縣、天祝藏族自治縣、夏河縣、瑪曲縣和舟曲縣的邊緣地帶,其面積占高寒草甸總面積的13.52%;高波動區域(Cv>0.30)主要分布在肅南裕固族自治縣的西北部,總面積1 095 km2(圖6)。

圖6 2000—2019年高寒草甸FVC變異系數分級Fig.6 Alpine meadow FVC coefficient of variation classification from 2000 to 2019

3 討論

3.1 高寒草甸FVC遙感反演模型

植被指數是將遙感數據相關的光譜信號綜合在一起,通過不同波段的組合,增強不同的植被信息,降低干擾[36]。植被指數受到植被自身、大氣、地物等因素的影響,能夠反映綠色植物的相對豐度,與FVC和葉面積指數(leaf area index,LAI)等指標具有較強的相關性[37]。然而,植被指數不具備物理意義,只有將其轉換為具體的生態參數(生物量、植被覆蓋度、葉面積指數等)后,植被指數才能夠定量分析生態環境的變化過程[38]。

本研究發現,NDVI和SAVI與FVC的相關性高于其他植被指數,而RVI與FVC間的相關性低于0.60,可能是由于高寒草甸地勢復雜、日照強烈、紫外線強、潮濕寒冷、晝夜溫差大、植被生長期短等因素導致[39]。利用像元二分法和經驗回歸模型法對高寒草甸植被覆蓋度構建的反演模型,其中經驗回歸模型法NDVI和SAVI兩個植被指數與植被覆蓋度的相關性較強,NDVI二次多項式回歸模型精度驗證擬合效果達到0.89,能夠較好地反演高寒草甸植被覆蓋度,這與巴音布魯克高寒草原遙感反演模型研究的結果相一致[2]。另外,陳建軍等[38]研究也發現,NDVI二次多項式模型對黃河源區高寒草地FVC具有較高的反演精度。盡管部分研究表明NDVI會出現飽和現象,但由于高寒草甸是以嵩草為主的草甸草原,單位面積內的生物量有限、草層較低導致NDVI的飽和現象不明顯[38]。像元二分模型的反演精度略低于NDVI二項式模型,這與李藝夢等[40]基于Landsat8影像植被覆蓋度估算的研究結果有相似之處。

3.2 高寒草甸FVC趨勢變化

本研究分析結果表明,2000-2019年高寒草甸年均FVC呈現波動上升趨勢,高寒草甸FVC呈現西北低、東南高的分布格局。另外,20年平均FVC顯示,高寒草甸FVC整體呈現較高水平,其中高FVC區域主要在甘南州地區,這與王偉軍等[41]在甘南州草地NDVI動態變化的研究結果相一致。周偉等[42]利用MODIS數據研究中國草地覆蓋度時空動態時也發現,甘南州整體FVC處于較高水平。趨勢分析結果表明,20年間高寒草甸FVC呈增加趨勢的面積高于降低趨勢的面積。其中呈減少趨勢的地區主要在瑪曲縣的東部和碌曲縣的中部,這與楊瑞瑞[43]在黃河源植被覆蓋度的研究結果相一致。其他地區均呈現增加趨勢。變異系數結果表明,高寒草甸FVC相對穩定的比例占86.57%,說明20年間高寒草甸FVC整體呈現穩定趨勢。

究其原因,是因為甘肅省高寒草甸大部分位于甘南州境內,而甘南州又地處青藏高原東北邊緣,是青藏高原與黃土高原的交匯區,在氣候上屬于亞洲季風邊緣區域和高寒氣候的過渡區,屬于華西秋雨區[44]。秋季連陰雨對甘南高寒草甸草原的返青和枯黃具有十分密切的關系[45]。秋季較高的降水量是該區域高寒草甸FVC高于其他區域FVC主要原因,適宜的水份直接影響植物的光合作用和呼吸作用,是決定植物生長狀況的關鍵因素。由于甘南州高寒草甸面積大、分布廣,其FVC的水平直接影響甘南高寒草甸的整體水平,是導致高寒草甸FVC較高的一個重要原因。加之上世紀60年代開始的全球氣候變暖為高寒草甸植被生長提供了良好的溫度條件,由于溫度和水分均處于較適宜的條件,所以高寒草甸FVC呈現逐年增加的變化趨勢。另外,高寒草甸FVC的增加與我國生態恢復工程的不斷實施也有密切的關系。如1999年和2003年相繼開始實施的退耕還林還草和退牧還草工程,對我國草地植被恢復和牧草生長具有重要的促進作用[46],使得2002—2007年間甘肅省高寒草甸FVC呈現上升趨勢。同樣,2010和2015年相繼實施的二輪生態補獎政策對草地生態保護起到了重要作用,也促進了草地FVC的穩定增長[46]。甘肅省高寒草甸區位于青藏高原東段,是兩輪生態補獎中的重點補償區域,2010年以后高寒草甸FVC出現拐點,呈上升趨勢,尤其是2015—2019年高寒草甸FVC呈現明顯的上升趨勢,這與生態工程的實施有著密不可分的關系。

但對高寒草甸這一特殊地理環境下的獨特植被研究較少,對于該草地類型的時空變化的具體原因,有待進一步研究。

4 結論

本研究基于MODIS數據和高寒草甸FVC實測數據,采用Person相關性分析后發現DVI,NDVI,EVI,SAVI和MSAVI與高寒草甸FVC具有較高相關性。NDVI的二次多項式模型:y=—0.65x2+1.97x—0.23能較好的反演高寒草甸FVC,模型精度較高。2000—2019年,甘肅省高寒草甸年均FVC時間上呈現波動上升的趨勢,2002—2007年、2008—2013年和2015—2019年這3個時間段內FVC均呈現上升趨勢,大約每5年會出現明顯波動。空間上高寒草甸FVC由東南向西北逐漸遞減,其中FVC較高區域主要分布在甘南藏族自治州和天祝藏族自治縣,且具有較強的穩定性;FVC較低的主要集中于研究區的西北部且穩定性較差。

猜你喜歡
模型研究
一半模型
FMS與YBT相關性的實證研究
2020年國內翻譯研究述評
遼代千人邑研究述論
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
視錯覺在平面設計中的應用與研究
科技傳播(2019年22期)2020-01-14 03:06:54
EMA伺服控制系統研究
新版C-NCAP側面碰撞假人損傷研究
3D打印中的模型分割與打包
主站蜘蛛池模板: 欧美日韩亚洲国产主播第一区| 中文字幕乱码二三区免费| 日韩精品亚洲人旧成在线| 67194亚洲无码| 日本黄色a视频| 日韩精品成人在线| 国产精品3p视频| 狠狠色丁香婷婷| 成人精品午夜福利在线播放| 日本亚洲成高清一区二区三区| 高潮爽到爆的喷水女主播视频| 久久99精品久久久久久不卡| 色天天综合久久久久综合片| 538国产在线| 伊人久久综在合线亚洲2019| 国产男女免费完整版视频| 亚洲视频四区| 福利在线不卡一区| 亚洲VA中文字幕| 国内熟女少妇一线天| 国产精品视频观看裸模| 一本大道东京热无码av | 亚洲黄色成人| 国产尹人香蕉综合在线电影| 久久午夜夜伦鲁鲁片不卡| 久久一本精品久久久ー99| 99热这里只有精品在线观看| 久久精品66| 国产精品任我爽爆在线播放6080| 国产精品香蕉在线| 久久性视频| 尤物精品视频一区二区三区| 美女视频黄又黄又免费高清| 日韩小视频在线播放| 国产午夜不卡| 日韩精品成人在线| 国产精品免费入口视频| 国产精品第三页在线看| 久久久久国产一级毛片高清板| 午夜不卡福利| 欧美 国产 人人视频| 欧美不卡视频一区发布| 免费va国产在线观看| 久久久久久久97| 人禽伦免费交视频网页播放| 久久人午夜亚洲精品无码区| 亚洲精品在线影院| 中文字幕亚洲第一| 萌白酱国产一区二区| 国产成人综合日韩精品无码不卡| 99在线观看视频免费| 久久综合色视频| 亚洲中文字幕日产无码2021| 国产va在线| 国产精品无码AⅤ在线观看播放| 亚洲侵犯无码网址在线观看| 成人午夜在线播放| 亚洲精品综合一二三区在线| 国产制服丝袜91在线| 久久综合伊人 六十路| 午夜在线不卡| 波多野结衣久久精品| 91丝袜在线观看| 好紧太爽了视频免费无码| 日本高清成本人视频一区| 国产高清不卡视频| 中文字幕永久视频| 青青草久久伊人| 亚洲欧美人成人让影院| 人妻出轨无码中文一区二区| 亚洲色图另类| 成人午夜亚洲影视在线观看| 中文字幕永久视频| 国产理论一区| 国产一区亚洲一区| 亚洲第一成年免费网站| 青青草一区二区免费精品| 综合久久五月天| 国产高清国内精品福利| 波多野结衣亚洲一区| 中文成人在线视频| 小13箩利洗澡无码视频免费网站|