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

魯中南山丘區耕地地力的遙感反演模型與應用

2020-03-03 01:00:46李因帥趙庚星李建偉竇家聰范瑞彬
農業工程學報 2020年23期
關鍵詞:耕地評價模型

李因帥,張 穎,趙庚星,李 濤,李建偉,竇家聰,范瑞彬

魯中南山丘區耕地地力的遙感反演模型與應用

李因帥1,張 穎1,趙庚星1※,李 濤2,李建偉2,竇家聰3,范瑞彬4

(1. 山東農業大學資源與環境學院土肥資源高效利用國家工程實驗室,泰安 271018;2. 山東省土壤肥料工作總站,濟南 250013;3. 山東省農業技術推廣總站,濟南 250013;4. 山東省招遠市自然資源和規劃局,招遠 265400)

耕地地力是耕地生產能力的重要表征,地力的遙感快速準確反演是耕地資源利用管理的客觀需求。該研究針對魯中南山丘區,選擇東平縣和滕州市2個代表性縣市,利用東平縣的TM影像構建與篩選光譜指標,通過經典統計分析(一元線性回歸、曲線回歸、多元逐步線性回歸)和機器學習(BP神經網絡、極限學習機)方法構建與優選反演模型,進而在滕州市進行模型的驗證和應用,同時對不同時相的反演模型進行了比較分析。結果表明:5類光譜指標與耕地地力綜合指數均有顯著的相關性,其中改進型光譜指數的均>0.684,能更好地反映耕地地力狀況;最佳反演模型為經典統計分析方法中的改進指數組-多元逐步線性回歸(IIG-MLSR)模型(2=0.684,RMSE=5.674)和機器學習算法的改進指數組-BP神經網絡(IIG-BPNN)模型(2=0.746,RMSE=5.089);模型在山丘區具有較好的普適性,地力反演與評價的結果具有相似的空間分布特征、相近的耕地面積比例和較高的空間契合度。其中2個最佳模型的高中低3級耕地的面積比例差普遍低于5.55個百分點,空間契合度分別為84.50%和88.76%;模型動態反演分析結果顯示,2007—2016年滕州市耕地地力不斷提升,高級地由67.30%增加至80.72%;不同時相模型比較結果顯示,多時相遙感反演耕地地力具備可行性,4月份冬小麥返青拔節期是反演的最佳時相,10月份裸土時相次之,8月份夏玉米時相最差。該研究提出了山丘區耕地地力快速定量遙感反演的有效方法,對完善遙感反演指標與模型,提高評價效率有參考價值。

遙感;反演;模型;耕地地力;BP神經網絡;極限學習機;魯中南山丘區

0 引 言

耕地是農業發展的重要生產資料,也是促進社會經濟可持續發展的重要資源。而耕地地力的優劣直接影響到糧食產量、質量和農業的可持續發展[1]。山地丘陵區受自然與人為因素影響,耕地質量差且退化風險大,歷來是學者研究的重點,相關地力研究多采用基于GIS的傳統評價方法[2-3],此方法雖然評價精度高、應用范圍廣,但需消耗大量的人力、物力、財力資源,難以實現地力動態監測和大區域尺度上的快速評價。而遙感技術的應用為地力的快速定量評價提供了新途徑,如何借助遙感提高地力評價效率,監測地力動態變化,充分發揮耕地生產潛力,促進耕地的可持續利用是當前耕地地力評價的重要方向[4]。

近年來,利用遙感進行耕地地力評價的研究大致可分為3類,一是以遙感影像為數據源編制土地利用現狀圖等基礎圖件,提取或更新耕地空間分布信息[5];二是利用遙感影像解譯相關評價因子,構建耕地地力評價模型[6-7];三是利用遙感影像建立反演模型來評估地力狀況。后者作為當前研究熱點受到了較多關注,例如武婕等[8]基于TM影像構建了植被指數—耕地地力的Quadratic模型,實現了縣域級別耕地地力的反演;官炎俊等[9]利用OLI影像NDVI和糧食單產,構建了耕地質量反演模型;Liu等[10]從GF-1中提取與篩選植被指數作為土壤肥力指標,進行了耕地質量反演,其中GA-BPNN模型效果最佳;Xia等[11]比較分析了水稻不同生長期的植被指數模型對耕地質量的預測精度。近年來機器學習算法在耕地地力評價中的應用趨于增多[12-13],其優勢在于較好體現了耕地地力的動態性、隨機性和非線性等特點[14]。總體看,目前遙感反演耕地地力的系統研究尚淺,一是多以植被指數為解釋變量,模型反演參量較少,缺乏對光譜指標的優化選擇;二是反演方法局限,多使用經典統計模型,缺乏對深層關系的挖掘,有必要進一步探索機器學習算法;三是反演模型的時空應用性有待提高。

本文在總結前人研究基礎上,從魯中南山丘區的東平縣的TM影像中提取耕地光譜信息,利用多種分析方法構建與優選反演模型,并選取地力相近的滕州市從多個角度進行模型的精度分析與應用,旨在探討山丘區耕地地力反演的新模式與新方法,完善地力遙感反演體系,提高模型的實用價值,更好為耕地地力評價服務。

1 材料與方法

1.1 研究區概況

魯中南山地丘陵區位于黃河以南,大運河以東,沭河以西,膠濟鐵路以南,山東半島的中南部地區,主要由泰沂蒙山脈構成,周圍環繞著地勢逐漸降低的丘陵。區域內地形復雜,水土易流失,地力狀況差異大且具有不穩定性,如何快速準確的獲取其地力狀況以實現耕地質量保護顯得尤為重要。本文選取小麥、玉米大面積種植且耕地地力狀況相似的東平縣和滕州市作為典型研究區,兩區的具體情況如下:

