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

基于多時相植被指數的云南高原山地冬小麥識別與研究

2024-01-15 11:00:18楊永明牛昱杰安衛國顏定飛
農業工程 2023年9期
關鍵詞:耕地特征模型

楊永明, 牛昱杰, 安衛國, 郭 鈺, 顏定飛

(1.滇西應用技術大學地球科學與工程學院,云南 大理 671000; 2.昆明理工大學國土資源工程學院,云南 昆明650031; 3.云南省高校山地實景點云數據處理及應用重點實驗室,云南 大理 671006; 4.滇西應用技術大學多源數據融合實景三維構建研究科技創新團隊,云南 大理 671006; 5.海南大學土木建筑工程學院,海南 ???570228)

0 引言

糧食的安全保障是國家的基礎,冬小麥是我國重要的糧食作物之一[1-2]。云南省位于我國西南冬麥區,其農業特色是以壩子農業為主。對云南省冬小麥種植面積的快速預測,是準確預測糧食產量的基礎。此外,云南省高原山地冬小麥的空間分布對糧食政策的制定具有重要作用。傳統農作物種植面積和空間分布信息需要人工現場實地考察等方式獲取,不但耗費大量人力、物力,同時無法及時快速地更新有關信息[3]。而隨著中高分辨率遙感衛星技術的發展,利用多時相遙感數據便可以對地形復雜、植被茂密的云南省高原山地農作物的相關信息進行快速地獲取,從而優化工作流程,因此利用遙感數據對云南省高原山地農作物種植面積的提取具有重要意義[4]。

隨著衛星遙感影像技術的發展,光學衛星影像作為數據源被廣泛應用在耕地信息提取的各個領域[5-6]。MODIS 數據具有重返周期短、光譜范圍廣、更新頻率高、數據易接受和相對完善的處理產品等特點,近年來在作物耕地信息提取研究中時序數據被廣泛應用,但受限于MODIS 的250 m 空間分辨率而無法對耕地作物精細化提取[7-9]。Landsat 8 衛星數據(30 m)重返周期為16 d,相較于前者有較高空間分辨率,有研究指出利用Landsat 8 OLI 數據獲取的多特征信息,通過隨機森林分類方式能夠有效地對耕地信息進行提取[10]。GF 系列光學衛星數據可利用其較高的空間分辨率,降低影像中像元的混合比例,在對耕地信息精細化的提取研究中被廣泛使用[11]。2017—2022 年,歐空局開始提供較高空間分辨率、波段信息豐富、雙星互補重返周期為5 d 的哨兵2 號(Sentinel-2)衛星數據,相較于上述衛星數據具有高空間分辨率、重復周期短等優勢,在對河北省的冬小麥精確提取研究中,利用歸一化差異植被指數(Normalize Difference Vegetation Index,以下簡稱植被NDVI)長時序影像提取的總體精度高達92.8%,Kappa系數在0.855 8 以上,試驗表明,Sentinel-2 為農作物耕地信息精準提取提供了數據支撐[4]。

耕地信息的提取主要基于分類算法進行研究,隨著更多新算法的出現及算法技術的改進,適宜的算法可以逐漸提高對作物的識別效率及精度[5]。有研究表明,基于多源衛星數據集,利用不同指數特征信息方法對我國北方產糧區的冬小麥種植面積提取的精度高達96.31%;在我國南方地區,利用分層提取原理對冬小麥兩個關鍵期提取的種植面積精度約90.2%,而利用多源時序數據的多特征組合分析等方法則將提取精度提高至90.6%[12-15]。隨著算法技術的不斷進步,農作物耕地面積提取精度也在不斷提高。在空間分辨率2.5 m的SPOTS 遙感數據對耕地信息提取試驗中,SEaTH 算法與CART 算法的總體精度均超過90%,但SEaTH 算法更適合快速提取耕地信息[16]。隨著Sentinel-2衛星數據對耕地信息提取研究的廣泛應用,有研究表明,通過結合10 m 分辨率的Sentinel-2 衛星數據和其光譜特性、指數特征等信息,對比4 種不同分類算法對耕地信息的提取,證明利用RF 算法耕地信息提取結果最佳精度最高達到88.52%,有效地優化了提取結果[17]。而在耕地面積稀少、壩子農業特征明顯、植被覆蓋率高及光學衛星影像覆蓋缺失等因素背景下,對我國云南省高原山地地區的農作物耕地信息提取的研究較少。近年來,利用多期影像融合來提取耕地信息的技術被逐漸應用于耕地信息提取研究中,并獲得高精度耕地信息,利用增強植被指數(Enhanced Vegetation Index,以下簡稱植被EVI)提取南方冬小麥種植面積總體精度約為87.1%,但利用多時序影像的不同植被指數對比方式提取冬小麥種植面積的研究較少[9,18-20]。

