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

基于Landsat 8 OLI影像的延河流域土壤線提取及其應用研究

2017-03-27 09:52:56劉詠梅李曉煥趙牡丹
水土保持通報 2017年1期

王 玲, 劉詠梅, 常 偉, 李曉煥, 凱 楠, 趙牡丹

(西北大學 城市與環境學院, 陜西 西安 710127)

基于Landsat 8 OLI影像的延河流域土壤線提取及其應用研究

王 玲, 劉詠梅, 常 偉, 李曉煥, 凱 楠, 趙牡丹

(西北大學 城市與環境學院, 陜西 西安 710127)

[目的] 以延河流域為研究區,利用自動提取算法建立典型黃土的土壤線方程,為土壤調節植被指數的計算提供基本參數。[方法] 基于Landsat 8 OLI遙感影像,采用自動提取算法獲取土壤線參數;分析歸一化植被指數NDVI,土壤調節植被指數PVI,TSAVI,ATSAVI與實測蓋度的相關性,探討所構建的土壤線方程在黃土高原地區植被指數提取中的適用性。[結果] 通過自動算法與常規方法對比發現兩者偏差較小,且自動提取算法具有較高的精度和穩定性;各植被指數與實測蓋度的相關性大小為:PVI>NDVI>TSAVI>ATSAVI,PVI為延河流域植被蓋度反演的最優植被指數,NDVI次之,TSAVI與ATSAVI較差;與NDVI指數相比,PVI指數能夠較好地抵抗土壤噪聲的影響,對不同植被類型的敏感性較高,更適用于植被覆蓋度較低的黃土高原。[結論] 自動提取算法對延河流域土壤線的提取較為適用,所得參數適合于計算土壤調節植被指數。

土壤; 植被; 土壤線; 植被指數

文獻參數: 王玲, 劉詠梅, 常偉, 等.基于Landsat 8 OLI影像的延河流域土壤線提取及其應用研究[J].水土保持通報,2017,37(1):161-165.DOI:10.13961/j.cnki.stbctb.2017.01.029; Wang Ling, Liu Yongmei, Chang Wei, et al. Extraction and application of soil line in Yanhe River basin based on Landsat 8 OLI[J]. Bulletin of Soil and Water Conservation, 2017,37(1):161-165.DOI:10.13961/j.cnki.stbctb.2017.01.029

土壤線是指土壤的反射率或亮度值在可見光紅、近紅外波段之間的線性關系[1-3]。土壤的光譜特征受土壤有機質、氧化鐵、粗糙度等眾多因素的影響,但是對于同一種土壤來說,其反射光譜在紅與近紅外波段形成了特定的土壤線[4-6]。土壤線不僅反映土壤的光學特性,有助于了解土壤和植被的理化性質與生態特征,并且參與植被指數計算,消除觀測場土壤背景的影響,如垂直植被指數(perpendicular vegetation index, PVI)[3]、土壤調整植被指數(soil-adjusted vegetation index, SAVI)[7]和轉換型土壤調整植被指數(transformed soil-adjusted vegetation index, TSAVI)[8]等均基于土壤線方程的參數進行計算。另外,土壤線還可用于表征土壤的不同濕度狀況,詹志明等[9]基于土壤線斜率提出的土壤干旱指數(perpendicular drought index, PDI),可以有效監測土壤水分。

常規提取土壤線的方法,是利用地面實測土壤光譜、影像解譯及其他方法從遙感影像上識別出裸土像元,進行土壤線參數的計算。劉煥軍等[5]利用室內測定的土壤光譜反射率計算得到土壤線,并在研究東北黑土帶土壤線變異規律時,使用比值植被指數確定裸土像元的方法提取黑土帶不同分區的土壤線參數[10];武婕等[4]采用遙感圖像解譯裁剪出裸土像元的方法提取山東典型土壤的土壤線。常規方法較為復雜,受人為因素影響較多,并且對于中低分辨率的遙感影像來說,裸土的識別難度較大而難以實行。因此,越來越多的學者將研究重點轉向土壤線的自動提取,土壤線自動提取是通過一系列算法提取圖像上的裸露像元,如Fox等[11]和秦其明等[12]提出的土壤線自動提取算法;徐丹丹等[13]在提取加拿大混合草地的土壤線時采用的(R,NIRmin)法和分數位回歸法;楊國范等[14]在秦其明的算法基礎上進行改進提出的修正土壤線自動提取算法等。自動算法在保證模型精度的情況下,使得土壤線的提取易于實現。