東平縣位于魯中南山丘區的西部(35°46'24"N~36°10'20"N,116°02'52"E~116°39'44"E)。西臨黃河,東望泰山,地勢北高南低,東高西低,山丘、平原、湖泊面積各占1/3。春秋兩季干燥少雨,夏季高溫多雨,四季分明,雨熱同期。土壤包含2個土類,6個亞類,12個土屬,36個土種。耕地利用方式以水澆地、旱地、灌溉水田、菜地為主,主要農作物為冬小麥和夏玉米,是全國糧油商品生產基地縣。

滕州市位于魯中南山丘區的西南部(34°49'32"N~35°17'21"N,116°48'27"E~117°24'26"E)。北、南、東三面環山,西臨南四湖,地勢由東北向西南傾斜。屬于暖溫帶大陸性季風氣候,光照充足,雨量集中。土壤以褐土和潮土為主,分為5個土類、12個亞類、22個土屬、90個土種。耕地利用以水澆地、旱地、灌溉水田、菜地為主,主要種植小麥、玉米、馬鈴薯等,農業現代化與機械化水平高,有“魯南糧倉”之稱(圖1)。

圖1 研究區位置與遙感影像圖(RGB假彩色合成)

1.2 數據來源與預處理

1.2.1 評價資料收集處理與耕地地力評價

相關評價資料主要來源于2007年份的山東省測土配方施肥項目。采用《耕地地力調查與質量評價技術規程》(NY/T 1634-2008)[15]推薦的方法,首先將土壤圖、土地利用現狀圖疊置,分別得到7 980(東平)和6 764(滕州)個評價單元。進而采用特爾斐法和系統聚類法從立地條件(灌溉保證率、坡度、地貌類型)、物理性狀(耕層質地、土體構型、土層厚度、障礙層次)、化學性狀(有機質、有效磷、速效鉀、有效鋅、有效硼)3個方面進行評價因子的篩選,利用層次分析法確定因子權重,模糊評價法構建隸屬函數,確定隸屬度[2]。最后通過指數和法計算耕地地力綜合指數(Integrated Fertility Index, IFI),為服務于反演建模,本文將其擴大100倍。其中耕地地力等級統一使用等間距法劃分為10個等級。

1.2.2 遙感影像收集處理與耕地信息提取

本文選取東平縣2007年4月3日(冬小麥返青拔節期)的Landsat-TM影像,用于耕地信息提取和反演模型構建;同時選擇2009年4月8日(冬小麥返青拔節期)、8月30日(夏玉米乳熟期)、10月17日(裸地期)的TM影像用于多時相耕地地力遙感反演比較研究。選取滕州市2007年4月3日(冬小麥返青拔節期)的TM影像用于耕地提取和反演模型精度驗證分析,同時選取該市相同作物物候期的TM(2011年3月29日)與OLI(2016年3月26日)影像用于地力動態分析。

利用ENVI5.3軟件進行輻射定標、大氣校正以獲取地表真實反射率,參照地形圖進行幾何精校正以消除幾何畸變,再對影像進行配準和裁剪,得到其遙感影像圖。

耕地信息的提取統一采用ISODATA—人工交互合并及修正法[16]。3月末至4月初研究區的冬小麥正處于返青拔節期,而其他植被尚未返青,是提取耕地光譜信息的最佳時相。利用該時相影像提取兩區的耕地信息進行模型的構建與驗證,進而基于不同年份的返青拔節期影像提取滕州市的耕地信息反演分析地力動態變化,最后使用東平縣小麥玉米輪作耕地(通過將冬小麥與夏玉米影像提取的“植被覆蓋區”疊置獲取)進行多時相的耕地地力反演對比分析。

1.3 研究方法

本研究主要通過TM影像提取耕地光譜信息,利用多種建模方法實現耕地地力的模擬和時空尺度的推廣應用,旨在實現魯中南山丘區高精度、實時定量化的地力評價。技術路線如圖2所示。

圖2 技術路線圖

1.3.1 地力光譜指標的構建與篩選

在總結相關研究成果基礎上,廣泛搜集可能與耕地地力相關的遙感光譜指標,并加以改進構建了改進型光譜指數,主要光譜指標如下:

1)單波段光譜指標

由Landsat-TM多光譜影像的各個波段組成,反映耕地不同譜段的光譜特征。

2)光譜變換指標

①K-L變換(Karhunen-Loeve Transform)

K-L變換可去除波段之間多余信息、將其壓縮到更有效的幾個主分量中。一般前3個主成分包含原始影像絕大部分的光譜信息,是耕地光譜特征信息的增強集中反映[17]。

②K-T變換(Kauth-Thomas Transformation)

K-T變換即纓帽變換,是根據土壤、植被等在多光譜空間中的信息分布結構對圖像做的經驗性正交變化。其變換第1分量為亮度(Brightness Index, BI),反映耕地總體反射率狀況;第2分量為綠度(Greenness Vegetation Index, GVI),常用于地表植被覆蓋度研究;第3分量為濕度(Wetness Index, WI),反映耕地土壤濕度狀況[18]。

3)植被光譜指標

農作物狀況是耕地生產力水平的良好表征,植被指數可以綜合反映與耕地質量有關的農作物長勢與覆蓋度、土壤墑情以及有機質、全氮等土壤養分信息。本研究選取包含相關因子影響的比值植被指數(Ratio Vegetation Index, RVI)、三角形植被指數(Triangular Vegetation Index, TVI),消除綜合因子影響的歸一化植被指數、增強型植被指數(Enhanced Vegetation Index, EVI),消除大氣影響的大氣阻抗植被指數(Atmospherically Resistant Vegetation Index, ARVI),消除土壤影響的修正型土壤調節植被指數(Modified Soil Adjusted Vegetation Index, MSAVI)等參與地力反演[19-24]。

4)水分光譜指標

地表濕度狀況間接反映了農田灌排能力的高低,是影響耕地地力水平的重要因素。有研究表明,相較于可見光與近紅外波段,中紅外波段建立的光譜指數在反映土壤含水量時,參數更穩定,效果更好[25]。因此選取了中紅外波段參與構成的地表含水量指數(Surface Water Capacity Index, SWCI)、修正型地表含水量指數(Modified Surface Water Capacity Index, MSWCI)、陸地表面水分指數(Land Surface Water Index, LSWI)、歸一化差異紅外指數(Normalized Difference Infrared Index, NDII)等水分光譜指數[26-29]。

5)改進型光譜指標