本研究區位于云南省西部橫斷山系縱谷區,地形復雜,植被覆蓋率高,具有我國西南高原地區特有的壩子農業,受“異物同譜”“異譜同物”的影響,單期的哨兵遙感影像數據無法對云南省高原地區農作物種植面積進行精確提取[21]?,F有對耕地信息提取的研究多針對冬小麥產糧大區,但針對西南冬小麥區的相關研究較少,或多利用中低分辨率遙感影像對冬小麥范圍進行粗提取,無法實現對西南高原山地冬小麥的精細化識別[22-25]。鑒于此本研究提出了利用Sentinel-2數據,基于多時相影像對2020—2021 年云南省高原山地農作物種植面積的預測方法進行研究,揭示多時相植被指數數據可以有效反映云南高原山地農作物生長變化的規律,為當地的農業生產提供數據支撐,明確了利用多時相數據和多特征信息可有效對農作物耕作面積進行預測。因此,對高原山地農作物種植面積識別方法的研究,為西南冬麥區的識別與提取提供參考,具有一定的研究價值。

1 研究區概況與數據分析

1.1 研究區

本研究區位于南澗彝族自治縣(以下簡稱南澗縣),隸屬于云南省大理白族自治州(100°06′~100°41′E,24°39′~25°10′N),總面積達1 738.82 km2,地處橫斷山系縱谷區,地形復雜,多以山地為主,并且耕地較為零散,種植農作物歷史悠久,以冬小麥和玉米為主,南澗縣具有我國西南冬麥區典型耕種特征,為大理州冬小麥種植面積較大的縣域。自然生態資源豐富,種植冬小麥作物歷史悠久,2021 年南澗縣縣域內冬小麥種植面積為2 694 hm2,冬小麥生長期是10 月下旬播種,次年5 月收割,冬小麥物候期如表1 所示[26]。

1.2 數據來源介紹

Sentinel-2 衛 星 數 據 由 A( Sentinel-2A) 與B(Sentinel-2B)雙星互補提供,重返周期5 d。本研究區域由編號為RPH 一景影像覆蓋,使用影像數據來源于https://scihub.copernicus.eu 平臺所提供的S2MSI 產品,包含1C 級數據和2A 級數據,其中1C 為經過正射校正、幾何校正的產品,2A 級數據為經過BOA 正射校正的影像。Sentinel-2 衛星影像覆蓋有13 個波段,包含10、20 和60 m 共3 種空間分辨率,選擇的時相是從2020 年10 月—2021 年6 月,在冬小麥生長期的每個階段獲取一景云量<30%研究區影像,影像時間間隔盡可能保持一致。Sentinel-2 衛星數據使用SNAP 軟件處理工具,利用Sen2Cor 插件對其大氣校正和輻射定標并重采樣獲得10 m 影像,通過構建的數據集,可對其中地類光譜信息、各種植被指數等信息進行分析。

植被指數(Vegetation Index,以下簡稱植被VI)通過植被的光譜特征有效地表現植被的遙感影像特征信息,是云南省高原山地農作物夏糧快速識別的關鍵,可以反映特定植被變化特征。選擇了光譜特征中與農作物密切的藍色波段、綠色波段、紅色波段和近紅外波段4 個波段,如表2 所示。

長時間序列影像構建的植被曲線圖可以有效反映研究區植被的生長狀況。基于GEE(Google Earth Engine)平臺,結合高清衛星影像選擇研究區地類樣本降噪重建,最終構建研究區地類的長時序曲線圖。

2 研究方法

本研究包括多時相數據集準備、農作物特征識別與比較、監督分類模型提取、精度評定及作物特征可分離性討論。技術路線如圖1 所示。首先對多期Sentinel-2 衛星數據預處理并提取多種植被指數,然后將每種植被指數多期影像進行組合,最后使用監督分類模型對影像中農作物指數特征分析其可分離性,最終根據研究結果討論不同特征信息對農作物識別的影響及不同植被指數模型之間的差異。

圖1 技術路線Fig.1 Technical route

