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

基于不同植被指數和分類回歸樹的關中地區土地利用分類

2018-04-11 07:12:29馬孝義邢旭光
水土保持研究 2018年3期
關鍵詞:耕地分類

張 盼, 馬孝義, 趙 龍, 邢旭光

(1.西北農林科技大學 水利與建筑工程學院, 陜西 楊陵712100; 2.西北農林科技大學 機械與電子工程學院, 陜西 楊陵712100)

土地利用是人類對土地進行的使用、保護和改造活動,反映了人類社會的發展,影響全球生態環境的變化。土地利用的分類研究對土地資源利用結構的規劃和調整、生態環境的平衡和保護有著重要的作用。近年來,遙感技術以其時效高、范圍廣和成本低等優點而廣泛地應用于對地觀測活動中,對及時準確地獲取土地利用信息具有重要的意義[1-3]。

目前遙感數據具有不同時間、空間和光譜分辨率的產品。但是空間高分辨率的數據,由于回訪周期長,對作物關鍵物候期的描述不足,且幅寬小,限制了數據在大范圍中的應用。MODIS數據為中空間分辨率,但是由于其回訪周期短,幅寬大,光譜波段范圍廣而廣泛運用于較大尺度地區的土地利用信息的獲取[4-5]。Muchoney等[6]和劉建光等[5]利用MODIS影像分別對美國中部和北京市進行土地覆被分類;許青云等[7]基于陜西省多年時序的MODIS數有效地識別大尺度農作物信息。植被指數可以反映植被的生長狀況,是植被生長狀態和空間分布密度的指示因子[8-9],具有時間序列的植被指數可以反映植被的物候特征,常常作為不同植被和農作物分類的依據。Wardlow等[4]和Damien等[10]分別MODIS/NDVI和MODIS/EVI時間序列進行作物分布信息的提取;何超英等[11]指出在單季相Landsat影像分類時引入MODIS/NDVI時間序列能有效地解決耕地與其他類型混分的問題;白文龍[12]基于逐旬的NDVI數據采用最大似然法實現了關中地區土地利用分類;周玉潔等[13]基于MODIS/EVI時間序列數據,基于諧波分析法和線性光譜模型對關中地區的耕地信息進行提??;賈明明等[14]綜合環境衛星和MODIS NDVI時序數據實現了雙臺子河口土地覆蓋分類。但目前不同植被指數對土地利用分類的影響研究較少,左麗君等[15]基于MODIS NDVI和EVI數據在河西走廊綠洲東部樣區進行耕地提取試驗,表明EVI較NDVI有更強的識別能力。

為了比較不同植被指數時間序列對土地利用分類精度的影響,以及組合植被指數和單一植被指數時間序列分類精度的差異,本研究以關中地區為研究區域,選取應用較多的歸一化植被指數(Normalized Difference Vegetation Index,NDVI)和增強型植被指數(Enhanced Vegetation Index,EVI),基于MODIS植被指數產品,采用迭代自組織分析算法與分類回歸樹算法結合的方法,實現關中地區土地利用分類,并進行空間精度驗證和定量精度驗證,為大尺度地區土地利用分類及農業信息提取提供一定的參考。

1 研究區概況及數據來源

1.1 研究區概況

關中地區地處33°35′—35°52′N,106°18′—110°38′E,位于陜西省中部,包括西安、咸陽、寶雞、渭南和銅川五市,東西長約390 km,南北寬約200 km,海拔227~3 772 m,關中北部位于渭北高原南部,南倚秦嶺山脈,渭河從中穿過,形成沖積平原。屬溫帶季風性氣候,年均溫6~13℃,年降水量500~800 mm,其中6—9月份占60%,多為短時暴雨,冬春降水較少,春旱、伏旱頻繁。研究區種植獼猴桃、蘋果、梨等果樹,小麥、玉米等作物,低海拔地區的渭河階地為一年兩熟制,北部黃土臺塬區和南部秦嶺山地由于海拔較高,灌溉條件差,為一年一熟制[12-13]。

圖1 研究區位置及高程

1.2 數據來源

1.2.1植被指數數據集本研究采用NASA USGS提供的植被指數產品MOD13Q1(https:∥lpdaac.usgs.gov),空間分辨率為250 m,時間分辨率為16 d。選取2014年10月1日至2016年3月15日共33個時相的遙感數據,構建NDVI和EVI時間序列曲線。

