基于時序Sentinel-2A影像的玉米秸稈覆蓋區(qū)智能識別研究

陶萬成, 張 穎, 謝茈萱, 王新盛, 董 鐿,張明政 , 蘇 偉*, 李佳雨, 軒 阜

1.中國農業(yè)大學土地科學與技術學院，北京 100083 2.農業(yè)農村部農業(yè)災害遙感重點實驗室，北京 100083

引 言

黑土是寶貴的自然資源，是肥力最高、最適宜農耕和最具生產潛力的土壤。黑土地保護性耕作模式是保護黑土地，促進農業(yè)可持續(xù)發(fā)展的重要措施。保護性耕作是指作物收獲后地面上秸稈覆蓋度大于30%的耕作和種植系統(tǒng)[1]，針對黑土地的保護我國東北地區(qū)正積極實施、推廣保護性耕作模式[2]。秸稈還田不僅有利于資源合理利用，還能避免因秸稈露天焚燒造成的環(huán)境污染[3]，與傳統(tǒng)耕作方式相比，長期實施秸稈還田有助于降低土壤容重，提高土壤孔隙度、土壤有機碳和土壤大團聚體百分比等[4]，對土壤保護與資源可持續(xù)利用起著重要作用。在亟需對黑土地實施保護措施的形勢下，究竟哪些地塊、有多少黑土地實施了秸稈覆蓋，對黑土地保護相關科學研究和當?shù)卣咧贫ㄅc實施具有重要意義。

傳統(tǒng)作物秸稈覆蓋區(qū)地面調查和識別的方法費時、費力、成本高，難以在區(qū)域范圍內實施。遙感是快速獲取區(qū)域范圍內地表信息的技術手段，特別是歐空局(European Space Agency，ESA)發(fā)布的Sentinel-2衛(wèi)星影像具有對秸稈纖維素和木質素敏感的短波紅外波段(SWIR)，是快速、準確監(jiān)測區(qū)域范圍內黑土地上作物秸稈覆蓋的有效數(shù)據(jù)源。已有許多研究者根據(jù)秸稈以2 100 nm處的強吸收谷與木質素、纖維素的高度相關性，利用多光譜影像、地面實測秸稈高光譜數(shù)據(jù)構建光譜指數(shù)展開秸稈覆蓋度研究[5-7]。秸稈覆蓋度研究一般僅利用作物收獲后單一時相影像估算，無法表征不同秸稈覆蓋類型的實際秸稈覆蓋量，例如，高留茬秸稈覆蓋從覆蓋度的角度分析其覆蓋度并不大，但實際上是4～50 cm的留茬直立在地表，其實際覆蓋量并不小。因此，本研究進行秸稈覆蓋區(qū)識別，是進一步開展地塊尺度秸稈覆蓋的基礎，也是對秸稈覆蓋類型精細分類(如覆蓋還田、粉碎還田、高留茬等)的前提，通過遙感影像直接獲取高精度的玉米秸稈覆蓋區(qū)更具有實際意義。目前，利用遙感技術的秸稈覆蓋區(qū)識別研究還較少，且沒有綜合考慮作物生長期和收獲后的狀態(tài)信息。本工作基于作物不同狀態(tài)信息對玉米秸稈覆蓋區(qū)開展精細、智能識別研究。

玉米是我國東北黑土地上主要種植的一種作物；吉林省四平市內有針對黑土地保護創(chuàng)建的“梨樹模式”區(qū)域，廣泛存在少耕、免耕等保護性耕作方式，又大量存在翻耕等傳統(tǒng)耕作方式，非常適合開展秸稈覆蓋區(qū)識別研究。由于玉米秸稈覆蓋地塊的反射率較高，且具有大面積連通的自然特性，在Sentinel-2A影像上表現(xiàn)為亮度偏高特征，呈現(xiàn)區(qū)域規(guī)則連續(xù)分布。根據(jù)秸稈覆蓋大于30%時為保護性耕作的規(guī)定，遙感影像可視為秸稈覆蓋區(qū)和背景兩類。隨機森林通過給定樣本的學習訓練形成分類規(guī)則，具有分析復雜地理信息系統(tǒng)分類特征的能力，適用于分類和高維數(shù)據(jù)的回歸分析，同時具備更高的精度和穩(wěn)健性。在作物倒伏識別[8]和土地覆蓋分類[9]等方面，隨機森林都具有很好的分類效果和處理速度，因此受到廣泛關注?；谙袼氐碾S機森林算法分類結果中普遍會出現(xiàn)細小圖斑。為克服不足，應用形態(tài)學濾波去除圖斑。

