基于多時相MODIS數據監測水、旱作物種植面積及空間分布

姜藍齊，王 萍，姜麗霞，宮麗娟，于成龍，李秀芬

（1黑龍江省氣象科學研究所，哈爾濱 150030；2中國氣象局東北地區生態氣象創新開放實驗室，哈爾濱 150030）

0 引言

農田作物信息的快速提取與解析，是精準農業實踐的基礎和前提。其中農作物種植面積遙感提取是作物估產的基本要素，作物空間分布是研究其時空動態變化的基礎，及時準確獲取上述信息不僅在農業生產管理、農業政策規劃等方面具有重要作用，也是理解農業對全球變化影響和歸因的重要內容[1-4]。黑龍江省是中國農業大省，農作物面積和產量作為制定糧食政策和經濟計劃的重要依據，受到政府部門的高度重視。近年來，在經濟和氣候變暖的雙重影響下，黑龍江省水田和旱田作物種植面積及其空間分布發生顯著變化[5-6]，而準確快速地掌握作物種植面積及其時空分布，可為輔助政府有關部門制定合理的糧食保障政策提供基礎依據。

隨著衛星遙感技術的迅速發展，多傳感器、多時間分辨率和多空間分辨率的遙感數據已經廣泛應用于農作物種植遙感提取中，并取得長足進展[7-10]。目前，中高分辨率衛星如Landsat、SPOT、HJ-1/2、GF-1等數據由于其較高的空間分辨率提高農作物分類識別的準確度而占有優勢[11-17]，但其回訪周期長（時間分辨率低）、覆蓋范圍小、成本高的特點又限制了在農業遙感中的應用。MODIS數據具有回訪周期短、更新快、易獲取、覆蓋范圍廣等優勢，在大范圍作物監測、分類提取、估產與預報研究中被國內外學者廣泛應用。如Lobell和Asner[18]、林文鵬等[19]、熊勤學和黃敬峰[20]利用時間序列MODIS NDVI數據分別通過混合像元分解法、Fuzzy ARTMAP神經網絡法和BP神經網絡法進行作物提取研究；閆峰等[21]構建MODIS地表溫度與增強植被指數時間序列數據，采用ISODATA非監督分類方法識別冬小麥；劉珺和田慶久[22]利用最大似然法對MODIS EVI時間序列數據進行分類，基于掩膜處理和決策樹規則提取夏玉米和夏棉花；平躍鵬和臧淑英[23]對MODIS NDVI時序數據進行支持向量機(SVM)分類，結合物候參數和歸一化水體指數(NDWI)提取玉米、水稻和大豆；Xiao[24-25]和Zhang等[26]利用MODIS計算歸一化雪指數(NDSI)、歸一化植被指數(NDVI)、增強植被指數(EVI)等建立掩膜規則提取水稻；Wardlow和Egbert[27]、Zhong等[28]、Wang等[29]分別用MODIS數據建立決策規則對美國大平原、巴西南部巴拉那州、美國密蘇里西部和堪薩斯等地區作物分布進行制圖。綜上可見，MODIS數據在大范圍作物監測、提取研究中發揮了重要作用，主要通過不同方法結合不同植被指數、作物物候等信息實現。但由于分類方法多樣、尺度敏感性等問題，使得目前作物遙感提取研究尚未系統地梳理，難以開展普適性研究和廣泛應用。

從農業遙感監測業務運行的角度出發，需要通過研究區作物發育期特征標準時間曲線，建立適合時間/時相影像的快速、直觀分類方法，以提高業務運行效率，實現大宗作物的準確識別問題。基于此，筆者在前人研究基礎上，選取適用性廣、代表性強的3個指標——歸一化植被指數(NDVI)、增強植被指數(EVI)、地表水指數(LSWI)用于識別水稻和旱田作物，通過構建作物發育期特征標準曲線結合黑龍江省作物積溫條件選取作物關鍵發育期影像，構建決策規則，最終實現黑龍江省水稻、旱田的分類監測提取。本文的特色之處在于引用作物生長發育所需積溫條件確定關鍵發育期影像，以期為解決大范圍作物識別中使用統一決策規則造成的局地識別精度低、差異性大的問題提供了一種可行思路。

1 研究區概況