1.2.2訓練及驗證數據集利用Google Earth數據,以目視解譯為主,結合實地考察的結果,選擇具有代表性的樣本作為感興趣區域,盡最大可能覆蓋不同果樹及不同作物熟制的耕地,將其中一部分作為土地利用分類的訓練數據集,一部分作為空間精度的驗證數據集。同時,選擇由陜西農業網(http:∥www.sxny.gov.cn/)農業統計數據庫提供的耕地和果園統計信息作為定量驗證數據集。

2 研究方法

2.1 數據預處理

利用MODIS產品批處理工具MRT軟件對MOD13Q1數據進行空間拼接、格式和投影轉換的處理,輸出坐標系為WGS_1984_UTM_Zone_49 N,并根據研究區矢量圖進行裁剪,提取歸一化植被指數NDVI和增強型植被指數EVI波段,構建植被指數時間序列。

MOD13Q1雖然糾正了分子和氣溶膠的影響,但由于云雪覆蓋等,仍存在一些噪聲,在ENVI/IDL平臺中采用Savitzky-Golay濾波法,對植被指數時間序列進行重構。經多次試驗,設置平滑窗口半徑為3,多項式次數為2效果較好,可準確描述波形和波峰的值,能夠滿足土地利用分類的要求。由于S-G濾波法存在邊緣效應,且考慮到作物的生育期,在2015全年數據的前后分別添加5個時相的數據參與平滑計算,實現時間序列數據的降噪、平滑和重構。

2.2 土地分類系統

根據《土地利用現狀分類》國家標準(GB/T21010—2007)兩級分類系統,結合本研究區的實際情況,確定研究區土地的5類分類系統,具體如表1所示。

2.3 典型地物植被指數時間序列曲線

根據本文確定的研究區分類系統,結合Google Earth 選擇典型地物的感興趣區,并提取植被指數,計算均值,構建典型地物植被指數時間序列曲線見圖2。

由圖2可以看出,建設用地及其他和水體植被指數年內變化不大,且水體的NDVI和EVI小于0.2,由于水體對近紅外波段的反射率很低,對可見光高反射率,植被指數存在負值;建設用地及其他的植被指數通常大于水體,但NDVI值不大于0.4,EVI值不大于0.3。果園和林草地的植被指數隨時間的變化趨勢相似,林草地在第15至20時相,即5月24日至8月12日,NDVI值無明顯變化,同期EVI依然隨時間變化達到峰值后下降;果園在第18至21時相,進入果實生長發育期,NDVI變幅很小,同期EVI值變化微弱。果園的NDVI時間序列和EVI時間序列的峰值都小于林草地,果園分別達0.74,0.52左右,而林草地可達0.86,0.62左右。一年兩熟制的耕地植被指數時間序列兩個峰值,夏糧作物的峰值在第12時相左右,即4月中旬,秋糧作物的峰值在第19時相左右,即7月下旬;一年一熟制的耕地植被指數只有一個峰值;不同作物由于物候期不同而峰值及其出現的時間有一定的偏差。由此可以看出不同地物的植被指數時間序列存在顯著差異。

表1 研究區分類系統

圖2 典型地物植被指數時間序列

為了比較單一NDVI,EVI,及組合植被指數用于土地分類的效果差異,分別采用NDVI,EVI,NDVI+EVI組合以及EVI+NDVI組合的時間序列,用于研究區土地利用分類中。

2.4 分類方法

2.4.1ISODATA算法ISODATA算法(Iterative Selforganizing Data Analysis Techniques Algorithm)即迭代自組織數據分析算法,是一種非監督分類方法,通過計算數據空間中均勻分布的類均值,基于輸入的閾值參數,用最小距離技術將剩余像元迭代聚集,對判別函數的不斷訓練和調整、自動進行類別的合并和分裂,實現地物的聚類,具有自組織、啟發式的特點[16]。

2.4.2CART算法CART(Classification and Regression Tree)是分類回歸決策樹構建算法,由Breiman等[17]于1984年提出并不斷改進。其基本原理是將訓練樣本分為預測變量和目標變量數據集,將基尼系數作為選擇最佳測試變量和分割閾值的準則,循環二分形成二叉樹式的決策樹結構,并采用交叉驗證的方法進行修剪,最終得到一個最優二叉樹,構成地物分類效果最佳的判斷條件組合。該算法對輸入數據沒有任何統計分布的要求,運行速度快、準確性高,且實現簡單,結構清晰易于理解[18-19]。該算法以其處理大量數據和高維數據的有效性,被廣泛用于土地利用分類及提取中[20],且已有研究表明,利用CART算法進行土地利用分類具有較高精度[5]。