通過對比以上光譜指標,發現植被指數大多由可見光-近紅外范圍的波段組合產生,水分指數則多在中紅外基礎上,引入其他波段通過歸一化組合產生。由此可見,尚缺乏可見光-近紅外-中紅外波段構成的各類不同組合方式的光譜指數。已有研究證明包含該波譜信息的改進指數可以有效降低波段間的信息冗余,提高模型反演效果[30-31]。因此,綜合分析TM波段本身特性和現有光譜指標中常用的歸一化、差值、比值波段組合運算方法構建了改進型光譜指數(表1),以更好反映地力狀況。

表1 改進型光譜指標

注:、、NIR、SWIR1、SWIR2分別為TM影像的B2、B3、B4、B5、B7。

Note:,, NIR, SWIR1and SWIR2are B2, B3, B4, B5 and B7 of TM image respectively.

使用SPSS分析光譜指標與IFI的相關性。通過光譜指標篩選得到變量組,作為模型的輸入變量。

1.3.2 耕地地力反演模型的構建與優選

為便于使用機器學習算法,保證不同建模方法間的可比性,將評價單元按照2:1的比例隨機分為建模集與驗證集。進而采用一元線性回歸(Simple Linear Regression, SLR)、曲線回歸(Curvilinear Regression, CR)、多元逐步線性回歸(Multiple Linear Stepwise Regression, MLSR)3種經典統計分析和BP神經網絡(Back Propagation Neural Networks, BPNN)[32]、極限學習機(Extreme Learning Machine, ELM)[33]2種機器學習算法,以光譜變量組為輸入變量,IFI為輸出變量進行建模。最后通過2與均方根誤差RMSE篩選得到最佳反演模型。

運用SPSS軟件實現統計分析模型的模擬,其中SLR與CR選取各變量組中最大的指數作為輸入變量。運用MATLAB軟件實現機器學習模型的建模與驗證,其中BPNN采用Kolmogorov經驗公式[34],ELM采用“試錯法”[35]輔助確定隱含層節點數。

1.3.3 耕地地力反演模型的精度分析

為驗證反演模型在魯中南等丘陵山地區的空間普適性,將基于東平縣數據構建的反演模型應用于滕州市。一方面,通過各等級耕地的面積比例差來評價反演精度。另一方面,將地力反演與地力評價結果疊加,分析反演結果的空間契合度。為使結果更為直觀,進一步將10級地合并為高(1~3級)、中(4~6級)和低(7~10級)地力等級耕地,并分析滕州市耕地地力的空間分布規律。

1.3.4 耕地地力的遙感監測應用

利用地力反演模型對滕州市不同年份的地力狀況進行預測,按高、中、低3級比較分析該市近10 a的耕地地力動態變化。

1.3.5 不同時相耕地地力反演模型的對比分析

使用東平縣不同作物物候期的TM影像與基于GIS的耕地地力評價結果,對各變量組的指標進行篩選,采用最優建模方法,依據、2與RMSE比較分析不同時相的反演模型,探索多時相遙感定量反演耕地地力的可行性并確定最佳時相。

2 結果與分析

2.1 敏感光譜指標的篩選結果

從相關性分析結果中選擇與耕地地力關系顯著且的絕對值大于0.5的光譜指標。從表2可知,單波段光譜指標中的近紅外波段B4的=0.702,篩選為單波段組敏感指標;光譜變換指標中的PC1、PC2、GVI與IFI相關性較強,篩選為光譜變換組敏感指標;其余光譜指標均與IFI表現出較高的相關性,分別篩選為植被指數組、水分指數組和改進指數組(Improved Index Group,IIG)的敏感指標,將各敏感指標組成的變量組作為輸入變量用于耕地地力反演模型的構建。

表2 光譜指標與耕地地力綜合指數IFI的相關系數

注:PC1、PC2、PC3分別為主成分分析的第1、2、3主分量,BI、GVI、WI分別為亮度、綠度和濕度;ARVI、RVI、EVI、MSAVI、TVI分別為大氣阻抗植被指數、比值植被指數、增強型植被指數、修正型土壤調節植被指數、歸一化植被指數和三角形植被指數;SWCI、MSWCI、LSWI、NDII分別為地表含水量指數、修正型地表含水量指數、陸地表面水分指數和歸一化差異紅外指數。下同。

Note: PC1, PC2 and PC3 are the first, second and third principal components of PCA, and BI, GVI and WI are Brightness index, Greenness vegetation index and Wetness index; ARVI, RVI, EVI, MSAVI and TVI are Atmospherically resistant vegetation index, Ratio vegetation index, Enhanced vegetation index, Modified soil adjusted vegetation index, and Triangular vegetation index; SWCI, MSWCI, LSWI and NDII are Surface water capacity index, Modified surface water capacity index, Land surface water index and Normalized difference infrared index. The same as below.

2.2 耕地地力反演模型

2.2.1 經典統計分析反演模型

表3為各變量組的經典統計分析模型,由于篇幅限制,CR模型僅列出了擬合效果最好的2組。

由表3可知,對比不同變量組,在SLR模型中,單波段組的2<0.5,RMSE>7,預測效果最差。整體上各變量組的反演效果從差至優依次為:單波段組、光譜變換組、植被指數組、水分指數組、改進指數組;在SPSS內嵌的11種CR模型中,各變量組均以Quadratic、Cubic的擬合效果最好,與武婕等[8]的研究結果一致。各變量組的反演效果為:單波段組、光譜變換組、水分指數組、植被指數組、改進指數組;在MLSR模型中,改進指數組的驗證集2=0.684,RMSE=5.674,精度最高。各變量組的反演效果從差至優依次為:單波段組、光譜變換組、植被指數組、水分指數組、改進指數組。對比不同建模方法,各變量組均以SLR模型的擬合效果最差,除植被指數組中Cubic的反演效果優于MLSR外,其他變量組均以MLSR模型的預測精度更高。

經過綜合分析,確定最佳經典統計分析模型為IIG-MLSR模型:

=272.858+25.184RPSI+1137.219DGSI?1114.913DRSI?

184.153RSSI?331.688NDGSI+152.518NDRSI。

2.2.2 機器學習反演模型

在MATLAB中以各變量組為輸入變量,分別構建BPNN和ELM模型,結果如表4所示。

