基于時序特征的玉麥輪作制耕地提取方法研究

李百紅 彭勃 董超

摘要 [目的]該研究對小麥、玉米輪作制下耕地的面積與分布有重要意義。[方法]基于目視判別收集樣本點和GEE平臺，分析地物的NDVI、EVI等指數時序特征，開發特征識別算法提取2018年玉麥輪作區面積分布，同時與CART算法監督分類結果進行精度分析。[結果]與基于多光譜的監督分類相比，基于時序特征的識別算法準確率、精準率、F1 Score和Kappa系數分別提高了0.121、0.110、0.136和0.246。[結論]基于時序特征的特征識別算法可有效識別玉麥輪作區，該算法可為輪作區農業生產提供基礎數據支持。

關鍵詞 GEE;時序;歸一化指數;提取算法

中圖分類號 S 126? 文獻標識碼 A? 文章編號 0517-6611（2021）19-0214-04

doi：10.3969/j.issn.0517-6611.2021.19.056

開放科學（資源服務）標識碼（OSID）：

Study on Extraction Method of Wheat-Maize Rotation System Based on Time Series Characteristics

LI Bai-hong1， PENG Bo2， DONG Chao2

（1.Taian Natural Resources and Planning Bureau of Shandong Province， Taian， Shandong 270000;2.College of Information Science and Engineering， Shandong Agricultural University， Taian， Shandong 271018）

Abstract [Objective] It is important to study the area and distribution of cultivated land under wheat and maize rotation system. [Method] Based on the collection of sample points and GEE platform， we analyzed the NDVI， EVI and other index time series features of ground objects， developed feature recognition algorithm to extract the 2018 wheat-maize rotation area distribution， and conducted precision analysis with the supervised classification results of CART algorithm. [Result] Compared with the multi-spectral supervised classification， the accuracy， F1 Score and Kappa coefficient of the time-series feature recognition algorithm were improved by 0.121， 0.110， 0.136 and 0.246， respectively. [Conclusion] The feature recognition algorithm based on time series features could effectively identify wheat-maize rotation area， and the algorithm could provide basic data support for agricultural production in rotation area.

Key words GEE;Time series;Normalized index;Extraction algorithm

基金項目 山東省重點研發項目（14032761，140380198）。

作者簡介 李百紅（1985—），女，山東青州人，工程師，碩士，從事土地資源利用與信息化工作。*通信作者，講師，博士，從事農業信息技術研究。

收稿日期 2021-02-26

快速準確獲取作物種植的分布信息有利于實現對耕地的精確管理、產量估算、制定科學的糧食管理政策[1]。當前，作物識別多基于遙感影像光譜差異，如熊勤學等[2-3]利用NDVI提取小麥的種植面積，路中等[4-6]提取玉米的種植面積，劉佳等[7-8]完成了中國冬小麥制圖和大規模春小麥提取研究。受限于傳統數據處理工具的低效，這些研究只使用了一種指數提取地物信息，難以發揮高分辨率多光譜遙感數據的優勢。多時相的遙感數據進行分類不僅可以提取地物物候特征，而且可有效避免單一時相存在的異物同譜和同譜異物狀況[9]。Google Earth Engine（GEE）平臺可對大尺度區域進行分析，還支持多源多時相海量數據的協同研究[10]。鑒于此，筆者選取華北平原小麥-玉米輪作區，基于GEE平臺，利用多源遙感數據分析地物樣本點光譜指數的時序特征，依據各類地物的時序特征設計有效的識別算法，通過精度評價評估識別算法的優劣。

1 研究區域與數據

1.1 研究區概況

華北平原小麥-玉米輪作區（32°16′～39°38′N，112°50′～118°57′E）總面積約221 838 km2，跨河北、河南、山東3省（圖1）。研究區大部分平原區為栽培植被，其中以冬小麥、夏玉米輪作最為常見，主要分布在河南北部和山東西北的引黃灌溉區、太行山前平原水源條件較好的區域。研究區內小麥10月中上旬播種，次年3月拔節，6月上旬收割;6月初播種玉米，9月下旬收割[11]。

1.2 數據及預處理

1.2.1 衛星影像選取。

研究所用數據包括：Landsat7 SR（Surface Reflectance，SR）數據集;Landsat8 SR數據集;MOD13Q1.006數據集的植被指數（MOD13Q1）第6版數據;Sentinel-1 SAR GRD數據集。研究影像數據來自GEE平臺（https：//earthengine.google.com/）。