為探究不同耕地特征對云南高原山地農作物識別的影響因素,本研究選取植被指數、光譜曲線、物候特征、紋理和長時序植被指數曲線特征,構建多期植被指數影像合成模型進行農作物耕地范圍的識別與提取。

基于Sentinel-2 時序遙感影像的高原山地地區耕地作物監測及信息提取,主要涉及兩個部分,一是利用哨兵影像對耕地信息的識別,二是耕地信息提取方法的選擇。在哨兵影像數據10 m 分辨率對云南高原山地研究區作物識別過程中,受到云、霧等自然條件的影響,通過單期影像對冬小麥耕地的識別具有一定的局限性,因此,需要基于面向對象的方法,結合冬小麥的物候譜、光譜曲線和長時序植被指數曲線重構等多種條件對冬小麥耕地信息識別。

2.1 時序植被指數合成

時序植被指數影像可以有效突出冬小麥的地類特征,更好展現出不同地物的可分離性。本研究根據植被指數特性,選用3 種植被指數,包括植被NDVI、植被EVI和植被RECI,以此表現出不同地類在時序指數影像種的特征。如圖2 和圖3 所示,圖中綠色(中間部分)地類表示以冬小麥為主的指數呈增長趨勢的作物。為進一步確定冬小麥所對應圖斑的顏色,在選冬小麥樣本的時候,遵循冬小麥在對應生長期的光譜波段特定數值區間的要求,以及對照高清影像等對地類樣本進行確認。

圖2 農作物樣本示意Fig.2 Schematic diagram of crop samples

圖3 植被指數融合影像示意Fig.3 Schematic diagram of vegetation index fusion image

時序植被指數模型是指將一定生長期間隔的植被指數,通過波段合成技術(Layer Stacking)將其融合為一幅影像,從而獲得一張最大程度反映不同地物特征差異的圖像,并對不同時段像元地類的植被指數的變化情況進行識別與分析[27]。將融合后模型中的不同地類圖斑與多期的原始影像變化特征、高清影像做分析比對,從而得出不同圖斑地類分別代表的地類信息。

本研究中涉及的植被NDVI、植被EVI和植被RECI相關特性:植被NDVI可適用于跟蹤植被發育狀態及植被覆蓋率,是一種可反應植被空間分布密度的指標;植被EVI則是用于植被茂盛地區的一種指數,用來增強對植被覆蓋區的敏感性,較好地處理植被NDVI的飽和問題;植被RECI是指當植被覆蓋度較高的時候,對植被異常敏感,并且對植物中氮含量敏感,可用于植物的生長發育期[28-33]。

式中NDVI——植被NDVI值

EVI——植被EVI值

RECI——植被RECI值

ρNIR、ρRED、ρBLUE——近紅外、紅色和藍色波段的光譜反射率,在Sentinel-2 影像數據中,表示波段8(842 nm)、波段4(655 nm)和波段2(490 nm)

生長關鍵期是指農作物區別于其他地類在植被指數變化中特有的時間段,因此選擇農作物生長期內的關鍵期能夠得到精度較高的分類結果。本研究選取冬小麥、同期作物、林地、水域和不透水層5 種地類樣本,獲取并構建的3 種植被指數的時序曲線,如圖4、圖5 和圖6 所示。研究表明,冬小麥的植被指數在每年11 月—次年1 月的時序曲線區別于其他地物植被指數在冬小麥收獲期(次年5 月)前后,植被指數會驟降,但在云南省高原地區每年的5—8 月受到云量增多等因素的影響,導致影像缺失。因此,主要選擇出苗-分蘗時間段(11 月上旬—次年1 月)為農作物生長關鍵期對冬小麥進行研究分析?;贕EE 平臺,獲取Sentinel-2A 數據影像,構建云南省高原山地各地類的不同植被指數長時間序列的變化規律。

圖4 不同土地利用類型植被NDVI 時序曲線Fig.4 Time series curves of vegetation index of different land use types

圖5 不同土地利用類型植被RECI 時序曲線Fig.5 Time series curves of vegetation index of different land use types

圖6 不同土地利用類型植被EVI 時序曲線Fig.6 Time series curves of vegetation index of different land use types

2.2 物候及光譜曲線特征

農作物的物候是反映作物受自然環境條件的影響產生的周期性的生長規律,有研究指出農作物的植被指數會隨著其所處生長期的不同而產生一定的數值變化,通過對比不同植被指數對農作物特征呈現情況,選擇最適宜的植被指數對農作物進行識別[34]。光譜數值隨著作物生長期發生改變,在對農作物識別時通過對波段閾值的選擇來識別。