延河流域作為黃土丘陵溝壑區的第二副區,其土壤類型在該地區具有代表性。準確提取典型黃土的土壤線,將為各種土壤調節植被指數提供基本參數,對于黃土高原地區的水土流失監測和植被生態建設有重要意義。由于土壤類型復雜多樣,目前對于土壤線的研究仍具有一定的地域限制,國內對于土壤線的研究集中于東北黑土區[5,10]、山東省[4]、寧夏[12]與遼寧[14]等省份,而對陜北黃土高原典型黃土的土壤線研究甚少。其次,目前自動提取算法的應用未得到廣泛推廣,準確提取土壤線參數仍具有一定難度。本文擬以延河流域為研究區,基于Landsat 8 OLI遙感影像,采用自動提取算法獲取典型黃土土壤線參數;通過歸一化植被指數NDVI,土壤調節植被指數PVI,TSAVI,ATSAVI與實測蓋度的相關性分析,探討所構建的土壤線方程在黃土高原地區植被指數提取中的適用性,以期為區域植被覆蓋度的遙感估算提供基本參數。

1 研究區與數據預處理

1.1 研究區概況

延河是黃河的一級支流,發源地為靖邊縣天賜灣鄉周山,主要支流有杏子河、西川、蟠龍川和南川等,由西北向東南注入黃河。延河流域位于陜北黃土高原丘陵溝壑區,東經108°45′—110°28′,北緯36°23′—37°17′,海拔495~1 795 m,流域面積7 725 km2,年平均氣溫9 ℃,年平均降水量520 mm,屬典型的溫暖半干旱大陸性氣候。流域內黃土丘陵溝壑占總面積的90%,主要土壤類型為粉沙質黃綿土,占土地面積的85%以上。由于自然條件和人為活動的原因,該地區土壤侵蝕十分嚴重。植被分布具有明顯的地帶性規律,從南向北依次為森林區、森林草原區和草原區[15]。

1.2 數據源及其預處理

1.2.1 遙感數據 黃土高原在冬季植被覆蓋度低,裸露地面較多,故研究中選取覆蓋延河流域的2015年1月22日127/34,127/35與2014年12月30日126/35軌道上的3景Landsat OLI遙感影像(含云量分別為0.11,0.07與6.28)提取土壤線。選取2014年8月15日127/34,127/35與2014年8月24日126/35軌道上的3景Landsat OLI遙感影像(含云量分別為0.08,0.54與2.77)用于計算研究區各種植被指數,所使用的投影類型為UTM,橢球體為WGS84。

遙感影像的預處理主要包括以下幾個部分:① 對冬、夏兩期遙感影像分別進行輻射定標,將DN值轉換為輻亮度值; ② 對定標后的影像進行Flaash大氣校正,得到地表真實反射率值; ③ 進行遙感影像拼接并利用延河流域邊界矢量文件裁剪出研究區;④ 對野外實測點與遙感影像進行配準,以確保實測樣點位置的準確性。

1.2.2 實測數據 為了研究各種土壤調節植被指數反演植被蓋度的適用性,2012—2014年8月在研究區進行了群落蓋度實地調查,選取典型植被群落樣方54個,主要包括刺槐林、楊樹林、沙棘灌叢、檸條灌叢、退耕草地等。調查內容包括主要植被類型、群落總蓋度,樣地經緯度坐標、坡度、坡向、海拔高度等。為保證群落總蓋度測量的準確性,由多名野外調查人員目測估算蓋度,取平均值作為樣地群落的總蓋度。

