基于R語言的作物分類研究

刁雨晴，冉谷林，查元源

（武漢大學水資源與水電工程科學國家重點實驗室，武漢 430072）

0 引 言

灌區的農業生產規劃情況影響著全國的糧食生產，對糧食生產具有關鍵的意義。作物的空間分布影響著農業生產管理方法以及農業政策的制定，快速、精準地提取灌區的作物種植信息對推進灌區生產的合理規劃及農業現代化具有重要意義。

遙感是一門高效、便捷的監測技術，其檢測范圍廣，信息多元。精準農業是農業發展的趨勢，信息采集是發展精準農業的基礎，將農業信息監測與遙感技術相結合，可以快速、準確的提取到所需的農業信息，目前，遙感技術已廣泛地應用于農業信息的監測[1,2]。

作物分類方法根據數據源的差異又可分為基于單一遙感數據的分類和基于多源數據融合的分類。單一的遙感數據源往往受時間分辨率或空間分辨率的限制，導致其不能準確表達作物的光譜特征，從而導致分類結果較差。基于多源數據融合的分類就是把各種遙感及非遙感數據融合為一個整體的數據集，這種方法可以有效地減少由于單一遙感數據分辨率引起的“同物異譜”及“同譜異物”的現象，從而提高分類的精度。根據宋茜[3]等統計，目前農作物遙感識別的研究中基于單一數據源的分類仍占主導，但基于多源遙感數據的分類越來越多地得到應用，單一數據源在單一農作物的識別和提取方面應用較為普遍，而多源數據在農作物分類方面的應用更多。

目前，國內外已有大量基于遙感的作物分類研究。使用遙感數據進行作物分類的方法主要分為統計分析的圖像分類、人工神經網絡分類?；诮y計分析的遙感圖像分類方法又分為監督分類和非監督分類。張健康[4]等使用TM/ETM+影像數據和MODIS EVI 影像數據，采用監督分類與決策樹分類相結合的人機交互解譯方法，建立決策樹識別模型，對黑龍港地區的主要作物進行了遙感解譯；王冬利[5]等使用GF-1 多光譜數據，以歸一化植被指數作為冬小麥信息提取的判別指標，運用非監督分類結合多尺度技術提取了辛集市冬小麥種植信息。張東彥[6]等基于Sentinel-2 數據，分析內蒙古自治區興安盟扎賚特旗現代農業示范園區作物的典型植被指數時序變化特征，采用隨機森林等機器學習算法對研究區主要作物進行識別。

在應用遙感數據進行作物分類的過程中，分析處理遙感數據的工具至關重要。目前的研究多使用arcGIS、ENVI 等軟件，這些軟件具有成熟的數據分析功能，但對于遙感數據的下載以及批量處理等方面有所欠缺。R語言是近些年來較為流行的一款開源性編程語言，其操作界面簡單，可視化功能完善，語法靈活[7-9]。此外，R 語言中的包具有眾多功能，用戶可根據需要在R 語言中直接下載安裝包，這些包涵蓋了遙感數據的下載、預處理及計算分析等各個方面。由于遙感數據的下載、格式等差異，使用集成化軟件處理遙感數據需要跨越多個平臺操作，步驟繁瑣，效率低下，而R 語言中只需安裝對應功能的包并輸入相應函數，即可批量完成遙感數據的處理過程，這種流程化的處理方式避免了跨平臺操作，極大地提高了遙感數據的處理效率。

因此，本研究基于R 語言，采用多源遙感數據融合技術，提出了基于機器學習的作物分類方法，為實現多源遙感數據融合，提高遙感數據處理效率，提升作物分類精度提供了新思路。

1 研究區及數據

1.1 研究區概況

研究區漳河灌區位于湖北省中部，江漢平原西部，地理坐標為111°54′～112°42′E，30°23′～31°34′N，總灌區面積約5 543 km2，地勢西北部高，東南部低，地形主要為平原與丘陵。灌區屬長江中游亞熱帶季風氣候區，日照充足，降水充沛，年平均氣溫15.6～16.3 ℃，年平均降水量804～1 067 mm。農作物以水稻、棉花、冬小麥、油菜為主，部分地區有早稻-晚稻連作、水稻-油菜輪作及小麥-棉花輪作的種植模式。漳河灌區是湖北省最重要的商品糧食基地之一，具有南方灌區的特點，以該區域作為研究區進行作物分類具有一定的代表性，見圖1。