鑒于此，本研究基于GEE云平臺[10-11]，綜合考慮玉米收獲前、后的狀態(tài)信息，結合研究區(qū)內2020年5月—11月份的時序Sentinel-2A影像創(chuàng)建的光譜指數(shù)和分位(quartile, QT)特征構建數(shù)據(jù)集，應用優(yōu)化后的隨機森林方法準確識別玉米秸稈覆蓋區(qū)域；根據(jù)秸稈覆蓋大范圍連通的特性，利用連通域標定法優(yōu)化分類結果，實現(xiàn)玉米秸稈覆蓋區(qū)制圖。

1 實驗部分

1.1 研究區(qū)概況

研究區(qū)為吉林省四平市，地理范圍為42°31′—44°09′N，123°17′—125°49′E，總面積約1.4萬km2，地理位置如圖1(a)所示。四平市處在黑土地和吉林省“黃金玉米帶”核心地區(qū)，從2007年開始探索并實施以“秸稈覆蓋、條帶休耕”為主要內容的保護性耕作。開展區(qū)域范圍內的玉米秸稈覆蓋區(qū)準確識別，可應用于評估保護性耕作模式的實施程度，為監(jiān)測保護性耕作推廣，相關政策擬定提供可靠的數(shù)據(jù)支持。

圖1 研究區(qū)位置與Sentinel-2A影像(日期：2020年11月4日)(a)與秸稈覆蓋區(qū)照片(b)

1.2 數(shù)據(jù)源

1.2.1 地面調查數(shù)據(jù)

2020年11月3日—11日，在研究區(qū)進行了田間玉米秸稈覆蓋實地調研。11月為研究區(qū)內玉米收獲后秸稈覆蓋的穩(wěn)定時期，存在大量的樣本地塊。根據(jù)Sentinel-2A影像的最高空間分辨率(10 m)，在樣本地塊中隨機選取10個分布均勻、大小為1 m×1 m的樣方，采用拉繩法測量。秸稈調查時使用兩根以0.1 m間隔作為標記的1 m刻度繩，垂直擺放在樣方中，計算標記點與地面秸稈相交點數(shù)量與總標記數(shù)比值，統(tǒng)計10個樣方的均值；依據(jù)保護性耕作的定義[1]對秸稈覆蓋區(qū)進行識別，圖1(b)中當統(tǒng)計的樣方均值大于30%時為秸稈覆蓋區(qū)，主要由機械收獲后產生，包括低和高留茬；低于30%時為非覆蓋區(qū)，主要是人工收獲后產生的根茬。為了識別玉米秸稈覆蓋區(qū)，采用實地調查和目視解譯相結合的手段獲取樣本數(shù)據(jù)集，在剔除異常值后保留了646個樣本：包括玉米秸稈覆蓋、水稻秸稈覆蓋、建筑用地、樹林和水體樣本，樣本數(shù)量分別為200，103，118，137和88。

1.2.2 Sentinel-2A遙感影像

從Google Earth Engine(GEE)云平臺獲取時序Sentinel-2A影像，時間分辨率為5 d，光譜包括可見光、近紅外、紅邊、水汽、卷云和短波紅外波段，各波段分辨率分別為10，20和60 m。為避免低分辨率波段對分類任務的影響，舍棄空間分辨率為60 m的水汽、卷云和可見光中的氣溶膠波段，選擇B2(Blue)，B3(Green)，B4(Red)，B5(Red-edge1)，B6(Red-edge2)，B7(Red-edge3)，B8(Nir)，B8A(Red-edge4)，B11(SWIR1)和B12(SWIR2)較高空間分辨率波段。