1.2.2 影像預處理。

影像預處理包括：①使用研究區對影像統一進行掩膜;②利用Landsat7和Landsat8影像中由CFMASK算法生成的質量評估波段去云[12];③對Landsat7和Landsat8計算歸一化光譜指數。NDVI 計算公式如下：

NDVI=（NIR-R）/（NIR+R）（1）

根據NIR、R與B波段計算EVI計算公式如下：

EVI=2.5×（NIR-R）/（NIR+6×R-7.5×B+1）（2）

根據G與MIR波段計算MNDWI，計算公式如式下：

MNDWI=（G-MIR）/（G+MIR）（3）

將Landsat7和Landsat8的影像合成指數的數據集，以彌補Landsat系列衛星回訪周期較長的缺點。

1.3 樣本點與驗證點

基于時間序列將研究區地物類別分為小麥-玉米輪作區、林地、城鎮建設用地、水域和其他耕地類型。利用在線高分辨影像及歷史時間軸，結合研究區內各區縣的相關文獻資料和統計報告，通過目視的方法選取地物樣本點和驗證點（表1）。

2 時序數據分類方法

2.1 時序特征的組合

基于指數數據集生成樣本點指數時序曲線，取同類指數的平均值生成各類別指數年積日（Day of year，DOY）時序曲線。指數時序曲線可反映類別特點，研究基于NDVI、EVI和MNDWI共3種指數的時序特征，篩選區分地物類別的時序參數。

2.2 輪作區耕地的識別

采用分層分類加決策樹的方法，基于指數的時序特征對研究區中的玉麥輪作區外進行分層掩膜[13]，掩膜順序為：水體、建筑、森林、輪作區。由于Landsat重訪周期較長，受天氣影響較大，在識別算法中使用MODIS的NDVI產品來對輪作耕地進行最后的掩膜。

2.3 監督分類

利用目視樣本點作為訓練數據集，通過GEE提供的云評分算法制作研究區最小云量合成影像，并基于影像的藍光、綠光、紅光、近紅外、短波紅外等7個波段的合成影像，使用分類回歸樹（Classification And Regression Tree，CART）算法完成輪作區耕地提取[14]。

2.4 精度評價

利用獨立的驗證點對結果進行分類評價，將非輪作區耕地的地類視為一類與輪作區耕地進行二分類評價，通過準確率、精準率、F1 Score以及Kappa系數對提取精度進行評價。

3 結果與分析

3.1 時序特征分析

由圖2a可知，水域MNDWI數值全年都為正值，且均值大于0.4明顯區別于其他地類，基于MNDWI提取水域效果最優。圖2b中森林和耕地在EVI上有較為明顯的差異，在DOY120～273 d，森林的EVI均值大于2，明顯區別與小麥收割后裸地。圖2c中建設用地的NDVI曲線較為平緩，全年的數值偏低，且年均值小于0.3，以此來區分建筑用地。圖2d中玉麥輪作區在NDVI時序上具有明顯雙峰特征，明顯區別于其他地類。儒略日第79～154天是第1個峰，第198～271天是第2個峰。

3.2 監督分類結果

為比較算法的精度，采用監督分類與識別算法進行比較。二分類結果如圖3所示，分類結果中輪作耕地分布過多，明顯與調查情況不一致。

二分類混淆矩陣如表2所示，可見分類結果中，有大量輪作耕地被誤分為非輪作耕地，同時也有大量非輪作耕地被誤分為輪作耕地，2種地類相互混淆，提取精度較低。

監督分類結果中主要城市、河流湖泊分類準確，但北部只有少量其他耕地被正確識別，且分類結果中林地分布偏少，面積為33 651.6 km2，輪作耕地分布過多，面積為121 661.7 km2，與調查情況不符。混淆矩陣如表3所示，林地、其他耕地與輪作區耕地之間錯分較多，有較大比例的林地和其他耕地被誤分為小麥-玉米輪作區。

3.3 時序分類結果

基于時序數據使用識別算法提取結果精度較高，輪作區耕地面積為68 032 km2，其中東北部（即黃河流域）、中部靠太行山一側和偏南部分布較多，分布情況基本與調查結果相符。二分類混淆矩陣如表4所示，只有少量輪作耕地和非輪作耕地發生錯分，該算法能較為精確地對輪作區進行識別。

