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

基于深度學習和多源遙感數據的玉米種植面積提取

2023-02-24 09:47:05呂偉宋軒楊歡
江蘇農業科學 2023年23期
關鍵詞:深度學習

呂偉 宋軒 楊歡

摘要:玉米作為我國主要糧食作物之一,及時準確監測其種植范圍及面積對農業產能評估、保障糧食安全具有重要意義。以華北平原典型農業區——原陽縣為例,基于歐空局Sentinel-1 SAR和Sentinel-2 MSI遙感影像數據,在谷歌地球引擎云平臺的支持下通過提取雷達后向散射系數時序曲線以及歸一化植被指數(NDVI)時序曲線,搭建卷積神經網絡(CNN)模型,并將時序數據輸入模型得到典型地物分類結果,提取了研究區玉米種植區域,利用野外調查數據進行精度驗證,并與隨機森林分類對種植區的提取結果進行對比。結果表明,基于光學和SAR融合遙感影像數據的識別效果最佳,總體精度達到93.33%,κ系數為0.911;與隨機森林分類法相比,卷積神經網絡分類的總體精度更高,分類效果更好。因此,采用卷積神經網絡以及多源遙感數據的融合能夠實現玉米種植面積的準確監測。

關鍵詞:深度學習;多源遙感;融合數據;卷積神經網絡;種植識別;時序數據;玉米

中圖分類號:TP79;S127? 文獻標志碼:A

文章編號:1002-1302(2023)23-0171-07

玉米作為我國主要糧食作物之一,及時準確監測玉米的種植面積對掌握玉米產量、保障糧食安全以及經濟和環境的可持續發展有至關重要的意義[1]。傳統玉米種植面積監測主要通過問卷調查和實地訪談統計信息并逐層上報,需要大量的人力物力[2],且耗時較長[3]。遙感分類法是基于1幅或多幅玉米種植影像進行面積提取的方法[4],衛星遙感技術以其覆蓋面大、重訪周期短、多時空分辨率等優勢,在植被面積監測方面取得了巨大成就[5]。遙感分類法主要包括2種,一種選取作物光譜特征差異明顯的1期影像來實現分類,另一種是基于多期遙感影像,提取作物生育期內光譜變化的時序特征進行面積提取。后者精度更高,是目前作物精細分類的主要方式[6],但是需要人工預處理并提取相關特征[7],數據重建工作復雜,不利于快速與自動化分類作物。近年來,深度學習作為機器學習和數據挖掘領域的突破性技術,廣泛應用于圖像分類、目標識別等領域。深度學習在特征表示方面具有特有的靈活性,基于原始輸入數據和輸入標簽,使用卷積運算自動提取特征,實現模型端到端訓練,不需要專家先驗知識即可實現高精度、自動化分類[8]。Kussul等基于多源多時序遙感數據對比分析后認為,卷積神經網絡(CNNs)要優于隨機森林方法,各種作物的目標分類精度均在85%以上[9]。Garnot等使用CNN、RNN、R-CNN對時序哨兵二號數據進行分類,發現所有深度學習模型均優于隨機森林模型,突出了深度學習模型在農業地塊分類方面的潛力[10]。目前國內外作物面積提取主要采用多光譜遙感影像,但是夏季作物生長關鍵期云雨天氣頻繁[11],難以獲得有效的多時相遙感數據[12]。相比之下,合成孔徑雷達能穿透云、霧,不受天氣影響,但雷達影像中存在很大的散斑噪聲[13]和幾何形變現象(如透視收縮、疊掩和陰影等)[14],會影響地物識別的精度。如何將多光譜和雷達圖像進行融合,從而提高作物種植面積提取精度顯得尤為重要。因此,本研究以華北平原典型玉米種植區的識別為例,基于Sentinel-1 SAR和Sentinel-2 MSI融合的時間序列遙感數據,采用卷積神經網絡算法(CNN),評估CNN在時間序列分類任務中的適用性,探索光學和雷達遙感融合數據在作物種植類型識別中的優勢。