2 研究方法

2.1 土壤線提取方法

根據前人研究成果,本文使用歸一化水體指數與閾值相結合的方法去除水體像元,再結合基于二維光譜特征空間的(R,NIRmin)法提取土壤線,其中R,NIRmin分別指紅光波段反射率值、近紅外波段反射率的最小值。(R,NIRmin)法的核心思想是:在NIR-RED空間中,由于水體與植被像元明顯地分列于土壤線兩側[16-17],且裸土像元與植被覆蓋像元整體呈典型三角形分布[12],剔除水體像元的影響后,土壤像元多分布于散點圖的下方,即橫坐標所對應縱坐標的最小值,據此找出土壤像元,計算土壤線參數。本文提取土壤線過程分為以下2個部分:

(1) 水體像元的剔除。(R,NIRmin)自動提取算法的精度受水體影響較大,選擇合適的水體提取方法尤為重要。常用的水體提取方法有:單波段閾值法、多波段譜間關系法和水體指數法等[18-20]。程磊等[19]指出各種水體識別方法均不能提取滿意的水體信息,反映了黃土高原地區水體識別的復雜性和下墊面的特殊性,給水體識別帶來了很大的不確定性,相比這3種方法,歸一化水體指數(normalized difference water index, NDWI)[21]和閾值法的綜合應用能夠有效去除山體陰影的影響,提取的水體信息更為真實準確。故本文選取NDWI指數與閾值相結合的方法提取水體并進行掩膜處理,NDWI的表達式如下:

NDWI=(ρGreen-ρNIR)/(ρGreen+ρNIR)

(1)

式中:ρGreen和ρNIR——影像綠波段與近紅外波段的反射率值。通過對水體像元的NDWI值進行統計,最終確定的NDWI閾值范圍大于-0.13。

(2) 土壤線參數的獲取。本文提取土壤線參數的獲取流程如圖1所示。

該過程主要在Envi 5.2環境中實現,首先,以剔除水體后的延河流域OLI影像波段4為橫軸,波段5為縱軸繪制二維散點圖;其次,使用ROI工具將散點圖中分布于較下方的散點篩選出來,并將其以感興趣區的形式輸出在影像上;然后,將ROI導出為CSV格式文件并用Excel等工具打開,該文件中含有每個ROI像元在遙感影像上對應的所有波段信息。由于遙感影像的像元數過多,尋找所有土壤點的效率較低,為了保證土壤點集的準確度,將NIR-RED空間的橫坐標按照反射率為0.002 5的間隔分為200組,分別求取每個間隔的反射率值所對應的波段5的最小值,得到初始土壤點集;接下來,為了數據的代表性與方程的穩定性,分別計算以橫坐標總區間范圍的0%~50%,0%~100%,25%~100%,50%~100%,0%~75%與25%~75%為區間子集的土壤點集的相關系數[12],對初始土壤點集進行篩選,取相關系數最大的為最佳橫坐標區間,得到最終的土壤點集;最后,采用最小二乘法擬合方程,得到土壤線參數。

圖1 土壤線提取流程

2.2 土壤調節植被指數

為了驗證所提取的土壤線參數應用于土壤調節植被指數計算的有效性,選用了以下植被指數:PVI,TSAVI與ATSAVI(改進轉換型土壤調整植被指數,adjusted transformed soil-adjusted vegetation index)[8],分析PVI,TSAVI,ATSAVI以及應用最為廣泛的植被指數NDVI與實測蓋度之間的相關性,并通過與NDVI的對比探討研究區最適宜的植被指數。其中,各植被指數的計算公式如下所示[13,22]。

(2)

(3)

(4)

(5)

式中:R,NIR——影像紅光波段、近紅外波段的反射率;a,b——土壤線的斜率、截距;ATSAVI計算公式中X=0.08。