本研究將ISODATA方法與CART方法結合對關中地區土地利用進行分類?;?種不同植被指數時間序列,根據分類系統和實際地物情形,設置ISODATA最小分類數為5,最大分類數為10,迭代20次,像元變化閾值5%,獲得研究區土地利用的聚類結果,并將聚類結果與對應的植被指數時間序列合并,結合Goole Earth提取的訓練樣本,利用CART算法構建分類決策樹,獲取分類規則,實現關中地區土地利用分類。

3 結果分析與討論

3.1 分類結果與分析

土地利用分類結果見圖3。

由于本研究采用數據空間分辨率為250 m,且考慮到研究區實際,確定的分類系統僅有5種地類,由于草地多夾雜在林地中,與部分林地難以區分,訓練樣本獲取難度大,將林地與草地劃分為一類。基于不同植被指數時間序列的土地利用分類結果見圖3。渭河從關中穿過,但分類結果并不存在明顯的連續水體,在寶雞東部和咸陽西部幾乎沒有水體顯示,可能是由于渭河河面寬度較小,且遙感數據為中等分辨率,混合像元的存在導致將其錯分為建設用地。分類結果中,建設用地的分布基本與行政中心、城鎮建設用地和農村居民點等分布一致。關中中部地區地勢平坦,土壤肥沃,農業種植業發展,分類結果顯示大面積耕地連續分布,與之相一致,而南部和北部由于海拔較高,地形比較復雜,南部秦嶺山區主要是以林地覆蓋,北部黃土梁峁和臺塬區,由于受地形的限制,耕地破碎,因地制宜地布置林地和果園。4種情形的分類結果基本相同,局部地區耕地和果園面積之間存在差異。在寶雞南部山區、咸陽西北部,渭南東北部以及銅川地區,基于NDVI時間序列的分類結果中,果園面積明顯大于其他幾種情形,主要是由于紅光通道易飽和,且沒有考慮背景土壤噪音的影響,NDVI數據在高植被覆蓋率區域具有較低靈敏度[15,21],造成林草地與果園混分;而寶雞南部山區,基于EVI時間序列的分類結果中,耕地面積明顯大于其他幾種情形,由于隨海拔的上升,植被由茂密的林地逐漸變化為稀疏的植被、草地等,而EVI通過藍光來修正大氣對紅光的影響,增加土壤調節參數,減弱了冠層背景與一年一熟制的耕地和果園混分?;诮M合和土壤變化對植被指數的影響,可以識別該區域植植被指數的分類結果與單一植被指數分布區域大致相似,在數量上存在一定的偏差,需進一步驗證。

圖3 基于不同植被指數時間序列土地分類

3.2 分類精度分析

為了更好地評價土地分類精度,從空間和定量兩個方面進行精度分析與評價,分別利用混淆矩陣和統計數據進行驗證。由于農業是國民經濟的基礎,農業用地的變化對于農業生產及社會穩定具有重要影響,因此對研究區內的果園和耕地進行重點分析。

3.2.1空間精度分析利用Google Earth數據提取感興趣區,作為空間精度分析的驗證樣本,基于ENVI平臺計算混淆矩陣,評價分類精度,具體見表2—3所示.

通過基于單一植被指數分類結果的空間精度比較,不難得出:基于NDVI時間序列分類結果的總體精度和Kappa系數均小于EVI的結果;前者耕地的制圖精度和用戶精度均表現不足,雖然果園的制圖精度為100.00%,事實上并不與實際情況完全相同,分類結果中的部分果園實則為林草地和耕地等其他地類,用戶精度仍然小于EVI的結果。總體來說,基于EVI的分類結果優于NDVI,總體精度較高,EVI對植被的識別能力高于NDVI,具有更高的提取精度??梢?,EVI在覆蓋度比較高的地區具有更高的靈敏度。