利用GEE獲取研究區(qū)內2020年的云占比小于10%的Sentinel-2A影像；同時選取20個實地調查的玉米秸稈覆蓋樣本點，統(tǒng)計年內歸一化差值植被指數(shù)(normalized difference vegetation index, NDVI)和歸一化差值秸稈指數(shù)(normalized difference residue index, NDRI)特征曲線和置信區(qū)間，如圖2所示。S1為玉米出苗期，H為玉米收獲期，S2為玉米秸稈覆蓋期，橫坐標為年積日(day of year, DOY)，圖2(a)縱坐標為NDVI均值(NDVIM)和NDVIM+/-標準差(Std)，圖2(b)縱坐標為NDRIM+/-Std。圖2(a)中DOY為204時達到NDVIM最大值0.864。圖2(b)中NDRI和相應時間區(qū)間的NDVI趨勢相似，但值范圍相差較大。為降低計算量，選取從5月13日到11月29日的182景云占比小于10%的Sentinel-2A作為影像源，充分合理地利用時間序列特征。根據(jù)研究需求，對獲取影像進行波段篩選、拼接和裁剪等預處理。

圖2 研究區(qū)2020年內玉米種植區(qū)時序NDVI和NDRI變化曲線

1.3 方法

基于時序Sentinel-2A的玉米秸稈覆蓋區(qū)識別包括3個步驟。首先，基于Sentinel-2A時序影像計算光譜和指數(shù)特征，同時利用分位法獲取分位QT特征，進而構建數(shù)據(jù)集；然后，對樣本數(shù)據(jù)按7∶3隨機劃分為訓練集和測試集，結合RF模型分類出玉米秸稈覆蓋區(qū)域，而基于像素的RF分類結果中極易產生細小圖斑，研究采用連通域標定法優(yōu)化全局；最后，通過評價指標定量分析結果。具體技術流程如圖3所示。

圖3 玉米秸稈覆蓋區(qū)識別技術流程圖

1.4 數(shù)據(jù)集構建

1.4.1 光譜特征

對研究區(qū)內2020年11月4日獲取的Sentinel-2A多光譜反射率影像進行統(tǒng)計分析得知，建筑物光譜特征較明顯、容易識別，因此重點分析玉米秸稈覆蓋區(qū)、水稻秸稈覆蓋區(qū)、樹林和水體樣本光譜特征。統(tǒng)計樣本內各地類波段反射率值并繪制箱線圖，結果如圖4(a—d)。

從圖4(a,b)可以看出，玉米秸稈和水稻秸稈的光譜反射率區(qū)間范圍為0.10～0.45，圖4(c,d)水體和樹林光譜反射率區(qū)間范圍分別為0.00～1.05和0.00～0.30。通過對比發(fā)現(xiàn)，水體和樹林相對秸稈反射率明顯較低，易于分類。水稻秸稈和玉米秸稈的光譜反射率變化趨勢和區(qū)間范圍基本相同，B2—B8波段反射率呈上升趨勢，B8A—B12波段呈下降趨勢，僅利用單期遙感影像很難進行地類判別。鑒于此，考慮到不同作物的物候和結構特性，本研究使用7月22日Sentinel-2A影像提高水稻和玉米地類的分離度，即在作物生長季中玉米種植區(qū)的NDVI值大于水稻。為了判斷該期影像是否有助于區(qū)分不同作物，利用樣本統(tǒng)計繪制出玉米和水稻的光譜反射率箱線圖，結果如圖5(a,b)。圖5(a)中玉米B2—B5波段反射率明顯低于圖5(b)水稻，B7—B8A反射率高于水稻，其他波段反射率相似，其原因可能與兩種農作物的冠層結構相關，這些特征有利于分類。因此，光譜特征集是結合7月22日和11月4日的Sentinel-2A影像建立。