3 結果與分析

3.1 土壤線提取結果

采用自動算法獲取土壤線,并將其加載在延河流域OLI影像所構成的NIR-RED光譜特征空間中,結果如圖2所示,該結果表明Landsat 8 OLI的第4,5波段反射率之間具有較好的線性關系,提取出的土壤線在NIR-RED光譜空間中的實際位置與理論上的相符。

圖2 自動算法提取土壤線擬合結果

目前,有學者在計算延河流域土壤調節植被指數時,采用純凈像元指數法確定土壤線參數[16]。為了檢驗自動提取算法的可靠性,本文在使用自動提取算法的同時,在Envi 5.2中對影像進行MNF變換,計算純凈像元指數,在N維可視化窗口中選取純凈裸土像元,并在NIR-RED空間中進行回歸分析擬合出土壤線,結果如圖3所示。將自動提取算法與常規方法得到的土壤線參數進行對比,兩者斜率相差0.067,截距相差0.061,由此可看出,2種方法提取的土壤線參數相差較小,且前者方程的擬合度較高,說明自動提取算法對延河流域土壤線的提取較為適用。

圖3 常規方法獲取土壤線結果

3.2 土壤線參數應用

3.2.1 植被指數與實測蓋度的相關性分析 利用自動提取和常規方法得到的土壤線參數計算垂直植被指數PVI,轉換型土壤調整植被指數TSAVI,改進轉換型土壤調整植被指數ATSAVI,歸一化差值植被指數NDVI,提取實測樣方蓋度對應的各植被指數值,計算各植被指數與實測蓋度之間的相關性,得到的結果如表1所示。

從2種不同方法所計算的相關性結果可以看出,采用自動算法計算的植被指數與實測蓋度的相關性在一定程度上有所提高,說明自動提取算法能夠提高基于土壤線參數的植被指數對植被覆蓋度的敏感性。這主要是因為自動算法能夠為植被指數提供更加準確和穩定的土壤線參數,消除不確定因素的影響。

表1 植被指數與實測蓋度的相關性

注:置信度為95%; 數據檢驗的有效性sig.<0.05; 上述回歸分析均顯著相關。

從各植被指數與實測蓋度的相關性結果可以看出,NDVI,PVI,TSAVI,ATSAVI與實測蓋度之間存在顯著相關性,均可以用于反演植被覆蓋度。各植被指數與實測蓋度的相關性大小依次為:PVI>NDVI>TSAVI>ATSAVI,說明PVI為該區域最優植被指數,其次是NDVI,TSAVI與ATSAVI最差。黃土高原地區植被覆蓋度較低,NDVI等指數受土壤背景的影響比較強烈,而PVI充分考慮研究區土壤背景的影響,引入土壤線參數提高了植被指數表征植被覆蓋度的精度。

3.2.2 PVI與NDVI指數的對比分析

(1) 不同植被指數抗土壤噪聲的能力對比。植被信號—土壤噪聲比(S/N)[23]用以表征不同植被指數對土壤噪聲的敏感程度,S/N值越大,代表植被指數削弱土壤背景影響的效果越好,植被信號越突出。S/N的計算公式為:

(6)

為了研究PVI與NDVI在不同蓋度下抵抗土壤噪聲的能力,分別以PVI和NDVI為基礎,采用像元二分模型估算研究區植被覆蓋度(見公式7),其中PVI,NDVI的置信區間取5%~95%;然后,按蓋度從低到高依次選取若干個以30×30像元大小窗口作為基本研究樣區,且要求每個樣區內的蓋度均勻;計算每個樣區內PVI與NDVI指數的S/N值。結果如表2所示。

(7)

式中:Fc——植被覆蓋度;NDVImax,NDVImin——NDVI指數的最大值、最小值;PVImax,PVImin——PVI指數的最大值、最小值。

表2 不同蓋度下植被指數的S/N值

