基于多時相Sentinel-2A的縣域農作物分類

吳 靜 呂玉娜 李純斌 李全紅

(1.甘肅農業大學資源與環境學院， 蘭州 730070； 2.甘肅農業大學管理學院， 蘭州 730070)

0 引言

農作物種植種類、面積和產量估算是國家糧食生產安全和經濟運行的重要監測指標。遙感技術由于能快速獲取作物類型和分布，已成為大尺度范圍農作物監測的主要手段之一[1]。農作物的遙感精細分類是精確估產的前提，是目前農業遙感的研究熱點之一[2-5]。

充分利用農作物的典型反射光譜特征和季相節律特征是區分不同作物的關鍵理論依據[6]，也是利用光學遙感手段進行農作物分類的主要著力點。國內外相關研究多采用農作物生長發育過程中多時相的光學遙感數據提取時間序列的植被指數曲線來分類。采用的遙感數據包括MODIS[7-11]、Landsat[2-3,12-14]、HJ衛星[15-17]、GF-1[18-20]等以及多元數據融合[21-24]，計算的植被指數包括NDVI、EVI等；提取的農作物信息包括單一種類農作物識別及分布信息[3-4,7,14,17]，以及多種農作物種植結構[18-27]。

圖1 景泰縣地理位置和行政區劃Fig.1 Location and administrative zoning map of Jingtai County

近年來，國內對多種農作物分類的研究主要集中于影響因素可控的小區域(如農場)[1,15-16,21]或地勢平坦的東北平原[7-10,19,22,24]，在地形較復雜的縣域尺度上進行農作物分類有一定的挑戰性，因為可能遇到地塊較為破碎，分布比較零散、作物種植結構相對復雜等問題。

Sentinel-2A衛星能夠提供10 d重訪周期和最高10 m空間分辨率的多光譜數據，尤其是其紅邊波段數據為農作物類型識別和產量監測提供了強有力的數據支持[26]。目前對于利用Sentinel-2A數據的時空優勢及其時序特征實現農作物提取的相關研究正在陸續展開[1,26-28]。然而農作物種植結構具有強烈的區域性，不同的區域地形地貌、氣候、水文特征不同，農作物的種類、分布、物候特征也不同，Sentinel-2A數據的全面應用需要以實際案例研究為基礎，探討其在具體區域農作物分類方面的能力和應用潛力。本文以西北內陸溫帶干旱區的景泰縣為研究區，在實地調查研究區主要作物及分布的基礎上，結合農作物物候信息、土地變更調查數據，利用2018年多時相Sentinel-2A數據，計算時序NDVI和RENDVI指數及其組合為特征參數，采用隨機森林法對研究區的主要農作物進行分類，并結合實地調查選取驗證數據進行結果驗證，對比分析不同特征參數提取作物信息的精度，以探討多時相Sentinel數據的不同特征參數應用于干旱區縣域農作物分類研究的可行性，以期為當地農作物種植結構信息的提取提供新思路。

1 研究區與數據源

1.1 研究區概況

景泰縣位于甘肅省中部(103°33′～104°43′E，36°43′～37°38′N)，河西走廊東端，甘、寧、蒙三省(區)交界處，是黃土高原與騰格里沙漠過渡地帶，如圖1所示。全縣面積5 483 km2，境內海拔1 276～3 321 m，地勢由西南向東北傾斜，地貌類型以傾斜平原為主，地勢相對平坦。景泰縣屬于溫帶干旱型大陸氣候，年均降水量185 mm，年均蒸發量3 038 mm；年日照時數2 652 h，日照百分率60%，年平均太陽輻射量619 kJ/cm2，年均溫度8.2℃，大于等于0℃年活動積溫3 614.8℃，大于等于10℃年有效積溫3 038℃，無霜期191 d，適宜一年一季的農業生產[29]。黃河流經縣境110 km，境內有被譽為“中華之最”的景電一、二期高揚程提灌工程，農業生產條件較為便利，是《全國新增1 000億斤糧食生產能力規劃(2009―2020年)》中確定的甘肅省7個產糧大縣之一。

本研究為減少其他地類的干擾、提高作物分類精度，首先提取出耕地信息，然后在耕地范圍內提取作物信息。農作物主要有春小麥、玉米、水稻、秋油葵、胡麻、洋蔥、大棚作物、馬鈴薯、砂田瓜果9種類型，其中，春小麥-秋油葵是景泰縣傳統輪作模式，除此之外，其他作物種植都屬于一年一季。各種作物物候期見表1。

表1 景泰縣主要農作物物候期Tab.1 Crops development period of Jingtai County

