基于Sentinel-2時序諧波特征的縣域農作物分類

孫慶松， 張曉楠， 陳利東， 王 坤， 宋宏利

(1.河北工程大學地球科學與工程學院，河北 邯鄲 056038； 2.河北工程大學礦業(yè)與測繪工程學院，河北 邯鄲 056038； 3.河北省地礦局第六地質大隊，河北 石家莊 050000)

農作物精細化分類是農作物長勢監(jiān)測、產量預測及災害評估的重要前提和基礎。及時準確地獲取農作物類型、面積和空間分布信息，對于優(yōu)化農作物種植結構、實現農業(yè)生產精準化管理及確保國家和地區(qū)糧食安全具有重要意義[1-3]。傳統(tǒng)農作物類型信息的獲取主要通過實地走訪調查，然而進行實地精細調查費時費力，且容易受到人為因素干擾，因此更新不夠及時。隨著對地觀測技術的發(fā)展，利用遙感數據提取農作物信息已經成為重要手段。

目前，越來越多的研究者利用植被光譜特征的季節(jié)性變化來提高作物分類的準確性，已有學者利用多時相光學遙感數據提取增強型植被指數(EVI)、歸一化植被指數(NDVI)、土壤調節(jié)植被指數(SAVI)等時序指數集開展農作物分類研究，并取得了較好的分類效果[4-8]。目前國內農作物信息提取主要是依據指數特征進行分類，利用諧波特征提取農作物信息的研究較少。諧波特征分類優(yōu)勢在于并不僅僅依靠某時相內指數值的差異，同時利用不同作物時序曲線中唯一對應的相位值和諧波余項值進行分類[9]。諧波分析[10-12]技術最早應用于電力信號分解和水文氣候預測。在作物信息提取方面最初由Jakubauskas等[13]提出，他們利用諧波分析技術對26個年內NOAA氣象衛(wèi)星傳感器(AVHRR)獲取的NDVI時間序列進行分析，并提出了一種基于NDVI值隨時間變化的作物識別方法。

目前常用的遙感數據源主要包括Landsat[14]、GF系列[15]、MODIS[16]及多源融合數據[17]。與Landsat相比，Sentinel-2衛(wèi)星具有高時間分辨率和高光譜分辨率，同時具有3個紅邊波段，近年來被廣泛應用于農作物分類研究[18-19]，并展現出其在農作物分類以及農情監(jiān)測應用中的優(yōu)勢[20-21],為縣域尺度范圍農作物信息提取提供了新的選擇。

以往研究主要集中于影響因素可控的小區(qū)域[22]和地勢平坦的西北及東北區(qū)域[23-24]，而在地塊破碎、種植結構復雜的華北縣域范圍內，關于農作物精細化分類的研究相對較少。因此本研究結合中國探索實行耕地輪作休耕制度的[25]背景，以中國地下水漏斗區(qū)黑龍港流域南宮市為研究對象，針對縣域尺度范圍內農作物種植信息提取精度低的問題，利用多時相Sentinel-2數據，通過構建指數時序集及其諧波特征，組合多種分類方案，對研究區(qū)內農作物信息進行提取，以期為縣域尺度農作物分類研究提供參考。

1 材料與方法

1.1 研究區(qū)概況

南宮市位于河北省東南部(37°05′～37°27′N，115°08′～115°45′E)，黑龍港流域南部(圖1)，屬黃河沖積平原，地面平均海拔28.6 m，地勢東南稍高西北低，坡降為1/7 000，地形相對平坦。屬暖溫帶亞濕潤季風型氣候，光照充足，雨熱同季，年均降水量476 mm。全市境內共有3條主干河道，24條分支，均屬于黑龍港水系，總長367 km。全市面積863.3 km2，總耕地面積597.3 km2，占全市面積的69%，主要農作物有冬小麥、棉花、夏玉米、谷子、紅薯、辣椒等，是華北地區(qū)重要產糧基地之一。