圖1 漳河灌區區位圖Fig.1 Zhanghe irrigation area location map

1.2 數據及預處理

本研究采用2017年12月至2018年12月的GF-1、Landsat 8 OLI 和Sentinel-2 數據。GF-1 衛星是我國于2013年4月26日發射的一顆高分辨率對地觀測衛星，其波段1-4分別對應可見光波段及近紅外波段，空間分辨率16 m；Landsat 8 衛星于2013年2月11日由美國航天航空局成功發射，衛星共有11 個波段，其中波段2-5為可見光波段及近紅外波段，其空間分辨率均為30 m[10]；Sentinel-2 衛星于2016年6月23日成功發射，其共包含12個波段，波段2-4對應可見光波段，波段8對應近紅外波段，可見光及近紅外波段的空間分辨率均為10 m。在研究區內根據云量小于10%篩選，GF-1全年可用數據共3期，Landsat 8 覆蓋研究區需2 景影像，全年可用數據共計5 期，Sentinel-2 數據覆蓋研究區需4-6 景影像，全年可用數據共4期。影像數據見表1。

表1 影像列表Tab.1 Image list

對于研究所選取的遙感數據的下載及預處理，R語言中有大量的處理遙感數據的包，這些包涵蓋了遙感數據的下載、預處理及使用、輸出等各個方面。getSpatialData包是一個用于下載遙感數據的包，研究所選的Landsat 8 數據及Sentinel-2 數據均可以通過該包進行下載，除此之外，getlandsat 包[11]和sen2r 包[12]也可以分別用于下載Landsat 數據和Sentinel 數據。除了下載，sen2r包還集成了針對Sentinel數據的校正、裁剪等各種構建完整的Sentinel-2 處理鏈所需的所有步驟，無須任何外部工具。GF-1 數據需要在中國資源衛星應用中心進行下載，但其處理可以在R 語言中完成。對于除Sentinel-2 數據外的遙感數據處理，RSToolbox 包[13]是一個專門處理遙感數據的包，可以對遙感數據進行導入、預處理、數據分析等多種處理；raster 包[14]是針對柵格對象處理的包，可以用于遙感圖像的拼接、裁剪、分析計算及保存等多個功能。使用上述R 包，對研究所選取的數據進行校正、拼接，將影像數據根據研究區范圍進行裁剪，導出藍、綠、紅、近紅外波段，并分別按日期保存得到研究所需的遙感數據。本研究使用的R 包及其功能匯總見表2。

表2 研究中使用的R包及其主要功能介紹Tab.2 The R package used in the study and its main functions

通過查閱相關文獻，訪問農業信息網，并結合Google Earth 高清歷史影像，可以基本掌握各種作物的種植區域、紋理和色相等特點。油菜廣泛分布于漳河灌區各個區域，其花期在10月下旬，影像上呈亮黃綠色；棉花主要產于馬良鎮附近區域，其花期在7月下旬，花期時棉花在影像中呈黃白色；除了油菜與棉花，各種作物都有比較獨特的特征，可以在影像中較好的分辨各種作物。在Google Earth 上根據目視解譯選取作物樣本點共3 092 個，按70%作為訓練樣本，30%作為測試樣本，見表3。

表3 地物樣本數量Tab.3 Number of sample points

2 研究方法

以遙感數據為基礎，分別計算各日期的歸一化差值植被指數(Normalized Difference Vegetation Index,NDVI)、增強型植被指數(Enhanced Vegetation Index,EVI)及比值植被指數(Ratio Vegetation Index,RVI)，并分別構造3種植被指數的時間序列數據。根據作物樣本點分別提取各作物的時間序列曲線，以此作為分類依據，分別使用CART 決策樹、隨機森林(Random Forest,RF)、樸素貝葉斯(Naive Bayes,NB)、支持向量機(Support Vector Machines,SVM)及K 近鄰(K Nearest Neighbor,KNN)算法進行作物分類，并對比分類效果。

2.1 特征選擇