通過基于組合植被指數分類結果的空間精度比較,可以看出:基于EVI+NDVI時間序列的分類結果,總體精度和Kappa系數大于NDVI+EVI的結果,整體效果較好,可能是由于組合植被指數增加了地類信息,不同類別之間的判別條件更為豐富,構建的決策樹結點更多,從根結點出發可以通過多條路徑達到達每一種地類的葉結點,實現分類,準確性更高;前者耕地的制圖精度和用戶精度均小于后者,果園的制圖精度小于后者,用戶精度卻大于后者,并不能體現對耕地和果園的提取效果差異,需要進一步利用統計數據進行定量精度評價。

表2 基于單一植被指數土地分類的空間精度比較

注:括號外為基于NDVI時間序列土地分類的混淆矩陣,括號內為基于EVI時間序列土地分類的混淆矩陣。

表3 基于組合植被指數土地分類空間精度比較

注:括號外為基于NDVI+EVI時間序列土地分類的混淆矩陣,括號內為基于EVI+NDVI時間序列土地分類的混淆矩陣。

綜合表2和表3的內容可知:在以上4種分類情形中,總體精度均超過96%,Kappa系數均超過0.94,分類結果可靠,本研究采用的分類方法可行,且精度存在EVI+NDVI>NDVI+EVI >EVI >NDVI的規律,組合植被指數的分類結果優于單一植被指數,但差別不大。

3.2.2定量精度分析利用混淆矩陣評價分類結果的空間精度,局限于感興趣區的選擇,易受人為因素干擾,并不能全面反映其精度,故利用統計數據再進行定量精度評價。由于農業用地是人類重要的食物來源,且受人類活動影響顯著,故對農業用地進行重點分析,提取各地市的農業用地,包括耕地和果園。為避免各個區域分類結果精度差異的影響,不僅比較關中地區提取面積與統計面積的大小,且根據行政區的劃分,比較各地市二者的差異,以提高精度評價的科學性。

表4 耕地提取定量精度評價

對4種情形的分類結果進行西安、銅川、寶雞、咸陽和渭南五市以及關中地區的耕地的提取(表4),在各個地市耕地的提取面積均大于統計數據,可能是混合像元及一年一熟耕地與一年生草本植物和果樹混分存在導致。利用EVI時間序列進行耕地提取時的表現優于NDVI,具有較高的精度,對耕地的具有更高的識別能力;組合植被指數的提取精度,較單一植被指數時間序列的分類結果更高,各個地市雖存在差異,但除銅川地區外,基于EVI+NDVI時間序列的提取精度均為最高。在整個關中地區,基于EVI+NDVI時間序列的提取精度為89.36%,較單一NDVI和EVI時間序列提取精度分別提高11.74%和10.16%,較NDVI+EVI時間序列提取精度提高8.65%。在銅川地區,基于EVI時間序列的提取結果與統計數據一致性最高,達90.36%,組合植被指數在銅川并沒有精度優勢,可能是由于銅川地區地形復雜,土地利用類型較為破碎,植被復雜影響分類精度。

對4種情形的分類結果進行關中地區以及西安、銅川、寶雞、咸陽和渭南五市的果園的提取(表5),各地市果園的提取面積基本上都大于統計數據,緣于果園與林地及秋糧作物在某些生長階段植被指數時間序列曲線具有相似性,造成混分,僅在銅川地區,基于組合植被指數時間序列的提取結果小于統計數據,基于組合植被指數的決策樹可能存在過度擬合,且連續面積小和混合像元的問題而造成漏分。關中地區基于EVI+NDVI時間序列的果園提取精度最高,達65.47%,較基于單一NDVI和EVI時間序列的提取精度分別提高24.69%和17.12%,較NDVI+EVI組合的時間序列提取精度提高0.75%;由于銅川地區的提取面積小于統計數據,剔除銅川果園面積后,基于NDVI+EVI時間序列的提取精度為67.87%,比EVI+NDVI時間序列的提取精度68.62%低0.75%,差異并不明顯。由此說明EVI對果園的識別能力優于NDVI,在進行分類時,組合植被指數的提取精度高于單一植被指數,且EVI優先時精度更高。

表5 果園提取定量精度評價

綜合空間和定量精度分析可知:在4種情形的分類結果中,基于單一EVI時間序列土地利用分類精度和農業用地提取精度優于NDVI;組合植被指數在土地利用分類和農業用地的提取中表現優于單一植被指數,且EVI優先時,精度最高。

3.3 討 論

