田 苗 單 捷 盧必慧 黃曉軍
(江蘇省農業科學院農業信息研究所, 南京 210014)
及時準確獲取水稻種植面積及空間分布,對水稻生產管理、農業政策制定及糧食安全分析等具有重要意義。目前水稻面積的遙感提取方法主要有常規分類法和基于MODIS的時間序列方法。常規分類法一般使用中、高分辨率的遙感數據[1-3],而此種數據重訪周期長,易受天氣影響,不易抓住水稻生長的物候信息,因而,難以進行大范圍的水稻面積業務化運行。MODIS數據具有較高的時間分辨率和較廣闊的空間覆蓋,而且可以免費獲取,所以適合大區域水稻面積提取和推廣應用[4-5]。基于MODIS時間序列的方法往往根據水稻的物候歷,確定水稻識別的移栽期、生長期和收獲期等關鍵時期,根據水稻在不同時期表現的遙感特征指數差異來識別水稻[6-11]。較少有學者考慮與水稻生產能力密切相關的季節積分和生長季振幅兩個物候指標。經過眾多學者對植被物候的研究,開發了許多判斷植被關鍵生育期的方法,主要包括植被指數閾值法[12]、植被指數中點法[13]、最大變化斜率法[14]、植被指數滑動移動平均法[15]、通過NDVI時間序列累積生長天數來計算物候信息[16]、植被物候期的頻率分布與遙感數據相結合法[17]等。目前圖像處理已經進入人工智能時代,機器學習回歸算法的高度適用性和通用性使其比傳統的擬合函數能更精確地估計物候趨勢。
本研究以江蘇省為研究區域,應用機器學習算法對MODIS-EVI時間序列進行重構,并提取水稻物候信息,選取季節積分和生長季振幅兩個指標,采用閾值法提取水稻種植面積,并應用統計數據和Landsat8 OLI影像對提取結果進行精度驗證。
江蘇省地處亞熱帶向暖溫帶過渡的區域,其自然條件對水稻生長有利,是全國水稻主產省之一。江蘇省南北縱跨5個緯度(30°45′~35°20′N),東西橫跨5個經度(116°18′~121°57′E),其中,北部地處我國秈稻種植最北邊緣,南部地處我國兩熟制粳稻種植最南邊緣,東部受海洋氣候影響,西部受陸地氣候影響,形成了境內多種稻作生態區、多種類型稻作熟制和稻作方式、秈粳交錯和多種熟制品質并存的水稻生產格局,在全國稻作生產中具有特殊意義。
本文采用2019—2020年地表反射率數據(MOD09A1)作為基礎數據,數據來源于美國國家航空航天局網站(https:∥modis.gsfc.nasa.gov/),MOD09A1數據的空間分辨率為500 m,時間分辨率為8 d。使用MODIS重投影工具(MODIS Re-projection tool, MRT)對數據進行重采樣、格式轉換、拼接和重投影等預處理。采用中高空間分辨率的Landsat8 OLI遙感影像數據進行水稻信息目視解譯,選取水稻建模樣點及精度評估樣點,建立水稻面積提取模型并對提取結果進行精度驗證。統計數據來自2019年和2020年江蘇省農村統計年鑒,包括2019年和2020年江蘇省各地級市水稻種植面積。
江蘇省所在的MODIS影像序列號為h27v05和h28v05。文獻[18]研究證明覆蓋作物關鍵生育期的時間序列長度即可滿足一定分類精度,因此,本研究選擇的影像時期為2019年和2020年6月2日—11月9日(DOY153~DOY313),包含水稻移栽期至成熟期。對得到的圖像進行波段運算得到EVI時間序列數據集。EVI計算公式為
(1)
式中ρnir——近紅外通道反射率
ρred——紅色通道反射率
ρblue——藍色通道反射率
L——土壤調節參數,取1
基于遙感手段監測植被物候涉及兩個重要步驟,一是植被指數的重構,二是植被物候參數提取方法的選擇。目前,圖像處理已進入了人工智能時代,機器學習算法已成為時間序列數據處理的有力工具[19-21]。機器學習回歸算法(Machine learning regression algorithm, MLRAs)的高度適應性和通用性使其比傳統的擬合函數能更精確地估計物候趨勢[22]。高斯過程回歸(Gaussian process regression,GPR)[23]作為一種最優機器學習回歸算法在成功重建植被指數和提取可靠的物候指標方面具有潛力。本研究應用GPR算法對時間序列進行重構,在此基礎上提取物候信息,提取方法為:生長季開始日期/生長季結束日期在曲線的左/右部分達到曲線上升/衰減部分的季節振幅的30%時被識別。物候指標包括:生長季振幅、最大值、最大值日期、生長季開始日期(Start of season, SOS)、生長季結束日期(End of season, EOS)、季節積分、生長季長度。
在生長發育期識別的各個指標中,季節積分和生長季振幅是兩個非常重要的指標,不僅顯示了水稻的生育特點,還能夠反映水稻的生長狀況[24]。經過反復實驗,本研究選用這兩個指標監測江蘇省水稻種植面積。
本研究采用閾值法開展水稻面積提取。閾值法是圖像分割中的經典方法,它利用圖像中要提取的目標與背景在灰度上的差異,通過設置閾值來把像素級分成若干類,從而實現目標的提取。最優閾值選擇方法通常有人工經驗選擇法、直方圖法和最大類間方差法[25]。本研究閾值確定方法選用直方圖法,結合Landsat8選取的水稻樣點和統計面積數據,分別確定江蘇省13個地級市季節積分和生長季振幅的閾值。確定步驟為:
(1)以高分辨率遙感影像(Landsat8 OLI)為基礎,以目視解譯的方式在每個地級市選擇若干個水稻樣點(圖1),并提取各個水稻樣點對應的MODIS-EVI時間序列。