在遙感技術發展初期，研究人員為了探索植被與遙感光譜波段之間的關系，將各波段以不同方式組合，分別探究各種組合與植被之間的相關關系，在此過程中發現了植被在紅光波段具有強吸收、在近紅外波段具有高反射的特點，提出了植被指數這一概念。隨著遙感技術的進一步發展，各種形式的植被指數被提出，植被指數逐步被應用在環境、生態、農業等領域[15]。在作物分類方面，因為各作物的生育期不同，光譜特征隨作物生育期發生變化，植被指數也會隨之變化，不同作物的植被指數變化趨勢不同，因此，可以以植被指數的時間序列為分類依據來對作物進行分類。

NDVI是在作物分類中應用最廣泛的植被指數，它是近紅外波段與紅光波段反射率之差與二者之和的比值。根據NDVI值，可以將植被從水體和土地中分類出來，NDVI值越大代表植被覆蓋度越高。但以NDVI時間序列作為作物分類依據有兩個缺陷:一是NDVI缺乏對大氣及土壤背景干擾的處理；二是NDVI在覆蓋度較好的地區飽和問題比較嚴重[16]。EVI和RVI可以改進NDVI的這兩個缺陷。針對大氣及土壤背景的干擾問題，EVI引入了藍光波段，減少了氣溶膠及土壤背景對植被指數的影響；在高植被覆蓋區，RVI具有比NDVI更高的靈敏度，且RVI可以趨于無窮大，不存在飽和問題。

基于以上分析，研究選用NDVI、EVI、RVI3 種植被指數為分類指標，各指標介紹如表4所示。

表4 植被指數介紹Tab.4 Calculation formula of vegetation index

2.2 分類方法

作物分類就是將各種作物的特征信息作為樣本集，通過分類器從樣本集中尋找到各類作物的共性，從而形成一種統一的分類規則。將測試樣本代入該分類規則中即可獲得分類的精度，將新數據代入該規則中即可得到預測分類結果。研究使用決策樹、隨機森林、樸素貝葉斯、支持向量機及K 近鄰算法5種分類方法進行分類。

決策樹是一種樹形的分類模型，它在樹枝節點處依據特征值對數據類別進行判斷，把數據集劃分為兩個分支，在分支的節點處再進行類別判斷，如此重復下去，將訓練的數據集按級依次分割，最終將具有類似特征條件的數據劃分在同一個集合內，完成對數據集的分類。本文選用決策樹中最基礎的CART決策樹進行分類。

隨機森林算法是一種基于決策樹算法而改進的集成學習方法，它組合了多棵決策樹，隨機選擇樣本并對樣本數據集進行重復預測，每棵決策樹都會有不同的分類結果，使用多數投票法對不同分類結果進行投票整合，最終得到對數據集分類的結果。

樸素貝葉斯分類方法是以貝葉斯定理為基礎理論的分類方法，是貝葉斯分類器中應用最為廣泛的模型之一。它計算給定樣本在各類別上的后驗概率，通過概率推理的方法，將樣本判定為最大后驗概率所對應的類別，從而完成分類[17,18]。

支持向量機算法是一種基于統計學習理論的學習方法[19,20]，該算法采用結構風險最小化的原則，求解一個能夠劃分訓練數據類別且能使各類別之間的幾何間隔最大的分類面，通過該分類面實現數據的分類。

K近鄰算法是一種非參數的分類技術，它通過在分類樣本集中選取與待分類樣本點最相近的K 個已知類別的樣本，從而將該待分類樣本點劃分為該類別，達到分類的效果。該算法簡單且易于實現，但計算量龐大，分類速度比較慢[21]。

2.3 精度評價

研究使用總體分類精度和Kappa系數評價研究區地物分類效果?？傮w分類精度為正確分類的像元數占分類總像元數的比例，即混淆矩陣的對角線之和與混淆矩陣之和的比值；Kappa系數是通過把所有地表真實分類中的像元總數乘以混淆矩陣對角線的和，再減去某一類地表真實像元總數與該類中被分類像元總數之積對所有類別求和的結果，再除以總像元數的平方減去某一類地表真實像元總數與該類中被分類像元總數之積對所有類別求和的結果所得到的[22]。

使用用戶精度(User’s Accuracy,UA)、生產者精度(Producer’s Accuracy,PA)評價各類地物分類效果。用戶精度表示某一類地物正確分類的像元數與該類別地物實際總像元數的比值；生產者精度表示某一類地物正確分類的像元數與整個分類過程中分類到該類別的像元數的比值[23]。