表2顯示,在蓋度大約為45%時,PVI與NDVI的S/N值大小關系開始發生變化。當蓋度小于45%時,PVI的S/N值大于NDVI的S/N值;當蓋度大于45%時,情況剛好相反。這一結果說明在蓋度小于45%時,PVI指數抗土壤噪聲的能力較NDVI的強,更能突出植被信號;而在蓋度大于45%時,NDVI指數所表征的植被信號比PVI指數強。上述結果說明,PVI指數對低植被覆蓋區植被信息的提取效果較好,NDVI指數則更好反映了高植被覆蓋區的植被信息,這一結論與池宏康[24]關于黃土高原地區植被信息提取方法的研究結果相一致。

(2) 變異系數對比。標準差與平均數的比值稱為變異系數,變異系數反映了數據的相對分散程度。變異系數越大,數據離散程度越大,對不同植被類型特征的敏感性越高。延河流域PVI,NDVI指數的統計特征值與變異系數如表3所示。

表3 各植被指數的特征值統計

表3顯示,PVI指數的變異系數比NDVI大,表明PVI對于不同植被類型的敏感度相對較高。結合黃土高原植被覆蓋低、大部分地區受土壤背景影響較大的特點來說,PVI比NDVI更適用于該區域的植被覆蓋度反演。

4 討論與結論

已發表的土壤線參數變化范圍:斜率為1.06~1.6,截距為-0.01~0.07(與斜率的變化趨勢相反)[6],本文得到的土壤線參數均在以上范圍內,且與常規方法相比,自動提取算法能夠在一定程度上提高植被指數對植被蓋度的敏感性,說明了利用自動提取算法提取延河流域土壤線方程的適用性。

將提取的土壤線參數應用于土壤調節植被指數的計算,各種植被指數與實測蓋度的相關性大小依次為:PVI>NDVI>TSAVI>ATSAVI,PVI為延河流域反演植被蓋度的最優植被指數,NDVI次之,TSAVI與ATSAVI較差。通過對NDVI與PVI指數的對比分析可知,PVI指數能夠更好地抵抗土壤噪聲的影響,對不同植被類型的敏感性更高,對低植被覆蓋區植被信息的提取效果較好。采用自動方法提取的典型黃土土壤線參數適用土壤調節植被指數的計算,因而具有實際的應用意義。

土壤線提取的關鍵是裸土像元的識別,在使用自動提取算法時,裸土識別的準確度受水體的影響較大,黃土高原水體識別的復雜性和不確定性,為準確提取土壤線參數帶來一定程度的影響,因此,后續研究中需進一步提高黃土高原水體識別的準確度,以提高自動提取算法的效果。

[1] Baret F, Jacqemoud S, Hanocq J B. About the soil line concept in remote sensing[J]. Advances in Space Research, 1993,13(5):281-284.

[2] Kauth R J, Thomas G S. The tasseled cap: A graphic description of the spectral-temporal development of agricultural crops as seen by LANDSAT[C]∥Proceeding of the Symposium on Machine Processing, Indiana, 1976:41-51.

[3] Richardson A J, Wiegand C L. Distinguishing Vegetation from Soil Background Information[J]. Photogrammetric Engineering & Remote Sensing, 1978,43(12):1541-1552.

[4] 武婕,李玉環,李增兵.基于遙感影像的土壤線空間變異規律及影響因素分析[J].自然資源學報,2014,29(4):702-708.

[5] 劉煥軍,張柏,宋開山,等.基于室內光譜反射率的土壤線影響因素分析[J].遙感學報,2008,12(1):119-127.

[6] Galv?o L S,caro V. Variability of laboratory measured soil lines of soils from Southeastern Brazil[J]. Remote Sensing of Environment, 1998,63(2):166-181.

[7] Huete A R. A soil adjusted vegetation index(SAVI)[J]. Remote Sensing of the Environment, 1988,25(3):295-309.