黑龍江省位于中國東北部地區，經緯度范圍為45°06′—52°36′N、123°42′—134°45′E（圖1）。地貌上西部為松嫩平原、東部為三江平原，地面起伏平緩，土層深厚，利于大型農業機械作業，為農業發展提供了良好的條件；研究區是世界三大黑土帶之一，土壤農業生產能力較高；農作物熟制以一年一熟為主，主要作物有玉米、水稻、大豆，是中國重要的商品糧基地。研究區屬于寒溫帶—溫帶濕潤—半濕潤季風氣候，南北溫差大，1月平均氣溫-31～15℃，7月為16～23℃，無霜期3～4個月，年平均降水量300～700 mm。

圖1 研究區示意圖

2 數據來源與處理

2.1 MODIS數據獲取及處理

本文采用MOD09Q1反射率數據，來源于美國地球數據中心(http://earthdata.nasa.gov/)，行列號為h25v03、h26v03、h26v04、h27v04，空間分辨率為250、500 m，該反射率數據經過8天最大值合成(Maxmum Value Composite,MVC)處理。研究時段為2017年第129天至273天中的17個時相，包括整個作物生長周期。利用MODIS產品批處理工具MRT(Modis Reprojection Tool)軟件對下載影像進行空間拼接、投影轉換及重采樣處理，統一輸出為Albers投影，空間分辨率為250 m的數據。

利用ENVI 5.5.2軟件計算歸一化植被指數(NDVI)、增強植被指數(EVI)和地表水指數(LSWI)，計算如式(1)～(3)所示，其中ρnir為近紅外波段、ρred為紅光波段、ρblue為藍光波段、ρswir為短波紅外波段。

2.2 高分數據獲取及處理

采用的高分一號(GF-1)數據來源于中國資源衛星應用中心，考慮研究區天氣狀況、影像質量以及作物物候特征，下載2017年作物成熟收獲時期覆蓋松嫩平原與三江平原主要種植區GF-1 WFV影像共16景，空間分辨率為16 m。利用ENVI 5.5.2軟件對L1級數據依次進行正射校正、輻射定標、FLAASH大氣校正、幾何校正，對相鄰影像進行融合，完成數據預處理。結合野外實地調查，對2017年GF-1 WFV影像進行目視解譯，獲取覆蓋黑龍江省各積溫帶森林、草地、濕地、水域、水稻、旱田樣本。

2.3 作物發育期和氣象數據

（1）逐日平均氣溫。黑龍江省2017年逐日平均氣溫數據來自于黑龍江省氣候中心，共包括基準氣象站（或基本氣象站）80個，統計無霜期內連續≥10℃平均氣溫累加結果，得到研究區≥10℃活動積溫。采用克里金插值方法，實現研究區≥10℃活動積溫的空間表達，并劃分為6個積溫帶（圖2）。

圖2 黑龍江省積溫帶劃分空間分布

（2）作物發育期數據。來源于黑龍江省氣象科學研究所農業試驗站，包括44個農業氣象站點2013—2017年玉米、大豆、水稻的全發育期時間。按照站點所處積溫帶及作物種類匯總黑龍江省水稻、旱田作物發育期時間（圖3）。

圖3 黑龍江省水稻、旱田作物發育期

2.4 作物樣點及精度驗證數據

（1）土地利用數據。數據來源有2個，一個是2015年中國土地利用現狀遙感監測數據，來源于中國科學院地理科學與資源研究所，該數據基于Landsat 8遙感影像，通過人工目視解譯生成，土地利用類型包括耕地、林地、草地、水域、居民地和未利用土地6個一級類型以及25個二級類型，其中耕地分為水田和旱田（圖4a）。另一個是全球30 m地表覆蓋數據集(GlobelLand30)，來源于國家基礎地理信息中心，該數據集包含10個主要的地表覆蓋類型，分別是耕地、森林、草地、灌木地、濕地、水體、苔原、人造地表、裸地、冰川和永久積雪（圖4b）。

圖4 黑龍江省土地利用空間分布圖

（2）精度驗證數據。野外實地調研數據獲取時間為2017年7月，包括松嫩平原、三江平原作物樣點，其中水稻田88個、旱田147個。農作物播種面積來源于2018年《黑龍江統計年鑒》[30]。

3 方法與結果

3.1 發育期下墊面植被指數特征曲線

水稻和旱地在移栽期的水分和生長期植被指數的變化（NDVI、EVI和LSWI）是水稻區別于旱地和其他植被的重要特征，本研究利用該特征實現水稻和旱地的識別。基于ArcGIS 10.2軟件，提取黑龍江省林地、草地、濕地、水域、水田、旱田樣點發育期（第129～273天中的17個時相）NDVI、EVI和LSWI，計算各類樣點不同時像NDVI、EVI和LSWI平均值，繪制特征曲線，如圖5所示。圖5a為各下墊面NDVI特征曲線，可以看出，在第121～169天林地和水體的NDVI值與草地、濕地、耕地有明顯差異，表現為林地NDVI值顯著高于其他下墊面，而水體NDVI值顯著低于其他下墊面，因此可以通過NDVI值構建公式去除林地和水體。圖5b為濕地、草地、旱地和水稻EVI特征曲線，可知在第121～169天濕地和草地的EVI值大于耕地，第193～233天濕地和草地的EVI值小于耕地，因此可以通過EVI值對濕地/草地和耕地進行區分。圖5c為水稻和旱地LSWI特征曲線，可知水稻發育期LSWI值大于旱地，其中第121～161天水稻與旱地LSWI差異最大，主要表現為水稻LSWI大于0，旱地小于0，可以通過LSWI區分水稻和旱地。

圖5 黑龍江省不同下墊面NDVI指數特征曲線

由于黑龍江省地跨6個積溫帶（圖2），黑龍江省西部松嫩平原處于1、2積溫帶，熱量條件最好；東部三江平原處以3、4積溫帶為主；北部黑河市、伊春市、大興安嶺地區處于5、6積溫帶，熱量條件較差。由于熱量條件的限制，全省三大主要作物（玉米、大豆、水稻）所處積溫帶不同，其作物發育期時間有所差異（圖3），以玉米和大豆為主的旱地在5、6積溫帶播種期晚于1～4積溫帶，而在1、2積溫帶出苗期早于3、4積溫帶；水稻主要種植在1～4積溫帶，播種期在1、2積溫帶為4月中下旬，3、4積溫帶為4月下旬。因此，同一時像不同積溫區植被指數將有所差異，在提取全省水稻、旱田作物分布時，充分考慮積溫條件分區提取，可以通過構建真實作物發育期特征曲線提高分類精度。綜上，在提取耕地與其他，水稻和旱田分類中，按樣點所處積溫帶統計作物發育期EVI（圖6）和LSWI指數特征曲線（圖7），結合積溫條件，選取不同時像和參數依次提取耕地、水稻和旱地。

如在1、2積溫帶（圖6a），第153天和161天草地和濕地EVI指數均大于0.2，明顯高于水稻和旱地，而在第193～241天草地和濕地EVI指數均小于0.5，明顯低于水稻和旱地，因此可通過構建決策樹規則EVI153d＞0.2且EVI209d＜0.5區分第1、2積溫帶濕地/草地和水稻/旱地。在3、4積溫帶（圖6b），由于第169天之前濕地與水稻/旱地EVI指數沒有顯著區分，第193～233天草地/濕地EVI指數均小于0.6，明顯低于水稻/旱地，因此以EVI209d＜0.6為指標區分第3、4積溫帶濕地/草地和水稻/旱地。在5、6積溫帶（圖6c），濕地/草地EVI指數在第129～169天均大于0.2，顯著高于水稻/旱地，因此以EVI153d＞0.2為指標區分濕地/草地與水稻/旱地。

圖6 不同積溫帶濕地/草地和耕地EVI特征曲線

在1、2積溫帶（圖7a），第129～169天水稻與旱地LSWI有顯著差異，表現為水稻LSWI指數小于0.1，旱地LSWI指數大于0.1，因此構建指標LSWI145d＞0.1判別第1、2積溫帶水稻和旱地。在3、4積溫帶（圖7b），水稻和旱地LSWI指數在第129～169天差異表現為是否大于0，因此以LSWI145d＞0判別第3、4積溫帶水稻和旱地。在5、6積溫帶，第153～169天水稻LSWI指數大于0而旱地LSWI指數小于0，因此選取LSWI153d＞0判別第5、6積溫帶水稻和旱地。

圖7 不同積溫帶水稻和旱地LSWI特征曲線

3.2 決策樹規則構建與提取結果

綜上，通過下墊面植被指數特征曲線的分析，結合作物關鍵發育期和積溫條件，選取相應時像，構建決策樹規則，在ENVI 5.5.2軟件中實現。具體操作步驟為ENVI Classic-Classification-Decision Tree-Build New Decision，流程詳見圖8，提取結果如圖9所示。

圖8 分類決策樹

圖9 基于MODIS數影像的分類結果

可以看出耕地集中分布在黑龍江省西部松嫩平原與東部三江平原，整體來看旱地種植面積大于水稻，松嫩平原以旱作農業為主，三江平原水稻種植面積較為廣泛，分類結果空間分布與黑龍江省土地利用圖（圖4）大體一致。基于ArcGIS 10.2統計水稻、旱地面積分別為31439.70、115452.39 km2。具體在空間分布上，中北部地區黑河市、伊春市，南部牡丹江市以旱地為主，本研究分類結果空間上少于土地利用分類結果；東部三江平原，本研究分類結果水稻面積整體略高于土地利用分類結果。

3.3 水田、旱田分類結果精度驗證

基于ENVI 5.5.2，以2017年野外實地調查獲取的地面驗證樣點對分類結果進行混淆矩陣精度驗證。具體操作流程為：ENVI Classic-Classification-Post Classification-Confusion Matrix-Using Ground Truth ROIs。表1為混淆矩陣驗證結果，分類結果的總體精度為90.68%，Kappa系數為0.81，其中旱地制圖精度(98.46%)高于水稻制圖精度；水稻的用戶精度(97.73%)高于旱地，分類結果精度均較高。

表1 分類結果混淆矩陣精度評價表

4 討論

在基于遙感影像識別作物類型的研究中，判識指標的選取和分類結果精度驗證的準確性是決定分類結果準確與否的2個關鍵點。這也是本研究過程中的難點問題，值得深入討論。

（1）在已有研究中，水稻、旱田判識指標多選取關鍵期影像，基于不同植被指數構建決策規則，該判識指標可以實現作物分類的快速提取，但在大范圍作物提取研究中會存在局地作物識別精度低、不夠準確的問題。基于此本研究加入了積溫條件來界定局地關鍵期影像的選取，以提高大范圍作物局地空間提取精度。為研究積溫條件的引入對分類結果精度提高的作用，筆者討論了不考慮積溫條件，構建決策樹規則提取的水稻和旱田與本研究結果的差異，進一步明確積溫條件對提取水稻、旱地分類結果精度提高的重要性。

積溫條件引入對分類結果精度的影響。本文重建不考慮積溫條件的決策樹規則提取水稻和旱田，與本研究結果進行比較分析。具體流程為：以圖4黑龍江省不同下墊面植被特征曲線為標準建立決策樹規則，通過0.2＜NDVI129d＜0.1去除林地和水域，再通過EVI153d＞0.2且EVI209d＜0.6區分濕地/草地與水稻/旱地，最后通過LSWI145d＞0區分水稻和旱地，提取結果如圖10a所示。與本研究考慮積溫條件提取結果相比（圖10b），可以看出耕地空間分布趨勢與實際情況大體吻合，均分布在西部松嫩平原與東部三江平原地區，其中旱地面積均大于水稻種植面積，但整體來看耕地面積明顯大于實際情況，如黑龍江北部黑河市、大興安嶺地區；且東部、西部地區耕地成集中連片趨勢，空間精細化程度不夠。

圖10 水稻、旱地提取結果

以2017年野外實地調查獲取的地面驗證樣點對分類結果進行混淆矩陣精度驗證，結果表明總體精度為77.91%，Kappa系數為0.60，其中水稻制圖精度為58.56%，用戶精度為90.28%；旱地制圖精度為92.52%，用戶精度為90.07%。與充分考慮積溫條件相比，總體分類結果精度減少了12.77%，水稻制圖精度減少了22.57%，旱地制圖精度減少了5.94%。但在不考慮積溫條件下，提取的水稻和旱地用戶精度偏高，這可能是提取結果集中連片、空間精細化程度不夠，覆蓋野外實地樣點所致。

（2）目前已有研究在作物提取精度驗證中多采用混淆矩陣進行精度評價，缺少作物種植面積的精度驗證，尤其是缺少實際作物面積數據的支持。基于此，本研究進行了分類結果面積的精度驗證，比較分析作物實際播種面積、基于遙感影像提取的土地利用類型分類面積與本文分類結果的差異。

作物提取面積一致性評定。本研究以《黑龍江省統計年鑒》[30]記載的播種面積作為實際作物面積；將基于Landsat 8提取的土地利用類型圖作為遙感影像提取作物面積，與本研究提取結果進行對比。基于ArcGIS 10.2，按照所轄市統計面積結果，如表2所示。整體來看，基于Landsat 8遙感影像提取的水稻和旱地面積高于本研究提取結果和統計年鑒記載資料，其中水稻面積為32034.19 km2，與本研究結果相對誤差為1.85%，與統計年鑒結果相對誤差為22.04%；旱地面積為133945.79 km2，與本研究結果相對誤差為13.80%，與統計年鑒結果相對誤差為13.68%。《黑龍江省統計年鑒》中記載的水稻、旱地播種面積低于基于Landsat 8遙感影像與本研究提取結果，其中水稻播種總計面積溫24972.88 km2，與本研究結果相對誤差為25.89%，與Landsat 8提取結果相對誤差為28.27%；旱地播種總面積為115614.90 km2，本研究結果相對誤差為0.14%，與Landsat 8提取結果相對誤差為15.85%。可見，本研究提取的水稻種植總面積與基于遙感影像提取的結果一致性較好，而提取的旱地種植總面積與基于統計年鑒記載的結果一致性較好。

表2 黑龍江省水稻、旱田面積一致性評定

圖11 為黑龍江省各轄區水稻、旱地面積趨勢對比，可以更清晰直觀地表達各轄區水稻和旱地作物面積差異。圖11a整體看來，統計年鑒、遙感影像和本研究結果各市轄區旱地面積趨勢大致相同，其中統計年鑒記載面積整體小于遙感影像和本研究提取結果，僅大興安嶺地區統計年鑒記載結果大于兩者；各市轄區來看，齊齊哈爾和綏化地區遙感影像與本研究提取面積一致性較好，均大于統計年鑒記載面積；雞西、牡丹江、黑河地區本研究結果與統計年鑒結果一致性較好，均小于遙感影像提取面積；哈爾濱、鶴崗、雙鴨山、佳木斯地區旱地面積提取結果均表現為遙感影像大于本研究結果大于統計年鑒記載。圖11b整體來看，水稻種植面積本文提取結果與遙感影像一致性較好，統計年鑒播種面積在雞西、鶴崗、雙鴨山、佳木斯地區明顯小于本文和遙感影像提取面積，可能是統計年鑒各轄區記載面積沒有包括農墾總局導致的。

圖11 黑龍江省各轄區水稻、旱地面積趨勢對比

綜上，目前在作物提取面積精度驗證中仍然存在很大的不確定，雖然提取結果與實際記載播種面積有很強的相關性，但現有研究基于遙感影像提取的分類結果面積普遍大于實際作物播種面積；且基于不同方法和不同分辨率的遙感影像提取的面積結果也有明顯差異。同時由于MODIS數據受分辨率的影像，在作物面積提取中如何考慮混合像元分解問題，提高面積提取精度，都是在以后的研究中需要進一步解決的問題。

（3）在本研究中還存在一些不確定性需要進一步討論和解決。如雖然筆者以實地調查樣點的驗證結果分類精度較高，但仍需要進一步討論基于本研究方法提取結果的準確性；以GF-1號數據為基準提取了森林、草地、濕地、水域、旱地和水稻樣點用于植被指數特征曲線的分析，可能會在制圖、目視解譯過程中產生一定誤差。

5 結論

使用多時像高時間分辨率MODIS09Q1數據，結合作物物候特征，在前人研究基礎上引入積溫條件，建立參考NDVI、EVI和LSWI時序植被指數曲線，構建決策樹規則提取水稻和旱地。研究結果表明：

（1）該方法可以實現全省范圍水稻和旱田的快速有效提取，識別作物的總體精度為90.68%；水稻、旱地制圖精度分別為81.13%、98.46%，用戶精度分別為97.73%、87.07%。

（2）分類結果顯示黑龍江省耕地集中分布在西部松嫩平原與東部三江平原，其中旱地種植面積大于水稻，松嫩平原以旱作農業為主，三江平原水稻種植面積較為廣泛，分類結果空間分布與黑龍江省土地利用分類結果空間一致性強。

（3）積溫條件的引入有效提高了水稻、旱田的識別精度，其中總體提取精度提高12.77%，水稻制圖精度提高22.57%，旱地制圖精度提高5.94%。