圖4 不同地物類型光譜反射率特征(2020年11月4日Sentinel-2A影像)

圖5 不同作物光譜反射率特征(2020年7月22日Sentinel-2A影像)

1.4.2 指數(shù)特征

為區(qū)分玉米秸稈覆蓋區(qū)、水稻秸稈覆蓋區(qū)、裸地和其他相似地類，分別構建了差值植被指數(shù)(NDVI)和秸桿指數(shù)(NDRI)。其中，NDVI可以較好地反映出不同農作物的生長狀況，因此用于區(qū)別玉米和水稻種植區(qū)域。NDRI光譜指數(shù)與NDVI相似，同樣是利用雙波段進行計算，該光譜指數(shù)與玉米秸稈覆蓋度顯著相關[12]，其計算公式如式(1)

(1)

式(1)中，ρB12是短波紅外2反射率，ρB4為紅波段反射率。在玉米秸稈中含有大量的纖維素和木質素，反射率在2 100和2 300 nm附近有強吸收，即Sentinel-2A影像的B12波段反射率。使用Sentinel-2A的B12和B4波段反射率構建的NDRI可以最大限度地減少綠色植被的影響。綜上分析，在光譜特征集的基礎上，利用NDVI和NDRI建立指數(shù)特征集。

1.4.3 分位特征

分位QT特征的建立主要基于四分位數(shù)的思想。在統(tǒng)計學中，將所有數(shù)值按大小排列并分成四等份，處于三個分割點位置的數(shù)值為四分位數(shù)。同理，四分位數(shù)法可應用于時序遙感影像以構建QT特征，流程如圖6所示。首先，通過GEE云平臺獲取時序影像集，圖6(a)中為11月4日影像中藍波段熱度圖，顏色越深反射率值越小，相反則反射率值越大；以中間紅框內的像素為例，抽取該位置時間維度上的數(shù)據(jù)，如圖6(b)；對獲取的數(shù)據(jù)按照從小到大進行排序，V1—V5分別對應0%，25%，50%，75%和100%，依次選擇對映分位上的反射率值，如圖6(c)；最后，按照上述方法對影像中像素遍歷處理獲取光譜分位集，如圖6(d)。

圖6 分位QT特征構建示意圖

QT特征還可以應用于時序指數(shù)特征集，與計算光譜分位集相似，首先獲取時序指數(shù)特征集，然后應用QT特征提取指數(shù)分位集。QT特征集構建包括光譜分位集和指數(shù)分位集，這種提取特征的方式利用時序特性的同時避免了數(shù)據(jù)冗余，極大地降低了計算量。

1.5 隨機森林秸稈覆蓋區(qū)識別

隨機森林監(jiān)督分類算法是一種基于多棵CART決策樹組合構成的集成學習方法。為了完全劃分變量空間，利用自助抽取輸入樣本和節(jié)點隨機分裂技術構建多棵決策樹。通過單個決策樹的多數(shù)投票策略獲得預測結果，有效避免單個模型或參數(shù)組引起的誤差，具有更高的精度和魯棒性[13]。由于其分類的優(yōu)越性，該分類器得到遙感界的廣泛認可。

針對所設計的數(shù)據(jù)集，具有高維度特性，隨機森林可以將每個維度的重要性輸出到分類中，幫助特征選擇以提高效率。因此，本研究使用嵌入在GEE中的隨機森林算法對玉米秸稈覆蓋區(qū)進行分類。根據(jù)分類結果的準確性調整決策樹數(shù)量以期生成最佳的分類模型。

1.6 連通域標定和評價指標