對比不同變量組,在BPNN模型中,單波段組的驗證精度最低,R2=0.526,RMSE=6.950,其余變量組模型的決定系數均大于0.683,其中改進指數組的R2=0.749,R2=0.746,驗證精度最高;在ELM模型中,除單波段組外,其余變量組的2均大于0.6,RMSE均小于6。各變量組的反演效果從差至優依次均為:單波段組、水分指數組、光譜變換組、植被指數組、改進指數組。對比不同建模方法,兩種算法的反演效果相近,除光譜變換組外,BPNN的反演效果均要略優于ELM。因此,機器學習中的最佳反演模型為IIG-BPNN模型。

表3 基于不同變量組的經典統計分析模型

注:SLR、Quadratic、Cubic、MLSR分別為一元線性回歸、二次曲線回歸、三次曲線回歸和多元逐步線性回歸模型。

Note: SLR, Quadric, Cubic and MLSR are simple linear regression, quadratic curve regression, cubic curve regression and multiple linear stepwise regression models.

表4 基于不同變量組的機器學習模型

注:BPNN、ELM分別為BP神經網絡和極限學習機模型。

Note: BPNN and ELM are back propagation neural networks and extreme learning machine models.

綜上,在5類變量組中,改進指數組的光譜信息量最豐富,對地力狀況的擬合效果也最好;在5類建模方法中,機器學習的擬合精度明顯高于經典統計分析。如圖3所示,在改進指數組模型中SLR的預測精度最低,Cubic優于Quadratic,MLSR的2=0.684,RMSE=5.674,反演效果較好,而BPNN與ELM模型的預測精度顯著優于上述模型。

圖3 改進指數組驗證集的預測值與真實值散點圖

Fig.3 Scatterplot of predicted values and true values of improved index group verification set

最終經過篩選得到的最佳地力反演模型分別為經典統計分析方法中的IIG-MLSR模型和機器學習算法中的IIG-BPNN模型。

2.3 耕地地力模型的反演精度分析

2.3.1 耕地地力的反演面積精度分析

表5為滕州市地力評價與地力反演結果,可見其耕地的地力等級集中于2級附近,整體地力水平較高。IIG-MLSR模型反演的5~10級地面積比例較接近,面積差均低于2.72個百分點,1~4級地的面積比例差別較大,面積差均大于5.06個百分點。在IIG-BPNN模型中,2、3、5級地的反演面積比例差別較大,面積差普遍大于4.19個百分點。將10個等級歸納劃分為高、中、低3個等級時,IIG-MLSR模型的面積差普遍低于5.55個百分點,IIG-BPNN模型的預測結果更佳,面積差均低于2.64個百分點。2種反演模型預測的地力等級面積比例與地力評價的結果基本一致,證明反演模型適用于山丘區耕地地力的預測。

2.3.2 耕地地力的空間分布反演分析

由圖4可得,地力評價與地力反演的各等級耕地在空間分布上具有較好的一致性,高等級耕地(1~3級)主要分布在市域中部和西南部,中等級耕地(4~6級)主要分布在西部和東南部,而低等級耕地(7~10級)則主要分布在北部及東南部地區。所得的地力空間分布規律亦與張立文等[36]對滕州市冬小麥種植適宜區的劃分結果基本一致,說明了模型具有較好的空間普適性,遙感定量反演耕地地力具備可行性。

表5 耕地地力的反演面積比例分析

圖4 地力評價與地力反演的耕地地力空間分布圖

2.3.3 耕地地力的空間契合度反演分析

分別將地力反演與地力評價的結果進行疊置,得到各耕地單元間的等級差并進行面積統計(圖5)。由圖可得,劃分為10等級時,IIG-MLSR與IIG-BPNN反演模型中等級一致的耕地面積分別占48.14%和53.88%,相差一級的分別占40.60%和39.29%,絕大多數耕地的地力反演與評價結果的等級差別較小。另外當合并為高、中、低3級時,IIG-MLSR模型中等級一致的占84.50%,IIG-BPNN模型為88.76%,證明2種模型在山丘區具有較高的空間普適性。

圖5 反演模型的空間契合度分析圖

2.4 耕地地力的遙感監測結果

采用IIG-MLSR模型反演分析滕州市2007-2016年的耕地地力狀況(圖6)。

圖6 滕州市耕地地力動態反演結果

該市的高級地面積比例由67.30%(2007年)增加至75.96%(2011年),最后增加到80.72%(2016年);中級地面積比例不斷縮減,由2007年的27.97%逐漸減少到2016年的16.86%;低級地的面積比例則先由2007年的4.73%減少到2011年的1.38%,而后在2016年又增加到2.42%。總體看,近10 a滕州市的耕地地力水平呈提高趨勢,其中高級地的面積比例不斷增加,中級地的面積比例不斷縮減,而低級地的面積比例則相對穩定。

2.5 不同時相的耕地地力反演模型對比結果

由冬小麥、夏玉米和裸地3個不同時相的反演模型比較(表6)可以看出,基于冬小麥時相(4月8日)建立的MLSR模型的驗證集的=0.814,2=0.663,RMSE=5.859,預測效果較好。BPNN模型的驗證集=0.843,2=0.711,RMSE=5.431,預測效果更佳;裸地時相(10月17日)的MLSR反演模型中R2=0.558,R2=0.535,BPNN模型R2=0.639,R2=0.617,反演精度均低于冬小麥時相。而夏玉米時相(8月30日)所得模型的預測效果最差,但相關系數仍可達0.724以上,亦具備一定的反演能力。因此,利用多時相遙感數據反演耕地地力具備可行性,其中4月份冬小麥返青拔節期為最佳時相,10月份裸土時相次之,8月份夏玉米時相最差。

表6 耕地地力的多時相遙感反演結果

3 討 論

光譜指標的選擇是進行耕地地力反演的前提,確定準確反映耕地地力的遙感指標是確保反演精度的關鍵。前人在耕地質量的研究中多選用植被指數[8-11]作為反演模型的解釋變量,僅以地表植被的長勢與覆蓋狀況估測耕地地力的高低,難以精準反映地力狀況。本研究在選取植被指數的同時,引入并構建了與土壤水肥相關的其它光譜指標,研究發現,所選指標不僅能用來反映耕地質量要素(方琳娜等[37]、Liu等[38]、楊建鋒等[6]),且能較好間接表征整體的耕地質量狀況。多類型的遙感光譜指標豐富了地力反演的可用參量,對耕地質量光譜估測有參考意義。