1 研究區域和數據處理

1.1 研究區域概況

研究區選擇河南省新鄉市原陽縣,該縣地處豫北平原,南臨黃河,位于34.55°~35.11°N、113.36°~114.15°E,東接封丘縣,西鄰武陟縣、獲嘉縣,背靠新鄉縣、延津縣,南與鄭州市隔河相望,總面積1 022 km2。地貌屬沖積平原,地勢西南高、東北低。溫帶大陸性季風氣候,年均氣溫約為14.3 ℃,平均降水量為556? mm,全年無霜期227 d。

原陽縣農業發達,是華北平原主要農作物種植區之一,以1年2熟的作物輪作模式為主,玉米為夏季主要作物之一,一般于6月中上旬播種,8月進入生長旺盛期,10月中下旬收獲,生育期為每年的6—10月(圖1)。

1.2 數據來源及預處理

1.2.1 遙感影像數據

選取哨兵一號雷達影像(Sentinel-1 SAR )和哨兵二號多光譜影像(Sentinel-2 MSI),數據訪問、處理及模型構建通過谷歌地球引擎(Google Earth Engine,GEE)云平臺進行。哨兵系列影像具有重訪周期短、空間分辨率高、影像覆蓋范圍大、數據免費等優點,是現今較有優勢的一種對地觀測數據。

1.2.1.1 Sentinel-1 SAR GRD數據

該數據通過GEE云平臺獲取,采用干涉測量寬幅模式(interferometric wide swath,IW)下的VH (垂直發射水平接收)極化方式和VV(垂直發射垂直接收)極化方式,利用Sentinel-1 Toolbox (S1TBX)工具進行預處理,包括熱噪聲去除、輻射定標、地形校正等,空間分辨率為 10 m,重訪周期為12 d。

1.2.1.2 Sentinel-2 MSI數據

Sentinel-2 MSI采用Level-1C產品,共13個波段,其中分辨率為10 m的波段有4個,分辨率為20 m的有6個,分辨率為60 m的有3個,光譜覆蓋范圍從可見光到近紅外、短波紅外波段[15],時間分辨率 5 d,其自帶質量評估(QA)波段識別云像元。本研究通過GEE云平臺篩選出作物生長期云量少于25%的8幅影像。

1.2.2 野外數據

為準確獲取樣本信息,2020年7月利用GPS采集作物分布信息,同時利用Sentinel-1/2 數據選取其他地物典型像元作為樣本,并將所選樣本導入 Google Earth,調用同年份高分辨率影像,檢查樣本正確性。最終獲得的樣本信息見表1。以6 ∶[KG-*3]2 ∶[KG-*3]2的比例從中選擇720個樣本進行訓練,240個樣本用于驗證分類精度,240個樣本用于測試模型準確性。

2 研究方法和技術路線

2.1 技術路線

本研究技術路線見圖2。第一,數據集的構建,利用GEE云計算平臺對Sentinel-1/2數據進行預處理,并根據樣本坐標點信息提取各種地物類型的不同遙感數據作物生長期時序曲線。第二,對各類時序數據進行組合,在GEE云平臺上基于不同組合數據集構建單變量或多變量卷積神經網絡,對數據進行作物分類,對比試驗采用隨機森林方法。第三,分析測試數據的分類結果,得到作物分類識別的最優數據組合。

2.2 時間序列NDVI/后向散射曲線提取

2.2.1 GEE提取時序數據

歸一化植被指數(normalized difference vegetation index,NDVI)能較好地反映植被綠度變化,是植被生長狀態及植被覆蓋度的最佳指示因子[16],其結果在[-1,1]之間。計算公式為

NDVI=(ρNIR-ρR)/(ρNIR+ρR)。(1)

式中:NDVI表示歸一化植被指數;ρNIR表示近紅外波段反射率;ρR表示紅外波段反射率。

將樣本坐標點處理成shp文件,上傳至GEE云平臺顯示標注地圖,提取每一地物樣本數據集(以實地考察樣點為主)相應的時序數據。并通過光學影像質量標識剔除因云雨天氣導致質量不高的數據[17],最終得到相應的NDVI、VV后向散射系數、VH后向散射系數時序數據。