[8] Baret F, Guyot G, Major D J. TSAVI: A vegetation index which minimizes soil brightness effects on LAI and APAR estimation[C]∥Geoscience and Remote Sensing Symposium, 1989. IGARSS’89. 12th Canadian Symposium on Remote Sensing, 1989 International. IEEE, 1989:1355-1358.

[9] 詹志明,秦其明,阿布都瓦斯提·吾拉木,等.基于NIR-Red光譜特征空間的土壤水分監測新方法[J].中國科學:D輯,2006,36(11):1020-1026.

[10] 劉煥軍,宇萬太,張新樂,等.中國東北黑土帶土壤線空間變異規律[J].農業工程學報,2009,25(10):166-170.

[11] Fox G A, Sabbagh G J, Searcy S W. An automated soil line identification routine for remotely sensed images[J]. Soil Science Society of America Journal, 2004,68(4):1326-1331.

[12] 秦其明,游林,趙越,等.基于二維光譜特征空間的土壤線自動提取算法[J].農業工程學報,2012,28(3):167-171.

[13] XU Dandan, GUO Xulin. A Study of Soil Line Simulation from Landsat Images in mixed grassland[J]. Remote Sensing, 2013,5(9):4533-4550.

[14] 楊國范,翟光耀,張婷婷,等.基于環境小衛星的土壤線提取研究[J].遙感技術與應用,2014,29(3):428-432.

[15] 朱樂天,焦峰,劉源鑫,等.黃土丘陵區不同土地利用類型土壤水分時空變異特征[J].水土保持研究,2011,18(6):115-118.

[16] 趙婷.陜北延河流域結構性植被蓋度遙感模型研究[D].西安:西北大學,2013.

[17] 徐彬彬.我國土壤光譜線之研究[J].環境遙感,1991,6(1):61-71.

[18] 畢海蕓,王思遠,曾江源,等.基于TM影像的幾種常用水體提取方法的比較和分析[J].遙感信息,2012,27(5):77-82.

[19] 程磊,徐宗學,左德鵬,等.基于Landsat TM數據的黃土高原區水體識別方法研究[J].北京師范大學學報:自然科學版,2010,46(3):424-430.

[20] 王偉武,朱霞,孫躍池,等.基于ETM圖像的山地水體提取方法研究[J].系統仿真學報,2013,25(9):2196-2200.

[21] Mcfeeters S K. The use of the normalized difference water index(NDWI)in the delineation of open water features[J]. International Journal of Remote Sensing, 1996,17(7):1425-1432.

[22] 田慶久,閔祥軍.植被指數研究進展[J].地球科學進展,1998,13(4):327-333.

[23] Qi J, Chehbouni A, Huete A R, et al. A modified soil adjusted vegetation index[J]. Remote Sensing of Environment, 1994,8(2):119-126.

[24] 池宏康.黃土高原地區提取植被信息方法的研究[J].植物學報,1996,38(1):40-44.

Extraction and Application of Soil Line in Yanhe River Basin Based on Landsat 8 OLI

WANG Ling, LIU Yongmei, CHANG Wei, LI Xiaohuan, KAI Nan, ZHAO Mudan

(DepartmentofUrbanandEnvironmentalScience,NorthwestUniversity,Xi’an,Shaanxi710127,China)

[Objective] Taking Yanhe River basin as a study area, typical loess soil line equation was established by automatic algorithm to provide basic parameters for the calculation of the soil adjusted vegetation index. [Methods] Based on Landsat 8 OLI images, soil line was extracted by automatic algorithm, then the correlation coefficients between normalized difference vegetation index(NDVI), perpendicular vegetation index(PVI), transformed soil-adjusted vegetation index(TSAVI), adjusted transformed soil-adjusted vegetation index(ATSAVI) and test coverage were calculated; and the applicability of soil line equation in vegetation index extraction was also discussed in the Loess Plateau. [Results] (1) Compared with conventional method, the difference was small, and the automatic algorithm had high accuracy and stability; (2) The sequence of the correlation between vegetation index and test coverage was: PVI>NDVI>TSAVI>ATSAVI, which showed that PVI was the optimal vegetation index, NDVI was the next, TSAVI and ATSAVI were the worst to the extraction of vegetation coverage in Yanhe river basin; (3) Compared with NDVI, PVI index could reduce the influence of soil noise, and its sensitivity to different vegetation types was higher. PVI was more suitable for the Loess Plateau with low vegetation coverage. [Conclusion] It is suitable to extract the soil line by automatic extraction algorithm in Yanhe river basin, and the parameters can be used to calculate the soil-adjusted vegetation index.