耕地地力具有動態性、隨機性、非線性、空間變異性等特點[14],難以用簡單的線性模型來描述。現有的研究中,經典統計方法仍占主導[8-9,11],缺乏對耕地地力與遙感光譜指標間關系的深入挖掘。本研究在經典統計分析模型研究的基礎上同時采用了機器學習算法(BPNN與ELM),結果各指標機器學習模型的反演效果普遍優于經典統計分析模型。這與Liu等[10]在耕地質量反演、張智韜等[39]在鹽分反演中的研究結論基本一致,說明機器學習算法相較于經典統計分析方法表現出了強大的非線性擬合能力和出色的數據挖掘能力,能更好地模擬耕地地力與遙感光譜信息之間的復雜多元非線性關系,實現區域耕地地力更高精度的反演。

快速、高效地實現耕地地力評價結果的更新,長期監測耕地地力的動態變化已成為當前耕地資源利用管理的客觀需求[40],因此要求所建立的反演模型具備一定的動態監測能力,而以往的研究[8-11]中較少顧及模型的動態性,使模型的普適性較差。本研究將反演模型應用于耕地地力的時空動態監測,發現滕州市的耕地質量為提升趨勢。據統計資料顯示,2007、2011與2016年該市的冬小麥單產分別為6 814.0、7 467.0和7 717.1 kg/hm2,呈增加趨勢,整體上與遙感反演的結果一致,間接證明模型適用于耕地地力動態變化監測,顯示了魯中南山丘區耕地地力的動態趨向,亦為區域的耕地地力動態監測管理提供參考。

本研究利用多時相遙感反演耕地地力,相較于單時相遙感影像的反演模型[8-9],探索了在不同地表覆被狀態下遙感反演耕地地力的可行性。由于不同作物在不同生長期長勢、覆蓋度的變化,造成其光譜響應的差異性[11,41],使得多時相的地力反演模型的預測精度存在差異。本研究顯示冬小麥是耕地地力的良好指示作物,而夏玉米乳熟期的模型反演效果不夠理想。后續可進一步細分農作物進行耕地地力反演,并利用不同作物生長期的衛星數據進行模型修正,從而進一步提升模型反演的準確性。

基于主要糧食作物小麥和玉米,選擇魯中南山丘區的代表縣區構建地力反演模型,經地力評價與反演結果的對比,二者具有較高的一致性,顯示了模型在山丘區較好的普適性。但所得模型可能難以準確反映其他地形區的耕地地力狀況。袁秀杰等[42]研究發現使用線性回歸模型可實現平原與丘陵區地力評價結果的銜接與轉換,因此,實現不同地形區間模型的轉換與修正,進而建立起一套囊括各類地形區、覆蓋各類作物的耕地地力遙感反演體系將是下一步的研究重點。有必要考慮結合實地調查數據對氣象與病蟲災害、作物品種及農田管理差異等外界因素的影響,進一步完善模型,提高反演精度。

4 結 論

本文以魯中南山丘區的東平縣與滕州市為研究區,通過經典統計分析和機器學習算法構建并篩選光譜指標變量組—耕地地力反演模型,并從時空尺度進行模型的精度驗證和應用分析,得到以下結論:

1)所選5類光譜指標均能反映耕地地力狀況,與耕地地力綜合指數IFI均具有顯著的相關性。其中改進型光譜指數可以更有效的提取與凝聚地力相關光譜信息,與IFI的相關系數均大于0.684。

2)耕地地力最佳反演模型為改進指數組-多元逐步線性回歸(IIG-MLSR)模型和改進指數組-BP神經網絡(IIG-BPNN)模型,模型驗證集的2分別為0.684和0.746,RMSE分別為5.674和5.089。以滕州為驗證區的精度分析結果顯示:地力反演與實際評價的各等級耕地的面積比例相近,十等級的面積差均合并為高中低3級時均低于5.55個百分點;MLSR與BPNN模型在劃分為高中低3級時的空間契合度分別為84.50%和88.76%,模型在魯中南山丘區具備較好的空間普適性。

3)根據耕地地力的動態反演分析,近10年滕州市耕地質量呈提升趨勢,高等級耕地由2007年的67.30%增加至2011年的75.96%,2016年增加到80.72%,同時中等級地的面積比例不斷縮減,低級地則相對穩定。

4)利用多時相遙感數據反演耕地地力具備可行性,其中4月份冬小麥返青拔節期為最佳時相,模型驗證集的均大于0.814,2均大于0.663,RMSE均小于5.859,10月份裸土時相次之,8月份夏玉米時相最差。

[1] 毛雪,孟源思,張東紅,等. 基于GIS的皖江流域耕地地力評價研究[J]. 中國農業資源與區劃,2019,40(7):110-118.

Mao Xue, Meng Yuansi, Zhang Donghong, et al. Evaluation of cultivated land capacity in the Yangtze River Basin in Anhui Province based on GIS[J]. Chinese Journal of Agricultural Resources and Regional Planning, 2019, 40(7): 110-118. (in Chinese with English abstract)

[2] 王瑞燕,趙庚星,李濤,等. GIS支持下的耕地地力等級評價[J]. 農業工程學報,2004,20(1):307-310.