對不同時序數據進行組合,共產生5組數據集:NDVI數據、交叉極化數據(VH數據)、單一極化數據(VV數據)、融合2類極化數據(VH+VV數據)、融合2類影像提取的NDVI和極化數據得到的融合數據(NDVI+VV+VH)。

2.2.2 Savitzky-Golay濾波平滑

由于云、氣溶膠、太陽高度角等因素,原始遙感數據集存在很多噪聲,在此情況下提取的時序曲線,波動性較大,不具有代表性。因此,需要對時間序列進行去噪處理。Savitzky-Golay(S-G)濾波被廣泛應用于數據流平滑去噪,可以在消除噪聲的同時確保信號的形狀、寬度不變,是遙感植被指數時間序列的主要濾波方法。它是一種時域低通濾波法,通過局部多項式回歸模型平滑時序數據,是移動窗口的加權平均算法,通過將給定高階多項式的最小二乘法在滑動窗口內擬合得到加權系數。基本思想是:基于多項式在濾波窗口內利用最小二乘法對數據進行擬合[18]。

式中:Y表示原始數值;Y′表示擬合值;Ci表示第i個點的權重;N=2n+1表示濾波窗口的大小。

進行數據清洗后,將每類地物所有樣本的時序數據在同一時間上分別求平均,通過 S-G 濾波器對平均時序數據進行濾波處理,得到較平滑的時序曲線[19]。為得到較好的平滑結果,經多次試驗后,濾波核左右各選取5點,設平滑多項式次數為 2,最終的標準時序曲線分別見圖3、圖4、圖5。

2.3 卷積神經網絡

卷積神經網絡(CNN)是深度學習方法中最成功的網絡架構之一,是建立在傳統人工神經網絡上的一種深度學習算法,也是第1個成功訓練多層網絡的算法[20]。CNN具有局部連接、權值共享和池化層降采樣的特點,與傳統神經網絡相比,可以減少參數量、降低模型復雜度,并賦予模型對平移和形變的容忍性[21]。卷積層通過計算卷積核與覆蓋區域信號值的點擊來確定神經元的輸出。

CNN一般由用于分層提取特征的卷積濾波層和用于計算輸出值的全連接層2個部分組成。網絡層數的選擇、卷積核的數量和大小、激活函數等,在CNN結構設計中尤為重要。經多次調參嘗試,本研究選擇4層Conv1D來提取特征值,每2層Conv1D后添加1層MaxPooling1D保留主要特征,減少計算量。每層卷積層使用修正線性單元ReLU函數作為激活函數來提高神經網絡對模型的表達能力。池化層被固定為最大池化層,窗口大小為2。后面接1個全連接層、1個Dropout層、1個Dense層進行Softmax分類。模型訓練選擇Adam優化器,損失函數采用交叉熵損失函數(CEloss),公式如下[22]。

式中:N表示每批訓練樣本的數量;C表示所有類別的集合;I(n)c表示當前批次第n個樣本獨熱標簽類;P(n)c表示模型預測樣本n為c類的概率,ln P(n)c用來表示懲罰模型對錯誤分類的樣本的預測。

Dropout是一種正則化技術,通過隨機失活一些神經元,防止數據集較小時容易造成過擬合[23-24],以提高神經網絡的性能,設置為0.5。全連接層經過Dropout層后,輸入Softmax層,最終輸出4種類別的概率結果。經過濾波處理的5組不同特征組合的時序數據集均輸入該模型,模型架構見圖6。

2.4 精度評價指標