注：表中數據來源于實地走訪調查及參考中國氣象數據網(http://data.cma.cn/)

1.2 遙感數據及預處理

Sentinel-2A攜帶一臺多光譜成像儀(MSI)，幅寬達290 km，重訪周期10 d。從可見光和近紅外到短波紅外可覆蓋13個光譜波段，最高空間分辨率為10 m，在光學遙感數據中，Sentinel-2A數據是唯一在紅邊范圍含3個波段(中心波長分別為560、665、705 nm)的數據，為細致地監測植被生長情況提供了更多可能性。

通過實地走訪調查和物候期分析，選取2018年4個時期的Sentinel-2A衛星遙感影像(成像時間為5月15日、7月24日、8月23日和9月22日，產品級別為Level-1C)。數據來源于Copernicus Open Access Hub(https:∥scihub.copernicus.eu)。覆蓋研究區的圖像包括4景，編號分別為48SUG、48SVG、48SUF、48SVF。為保證數據質量，在進行作物精細分類前，利用Sen2cor軟件包對Sentinel-2A影像進行大氣校正；利用SNAP軟件將校正結果重采樣為ENVI格式，在ENVI軟件中選取藍綠紅3個可見光波段、紅邊波段1和近紅外波段(編號分別為2、3、4、5、8波段，其中第5波段空間分辨率為20 m，其他波段空間分辨率為10 m)進行波段疊加、影像拼接及裁剪，得到包含5個波段的研究區影像圖。根據土地利用變更數據，掩膜獲取耕地分布信息。

1.3 農作物訓練樣本及驗證數據集

2018年7月底對研究區開展野外實地考察，并記錄所到之處的地理位置和作物種類。利用實地調查數據，結合農事歷，對影像進行目視解譯，選取8 760個樣本點，包括水稻、玉米、春小麥、春小麥+秋油葵、胡麻、洋蔥、大棚作物、馬鈴薯和砂田瓜果9種作物類型。樣本點分布較為均勻，其中2/3作為訓練樣本，1/3作為驗證數據。

2 研究方法

2.1 研究思路

影像預處理之后，首先計算各個時期遙感影像的植被指數，構建特征參數數據集，分析不同作物不同指數的時序變化曲線；然后采用隨機森林分類法進行農作物識別得到分類結果，最后根據驗證樣本進行精度評價。具體流程如圖2所示。

圖2 技術流程圖Fig.2 Flowchart of technical route

2.2 NDVI和RENDVI計算及時序變化曲線提取

NDVI和RENDVI是農作物監測中常用的光譜指數，計算式為

(1)

(2)

式中NIR、R、VRE——Sentinel-2A數據的波段8、4、5的反射率

分別計算4個時期的兩種指數圖像，再進行組合。

由于不同指數特征或其組合對分類結果的精度影響及貢獻程度不同[30]，本文設計了3種特征指數構建方案(表2)，共提取了5種特征指數作為分類特征：NDVI、RENDVI、NDVI+RENDVI、NDVI-RENDVI和NDVI&RENDVI(即NDVI和RENDVI通過layer stack組合)。

利用不同時期指數圖，統計分析不同作物的特征指數時序變化曲線，反映不同作物類型的物候差異，分析不同作物類別的可分性以及時相選擇的合理性。

2.3 多種特征的隨機森林分類

隨機森林分類是組合多棵決策樹的預算結果，然后進行投票判斷，預測準確率高，是多種遙感分類方法中精度較高的一種方法[31]。

本文將5種特征指數圖像與訓練樣本一起輸入到隨機森林分類器，得到5種分類結果。

2.4 精度驗證

本文利用混淆矩陣進行分類結果的精度驗證。基于混淆矩陣，可以計算總體分類精度、Kappa系數、每一類別的制圖精度和用戶精度。

表2 特征指數構建方案Tab.2 Scheme of feature parameters

總體分類精度是指所有被正確分類的像元數量與總像元數量的比值，計算式為

(3)

Kappa系數計算式為

(4)

式中k——混淆矩陣行列數，代表分類的類別數量

xii——混淆矩陣中對角線上的值，代表正確分類的像元數

N——驗證的像元總數

xi+——混淆矩陣第i行元素相加之和

x+i——混淆矩陣第i列元素相加之和

3 結果與分析

3.1 不同作物指數特征的時間變化分析

基于訓練樣本數據，統計各種農作物的時序植被指數特征(NDVI、RENDVI)均值，時序植被指數變化曲線如圖3所示。

圖3 不同農作物的特征指數時序變化曲線Fig.3 VI curves of different crops

從圖3中可以看出，春小麥-秋油葵的特征曲線有2個峰，為一年兩熟的種植模式，而其他作物的特征曲線都只有1個峰，為一年一熟的種植模式；在各種特征曲線上，洋蔥、大棚作物、砂田瓜果3種作物都表現出較低值，并且生長期內起伏不大，尤其是砂田瓜果能與其他作物明顯區分。

根據圖3所示特征，9種作物大致可分為3種類型：高值型(水稻、玉米、胡麻、馬鈴薯)，曲線形態表現為鐘型，5月特征值較低且比較集中，7月達到峰值，8月略有下降，9月特征值分散下降。低值型(洋蔥、大棚作物、砂田瓜果)，曲線形態表現為比較平緩，全年特征值較低，7月有一小峰值。開口型(春小麥、春小麥-秋油葵)，曲線形態表現為開口型，即5月已有較高的特征值，9月仍處于較高值。

可以看出，5月是區分開口型和其他類型的最佳時期，7、8月是區分高值型和低值型的最佳時期。開口型春小麥、春小麥-秋油葵在5月中旬表現出較高的特征值，與其他類型作物的特征參數差異較大，可以明顯區分；高值型、低值型作物在5月下旬的特征參數比較接近，難以區分。從農事歷來看，大部分作物在5月處于苗期，植被指數值較低；到7月下旬，兩種類型的差異明顯分化，并保持到8月下旬；到9月下旬，大部分作物已成熟收割，高值型和低值型作物的植被指數值降低，兩種類型的分異趨緩。

高值型的4種作物在5、7月影像上指數值都比較接近，難以區分，根據農事歷，利用表現出成熟期時差的8、9月影像可以區分：馬鈴薯與胡麻成熟得較早，而水稻和玉米還處于生育期末期，因此馬鈴薯與胡麻的特征值低于水稻和玉米；低值型的砂田瓜果全年的特征值都低于其他作物，洋蔥和大棚作物彼此之間的差異也是顯著可見，洋蔥的特征值全年高于大棚作物；開口型的春小麥和春小麥-秋油葵的特征參數到8月下旬秋油葵長勢旺盛期可以明顯區分。

各種作物的RENDVI曲線總體形態與NDVI曲線相似，但曲線起伏趨緩，高值型的作物在生長旺盛期(7月)的特征值分異比NDVI更加明顯。

3.2 不同特征指數分類的精度對比

為了評估不同指數及其組合特征對分類精度的影響，本文采用隨機森林算法對NDVI、RENDVI、NDVI+RENDVI、NDVI-RENDVI、NDVI&RENDVI 5種特征指數分別進行分類，并利用混淆矩陣進行分類精度驗證，結果如表3所示。

5種分類結果中，采用NDVI&RENDVI特征組合分類的總體精度最高，其次是NDVI，最后是NDVI-RENDVI，NDVI&RENDVI分類精度比NDVI分類精度提高了約4個百分點，比NDVI-RENDVI分類精度高出約10個百分點。并且NDVI&RENDVI分類結果與驗證數據取得了較好的一致性，Kappa系數為0.83，NDVI次之，Kappa系數為0.78，最低的是NDVI-RENDVI，Kappa系數為0.71。

表3 各種分類結果總體精度及Kappa系數Tab.3 Overall accuracy and Kappa coefficient of different classifications

NDVI-RENDVI分類結果精度最低。從算法原理上分析，NDVI和RENDVI兩值相減將共同的紅波段消除而用近紅外與紅邊波段相減，削弱了植被由吸收谷紅光波段到反射峰近紅外的反射陡坎效應，導致NDVI-RENDVI分類的精度最低。

單獨利用RENDVI分類的精度較低，但是將RENDVI特征與NDVI特征組合到一起的精度有明顯提高。

從NDVI&RENDVI分類結果來看，春小麥、春小麥-秋油葵、砂田瓜果、馬鈴薯的分類效果最好(用戶精度大于90%)，胡麻、洋蔥、大棚作物分類效果相對較差(用戶精度小于80%)，被誤分為馬鈴薯、砂田瓜果的概率較高；在NDVI特征分類結果中，玉米、水稻、大棚作物3類的精度較低，而在NDVI&RENDVI特征組合分類結果中，這3類作物的用戶精度有了顯著提高，從而提高了總體的分類精度(提高了4.06個百分點)。

3.3 研究區作物種植結構與分布

景泰縣2018年農作物NDVI&RENDVI分類結果如圖4所示。

圖4 2018年景泰縣作物類型空間分布圖Fig.4 Spatial distribution of crops in Jingtai County in 2018

由圖4可以看出，具有干旱區特色的砂田瓜果(包括小金瓜、籽瓜等)播種面積最大，主要分布在海拔較高、氣溫偏低的正路鎮、寺灘鄉和喜泉鎮；大棚種植方式在景泰縣比較普遍，廣泛分布于全縣；玉米、洋蔥、馬鈴薯、胡麻是景泰縣重要的農作物，主要分布在北部和中部比較平坦的地區，紅水鎮、漫水灘鎮、條山鎮、草窩灘鎮等地；水稻面積較小，主要分布在水源條件較好的五佛鄉。 2018年景泰縣作物種植結構如表4所示。

4 結論

(1)特征曲線分析表明，根據作物物候期選擇的4個時期圖像可以較好地表現研究區作物的生長期差異，能有效區分不同作物類型。

(2)特征選取對分類精度有明顯的影響，NDVI-RENDVI分類精度較差，采用NDVI&RENDVI特征組合分類精度最高，較NDVI特征分類的總體精度高4.06個百分點，比RNDVI-RENDVI特征分類精度高10.56個百分點，說明這種特征組合的方式能有效提高分類精度。

表4 2018年景泰縣作物種植結構Tab.4 Crops structure of Jingtai County in 2018

(3)RENDVI特征輔助NDVI可以提高分類精度，單獨利用RENDVI分類精度不高，將NDVI和RENDVI組合在一起能夠明顯提高分類精度，說明Sentinel-2A特有的紅邊波段數據及其較高的空間分辨率在農作物精細分類上具有很大的潛力。