研究區(qū)中秸稈覆蓋地塊一般具有大面積連通的自然特征，而本研究基于像素級別對玉米秸稈覆蓋區(qū)進行分類，結果中會出現(xiàn)細碎的圖斑，因此有必要借助連通域標定法對結果全局優(yōu)化。連通域標定法示意圖見圖7，假定圖7(a)為分類結果中某一區(qū)域，標簽1表示為玉米秸稈覆蓋區(qū)，標簽0表示背景類，綠色區(qū)域為設定的4連通域滑動窗口。通過滑動窗口尋找標簽為1的連通域，保留圖7(b)左上角最大連通域(黃色區(qū)域)，刪除右下角最小連通域(黃色區(qū)域)。如果圖像中存在大量目標區(qū)域，則需標定所有的連通域，通過設置閾值對全局處理。

圖7 連通域標定示意圖

Kappa系數(shù)是衡量分類結果空間一致性的指標，明確揭示了分類結果的空間變化?？傮w精度(overall accuracy, OA)以預測正確與總體數(shù)量之間的比值反映分類結果的質量[14]。本研究選擇上述指標來定量定性評價分類結果。

2 結果與討論

基于時序Sentinel-2A影像設計的數(shù)據(jù)集，采用隨機森林算法識別玉米秸稈覆蓋區(qū)。為確定算法中最佳決策樹數(shù)量，研究對比5，10，20，30和40決策樹對分類模型的影響，具體見表1。從5到30整體呈上升趨勢，30達到最大值并趨于穩(wěn)定，當決策樹數(shù)量設置為40時，評價精度略有下降。因此，實驗統(tǒng)一標準，將決策樹數(shù)量均設置為30。

表1 不同決策樹的分類模型定量評價

2.1 基于不同特征集組合分類

在相同訓練樣區(qū)條件下，通過不同的特征集組合識別玉米秸稈覆蓋區(qū)，定量評價結果見表2。以光譜特征集為基礎，模型M5精度最高，Kappa/OA為97.41%/97.91%，比M1，M2和M4模型的Kappa/OA分別提高了5.84%/4.69%，4.52%/3.64%和1.29%/1.04%，表明相比于指數(shù)特征集，加入QT特征集可以明顯改善模型精度；以QT特征集為基礎，M4比M3模型提高1.95%/1.56%，表明與指數(shù)特征集相比，加入光譜特征集更有利于分類；以指數(shù)特征集為基礎，M3比M2模型提高1.28%/1.04%，表明加入QT特征集優(yōu)于光譜特征集。

表2 基于不同特征集的分類定量評價結果

本研究使用五種不同特征集組合模型，測試結果見圖8，圖8(a)和圖9(a)分別為2020年11月4日和7月22日真彩色合成圖，分類結果中藍色和白色區(qū)域分別表示玉米秸稈覆蓋區(qū)和背景。根據(jù)圖8(b)—(f)對比發(fā)現(xiàn)，M5模型一定程度上降低細小圖斑的產生，同時保留道路、農田等邊緣細節(jié)信息，效果最好；M1和M2模型分類結果較差，易產生噪聲，將道路等地類分為秸稈覆蓋區(qū)；M3和M4模型分類結果相對較好，保留了部分邊緣細節(jié)信息。實驗結果與定量評價結果基本一致，QT特征的加入極大改善了分類可視化結果。

圖8 基于不同特征集的分類可視化結果

2.2 基于不同時序尺度的特征分類

通過上述實驗證明了QT特征的重要性，為了評估不同時序尺度生成的QT特征集對識別秸稈覆蓋區(qū)的影響，本研究基于M5模型，對利用不同QT特征集的M5_1—M5_6模型定量分析，具體信息見表3。

根據(jù)表3，利用5月—11月影像創(chuàng)建的QT特征M5_6/M5模型分類效果最好。其中，M5_2比M5_1模型的kappa/OA低0.97%/0.77%，可能是10月為玉米收獲期，秸稈地塊中包括收獲與未收獲，收獲的地塊又包括秸稈成捆或散放在地里等，信息多樣，識別困難；M5_3相比于M5_2模型精度低，分析認為9月份研究區(qū)受到臺風影響，缺少質量好的影像。而M5_4，M5_5和M5_6精度相繼提升表明長時序影像創(chuàng)建的QT特征集有利于識別秸稈覆蓋。

表3 基于不同時間尺度特征分類定量評價結果