圖1 2019年水稻樣點分布圖Fig.1 Distribution map of rice sampling points in 2019
(2)應用MODIS-EVI時間序列提取每個水稻樣點的物候信息,得到每個水稻樣點的季節積分和生長季振幅指標。
(3)以地級市作為分區,應用每個區域的若干個樣點,做該區域水稻的季節積分頻數分布圖和生長季振幅頻數分布圖。
(4)根據頻數分布圖,初步確定季節積分和生長季振幅的閾值,并將該值用于提取該區域水稻種植面積。
(5)將步驟(4)提取的水稻種植面積與該區域水稻種植面積的統計數據對比分析,根據頻數分布重新調整兩個指標的閾值,使水稻種植面積的總體誤差控制在±10%以內,并將該閾值作為該區域的水稻提取閾值。
(6)應用2019年的數據,經過步驟(1)~(5),確定各個地級市水稻提取的閾值,應用該閾值提取2020年的水稻種植面積,并對2020年的提取結果進行精度評價。
水稻面積提取精度評價以兩種方式進行:與統計面積數據進行對比分析;通過高分辨率影像選取水稻和非水稻樣點,建立混合矩陣對提取結果進行定量分析。
本研究以2019年的數據為建模數據,2020年的數據為精度驗證數據,采用2019年中、高空間分辨率的Landsat8 OLI遙感影像數據進行水稻信息目視解譯,在江蘇省內選取了432個水稻樣點(圖1),應用GPR方法將2019年水稻樣點時間序列(6月2日—11月9日)進行重構,重構后的時間序列時間間隔為1 d(圖2a)。圖2a顯示GPR模型可以很好地將時間序列中的異常值和缺失值進行重構,使曲線更加平滑合理,在此基礎上提取水稻物候信息,物候提取參數生長季振幅設置為30%(圖2b)。

圖2 單點時間序列重構及水稻生育期提取Fig.2 Single-point time series reconstruction and extraction of rice growth period
圖2b顯示了水稻物候提取參數,其中橫坐標為日期,第1天為時間序列開始日期,即6月2日,通過與江蘇省水稻物候數據對比,該模型提取的SOS和EOS分別對應江蘇水稻的分蘗期和成熟期。按以上方法提取所有水稻樣點的物候信息,并分區域整理所有樣點的季節積分和生長季振幅。
在水稻樣點物候信息提取的基礎上,逐像素提取2019年和2020年水稻物候信息,得到2019年和2020年生長季振幅和季節積分分布(圖3),從圖3中可以看出,生長季振幅和季節積分的分布具有很好的區域性,水體和城市所在地區的值都較低,能夠很好的與周圍地物分開,植被覆蓋區顯示也較清晰,有利于水稻物候信息提取。

圖3 江蘇省2019年和2020年生長季振幅和季節積分Fig.3 Amplitude and area total of Jiangsu Province in 2019 and 2020
分析2019年432個水稻樣點的生長季振幅和季節積分,通過直方圖和水稻面積的統計數據確定13個市兩個指標的閾值(圖4),圖4顯示,生長季振幅和季節積分閾值的分布也具有區域相似性,如南京市、蘇州市、無錫市、鎮江市、常州市,兩個指標的閾值均相似,淮安市和宿遷市也具有相似性,連云港市、南通市和鹽城市具有相似性等。

圖4 江蘇省13個市的生長季振幅和季節積分閾值Fig.4 Amplitude and area total thresholds of 13 cities in Jiangsu Province
將閾值相似的區域合并(圖5),圖5充分體現了江蘇省水稻種植的特點:東部受海洋氣候的影響(連云港市、鹽城市和南通市),其水稻種植具有一定的相似性,西部受陸地氣候的影響,形成了多種稻作生態區,南北縱跨5個緯度,南北水稻種植方式和格局有較大差異。徐州市作為江蘇省最北端的地級市與蘇南5市的差異較大,蘇中水稻種植情況也具有一定差異性,徐州市、連云港市和鹽城市等蘇北地區也具有一定的相似性(圖4)。因此,圖4和圖5從側面反映了模型建模的合理性。