構建混淆矩陣時選取用戶精度(user accuracy,UA)、生產者精度(producer accuracy,PA)、總體精度(overall accuracy,OA)、κ系數、F1分數(F1-score,

別稱平衡F分數)5個評價指標。其中,OA用于描述驗證樣本與分類結果一致的概率;PA表示地面的某類別被正確分類的概率;UA表示正確分類某類別的概率,主要用來評價分類結果的可信度;F1分數指標為綜合考慮PA和UA的調和值[24];κ系數用于描述混淆矩陣的一致性,表征分類結果的可信度,當κ>80%時,表示分類精度最高;當 40%<κ<80%時,表示分類精度中等;當κ<40%時,表示分類精度最差[25]。

3 結果與分析

整體而言,采用1D-CNN分類方法的各種組合數據集均獲得較好的分類效果,總體分類精度均高于85%,κ系數高于80%,所有試驗結果中,NDVI數據與SAR極化數據融合分類效果最好。分類結果分別見表2、表3。隨機森林分類效果較差,其總體分類精度及κ系數結果見表3。

3.1 基于時序NDVI數據的典型地物分類

由表3可知,基于時序NDVI數據的典型地物分類總精度為92.08%,κ系數為0.894,高于SAR極化數據分類精度。從單個類別的分類效果看,采用NDVI數據分類在水稻、玉米、建筑等3類地物的提取中有較好的分類效果,僅低于NDVI數據和極化數據融合的分類效果,其中玉米提取的F1分數為87.3%,但在水體提取上效果最差。

3.2 基于時序SAR后向散射數據的典型地物分類

首先比較2種極化方式分類結果,發現VV極化數據分類效果最差,基于時序VV極化數據的識別總精度為86.25%,κ系數0.817;VH極化數據分類效果較好,優于VV極化效果,識別總精度為89.17%, κ系數0.856; 與單一VH極化數據相比,VH+VV的組合數據分類效果更好,總精度為91.21%,κ系數為0.883,相比于單獨利用VH或VV極化數據,VH+VV的組合分類總精度分別提高4.96、2.04百分點,κ系數分別提高0.066、0.027。從單個類別的分類效果來看,后向散射系數對水體的分類效果較好,因為水體與其他地物后向散射系數差異明顯,單獨使用VV極化數據和單獨使用VH極化數據,對水體分類的F1分數均能達到99%以上,但是對玉米這類旱地作物分類效果較差,尤其是VV極化數據,對玉米分類的F1分數僅能達到76.34%,但VH極化數據對玉米分類的F1分數能達到86.21%,說明在提取玉米面積時,雷達圖像具有可分性。

3.3 基于融合時序NDVI和SAR后向散射數據的典型地物分類

由表3可知,融合數據分類總精度為93.33%,κ系數為0.911。相比NDVI數據總精度提高1.25百分點,κ系數提高約0.017;相對于單一極化數據,融合數據總精度分別提高7.08、4.86百分點,κ系數分別提高0.094、0.055;相對于VV+VH組合數據,融合數據總精度提高了2.12百分點,κ系數提高0.028。融合數據的玉米提取的F1分數為88.71%,較其他數據方案均有明顯提升。綜上,融合2類遙感影像提取出的數據用于作物分類效果更佳,對玉米提取精度的提高有明顯效果。

3.4 CNN在玉米面積提取中的表現

整體而言,利用CNN對5類組合數據進行地物分類,均得到較好的分類效果,總體分類精度均高于85%,高于隨機森林分類方法,κ系數均高于80%,說明使用卷積神經網絡能夠有效地應用于典型地物識別。各數據組合方案中,除單一VV極化數據外,玉米提取F1分數均在85%以上,說明卷積神經網絡能有效區分玉米,能達到玉米識別精度要求。不同數據集整體分類結果較接近,但NDVI數據分類結果有一定的錯分現象,含有后向散射系數數據的分類結果椒鹽現象較明顯,這是因為雷達影像存在斑點噪聲,單個像元受噪聲影響較大,影響了時序曲線的構建。結合NDVI時序數據以及后向散射系數數據的分類精度更高,椒鹽現象及錯分現象都得到明顯改善,減少了將建筑錯分為玉米的現象,有利于玉米的遙感提取與制圖。

3.5 玉米面積提取