3.4 精度對比

基于時序特征的識別算法準確率、精準率、F1 Score和Kappa系數分別為0.987、0.970、0.975和0.956，分別比監督分類高出0.121、0.110、0.136和0.246。準確率、精準率、召回率和F1 Score越高，分類精度越高，這4個指標都是時序特征識別算法占優。Kappa系數表示混淆矩陣的不平衡度，監督分類結果的Kappa系數遠低于時序識別算法。

4 結論

該研究基于GEE平臺，分別利用識別算法與監督分類，對研究區小麥、玉米輪作區耕地進行提取，并進行精度對比分析，得到以下結論：

（1）該研究實現的提取算法可以快速有效地實現大區域小麥-玉米輪作耕地的提取，其二分類準確率、精準率、F1 Score和Kappa系數分別為0.987、0.970、0.975和0.956，分別比基于多光譜的分類高出0.121、0.110、0.136和0.246。

（2）某些特殊時期，林地、輪作耕地及其他耕地具有相似的光譜特征，傳統基于多光譜數據的監督分類難以區分，導致其難以準確地提取輪作耕地。

（3）研究區2018年約有小麥-玉米輪作耕地68 032 km2，約占研究區總面積的30.7%，主要分布在河北東南部、河南東北部和山東西北部區域。

利用遙感數據的各類歸一化指數的時序特征，能夠比利用多光譜數據的監督分類更準確地區分具有類似光譜特征的地物，如輪作耕地與其他耕地，而在GEE平臺的支持下，能快速獲取并處理大區域、長時間的遙感數據，從而快速通過時序特征來區分不同地物。

參考文獻

[1] 史舟，梁宗正，楊媛媛，等.農業遙感研究現狀與展望[J].農業機械學報，2015，46（2）：247-260.

[2] 熊勤學，黃敬峰.利用NDVI指數時序特征監測秋收作物種植面積[J].農業工程學報，2009，25（1）：144-148.

[3] 王菲，李紅華，蘇潔.基于MODIS-NDVI的邯鄲市冬小麥和夏玉米面積提取[J].現代農業科技，2018（24）：3-5.

[4] 路中，雷國平，馬泉來，等.基于重構的Landsat 8時間序列數據和溫度植被指數的區域旱情監測[J].水土保持研究，2018，25（5）：371-377，384.

[5] 蘇偉，姜方方，朱德海，等.基于決策樹和混合像元分解的玉米種植面積提取方法[J].農業機械學報，2015，46（9）：289-295，301.

[6] 翟世常，王學強，頡耀文.黑河流域中游盆地玉米作物遙感估產研究[J].安徽農業科學，2015，43（31）：350-353.

[7] 劉佳，王利民，楊福剛，等.基于GF-1衛星影像的中國冬小麥制圖研究[J].中國農學通報，2019，35（9）：155-164.

[8] 劉沼輝，柳林，郭慧，等.基于GF1-NDVI時序影像對春小麥進行提取研究[J].北京測繪，2018，32（6）：643-646.

[9] 趙麗花，李衛國，杜培軍.基于多時相HJ衛星的冬小麥面積提取[J].遙感信息，2011，26（2）：41-45，50.

[10] 譚深，吳炳方，張鑫.基于Google Earth Engine與多源遙感數據的海南水稻分類研究[J].地球信息科學學報，2019，21（6）：937-947.

[11] 侯學會，隋學艷，姚慧敏，等.中國北方麥區冬小麥物候期對氣候變化的響應[J].麥類作物學報，2019，39（2）：202-209.

[12] FOGA S，SCARAMUZZA P L，GUO S，et al.Cloud detection algorithm comparison and validation for operational Landsat data products[J].Remote sensing of environment， 2017，194：379-390.

[13] 向海燕，羅紅霞，劉光鵬，等.基于 Sentinel-1A 極化 SAR 數據與面向對象方法的山區地表覆被分類[J].自然資源學報，2017，32（12）：2136-2148.

[14] DIDAN K.MOD13Q1 MODIS/Terra Vegetation Indices 16-Day L3 Global 250m SIN Grid V006 [EB/OL].[2020-11-05].https：//doi.org/10.5067/MODIS/MOD13Q1.006.