圖5 閾值合并后的分區Fig.5 Partition map after threshold merging
根據2019年數據確定13個市的生長季振幅和季節積分的閾值,分別提取2020年13個市的水稻種植面積及空間分布,合并得到江蘇省2020年水稻種植面積及空間分布(圖6)。從圖中可以看出,水稻的分布具有很好的區域性,能夠較好地將水體、城市和山區等剔除,符合江蘇省水稻種植的分布情況。
首先,應用2020年江蘇省13個地級市水稻種植面積的統計數據與MODIS提取數據進行對比分析(表1)。表1顯示2020年提取的全省水稻總面積的相對誤差為-6.10%。在13個市的面積統計中,連云港、淮安、宿遷、徐州、鹽城、揚州、泰州、南通和蘇州9個市的相對誤差都在±10%以內,南京、鎮江、常州、無錫4個市的相對誤差較大。經分析,可能原因有:①在建模過程中,蘇南4市選取的樣點較少,閾值的確定不夠準確。②蘇南4市水稻種植面積較少,在絕對誤差相同的情況下,相對誤差較高。③蘇南地形復雜,水稻種植情況不夠穩定,每年的種植情況變化較大(從統計面積上看,2020年蘇南4市的水稻種植面積較2019年減少了近10%),且水稻種植的空間分布情況變化也較大,因此,應用2019年的數據建立的模型,不適用于2020年的水稻面積提取。④在本研究中使用的MODIS影像空間分辨率為500 m,即一個像元內存在像元“同物異譜”和“同譜異物”的現象,對面積較小、地塊破碎的地區提取精度不夠,從而影響水稻的提取精度。除蘇南4市外,江蘇省大部分區域水稻種植面積隨時間變化具有一定的穩定性,可以用該模型進行水稻面積的提取。在后續研究過程中,還可以不斷地加入歷史數據或最新數據,不斷更新調整閾值,從而更準確地把握每個區域的水稻種植特點,使參數具有穩定性、代表性和可推廣性。

表1 2019年和2020年江蘇省水稻面積提取結果與農業統計數據對比Tab.1 Comparison of rice area extraction results and agricultural statistics in Jiangsu Province in 2019 and 2020
其次,僅僅依靠統計面積的比較進行精度驗證,在理論上存在一定程度的不可靠性,因為相同的面積可能會有不同的形式組合,用高分辨率的遙感影像驗證大尺度(低分辨率)提取精度,這種方法的優點在于省時省力效率高,彌補缺乏野外驗證的不足。本研究利用Landsat8和目視解譯方法選取各個地類的樣點共3 423個,通過混淆矩陣對提取結果進行定量分析(表2),分析結果表明,水稻的制圖精度為92.90%,用戶精度為89.09%。計算得水稻提取的總體精度為92.55%,Kappa系數為0.846 3。

表2 混淆矩陣精度評價Tab.2 Confusion matrix accuracy evaluation
將非水稻樣點分為水體、建筑、其他植被及其他4類,分析各個類別的提取精度(表3)。其中水體、建筑及其他類型由于與水稻差異較大,精度分別為99.52%、98.11%和92.33%,而其他植被由于與水稻具有相似性,因此錯分為水稻的概率較大,其精度為87.73%。

表3 非水稻類型精度評價Tab.3 Accuracy evaluation for non-rice types
在精度評價方面,盡管已經開展了大量工作,但是對這種光譜指數時間序列水稻面積提取方法的精度評價還不夠充分,本研究主要是針對統計年鑒統計數據進行的行政區尺度水稻面積評估,雖然本研究也應用高空間分辨率影像,提取樣點進行了精度驗證,但仍不能全面完整地體現面積提取方法的空間分布精度,對于水稻面積總量漏分和誤分的空間分布的精度評價工作還有待進一步深入研究。
(1)基于機器學習的GPR模型能夠較好地重構MODIS-EVI時間序列數據集的異常點和缺失值,為后續水稻物候期的提取提供了較好的基礎。根據提取的水稻物候信息,提取水稻種植面積,精度較好,與統計數據相比,大部分地級市的水稻提取精度都在90%以上。與高分辨率影像數據相比,水稻提取的總體精度為92.55%,Kappa系數為0.846 3,水稻的制圖精度為92.90%,水稻的用戶精度為89.09%。由于其他植被與水稻生長具有一定的相似性,其他植被類型的錯分率較高,其精度為87.73%,水體、建筑和其他類型的錯分率都較低,精度都在90%以上。
(2)江蘇省水稻種植方式、種植格局有較大差異,本研究分13個地級市確定水稻物候指標的閾值,通過閾值相似性合并,很好地體現了江蘇省水稻種植的特點。該結論一方面說明了本研究的合理性,另一方面說明對大區域水稻種植面積的分析與建模,需要分區域進行。
(3)本研究僅應用2019年的數據作為建模數據,將模型用于2020年水稻面積提取,并取得了較好的結果。