在監督分類過程中可視化每個地物的每個波段的特征,有助于確認地物的可分離性,此類圖表為光譜特征。分別選取了2020 年11 月20 日、2020 年12 月30 日和2021 年1 月29 日Sentinel-2A 遙感影像,并基于高清影像的判斷及現場識別等方式,分別選取地物樣本有冬小麥、水域、不透水層、林地和其他作物,獲得各日期的光譜曲線特征,如圖7 所示。

圖7 作物生長關鍵期光譜曲線Fig.7 Spectral curve of critical crop growth period

在對研究區冬小麥作物的影像光譜特征曲線統計分析中,農作物的生長關鍵期冬小麥的光譜曲線數值在第8 波段(842 nm)上升明顯,并且12 月的數值為3 500~4 000 nm,而到1 月該數值則會超過4 000 nm;第4 波段(665 nm)在關鍵期處于水平位,為500~1 000 nm,沒有明顯變化趨勢;第2、3 波段(490 nm、560 nm)由11 月—次年1 月先增加后沒有明顯變化。根據這一特性,在樣本選取時,對應的日期及波段值與這一特性相符,可以有效地提高樣本選擇精度。

2.3 紋理特征

圖像的紋理特征可以呈現對應地物的表面性質及空間領域灰度值分布情況[35]。而影像的紋理特征是一種結構化的表現形式,在Sentinel-2 遙感影像中可選擇紅色波段、近紅外波段和綠色波段對農作物有效識別。

3 結果與分析

為了尋找一種快速提取云南高原山地農作物耕地信息的方法,本研究使用植被指數、光譜曲線、物候特征和長時序植被指數曲線特征4 類特征研究對比,結合支持向量機(Support Vector Machine,SVM)算法進行識別分類試驗,分析不同特征條件下的影響因素對多期植被指數模型的分類差異。

3.1 不同植被指數的種植面積提取結果

本研究以2020—2021 年的數據進行數據提取,采用Sentinel-2 高分辨率多光譜衛星數據,進行冬小麥信息的目視解譯,使用多期植被合成的方法影像合成,基于規則在南澗縣內選取229 個樣本。

在對多期植被指數模型使用SVM 的監督分類后得到的分類結果如表3 所示,當年冬小麥耕地公布面積約2 694 hm2。

表3 不同植被指數3 期影像合成冬小麥提取結果Tab.3 Extraction results of three-stage image synthesis of winter wheat with different vegetation indexes

對比3 種不同植被指數的提取結果,多期植被NDVI模型提取結果優于另兩種植被指數,對細碎的冬小麥分布地區識別的用戶精度為93.28%,而植被EVI多時相融合模型和植被RECI多時相融合模型提取結果,無法滿足實際要求,在河流和道路邊界處存在錯提現象,細碎地塊的冬小麥受到周圍環境及混合像元的干擾無法識別提取,圖像像元存在同質化現象,因此造成農作物漏提現象的產生。對比3 者提取,結果表明,使用植被NDVI的多時相融合模型在SVM 監督分類下,該模型與云南省高原山地的冬小麥耕地信息有較高相關性,可以有效地對冬小麥面積進行預測。

3.2 不同時相組合的種植面積提取

表4 列出了基于本研究區內兩種日期組合方式下,多期植被NDVI模型的提取結果。本研究中對比植被NDVI在使用2 期指數模型和3 期指數模型下分類的差異,以及這兩種模型的提取結果。2021 年南澗縣統計年鑒給出的冬小麥種植面積為2 694 hm2,3 期植被NDVI實際預測種植面積為2 726.8 hm2,Kappa 系數為0.91,2 期植被NDVI實際預測種植面積為3 170.1 hm2,Kappa 系數為0.89[26]。對比植被NDVI的兩種合成期數在影像實際提取結果可以發現,3 期影像合成模型對冬小麥分類提取結果精度更高,2 期植被指數模型依然會存在一定的地物錯提現象。

表4 不同合成期數提取結果Tab.4 Extraction results of different synthesis periods

對比植被NDVI的2 種日期組合方式下在南澗縣的冬小麥面積預測結果表明,2 種多時相的組合方式均滿足對耕地面積的提取,其中3 期的組合結果更優。

3.3 精度驗證