本研究基于不同植被指數,采用ISODATA與CART結合的方法,實現了關中地區的土地利用分類,結果精度較好,與白文龍[12]采用的最大似然法相比,總體精度提高9%以上,Kappa系數提高0.10以上。同時建立的分類決策樹可以直接用于該地區其他氣候未發生特別重大變化、各地類植被指數變化不顯著的年份的土地利用分類中,為研究年際間土地利用變化提供快速、便捷、有效的分類方法。

基于不同的植被指數時間序列的分類結果精度均較高,但是水體的制圖精度和用戶精度較低,與建設用地及其他存在混分問題,可以考慮結合歸一化水體指數和歸一化建筑指數,以通過增加分類依據來提高其分類精度。雖然MODIS數據幅寬大,可以用于大尺度地區土地利用分類中,具有較高的時間分辨率,但是空間分辨率為250 m,屬中等分辨率遙感數據,混合像元現象顯著,在今后的研究中,可進行混合像元分解提高分類精度,或與其他類型的高空間分辨率遙感數據融合,兼顧幅寬問題的同時又提高空間分辨率,可提高分類的精度。本文采用的MODIS陸地3級產品,時間分辨率為16 d,若采用更高時間分辨率的產品,可獲取更多的物候信息,在進行土地覆被分類或者農作物提取時,提供更多的分類依據。對比關中地區五市的分類精度,銅川地區分類精度較低,可能與土地利用類型的破碎程度有關,在今后的工作中可以進一步探究。

4 結 論

(1) 4種情形分類結果與實際情況較為一致,精度較好,本研究方法適用于關中地區土地利用分類;

(2) EVI引入藍光修正大氣對紅光的影響,解決了紅光通道易飽和的問題,且考慮冠層背景和土壤變化對植被的影響,在進行土地分類時,對植被的識別能力大于NDVI,能較好地反映植被情況;

(3) 利用組合植被指數進行土地利用分類,與單一植被指數相比,構建決策樹時不同地類之間判別依據更多,分類規則準確性更高,且EVI優先時能獲得更高的精度。

參考文獻:

[1]馬玥,姜琦剛,孟治國,等.基于隨機森林算法的農耕區土地利用分類研究[J].農業機械學報,2016,47(1):297-303.

[2]林楠,姜琦剛,楊佳佳,等.基于資源一號02C高分辨率數據的農業區土地利用分類[J].農業機械學報,2015,46(1):278-284.

[3]宋富強,康慕誼,鄭壯麗,等.陜北黃土高原地區土地利用/覆被分類及驗證[J].農業工程學報,2011,27(3):316-324.

[4]Wardlow B D, Egbert S L. Large-area crop mapping using time-series MODIS 250 m NDVI data:An assessment for the U. S. Central Great Plains[J]. Remote Sensing of Environment, 2008,112(3):1096-1116.

[5]劉建光,李紅,孫丹峰,等. MODIS土地利用/覆被多時相多光譜決策樹分類[J].農業工程學報,2010,26(10):312-318.

[6]Muchoney D, Borak J, Chi H, et al. Application of the MODIS global supervised classification model to vegetation and land cover mapping of Central America. [J]. International Journal of Remote Sensing, 2000,21(6/7):1115-1138.

[7]許青云,楊貴軍,龍慧靈,等.基于MODISNDVI多年時序數據的農作物種植識別[J].農業工程學報,2014,30(11):134-144.

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

[9]李鑫川,徐新剛,王紀華,等.基于時間序列環境衛星影像的作物分類識別[J].農業工程學報,2013,29(2):169-176.

[10]Damien Arvor, Milton Jonathan, Vincent Dubreuil, et al. Classification of MODIS EVI time series for crop mapping in the state of Mato Grosso, Brazil[J]. International Journal of Remote Sensing, 2011,32(22):7847-7871.

[11]何超英,廖安平,陳志剛,等. NDVI時間序列在全球耕地提取中的應用[J].地理信息世界,2013(2):66-69.

[12]白文龍.關中地區植被覆蓋變化及其主要驅動因子分析[D].西安:陜西師范大學,2013.

[13]周玉潔,王卷樂,郭海會.基于諧波分析法和線性光譜模型的耕地信息提取[J].遙感技術與應用,2015,30(4):706-713.

[14]賈明明,王宗明,張柏,等.綜合環境衛星與MODIS數據的面向對象土地覆蓋分類方法[J].武漢大學學報:信息科學版,2014,39(3):305-310.