soil; vegetation; soil line; vegetation index

2016-06-21

2016-07-14

國家自然科學基金項目“面向土壤侵蝕評價的結構性植被蓋度遙感模型研究”(41171225), “梯田對坡度坡長因子的擾動特征研究”(41271284)

王玲(1989—),女(漢族),陜西省商洛市人,碩士研究生,主要從事遙感技術與應用方面的研究。E-mail:qqzjwl123@163.com。

劉詠梅(1970—),女(漢族),陜西省西安市人,博士,副教授,主要從事區域水土流失遙感測方面的研究。E-mail:liuym@nwu.edu.cn。

B

1000-288X(2017)01-0161-05

TP79, P641.6

主站蜘蛛池模板: 亚洲色图欧美激情| 亚洲国产中文在线二区三区免| 亚洲人成网站色7777| 91在线一9|永久视频在线| 尤物国产在线| 亚洲熟女中文字幕男人总站 | 久久久久人妻精品一区三寸蜜桃| 曰韩免费无码AV一区二区| 亚洲精品无码久久久久苍井空| 人人爱天天做夜夜爽| 91国内视频在线观看| 色综合久久久久8天国| av尤物免费在线观看| 婷婷五月在线| 色香蕉网站| 亚洲无线国产观看| 久久综合五月婷婷| 91欧美在线| 波多野结衣一二三| 国产成年无码AⅤ片在线| 亚洲一级毛片免费观看| www.日韩三级| 国产三级国产精品国产普男人 | 免费a级毛片18以上观看精品| 亚洲中文在线看视频一区| 欧美成在线视频| 国产一级妓女av网站| 91免费国产在线观看尤物| 1级黄色毛片| 国产女人在线观看| 伊人狠狠丁香婷婷综合色| 97综合久久| 国产欧美网站| 亚洲无码免费黄色网址| 久久综合丝袜日本网| 青青久在线视频免费观看| 国产理论精品| 亚洲品质国产精品无码| 欧美午夜视频在线| 亚洲国产中文精品va在线播放| www.91中文字幕| 久久久久久久97| 亚洲天堂区| 日韩国产欧美精品在线| 久久综合AV免费观看| 欧美日韩国产高清一区二区三区| 97国产成人无码精品久久久| 日本在线国产| 国产午夜一级毛片| 丁香婷婷久久| 三上悠亚精品二区在线观看| a级毛片在线免费| 亚洲av无码人妻| 国产成人一区在线播放| 性欧美在线| 国产成人禁片在线观看| 色悠久久久久久久综合网伊人| 999精品在线视频| 欧美97欧美综合色伦图| 国产国产人成免费视频77777| 最新加勒比隔壁人妻| 欧美黄网站免费观看| 日本道中文字幕久久一区| 人妻精品久久无码区| 国产福利影院在线观看| 精品黑人一区二区三区| 亚洲欧美日韩成人高清在线一区| 少妇人妻无码首页| 国产乱论视频| 国产一级视频在线观看网站| 日本欧美成人免费| 国产在线精品99一区不卡| 久久伊人操| 中文字幕av无码不卡免费| 国产 日韩 欧美 第二页| 激情成人综合网| 国产成人av大片在线播放| 韩国v欧美v亚洲v日本v| 性色生活片在线观看| 亚洲福利片无码最新在线播放| a亚洲天堂| 国产福利在线免费|