圖1 研究區(qū)及樣點分布Fig.1 Distribution of the research area and sample points

南宮市主要農作物是冬小麥、棉花、夏玉米、谷子、辣椒、紅薯、花生，其中，冬小麥-夏玉米是一年兩季傳統(tǒng)輪作模式，其他農作物均屬于一年一季。根據實地走訪調查及邢臺市統(tǒng)計局提供的農耕數據，總結出當地各類農作物物候歷(表1)。

表1 南宮市主要農作物物候期

1.2 數據源及預處理

本研究使用的遙感數據是歐洲航天局(ESA)(https：//scihub. copernicus. eu/dhus/#/home)免費提供的Sentinel-2多光譜數據。Sentinel-2包含A和B兩顆衛(wèi)星，單顆衛(wèi)星重訪周期為10 d，兩顆衛(wèi)星同時運行重訪周期為5 d，Sentinel-2衛(wèi)星最高空間分辨率為10 m，具有13個波段數據。為保證能夠完全覆蓋研究區(qū)內各類作物生長周期，選取2019年10月至2020年10月共63幅云量低于20%的光學影像，影像獲取時間信息如表2所示，其中影像最大時間間隔為30 d，最小為2 d。最大時間間隔在1月份，處于冬小麥的生長周期內，由于冬小麥在2月至6月期間，其植被指數值與其他作物有明顯差異，因此在1月份時，影像間隔對分類結果的影響較小。

表2 遙感影像數據

采用經幾何精校正的L1C級產品數據，根據植被指數計算公式，研究使用藍、綠、紅、2個紅邊波段及近紅外波段數據(其中2個紅邊波段空間分辨率為20 m，可見光及近紅外波段空間分辨率為10 m)。利用ESA提供的Sen2cor插件對L1C級數據進行大氣校正，并將所需波段的空間分辨率和反射率利用三次卷積法重采樣為10 m的分辨率和反射率。為便于后續(xù)調用各波段數據計算植被指數，在SNAP軟件中將同一影像中的波段數據合成為一個文件，并將文件轉換為ENVI可讀取的TIFF格式，之后在ENVI軟件中進行裁剪、鑲嵌、植被指數計算、最大值合成等操作，進而構建植被指數時間序列。利用MATLAB及IDL工具包提取時序數據中的諧波特征(對應各階振幅值圖像)，之后利用隨機森林分類器進行圖像解譯。

1.3 實地調查數據

2019年10月至2020年9月在南宮市進行多次實地調查。為避免手持全球定位系統(tǒng)(GPS)定位精度問題導致的作物種類邊緣誤差，在樣本選擇時，盡可能選擇地塊面積在30 m×30 m以上區(qū)域，在地塊中心用手持GPS接收機標記地塊經緯度坐標及地塊類型，樣本點均勻分布于全市(圖1)。共采集實地樣本1 090個，其中冬小麥547個，棉花282個，谷子118個，花生36個，辣椒12個，紅薯10個，建筑用地及裸地、水體和其他地類共85個。由于研究區(qū)內花生、辣椒、紅薯種植數量較少，為減少樣本數量懸殊導致的分類誤差，故合并為一類“其他農作物”。將2/3的實地采集數據作為訓練樣本，1/3的實地采集數據作為驗證樣本。

1.4 構建植被指數時序曲線

根據作物光譜反射特性，結合Sentinel-2衛(wèi)星中心波長及波段范圍，為獲取最佳組合分類方案，構建5種植被指數，分別為EVI、NDVI、NDVI705[26]、EVI+NDVI、EVI+NDVI+NDVI705，在構建2種組合指數方案時，采用指數值算數相加原理，其計算公式及構建規(guī)則見表3。計算每幅影像的植被指數后，采用最大值合成法(從相鄰影像中提取最大像元值重構一幅新的影像)，將全年63幅影像合成為30幅，確保每月有2～3幅影像，在保證時序數據集質量的同時，能最大程度減少云霧干擾，從而構建作物生長周期變化規(guī)律曲線。

3種植被指數(EVI、NDVI、NDVI705)被廣泛應用于作物估產、環(huán)境監(jiān)測、火災評估等領域。NDVI指數對于波段范圍限制較小，應用廣泛，能夠將植被與其他地物明顯區(qū)分，由于其本質是近紅外波段和紅外波段的非線性拉伸，數值容易飽和，導致在高密度植被區(qū)靈敏度降低。與傳統(tǒng)NDVI不同，NDVI705指數計算選用了葉綠素吸收特征較窄的波段，更容易受到植被內葉綠素含量的影響，在植被生長變化過程中，數值變化更加敏感[26]。EVI指數的提出是為了改進NDVI在植被高密度區(qū)靈敏度低的問題，加入藍光波段，以減少大氣折射造成的影響，其時序數據被廣泛應用于林地變化監(jiān)測研究。

表3 植被指數計算公式及構建規(guī)則

1.5 諧波特征

諧波分析又稱傅里葉分析，是指將一個隨時間變化的周期性函數通過傅里葉變換分解為無窮多個諧波分量，每個分量由頻率相同的正弦波和余弦波組成。在農業(yè)遙感領域，諧波分析的本質是將單維空間域影像轉換為多維頻率域數據，對頻率域內的波段數據進行級數分解，將不同植被物候變化相關的信號分隔開，最后再將頻域影像轉換為空間域影像。

諧波主要特征在于前三階諧波的振幅值和諧波余項，三階之后的諧波特征多是各類噪音誤差[10](圖2)。圖2a中關于時間的任意曲線(時間函數)都可以分解為圖2b中的正弦和余弦波，圖2c中冬小麥-夏玉米在NDVI時序中的函數可以分解為圖2d中的諧波特征分量，圖2b和圖2d諧波曲線的區(qū)別在于振幅值和諧波余項的不同，即諧波特征不同。Jakubauskas等[13]在2001年將諧波分析法應用于NDVI時間序列分解，并根據每個諧波分量在原始數據中占總方差的百分比進行農作物分類。諧波分析數學分解公式為：

(1)

式中，f(x)為時間序列函數，a0是諧波余項，L是周期函數的時間長度，x表示時間，an是第n次諧波余弦函數的振幅，bn是第n次諧波正弦函數的振幅。各類作物在NDVI時間序列中的諧波余項和振幅值如表4所示。

圖2 諧波特征示意Fig.2 Schematic diagram of harmonic characteristics

表4 各類作物諧波特征值

1.6 分類方法

2 結果與分析

2.1 不同作物時序曲線變化分析

為探究各類農作物植被指數值在生長周期內的變化趨勢，根據實地調查樣本數據，多次提取各類農作物植被指數值并取平均值，構建農作物生長變化曲線。

圖3表明，冬小麥-夏玉米的時序曲線在3個指數波段范圍內與其他地類均有明顯差異。冬小麥在10月上旬播種后，育苗期一直持續(xù)到12 月上旬，同時植被指數值升高；12月中旬之后進入越冬期，氣溫降低，光合作用緩慢從而導致植被指數值降低；次年1月中旬至3月中旬冬小麥進入返青發(fā)育期，植被指數值迅速增長；在3月中下旬時期，該地區(qū)受到北方寒潮影響，溫度驟降，植被指數值有小幅度下降；在4月至5月期間，冬小麥進入孕穗抽穗期，植被指數值迅速升高，達到生長發(fā)育期內的植被指數最大值，該時期冬小麥與其他作物物候差異最大，可以有效將冬小麥與其他作物區(qū)分開。夏玉米在6月份播種，之后進入育苗期，各類指數值呈上升趨勢，7月至8月中旬進入抽穗灌漿期，指數值到達第2個峰值，9月份以后進入成熟收獲期，指數值緩慢下降。

a:2019-10-05; b:2019-10-15; c:2019-10-30; d:2019-11-04; e:2019-11-14; f:2019-11-19; g:2019-12-04; h:2019-12-19; i:2020-01-03; j:2020-01-30; k:2020-02-17; l:2020-02-22; m:2020-03-10; n:2020-03-13; o:2020-03-18; p:2020-03-23; q:2020-04-12; r:2020-04-17; s:2020-04-22; t:2020-04-27; u:2020-05-02; v:2020-05-22; w:2020-06-03; x:2020-06-21; y:2020-07-06; z:2020-07-16; A:2020-08-22; B:2020-08-30; C:2020-09-04; D:2020-09-19。圖3 各類農作物時序曲線Fig.3 Time series curves of different crops

棉花在4月中下旬播種，5月中下旬開始出苗，在此期間，植被指數值沒有明顯變化。之后經歷現蕾、開花階段，植被指數值急速上升。7月份開花以后，其植被指數值基本維持穩(wěn)定，與其他非輪作農作物時序曲線在3種指數波段范圍內均有明顯差異，開花之后的指數值可以作為識別棉花的重要特征。9月份進入成熟期，植被指數值下降。谷子在5月中下旬播種，在出苗、拔節(jié)、抽穗時期，植被指數值呈上升趨勢，8月份成熟，此時植被指數值達到峰值。

辣椒和紅薯均在4月下旬播種，到6月上旬，2種作物植被指數值無明顯變化；6月下旬辣椒進入分蘗期，紅薯進入現蕾期，植被指數值上升；到8月份，植被指數到達峰值，成熟以后其植被指數值緩慢下降?；ㄉ?月中旬播種，進入出苗階段，植被指數值開始增長；6月份至7月份為幼苗期，植被指數值迅速升高；到8月份結莢時期，植被指數值保持相對穩(wěn)定；9月份成熟時收獲，植被指數值緩慢下降。由于這3類作物種植面積較小，因此在分類時，將其合并為一類。由圖3可以看出，依靠單時相影像數據，并不能將生長周期相似的農作物明顯區(qū)分，而時間序列數據可以凸顯物候差異特征，從而提高分類準確率。

2.2 不同指數組合的農作物分類方案精度對比

為探討不同指數組合及諧波特征對農作物分類精度的影響，選用EVI、NDVI、NDVI705、EVI+NDVI、EVI+NDVI+NDVI7055種植被指數，構建指數時序曲線，并提取作物生長曲線諧波特征。將分類特征和樣本數據輸入到隨機森林分類器中，結合實地采集的驗證樣本，采用混淆矩陣驗證分類結果，評價指標包括總體精度(Overall accuracy, OA)、用戶精度(User accuracy, UA)、生產者精度(Producer accuracy, PA)、Kappa系數，結果如表5、表6、表7所示。

表5 各種分類方案總體分類精度及Kappa 系數

從原始分類結果中可以看出，當指數特征和諧波特征共同作為分類依據時，相比于僅利用指數特征進行分類，5種時序集的總體分類精度均有明顯提高(表5)。其中EVI+NDVI時序集提升最為明顯，提高了9.21個百分點，其次是EVI+NDVI+NDVI705時序集，提高了8.75個百分點，NDVI705時序集提升最不明顯，提高了8.14個百分點，表明諧波特征的加入可以有效提高分類精度。

為進一步驗證不同植被指數對于農作物分類精度的影響，選取指數特征和諧波特征共同作為分類依據的5種指數組合方案，并對其精度評價指標進行詳細分析，結果見表6、表7，5種原始分類結果如圖4所示。

表6 3種單一指數的農作物分類方案精度評價結果

表7 2種組合指數的農作物分類方案精度評價結果

表6表明，以3種單一植被指數作為特征波段時，采用EVI指數特征作為分類依據的總體精度最低，為83.82%，利用NDVI705指數分類時總體精度最高，其值為90.00%，NDVI分類精度比EVI分類精度提高了約5個百分點。NDVI705分類結果說明，實際驗證樣本與模型分類結果具有高度一致性，Kappa系數為0.88；NDVI其次，Kappa系數為0.86；EVI最低，Kappa系數為0.80。在3種分類方案中，農作物中冬小麥-夏玉米的用戶精度和生產者精度均為最高，從種植模式上分析，研究區(qū)內在1月份至6月份期間只有冬小麥種植(圖3)，因此物候特征明顯，分類精度較高；棉花和谷子的分類精度略低，主要是由于其生長周期相似，且兩類作物種植結構呈嵌套式分布(圖4)，容易造成錯分情況；建筑用地、裸地和水體的光譜反射率并無明顯季節(jié)變換規(guī)律，與農作物物候特征有明顯差異，因此分類精度較高。NDVI705指數比NDVI指數分類精度高，從算法原理上分析，較窄的紅邊波段相比于紅光波段對于植被生長變化程度更加敏感，說明Sentinel-2特有的紅邊波段數據在農作物精細化分類上具有較大的應用潛力。

2種組合指數分類方案中，EVI+NDVI+NDVI705作為特征組合時，總體分類精度最高,為94.95%，比EVI+NDVI分類精度提高了約2個百分點(表7)，說明紅邊波段指數NDVI705的加入能夠明顯提高分類精度；EVI+NDVI組合指數比NDVI指數分類精度提高了約4個百分點，表明具有藍光波段的EVI指數的加入能夠提高分類精度。

EVI:增強型植被指數；NDVI：歸一化植被指數；NDVI705：紅邊歸一化植被指數。圖4 不同指數組合的農作物分類效果Fig.4 Crop classification effects of different index combinations

2.3 南宮市作物種植空間分布

南宮市農作物分類總體精度最高的指數組合方案是EVI+NDVI+NDVI705，其分類效果見圖4。南宮市2020年冬小麥-夏玉米、棉花、谷子、其他農作物的種植面積分別為21 580.77 hm2、16 968.16 hm2、3 000.81 hm2、403.36 hm2。種植面積最大的作物是冬小麥-夏玉米，主要分布在西北部和東南部；其次是棉花，主要位于中部和南部。冬小麥-夏玉米和棉花是當地重要的經濟作物，因此呈現大面積塊狀分布。谷子的種植分布比較零散，且地塊面積較小，多聚集在居民地附近，主要與當地的飲食習慣有關。

3 結論與討論

本研究基于Sentinel-2光譜數據，通過構建多種植被指數及其諧波特征，采用隨機森林算法，對南宮市內的主要農作物進行分類提取，結果表明：1)當指數特征和諧波特征共同作為分類依據時，相比于利用單一指數特征進行分類，5種時序集的總體分類精度均有明顯提高，EVI+NDVI時序集提升最為明顯，提高了9.21個百分點，NDVI705時序集提升最不明顯，提高了8.14個百分點，表明時序集中諧波特征的加入可以有效提高分類精度。2)在指數特征和諧波特征共同作為識別依據的分類方案中，基于EVI+NDVI+NDVI705組合指數的分類方案精度最高，比EVI+NDVI分類精度提高了2.57個百分點，而EVI+NDVI組合分類精度比NDVI指數分類精度提高了3.53個百分點。以上結果表明，紅邊NDVI705和EVI的加入能夠提高作物識別精度。3)通過對3種單一指數時序曲線特征的分析，可以明顯看出冬小麥-夏玉米與其他農作物生長周期不一致，而其他農作物物候周期相似，且指數值峰值無明顯差異，同時分類結果表明，利用單一指數進行分類并不能獲得最佳的分類效果，驗證了多特征組合分類方法在農作物精細化分類上的可行性。

盡管本研究取得了較好的分類效果，但仍存在一些問題。在作物生長關鍵期，單一數據源的影像獲取頻率較低，本研究沒有采用多源遙感來提升遙感數據的獲取概率。且本研究僅考慮了EVI和NDVI等植被指數，其他植被指數沒有納入考慮范圍。因此在今后的工作中可以結合高分系列遙感數據，融合其他植被指數特征開展研究。