[15]左麗君,張增祥,董婷婷,等. MODIS/NDVI和MODIS/EVI在耕地信息提取中的應用及對比分析[J].農業工程學報,2008,24(3):167-172.

[16]王瀚征.基于非監督分類與決策樹相結合的30 m分辨率土地利用遙感反演研究[D].石家莊:河北師范大學,2016.

[17]Breiman L, Friedman J H, Olshen R A, et al. Classification and Regression Trees[M]. California:Wadsworth International Group, 1984:1-358.

[18]趙萍,傅云飛,鄭劉根,等.基于分類回歸樹分析的遙感影像土地利用/覆被分類研究[J].遙感學報,2005,9(6):708-716.

[19]Yohannes Y, Hoddinott J. Classification and Regression Trees:an Introduction[M]. Washington D C:International Food Policy Research Institute, 1999.

[20]于文婧,劉曉娜,孫丹峰,等.基于HJ-CCD數據和決策樹法的干旱半干旱灌區土地利用分類[J].農業工程學報,2016,32(2):212-219.

[21]王正興,劉闖,HUETE Alfredo.植被指數研究進展:從AVHRR-NDVI到MODIS-EVI[J].生態學報,2003,23(5):979-987.

猜你喜歡
耕地分類
我國將加快制定耕地保護法
今日農業(2022年13期)2022-11-10 01:05:49
保護耕地
北京測繪(2021年12期)2022-01-22 03:33:36
新增200億元列入耕地地力保護補貼支出
今日農業(2021年14期)2021-11-25 23:57:29
分類算一算
垃圾分類的困惑你有嗎
大眾健康(2021年6期)2021-06-08 19:30:06
分類討論求坐標
耕地時節
數據分析中的分類討論
教你一招:數的分類
給塑料分分類吧
主站蜘蛛池模板: 欧美成人在线免费| 真人高潮娇喘嗯啊在线观看| 99精品免费在线| 国产亚洲视频播放9000| 狠狠亚洲婷婷综合色香| 精品国产一区二区三区在线观看 | 一本大道香蕉中文日本不卡高清二区| 日韩av电影一区二区三区四区| 国产午夜精品一区二区三| 丁香综合在线| 久久亚洲天堂| 青青久在线视频免费观看| 免费xxxxx在线观看网站| 久久久久久久97| 99视频在线看| 99久久国产精品无码| 日韩激情成人| 国产特一级毛片| 丰满人妻久久中文字幕| av一区二区三区在线观看| 伊人色天堂| 人妻熟妇日韩AV在线播放| 中文字幕在线日韩91| 国产精品无码AⅤ在线观看播放| 国产人前露出系列视频| 99re免费视频| 99人体免费视频| 国产第一页免费浮力影院| 亚洲中文字幕国产av| 怡红院美国分院一区二区| 日本在线亚洲| 国产第八页| 国产天天色| 欧美另类第一页| 99re在线视频观看| 亚洲综合精品第一页| 在线观看亚洲精品福利片| 人妻精品久久无码区| 国产欧美视频一区二区三区| 人妻精品久久无码区| 国产精品视频a| 99国产精品一区二区| 国产在线精品香蕉麻豆| 亚洲国产成人精品无码区性色| 国产小视频a在线观看| 国产精品成| a亚洲视频| 中文字幕 91| 一本大道无码日韩精品影视| 国产高清免费午夜在线视频| 91人妻日韩人妻无码专区精品| 国产一级视频久久| 亚洲性视频网站| 天堂中文在线资源| 无码内射中文字幕岛国片 | 操国产美女| 欧美国产精品不卡在线观看| 精品人妻系列无码专区久久| 91口爆吞精国产对白第三集| 国产精品综合久久久| 一区二区三区国产精品视频| 国产h视频在线观看视频| 久操中文在线| 欧美无专区| 国产福利拍拍拍| 日本影院一区| 久久香蕉国产线| 视频国产精品丝袜第一页| 国产成人综合亚洲欧美在| 九月婷婷亚洲综合在线| 国产喷水视频| 三级欧美在线| 亚洲 欧美 偷自乱 图片| 香蕉久久国产超碰青草| 欧美午夜小视频| 亚洲高清日韩heyzo| 日韩黄色在线| 国产色图在线观看| 精品国产欧美精品v| 波多野结衣一二三| 欧美亚洲一区二区三区导航| 中文字幕乱码中文乱码51精品|