2.3 結合連通域標定法分類

在M5模型基礎上，結合連通域標定法構建模型M6對全局優(yōu)化，設置連通域閾值為30，M6分類結果的Kappa/OA為96.76%/97.36%。相比于M5模型精度有所降低，但仍高于其他模型，同樣適用于玉米秸稈覆蓋區(qū)識別。選取兩組測試結果對比如圖9和圖10。

根據(jù)圖9(c)，M5[圖9(d)]和M6[圖9(b)]相比于M5_1模型分類結果可以很好的保留細節(jié)信息，進一步證明QT特征設計的有效性；M6相比于M5模型在保留絕大部分細節(jié)信息的同時，抑制了細碎圖斑的產生。圖10(a)是夏季的子圖塊，其中淺綠色為水稻地塊，圖10(b)是11月初秸桿覆蓋的子圖塊通過對比圖10(c)M5_1，圖10(d)M5和圖10(e)M6分類結果發(fā)現(xiàn)，在僅利用11月影像的M5_1模型容易將水稻秸稈覆蓋區(qū)分為玉米秸稈覆蓋區(qū)，而利用5月—11月份影像的M5模型可很好地將兩者區(qū)分開。此外，M5模型結果中有相對較少的細碎圖斑，結合連通域標定的M6模型在優(yōu)化細碎圖斑后得到更理想的分類結果，見圖10(e)。

圖9 結合連通域標定分類結果1

圖10 結合連通域標定分類結果2

通過實驗對比M6模型可視化效果最好，定量評價結果僅次于M5模型，因此研究使用M6模型對整個研究區(qū)進行預測，得到的分類結果見圖11。其中，左側為整個研究區(qū)識別結果，右側為研究區(qū)框選出的子區(qū)域識別結果，圖11(a)主要包括秸稈和建筑，圖11(b)主要包括秸稈，建筑，水域和樹林。從不同區(qū)域識別結果可以看出秸稈覆蓋區(qū)邊緣信息保持完好，幾乎沒有細碎圖斑，證明了所提出方法的有效性，適用于玉米秸稈覆蓋區(qū)識別。

3 結 論

以吉林省四平市為例，基于GEE云平臺，首先對時序Sentinel-2A影像預處理，利用光譜，指數(shù)和QT特征構建數(shù)據(jù)集。然后應用最佳決策樹個數(shù)的隨機森林分類器分析不同特征組合，不同時序尺度的QT特征和是否結合連通域標定對識別玉米秸稈覆蓋區(qū)域的影響。最終采用可視化效果最好，精度較高的M6模型對研究區(qū)影像預測，主要研究結論如下：

(1)利用不同特征集組合的分類模型所設計的分類模型M5定量評價結果最優(yōu)，通過模型對比表明加入QT特征能夠有效提升精度，保留邊緣細節(jié)信息，在一定程度上降低了噪聲出現(xiàn)的概率，優(yōu)于未加入QT特征的分類結果。

(2)基于不同時序尺度的QT和光譜特征集的分類模型，實驗利用5月—11月時序影像內QT特征比較短時序內影像的QT特征分類效果好，表明長時序特征有利于玉米秸稈覆蓋區(qū)識別，同時可以避免其他作物秸稈覆蓋區(qū)域的干擾。

(3)在模型M5分類結果的基礎上，結合連通域標定法的M6模型對全局連通域進行識別，通過設置閾值優(yōu)化結果，實驗結果表明該方法在保證較高玉米秸稈覆蓋區(qū)域識別精度的同時去除了細碎圖斑。

通過本研究提出M6模型可以準確的從Sentinel-2A影像中識別出玉米秸稈覆蓋區(qū)，對當?shù)乇O(jiān)測保護性耕作實施和相關服務政策的擬定提供便利。后續(xù)研究主要有兩點：① 將著重挖掘影像空間域紋理特征在玉米秸稈覆蓋識別中的應用潛力。②嘗試多源遙感影像融合進一步優(yōu)化結果。