在對冬小麥耕地信息提取結果的精度進行定量評定時,根據統計年鑒冬小麥種植面積數據對提取結果檢驗,隨機選取一定比例的冬小麥圖斑利用Google 地圖及World Imagery Wayback 的高清歷史影像來確認提取冬小麥的準確率,再使用混淆矩陣(Confusion Matrix)來分析誤差,通過上述研究可知利用3 期的植被NDVI構建的多時相植被指數融合模型提取的冬小麥種植面積較為精準,并且圖斑結果較為準確,但存在少量被其他作物與歷史影像對比后確認為被錯誤識別提取,及部分漏提現象的產生如表5 所示。

表5 2020—2021 年南澗縣樣本點個數及正確分類個數Tab.5 Number of sample points and correct classification in Nanjian County from 2020 to 2021

分別選取167 個為冬小麥樣本及33 個為非冬小麥樣本,經檢驗其中153 個樣本點正確,14 個樣本點錯誤,分類精度91.6%。分析提取誤差的原因主要是由于Sentinel-2 衛星數據為10 m 分辨率,在冬小麥區域邊界地帶存在混合像元,導致選擇樣本存在一定差異。此外,人工在選擇冬小麥樣本的標注過程中,樣本標注不準確也對結果存在一定的影響。

4 討論

受到云南省自然環境的影響,對高原地區作物的提取研究多是采用監督分類算法或植被NDVI加權指數(WNDVI)分類算法,在通過多時相的植被指數合成模型提取的研究較少[4]。恰當的多時相植被指數合成模型特征可以影響農作物提取效果及精度,在本研究的模型中涉及的農作物特征信息包含有光譜特征、影像紋理特征、空間關系特征及植被指數時序曲線特征參與到農作物信息的提取,其中光譜特征是提取農作物最主要的,再結合其他3 種特征可以有效地提高農作物提取精度[19]。在比較了3 種植被指數構建模型的影像提取精度,發現并非所有的植被指數模型都可以適用于云南省高原山地農作物的識別提取。因此,在選擇多時相植被指數合成模型時候,需要在長時間序列下,結合光譜曲線的變化,以及指數變化等特征,選擇適用于農作物變化規律的植被指數構建模型進行分類,以提高作物的識別精度。

本研究中,使用10 m 空間分辨率的Sentinel-2 衛星數據,構建時序指數模型,利用監督分類(支持向量機)算法,分別得到了植被NDVI總體分類精度93.97%(Kappa=0.91),植被EVI總體分類精度79.02%(Kappa=0.74),植被RECI總體分類精度84.82%(Kappa=0.73)。通過對3 種指數構建的模型對比高清歷史影像可以發現,植被NDVI構建的模型,對冬小麥提取效果較好,生產者精度96.39%。植被EVI算法中加入了藍波段,在觀察冬小麥關鍵期內各波段變化情況可知藍波段暫無明顯變化,故而造成其他植被變成干擾影像。植被EVI受植被氮元素影響較明顯,所以對生長期的植被都會敏感,因此錯提了很多河流及道路兩旁的其他生長中的植被。在對多時相植被NDVI模型提取結果進行隨機抽樣200 個樣本點,其中冬小麥分類精度91.6%,滿足對冬小麥種植面積的提取要求。在模型的實用性方面,該模型在對我國其他冬小麥產糧區精細化提取時,利用多時相的植被NDVI合成模型可較好地適用于我國大部分冬小麥產區的種植面積預測??傮w來看,多時相的植被NDVI合成模型可以符合云南省高原山地農作物種植面積的預測。

在研究中,由于受到光學衛星分辨率限制,在多期植被指數合成得到的影像中,地類邊界存在像元內地物混合現象,進而導致在選擇樣本的初期要通過不斷重復篩選冬小麥樣本以盡可能地提取準確的冬小麥樣本,避免錯提干擾地物。本研究通過固定冬小麥的單期植被指數閾值的方式,利用冬小麥每個生長期植被指數在一個閾值區間波段的特性,篩選混合像元內以冬小麥為主的樣本,基于上述樣本選取規則可用于對冬小麥種植面積的精準提取。

5 結束語

利用云南省南澗縣2020—2021 年Sentinel-2 時序數據集,基于ENVI5 平臺及GEE 平臺提取植被指數、光譜曲線、農作物物候和長時序植被指數曲線,將其特征規則錄入多期植被指數模型中,最后利用SVM 監督分類模型提取云南高原山地農作物耕地信息,獲取研究區的耕地面積,得出以下研究結論。