3 結果與分析

3.1 光譜時序特征分析

分別計算漳河灌區的NDVI、EVI、RVI，提取各作物樣本點的植被指數特征值，剔除異常值點，并按類別求取均值，得到漳河灌區主要地物的NDVI、EVI、RVI時間序列曲線見圖2～圖4。

圖2 各作物NDVI時序變化曲線Fig.2 Temporal changes of NDVI of all crops

圖3 各作物RVI時序變化曲線Fig.3 Temporal changes of RVI of all crops

圖4 各作物EVI時序變化Fig.4 Temporal changes of EVI of all crops

由圖2～圖4可以看出，建筑用地的NDVI、EVI、RVI值在全年都處于很低的水平，波動不大。森林的年內變化不大，其NDVI值全年都處于較高水平，均大于0.5，而RVI在高覆蓋區具有比NDVI更高的靈敏度，因此在夏季茂密的森林地區，RVI值呈上升趨勢，并在8月達到峰值。

棉花-小麥輪作區、早稻-晚稻間作區、水稻-油菜輪作區的植被指數時序曲線均具有明顯的雙峰特征，符合輪作區的植被指數變化規律。其中，棉花-小麥輪作區的NDVI曲線的第一個峰值在3月，對應的是小麥的返青及拔節，第一個谷值在6月，對應的是小麥的收獲，第二個峰值在6月，對應的是棉花開花；早稻-晚稻間作區的NDVI曲線自3月下旬早稻播種后緩慢降低，5月中旬早稻分蘗期NDVI值達到第一個谷值，然后NDVI值開始增加，在6月下旬早稻抽穗期達到第一個峰值，7月下旬早稻收獲，在8月中旬晚稻分蘗期時NDVI曲線達到第二個谷值，在9月中旬晚稻抽穗期達到第二個峰值。水稻-油菜輪作區的NDVI曲線第一個峰值是在3月油菜開花期，第二個峰值在7月中稻抽穗期。各輪作區的EVI、RVI時序曲線的趨勢均與NDVI時序曲線的趨勢基本相同，但RVI的變化幅度更劇烈。

3.2 分類方案及分類過程

將地物分類過程根據分類數據分成4 個方案，方案1、2、3分別以NDVI、RVI、EVI時間序列為分類數據，方案4采用按順序合成的3種植被指數時間序列為分類數據。每種方案均使用CART 決策樹、樸素貝葉斯、支持向量機、K 近鄰、隨機森林五種不同的算法為分類方法，對漳河灌區地物進行分類預測。

R 語言中有可以完成作物分類及驗證的包。例如，caret包[24]可以訓練決策樹、樸素貝葉斯分類模型，randomForest 包可以訓練隨機森林算法分類模型，e1071 包可以訓練支持向量機算法分類模型，class包[25]可以訓練K近鄰算法分類模型。對于訓練好的模型，將測試集數據去除類別信息，使用predict函數進行預測，并調用confusionMatrix 函數計算分類的混淆矩陣、總體精度值、Kappa值、用戶精度值及生產者精度值等指標。將整理好的植被指數時序圖像作為預測數據，使用predict函數，即可得到研究區域的作物分類圖見圖5。

圖5 漳河灌區地物分類預測圖Fig.5 Classification forecast map of Zhanghe Irrigation area

3.3 分類結果與精度評價

圖5為使用隨機森林算法對NDVI、RVI、EVI3 種植被指數混合數據進行分類得到的漳河灌區地物分類預測圖。從圖5可以看出，漳河灌區主要種植作物有水稻、油菜、棉花、小麥。其中水稻-油菜輪作的種植方式在漳河灌區最為普遍，分布較廣；而棉花-小麥輪作的種植方式主要集中在灌區東部臨近漢江的地帶；早稻-晚稻間作的種植模式在灌區分布較少，主要集中在灌區東南部靠近長湖的地帶。

使用訓練好的模型對測試數據進行分類預測，并調用confusionMatrix 函數，整理運行后的結果，得到分類總體精度值、Kappa 系數值見表5，各地物的用戶精度及生產者精度見圖6～圖7。

表5 分類精度統計Tab.5 Classification accuracy statistics