綜上,本研究采用卷積神經網絡和多源遙感時序數據,獲取原陽縣的玉米種植面積及分布(圖7),可見玉米為原陽縣主要夏季農作物,主要分布于除太平鎮和葛埠口鄉以外的大部分地區。《2021年新鄉統計年鑒》記錄原陽縣2020年的玉米種植面積為477.2 km2,分類結果中玉米種植面積為 446.1 km2,提取玉米面積與統計數據的一致性為93.5%。

4 結論與討論

4.1 結論

本研究基于Sentinel-1/2數據,采用5種組合數據方案,結合GEE云平臺以及1D CNN分類器,以華北平原典型農作物種植區為例,通過多光譜和雷達影像的數據融合,基于深度學習對玉米種植面積進行識別與提取,主要結論如下:(1)不同數據組合方式對分類識別結果存在差異,融合數據效果最好,總精度為93.33%,κ系數為0.911。SAR極化數據識別總精度均在85%以上,VV+VH的組合識別總精度最高,為91.21%,κ系數為0.883;NDVI數據分類精度高于SAR極化數據,總精度為92.08%,κ系數為0.894。(2)針對玉米提取,除VV極化數據效果較差,另外4種方案的F1分數均能達到85%以上。(3)采用卷積神經網絡,總體精度均在85%以上,κ系數均在0.8以上,分類效果均高于隨機森林的分類效果。

4.2 討論

本研究充分利用GEE云平臺強大的計算能力,通過建立研究區遙感地物分類模型,分析采用CNN算法使用Sentinel-1/2數據對華北平原玉米種植面積進行提取的效果。從方法上看,采用卷積神經網絡,使用光學遙感數據或SAR后向散射系數數據,均能有效提取玉米種植面積,相較于隨機森林算法而言,卷積神經網絡能更好地提取作物時序信息,5種數據組合分類效果均高于隨機森林的分類效果,總體分類精度均高于85%,κ系數均高于80%,說明使用卷積神經網絡能夠有效地應用于典型地物識別。玉米提取的F1分數除采用VV數據的分類方式外,均能達到85%以上,說明卷積神經網絡能有效提取玉米種植面積,滿足精細農業種植面積提取的需求,證明卷積神經網絡在作物分類等應用中具有巨大的潛力。從數據上看,光學遙感數據分類效果比SAR后向散射系數數據分類效果更好。光學遙感數據對作物提取更有效,水體的后向散射系數與其他地物后向散射系數差異明顯,故雷達影像數據對水體的提取效果更好。另外,在光學數據缺失的情況下,采用雷達后向散射時序數據也能達到較好的識別精度,對玉米的提取效果滿足精細農業種植面積提取的需求,說明雷達影像在玉米面積提取上具備可分性。融合數據分類精度較單一數據分類精度更高,NDVI數據分類結果有一定的錯分現象,含有后向散射系數數據的分類結果椒鹽現象較明顯,這是因為雷達影像存在斑點噪聲,單個像元受噪聲影響較大,會影響時序曲線的構建。但是,SAR數據的地物后向散射特征不同于光學影像,其穿透性不僅能獲取植被的表面信息,對植被的莖、枝、葉等信息也有一定的反應,結合光學影像的光譜特性,能增加作物種植區的識別精度。結合NDVI時序數據以及后向散射系數數據椒鹽現象及錯分現象都得到明顯的改善,降低將建筑錯分為玉米的現象,有利于玉米的遙感提取與制圖。綜上,數據融合在一定程度上可以提高作物可分性,對提高作物分類精度有一定的效果。

參考文獻:

[1]李 俐,孔慶玲,王鵬新,等. 基于時間序列Sentinel-1A數據的玉米種植面積監測研究[J]. 資源科學,2018,40(8):1608-1621.

[2]Sonobe R,Yamaya Y,Tani H,et al. Assessing the suitability of data from Sentinel-1A and 2A for crop classification[J]. GIScience & Remote Sensing,2017,54(6):918-938.

[4]潘 力,夏浩銘,王瑞萌,等. 基于Google Earth Engine的淮河流域越冬作物種植面積制圖[J]. 農業工程學報,2021,37(18):211-218.