(1)基于Sentinel-2 衛星影像的多期植被指數合成模型,分別對比植被NDVI、植被EVI和植被RECI構建的模型與耕地信息識別的相關性,其中多期NDVI指數融合模型,提取研究區冬小麥耕地信息的用戶精度達到93.28%,Kappa 系數0.91,多期植被EVI合成模型的用戶精度63.8%,Kappa 系數0.74,多期植被RECI合成模型的用戶精度66.36%,Kappa 系數0.73。證明多期植被NDVI合成模型與冬小麥耕地信息具有較高正相關性,結果滿足對農作物種植面積預測的要求。

(2)在不同數量的日期構建的多期植被NDVI合成模型與耕地信息提取結果的相關性研究中,分別對本研究區域內利用3 期(2020 年11 月、2020 年12 月、2021 年1 月)和2 期(2020 年11 月、2021 年1 月)構建的植被NDVI合成模型的總精度和Kappa 系數分別為93.97%、90.42%、0.91、0.89。結果表明,在對不同期數組合下,3 期植被指數構建的模型對提取結果改進,精度得到提升。此外也驗證了,在部分影像缺失地區,利用2 期植被指數模型對耕地信息的提取結果依然符合要求。

(3)本研究模型對南澗縣冬小麥耕地識別與提取,利用高清影像判讀得到的樣本點提取精度進行驗證,167 個冬小麥樣本中有153 為正確分類,提取的分類精度91.6%,提取的冬小麥預測面積2 726.8 hm2,本年度冬小麥實際統計面積2 694 hm2。通過對提取結果分析可知,受樣本提取規則、混合像元以及同期作物的影響,多期植被指數指數模型對冬小麥耕地提取中,存在將道路、河流兩側的綠植錯提、漏提細小地塊。

猜你喜歡
耕地特征模型
一半模型
自然資源部:加強黑土耕地保護
我國將加快制定耕地保護法
今日農業(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
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
如何表達“特征”
不忠誠的四個特征
當代陜西(2019年10期)2019-06-03 10:12:04
抓住特征巧觀察
主站蜘蛛池模板: 91尤物国产尤物福利在线| 波多野结衣一区二区三区AV| 国产特一级毛片| 狠狠色香婷婷久久亚洲精品| 国产精品成人久久| www.国产福利| 无码一区18禁| 2021最新国产精品网站| 在线毛片网站| 亚洲精品777| 精品一区二区三区四区五区| 无码专区第一页| Jizz国产色系免费| 91亚洲国产视频| 亚洲日本一本dvd高清| 国产91线观看| 毛片网站在线播放| AV熟女乱| 亚洲天堂在线免费| 欧美日韩一区二区在线免费观看| 亚洲精品无码av中文字幕| 免费看久久精品99| 欧美日韩国产在线人| 无码福利视频| 日韩欧美色综合| 色香蕉影院| 自慰高潮喷白浆在线观看| 欧美色99| 亚洲精品国产精品乱码不卞| 91午夜福利在线观看| 欧美专区日韩专区| 免费一级毛片在线观看| 国禁国产you女视频网站| 亚洲天堂区| 国产精品一区不卡| 国产精品男人的天堂| 97人妻精品专区久久久久| 97精品国产高清久久久久蜜芽| 欧洲精品视频在线观看| 超清无码一区二区三区| 好紧太爽了视频免费无码| 日本黄色不卡视频| 九九线精品视频在线观看| 国产在线高清一级毛片| 成人午夜久久| 国产成人在线无码免费视频| 国产美女主播一级成人毛片| 国产原创演绎剧情有字幕的| 国产高清不卡视频| 999精品色在线观看| 亚洲国模精品一区| 欧洲高清无码在线| 亚洲国模精品一区| 亚洲资源站av无码网址| 国产清纯在线一区二区WWW| 欧美成人日韩| 国产探花在线视频| 麻豆国产精品视频| 国产白浆一区二区三区视频在线| 亚洲第七页| 日韩天堂视频| 国产精女同一区二区三区久| 国产精品成人久久| 手机永久AV在线播放| 国产永久在线视频| 国产精品免费露脸视频| 国产成人精品一区二区三区| 免费人成视网站在线不卡| 青青网在线国产| 免费又黄又爽又猛大片午夜| 国模沟沟一区二区三区| 亚洲一区二区三区香蕉| 色天天综合| 欧美精品不卡| 免费视频在线2021入口| 91福利一区二区三区| 欧美日韩亚洲国产主播第一区| 免费毛片全部不收费的| 黄色网页在线播放| 久久91精品牛牛| 亚洲国产成熟视频在线多多| 久久久久国色AV免费观看性色|