根據表5，對方案4 采用隨機森林算法對灌區進行地物分類的分類效果最好，總精度達96.96%，Kappa 系數為0.948；對方案3 采用CART 決策樹對灌區進行分類的分類效果最差，總精度僅為81.52%，Kappa 系數僅為0.662。對比同一方案中不同分類器的分類效果，隨機森林算法的分類效果最好，K近鄰算法、支持向量機算法及樸素貝葉斯算法次之，CART 決策樹分類器的效果最差，其中，隨機森林算法的分類平均總精度達95.6%，K 近鄰、支持向量機及樸素貝葉斯算法的分類平均總精度值在93.2%～94.5%之間，CART 決策樹分類的平均總精度僅有83.2%。對比同一分類器下不同方案的分類效果，方案4 效果最好，方案1 及方案2 的平均總精度值及Kappa 系數略低于方案4，但差別不大，方案3 的分類效果最差，平均總精度值比其他方案低2%左右。

根據圖6～圖7對比各方案中地物的用戶精度值及生產者精度值，水體在各種分類方案下的用戶精度值均為100.0%，生產者精度均大于97.0%，其分類結果最可靠，分類效果最好；油菜-水稻輪作區、棉花-小麥輪作區、其他植被區及建筑用地的分類可靠性及效果比水體略差；早稻-晚稻間作區的分類結果最不可靠，分類效果最差。結合圖2～圖4，水體、建筑及其他植被區的植被指數曲線具有明顯特征，易于從其他地物中區分出來；油稻輪作區及棉麥輪作區的植被指數時序曲線雖有相同的雙峰趨勢，但達到峰值的時間不同，因此也利于區分；而早晚稻間作區由于面積較小，樣本點數目少，且植被指數曲線在夏季達到峰值的時間與油稻輪作區基本相同，容易錯分進油稻輪作區，因此分類效果不好。

圖6 各地物分類用戶精度Fig.6 User accuracy of local object classification

圖7 各地物分類生產者精度Fig.7 Producer accuracy of local object classification

綜上所述，對于分類器，使用隨機森林分類器的分類效果最好，使用CART 決策樹的分類效果最差；對于分類數據，方案4 的數據分類效果最好，方案3 的數據分類效果最差，但綜合考慮工作量與精度提升之間的關系，使用方案1的數據進行分類與方案4差別不大；對于地物分類結果，水體的分類結果最可靠，雙季稻區域的分類結果最差。

4 結 論

本研究以R 語言為工具，分別基于GF-1數據、Sentinel-2數據及Landsat 8 數據構建的湖北省漳河灌區NDVI、EVI、RVI時間序列，采用決策樹、支持向量機、K近鄰算法、樸素貝葉斯算法及隨機森林5種分類器對研究區地物進行分類，主要得到以下結論：

（1）以R 語言為工具，實現了遙感數據的下載、預處理、分析及可視化過程，避免了跨平臺數據處理，簡化了遙感數據的處理過程，驗證了R 語言在遙感數據處理過程中的適用性與優越性。

（2）采用CART 決策樹、樸素貝葉斯、支持向量機、K 近鄰算法及隨機森林模型對研究區地物分類。結果表明，隨機森林算法分類效果最好，平均總精度達到95.60%；CART決策樹算法分類效果最差，平均總精度僅83.15%。

（3）分別使用NDVI、RVI、EVI及三者組合的時間序列作為分類數據，對比4種方案的分類效果，結果表明，使用組合數據分類的效果最好，單獨使用EVI時間序列的分類效果最差，單獨使用NDVI作為分類數據進行分類與使用組合數據分類的精度差別不大。綜合考慮工作量與精度提升之間的關系，使用NDVI時間序列為分類數據進行分類最佳。

（4）使用用戶精度、生產者精度對各地物分類效果進行評價，結果表明，水體分類效果最好，雙季稻分類效果最差。

本研究的主要優勢在于使用R 語言簡化了遙感數據的處理及分析過程，為提高遙感數據處理效率提供了新思路。同時，對湖北省漳河灌區的作物進行分類，為當地農業布局提供了數據支撐。但研究還存在一些不足之處需要改進，目前研究表明，同時考慮地物的光譜特征和紋理特征可以極大提高分類的精度，后續研究應在本研究的基礎上考慮使用面向對象的分類方法對本研究進行改進。