[4]魏鵬飛,徐新剛,楊貴軍,等. 基于多時相影像植被指數變化特征的作物遙感分類[J]. 中國農業科技導報,2019,21(2):54-61.

[5]黃啟廳,曾志康,謝國雪,等. 基于高時空分辨率遙感數據協同的作物種植結構調查[J]. 南方農業學報,2017,48(3):552-560.

[6]羅 明,陸 洲,徐飛飛,等. 基于快速設定決策閾值的大范圍作物種植分布的遙感監測研究[J]. 中國農業資源與區劃,2019,40(6):27-33.

[7]王庚澤,靳海亮,顧曉鶴,等. 基于改進分離閾值特征優選的秋季作物遙感分類[J]. 農業機械學報,2021,52(2):199-210.

[8]Zou Q,Ni L H,Zhang T,et al. Deep learning based feature selection for remote sensing scene classification[J]. IEEE Geoscience and Remote Sensing Letters,2015,12(11):2321-2325.

[9]Kussul N,Lavreniuk M,Skakun S,et al. Deep learning classification of land cover and crop types using remote sensing data[J]. IEEE Geoscience and Remote Sensing Letters,2017,14(5):778-782.

[10]Garnot V S F,Landrieu L,Giordano S,et al. Time-space tradeoff in deep learning models for crop classification on satellite multi-spectral image time series[C]//2019 IEEE International Geoscience and Remote Sensing Symposium.Yokohama,2019:6247-6250.

[11]朱鳳敏,吳 迪,楊佳琪. 基于Sentinel-1B SAR數據的農作物分類方法研究[J]. 測繪與空間地理信息,2020,43(5):105-108.

[12]謝新喬,楊繼周,鄧邵文,等. 多時相Sentinel-1影像反演玉溪典型煙區烤煙種植分布的方法[J]. 農業資源與環境學報,2023,40(1):188-195.

[13]Ndikumana E,Minh D H T,Baghdadi N,et al. Deep recurrent neural network for agricultural classification using multitemporal SAR sentinel-1 for Camargue,France[J]. Remote Sensing,2018,10(8):1217.

[14]郭 交,朱 琳,靳 標. 基于Sentinel-1和Sentinel-2數據融合的農作物分類[J]. 農業機械學報,2018,49(4):192-198. [HJ2mm]

[15]成科揚,榮 蘭,蔣森林,等. 基于深度學習的遙感圖像超分辨率重建方法綜述[J]. 鄭州大學學報(工學版),2022,43(5):8-16.

[16]姜伊蘭,陳保旺,黃玉芳,等. 基于Google Earth Engine和NDVI時序差異指數的作物種植區提取[J]. 地球信息科學學報,2021,23(5):938-947.

[17]張 淼,吳炳方,于名召,等. 未種植耕地動態變化遙感識別——以阿根廷為例[J]. 遙感學報,2015,19(4):550-559.

[18]楊澤航,王 文,鮑健雄. 融合多源遙感數據的黑河中游地區生長季早期作物識別[J]. 地球信息科學學報,2022,24(5):996-1008.

[19]陳思寧,趙艷霞,申雙和. 基于波譜分析技術的遙感作物分類方法[J]. 農業工程學報,2012,28(5):154-160.

[20]趙子娟,劉 東,杭中橋. 作物遙感識別方法研究現狀及展望[J]. 江蘇農業科學,2019,47(16):45-51.

[21]羅榮輝,袁 航,鐘發海,等. 基于卷積神經網絡的道路擁堵識別研究[J]. 鄭州大學學報(工學版),2019,40(2):18-22.

[22]Kampffmeyer M,Salberg A B,Jenssen R. Semantic segmentation of small objects and modeling of uncertainty in urban remote sensing images using deep convolutional neural networks[C]//2016 IEEE Conference on Computer Vision and Pattern Recognition Workshops.Las Vegas,2016:680-688.