Wang Ruiyan, Zhao Gengxing, Li Tao, et al. GIS supported quantitative evaluation of cultivated land fertility[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2004, 20(1): 307-310. (in Chinese with English abstract)

[3] 周俊,楊子凡. 高臺縣耕地地力評價[J]. 中國農業資源與區劃,2018,39(6):74-78.

Zhou Jun, Yang Zifan. Evaluation of cultivated land force in Gaotai County[J]. Chinese Journal of Agricultural Resources and Regional Planning, 2018, 39(6): 74-78. (in Chinese with English abstract)

[4] 李婷,吳克寧. 基于遙感技術的耕地質量評價研究進展與展望[J]. 江蘇農業科學,2018,46(15):5-9.

Li Ting, Wu Kening. Research progress and prospect of cultivated land quality evaluation based on remote sensing technology[J]. Jiangsu Agricultural Sciences, 2018, 46(15): 5-9. (in Chinese with English abstract)

[5] 張海濤,周勇,汪善勤,等. 利用GIS和RS資料及層次分析法綜合評價江漢平原后湖地區耕地自然地力[J]. 農業工程學報,2003,19(2):219-223.

Zhang Haitao, Zhou Yong, Wang Shanqin, et al. Natural productivity evaluation of cultivated land based on GIS and RS data in Houhu Farm of Jianghan Plain[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2003, 19(2): 219-223. (in Chinese with English abstract)

[6] 楊建鋒,馬軍成,王令超. 基于多光譜遙感的耕地等別識別評價因素研究[J]. 農業工程學報,2012,28(17):230-236.

Yang Jianfeng, Ma Juncheng, Wang Lingchao. Evaluation factors for cultivated land grade identification based on multi-spectral remote sensing[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2012, 28(17): 230-236. (in Chinese with English abstract)

[7] Wang Z, Wang L M, Xu R N, et al. GIS and RS based assessment of cultivated land quality of Shandong Province[J]. Procedia Environmental Sciences, 2012, 12(4): 823-830.

[8] 武婕,李增兵,李玉環,等. 基于TM影像植被指數的耕地地力狀況反演研究[J]. 自然資源學報,2015,30(6):1035-1046.

Wu Jie, Li Zengbing, Li Yuhuan, et al. Arable land fertility inversion based on vegetation index from TM image[J]. Journal of Natural Resources, 2015, 30(6): 1035-1046. (in Chinese with English abstract)

[9] 官炎俊,鄒自力,張曉平,等. 基于歸一化植被指數的耕地質量反演模型研究[J]. 土壤通報,2018,49(4):779-787.

Guan Yanjun, Zou Zili, Zhang Xiaoping, et al. Research on the inversion model of cultivated land quality based on normalized difference vegetation index[J]. Chinese Journal of Soil Science, 2018, 49(4): 779-787. (in Chinese with English abstract)

[10] Liu S S, Peng Y P, Xia Z Q, et al. The GA-BPNN-based evaluation of cultivated land quality in the PSR framework using Gaofen-1 satellite data[J]. Sensors, 2019, 19(23): 5127-5140.

[11] Xia Z P, Peng Y P, Liu S S, et al. The optimal image date selection for evaluating cultivated land quality based on Gaofen-1 images[J]. Sensors, 2019, 19(22): 4937-4951.

[12] 王瑞燕,趙庚星,陳麗麗. 基于ANN—產量的耕地地力定量評價模型及其應用[J]. 農業工程學報,2008,24(1):113-118.

Wang Ruiyan, Zhao Gengxing, Chen Lili. Evaluation model of cultivated land productivity using artificial neural network and productivity and its application[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2008, 24(1): 113-118. (in Chinese with English abstract)

[13] Liu Y, Wang H F, Zhang H, et al. A comprehensive support vector machine-based classification model for soil quality assessment[J]. Soil and Tillage Research, 2016(155): 19-26.

[14] 閆一凡,劉建立,張佳寶. 耕地地力評價方法及模型分析[J]. 農業工程學報,2014,30(5):204-210.

Yan Yifan, Liu Jianli, Zhang Jiabao. Evaluation method and model analysis for productivity of cultivated land[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2014, 30(5): 204-210. (in Chinese with English abstract)

[15] 中華人民共和國農業部. 耕地地力調查與質量評價技術規程(NY/T1634-2008)[S]. 北京:中國農業出版社,2008.

[16] 張銀輝,趙庚星,趙文武. 縣級耕地遙感動態監測方法研究[J]. 測繪通報,2002(3):25-27.

Zhang Yinhui, Zhao Gengxing, Zhao Wenwu. Methods of cultivated land dynamic monitoring by remote sensing technology in counties[J]. Bulletin of Surveying and Mapping, 2002(3): 25-27. (in Chinese with English abstract)

[17] 鄧勁松,李君,王珂. 基于多時相PCA光譜增強和多源光譜分類器的SPOT影像土地利用變化檢測[J]. 光譜學與光譜分析,2009,29(6):1627-1631.

Deng Jinsong, Li Jun, Wang Ke. Detecting land use change using PCA-enhancement and multi-source classifier from SPOT images[J]. Spectroscopy and Spectral Analysis, 2009, 29(6): 1627-1631. (in Chinese with English abstract)

[18] 劉煥軍,楊昊軒,徐夢園,等. 基于裸土期多時相遙感影像特征及最大似然法的土壤分類[J]. 農業工程學報,2018,34(14):132-139.

Liu Huanjun, Yang Haoxuan, Xu Mengyuan, et al. Soil classification based on maximum likelihood method and features of multi-temporal remote sensing images in bare soil period[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2018, 34(14): 132-139. (in Chinese with English abstract)

[19] Jordan C F. Derivation of leaf-area index from quality of light on the forest floor[J]. Ecology, 1969, 50(4): 663-666.

[20] Broge N H, Mortensen J V. Deriving green crop area index and canopy chlorophyll density of winter wheat from spectral reflectance data[J]. Remote Sensing of Environment, 2002, 81(1): 45-57.

[21] Rouse J W, Haas R H, Schell J A, et al. Monitoring the vernal advancements and retrogradation of natural vegetation. NASA/GSFC, Type III, Final Report[R]. Greenbelt, MD, USA, 1974.

[22] Huete A R. A Soil-Adjusted Vegetation Index (SAVI)[J]. Remote Sensing of Environment, 1988, 25(3): 295-309.

[23] Kaufman Y J, Tanre D. Atmospherically Resistant Vegetation Index (ARVI) for EOS-MODIS[J]. IEEE Transactions on Geoscience and Remote Sensing, 1992, 30(2): 261-270.

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

[25] Lobell D B, Asner G P. Moisture effects on soil reflectance[J]. Soil Science Society of America Journal, 2002, 66(3): 722-727.

[26] 杜曉,王世新,周藝,等. 一種新的基于MODIS的地表含水量模型構造與驗證[J]. 武漢大學學報:信息科學版,2007,32(3):205-207.

Du Xiao, Wang Shixin, Zhou Yi, et al. Construction and validation of a new model for unified surface water capacity based on MODIS data[J]. Geomatics and Information Science of Wuhan University, 2007, 32(3): 205-207. (in Chinese with English abstract)

[27] 程若曦. 基于改進SWCI指數的農田土壤水分遙感估算及在干旱監測中的應用:以江漢平原為例[D].武漢:湖北大學,2016.

Cheng Ruoxi. Soil Moisture Monitoring of Farmland Based on Modified SWCI Index: A Case Study on Jianghan Plain[D]. Wuhan: Hubei University, 2016. (in Chinese with English abstract)

[28] Xiao X M, Hollinger D, Aber J, et al. Satellite-based modeling of gross primary production in an evergreen needleleaf forest[J]. Remote Sensing of Environment, 2004, 89(4): 519-534.

[29] Van Wagtendonk J W, Root R R, Key C H. Comparison of AVIRIS and Landsat ETM+ detection capabilities for burn severity[J]. Remote Sensing of Environment, 2004, 92(3): 397-408.

[30] 陳紅艷,趙庚星,陳敬春,等. 基于改進植被指數的黃河口區鹽漬土鹽分遙感反演[J]. 農業工程學報,2015,31(5):107-114.

Chen Hongyan, Zhao Gengxing, Chen Jingchun, et al. Remote sensing inversion of saline soil salinity based on modified vegetation index in estuary area of Yellow River[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2015, 31(5): 107-114. (in Chinese with English abstract)

[31] 蔡亮紅,丁建麗. 基于改進植被指數土壤水分遙感反演[J]. 干旱區地理,2017,40(6):1248-1255.

Cai Lianghong, Ding Jianli. Remote sensing inversion of soil moisture based on modified vegetation index[J]. Arid Land Geography, 2017, 40(6): 1248-1255. (in Chinese with English abstract)

[32] Mouazen A M, Kuang B, Baerdemaeker J D, et al. Comparison among principal component, partial least squares and back propagation neural network analyses for accuracy of measurement of selected soil properties with visible and near infrared spectroscopy[J]. Geoderma, 2010, 158(1/2): 23-31.

[33] Huang G B, Zhu Q Y, Siew C K. Extreme learning machine: Theory and applications[J]. Neurocomputing, 2006, 70(1/2/3): 489-501.

[34] 張佩,陳鄭盟,劉春偉,等. 冬小麥產量結構要素預報方法[J]. 農業工程學報,2020,36(8):78-87.

Zhang Pei, Chen Zhengmeng, Liu Chunwei, et al. Method for the prediction of wheat yield components[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2020, 36(8): 78-87. (in Chinese with English abstract)

[35] 楊一,張淑娟,何勇. 基于ELM和可見/近紅外光譜的鮮棗動態分類檢測[J]. 光譜學與光譜分析,2015,35(7):1870-1874.

Yang Yi, Zhang Shujuan, He Yong. Dynamic detection of fresh jujube based on ELM and visible/near infrared spectra[J]. Spectroscopy and Spectral Analysis, 2015, 35(7): 1870-1874. (in Chinese with English abstract)

[36] 張立文,黃燕玲,趙淑芳. 滕州市冬小麥精細化農業區劃研究[J]. 中國農業信息,2017(3):26-30.

Zhang Liwen, Huang Yanling, Zhao Shufang. Study on the regionalization of fine winter wheat agriculture in Tengzhou City[J]. China Agricultural Informatics, 2017(3): 26-30. (in Chinese with English abstract)

[37] 方琳娜,宋金平. 基于SPOT 多光譜影像的耕地質量評價:以山東省即墨市為例[J]. 地理科學進展,2008,27(5):71-78.

Fang Linna, Song Jinping. Cultivated land quality assessment based on SPOT multispectral remote sensing image: A case study in Jimo City of Shandong Province[J]. Progress in Geography, 2008, 27(5): 71-78. (in Chinese with English abstract)

[38] Liu Y S, Zhang Y Y, Guo L Y. Towards realistic assessment of cultivated land quality in an ecologically fragile environment: A satellite imagery-based approach[J]. Applied Geography. 2010, 30(2): 271-281.

[39] 張智韜,魏廣飛,姚志華,等. 基于無人機多光譜遙感的土壤含鹽量反演模型研究[J]. 農業機械學報,2019,50(12):151-160.

Zhang Zhitao, Wei Guangfei, Yao Zhihua, et al. Soil salt inversion model based on UAV multispectral remote sensing[J]. Transactions of the Chinese Society for Agricultural Machinery, 2019, 50(12): 151-160. (in Chinese with English abstract)

[40] 陳延華,王樂,張淑香,等. 我國褐土耕地質量的演變及對生產力的影響[J]. 中國農業科學,2019,52(24):4540-4554.

Chen Yanhua, Wang Le, Zhang Shuxiang, et al. Quality change of cinnamon soil cultivated land and its effect on soil productivity[J]. Scientia Agricultura Sinica, 2019, 52(24): 4540-4554. (in Chinese with English abstract)

[41] Kattenborn T, Fassnacht F E, Schmidtlein S. Differentiating plant functional types using re?ectance: Which traits make the difference?[J]. Remote Sensing in Ecology and Conservation. 2019, 5(1): 5-19.

[42] 袁秀杰,趙庚星,朱雪欣. 平原和丘陵區耕地地力評價及其指標體系銜接研究[J]. 農業工程學報,2008,24(7):65-71.

Yuan Xiujie, Zhao Gengxing, Zhu Xuexin. Linkage of evaluation index system for cultivated land fertility evaluation in plain and hill regions[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2008, 24(7): 65-71. (in Chinese with English abstract)

Remote sensing inversion and application for soil fertility of cultivated land in the hilly areas of central-south Shandong of China

Li Yinshuai1, Zhang Ying1, Zhao Gengxing1※, Li Tao2, Li Jianwei2, Dou Jiacong3, Fan Ruibin4

(1.,271018,; 2.250013,; 3.250013,; 4.265400,)

Soil fertility of a cultivated land is an important indicator of cultivated land productivity. It is necessary to obtain the rapid and accurate inversion of cultivated land fertility via remote sensing for the better utilization and management of land resource. In this study, a new inversion model was constructed and optimized using the classical statistical analysis (SLR, CR, and MLSR), and machine learning (BPNN and ELM). An effective way was also proposed for rapid quantitative remote sensing inversion of cultivated land fertility in hilly areas. Dongping County and Tengzhou City were selected as two representative counties and cities in the hilly area of center southern Shandong Province, China. In Dongping County, the TM image during the turning green and jointing stage was used to construct and screen spectral indexes of cultivated land fertility. Tengzhou City was selected to verify the spatial universality of inversion model for the soil fertility of a cultivated land. Furthermore, the remote sensing inversion model was used to quantitatively monitor the spatial-temporal dynamic status of cultivated land fertility in Tengzhou City in 2007, 2011, and 2016. The prediction accuracy of inversion models was compared in different periods. The results showed that there were significant correlations between the five kinds of spectral indexes in a remote sensing and the Integrated Fertility Index (IFI), among which the correlation coefficients of improved spectral index were greater than 0.684, indicating better reflecting the status of cultivated land fertility. The best inversion model was the IIG-MLSR model (R=0.684, RMSE=5.674) in the classical statistical analysis, while, theIIG-BPNN model (R=0.746, RMSE=5.089) in the machine learning. The obtained model demonstrated excellent universal applicability in hilly areas, where there were similar spatial distribution characteristics between the inversion and evaluation on the cultivated land fertility, and the similar proportion of cultivated land and high spatial compatibility. In the two best models, the difference in the area ratio of the high, middle, and low levels of cultivated land fertility inversion and cultivated land fertility evaluation was generally less than 5.55 percentage point, where the spatial fit was 84.50% and 88.76%, respectively. The dynamic inversion analysis showed that the cultivated land fertility of Tengzhou City increased continuously in recent 10 years (from 2007 to 2016). The area proportion of high-level land increased from 67.30% to 80.72%, whereas, that of middle-level and low-level land decreased. The multi-temporal remote sensing inversion of cultivated land fertility was feasible, compared with the remote sensing inversion models in different time periods. The optimal time phase to invert the cultivated land fertility was the turning green and jointing stage of winter wheat in April, followed by bare soil in October, and the worst in summer maize in August. The remote sensing inversion index and model can be used to effectively increase the evaluation efficiency of cultivated land fertility. At the same time, this finding can provide a positive reference for the related research of cultivated land quality.

remote sensing; inversion; models; cultivated land fertility; back propagation neural networks; extreme learning machine; hilly area of center-south Shandong Province

李因帥,張穎,趙庚星,等. 魯中南山丘區耕地地力的遙感反演模型與應用[J]. 農業工程學報,2020,36(23):269-278.doi:10.11975/j.issn.1002-6819.2020.23.031 http://www.tcsae.org

Li Yinshuai, Zhang Ying, Zhao Gengxing, et al. Remote sensing inversion and application for soil fertility of cultivated land in the hilly areas of central-south Shandong of China[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2020, 36(23): 269-278. (in Chinese with English abstract) doi:10.11975/j.issn.1002-6819.2020.23.031 http://www.tcsae.org

2020-09-29

2020-11-08

國家自然科學基金(41877003);山東省重大科技創新工程項目(2019JZZY010724);山東省“雙一流”獎補資金(SYL2017XTTD02)

李因帥,主要研究方向為土地資源與信息。Email:sdauzhlys@163.com

趙庚星,教授,博士生導師。主要研究方向為土地資源、遙感及信息技術應用。Email:zhaogx@sdau.edu.cn

10.11975/j.issn.1002-6819.2020.23.031

S127

A

1002-6819(2020)-23-0269-10

猜你喜歡
耕地評價模型
一半模型
自然資源部:加強黑土耕地保護
我國將加快制定耕地保護法
今日農業(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
SBR改性瀝青的穩定性評價
石油瀝青(2021年4期)2021-10-14 08:50:44
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
3D打印中的模型分割與打包
基于Moodle的學習評價
主站蜘蛛池模板: 免费在线成人网| 九九精品在线观看| 亚洲精品成人片在线观看 | 中文无码精品a∨在线观看| 亚洲欧美日韩另类| 女同久久精品国产99国| 国产高清国内精品福利| 国产99在线观看| 国产高清国内精品福利| 999在线免费视频| 国产微拍一区| 亚洲成人精品久久| 丁香五月激情图片| 色久综合在线| 欧美国产日产一区二区| 亚洲AV电影不卡在线观看| jizz国产视频| 国产欧美视频在线| 国产亚洲精品yxsp| 国产AV毛片| 高清码无在线看| 91丝袜乱伦| 在线观看国产精美视频| 国产精品毛片一区| 手机在线免费不卡一区二| 欧美激情伊人| 尤物成AV人片在线观看| 国产成人喷潮在线观看| 波多野吉衣一区二区三区av| 在线精品自拍| 亚洲熟女中文字幕男人总站| 一级高清毛片免费a级高清毛片| 一级毛片基地| 亚洲国模精品一区| 精品成人一区二区三区电影| 亚洲中文无码av永久伊人| 国产精欧美一区二区三区| 国产精品久久久免费视频| 亚洲国产精品一区二区第一页免| 国产自无码视频在线观看| 99国产精品免费观看视频| 91麻豆精品视频| 免费女人18毛片a级毛片视频| 熟女视频91| 国产欧美日韩精品综合在线| 9久久伊人精品综合| 国产欧美在线视频免费| 国产精品制服| 无码精油按摩潮喷在线播放 | 欧美成人第一页| 高清精品美女在线播放| 国产在线观看成人91| 2021天堂在线亚洲精品专区| 欧美成人免费午夜全| 欧美日韩第三页| 国产成人三级| 人妻中文久热无码丝袜| 国产另类乱子伦精品免费女| 综合五月天网| 欧美成人在线免费| 日韩久久精品无码aV| 色综合成人| 欧美激情视频一区| 一区二区三区成人| 日韩欧美国产中文| 国产产在线精品亚洲aavv| 91视频首页| 久久大香伊蕉在人线观看热2| 国产精品成人第一区| 亚洲熟女偷拍| 国产免费网址| 99精品一区二区免费视频| 亚洲欧美一区二区三区图片| 国产剧情无码视频在线观看| 午夜影院a级片| 国产成人精品亚洲日本对白优播| 热思思久久免费视频| 国产不卡在线看| 中文纯内无码H| 天天躁夜夜躁狠狠躁躁88| 青草视频久久| 欧美性爱精品一区二区三区|