[23]屈 煬,袁占良,趙文智,等. 基于多時序特征和卷積神經網絡的農作物分類[J]. 遙感技術與應用,2021,36(2):304-313.

[24]Zhong L H,Hu L N,Zhou H. Deep learning based multi-temporal crop classification[J]. Remote Sensing of Environment,2019,221:430-443.

[25]王小慧,姜雨林,傅漫琪,等. 海河低平原典型縣種植制度與農田景觀格局變化遙感監測[J]. 農業工程學報,2022,38(1):297-304.

猜你喜歡
深度學習
從合坐走向合學:淺議新學習模式的構建
面向大數據遠程開放實驗平臺構建研究
基于自動智能分類器的圖書館亂架圖書檢測
搭建深度學習的三級階梯
有體驗的學習才是有意義的學習
電子商務中基于深度學習的虛假交易識別研究
現代情報(2016年10期)2016-12-15 11:50:53
利用網絡技術促進學生深度學習的幾大策略
考試周刊(2016年94期)2016-12-12 12:15:04
MOOC與翻轉課堂融合的深度學習場域建構
大數據技術在反恐怖主義中的應用展望
深度學習算法應用于巖石圖像處理的可行性研究
軟件導刊(2016年9期)2016-11-07 22:20:49
主站蜘蛛池模板: 日本免费新一区视频| 九色在线视频导航91| 自拍偷拍欧美日韩| 一级一级一片免费| 99无码中文字幕视频| 一级毛片免费高清视频| 欧美一区日韩一区中文字幕页| 高清国产在线| 午夜啪啪福利| 久久久久无码精品国产免费| jizz在线免费播放| 伊人久久大香线蕉影院| 亚洲欧美成人在线视频| 亚洲永久色| 国产成人免费高清AⅤ| 粗大猛烈进出高潮视频无码| 最新国产麻豆aⅴ精品无| 色哟哟精品无码网站在线播放视频| 国产成人综合日韩精品无码不卡| 日韩无码一二三区| 国产综合色在线视频播放线视| 国产人前露出系列视频| 国产在线观看第二页| 国产欧美日韩视频怡春院| 91国语视频| 亚洲国产成人自拍| 亚洲无码精彩视频在线观看| 亚洲色成人www在线观看| 久久综合色播五月男人的天堂| 人妻一区二区三区无码精品一区| 欧美成人在线免费| 狠狠干综合| 91久久夜色精品| 亚洲一区二区三区国产精品| vvvv98国产成人综合青青| 国产精品99r8在线观看| 蜜桃视频一区二区| 在线免费观看AV| 中文无码精品A∨在线观看不卡 | 久久动漫精品| 免费一看一级毛片| 少妇极品熟妇人妻专区视频| 亚洲中字无码AV电影在线观看| 2021国产v亚洲v天堂无码| 亚洲国产日韩在线观看| 天天躁日日躁狠狠躁中文字幕| 午夜国产理论| 亚洲第一精品福利| 爱色欧美亚洲综合图区| 日本高清视频在线www色| 无码国内精品人妻少妇蜜桃视频| 欧美国产日本高清不卡| 91麻豆久久久| 性做久久久久久久免费看| 狠狠亚洲五月天| 欧美在线精品怡红院| 亚洲综合专区| 青青青国产免费线在| 嫩草国产在线| 国产精品香蕉在线观看不卡| 午夜精品国产自在| 亚洲综合经典在线一区二区| 日韩国产黄色网站| 岛国精品一区免费视频在线观看| 天天摸夜夜操| 69av在线| 91无码视频在线观看| 亚洲欧美成人综合| 无码高潮喷水在线观看| 中文字幕啪啪| 青青国产成人免费精品视频| 精品视频一区二区三区在线播| 天天综合亚洲| 亚洲国产综合精品一区| 成人毛片在线播放| 老司机久久精品视频| 乱人伦中文视频在线观看免费| AV天堂资源福利在线观看| 国产精品亚洲精品爽爽 | 国产精品观看视频免费完整版| 日韩欧美中文在线| а∨天堂一区中文字幕|