基于長時間序列MODIS影像的黑龍江省主要農(nóng)作物種植結(jié)構(gòu)變化研究

李 暉，陳雨琪，張紫滟，武洪峰，劉占宇

(1. 杭州師范大學(xué)遙感與地球科學(xué)研究院，浙江 杭州 311121； 2. 黑龍江省農(nóng)墾科學(xué)院情報研究所，黑龍江 哈爾濱 154005)

0 引言

隨著人口的快速增長和城市化趨勢的不斷增強(qiáng)，世界的農(nóng)作物分布格局都在發(fā)生大小不一的變化.農(nóng)作物的種植結(jié)構(gòu)信息、空間分布格局和農(nóng)業(yè)精細(xì)化管理對于保障國家糧食安全、應(yīng)對全球農(nóng)情變化以及農(nóng)業(yè)可持續(xù)發(fā)展等具有重要意義.農(nóng)作物種植結(jié)構(gòu)是空間分布格局的重要組成部分，它用來描述一塊區(qū)域或生產(chǎn)組織內(nèi)農(nóng)作物的組成和布局，即主要農(nóng)作物類型(種什么)和空間分布(種在哪里)[1-3].快速、準(zhǔn)確地獲取耕作區(qū)域農(nóng)作物的空間分布信息和時空動態(tài)變化信息不僅是監(jiān)測農(nóng)作物種植面積、預(yù)測區(qū)域農(nóng)作物產(chǎn)量、保障區(qū)域糧食安全平衡的主要數(shù)據(jù)來源，也是進(jìn)行農(nóng)作物種植結(jié)構(gòu)布局優(yōu)化的重要依據(jù)[4].同時，作為可持續(xù)化發(fā)展的重要形式，農(nóng)業(yè)土地的開墾利用、作物種植更替以及與其他地物的相互轉(zhuǎn)換，必然會對整個生態(tài)環(huán)境包括大氣圈、水圈和生物圈產(chǎn)生不同程度的影響[5].遙感技術(shù)具有覆蓋范圍大、重訪周期短、現(xiàn)勢性強(qiáng)和效率高等優(yōu)點，比傳統(tǒng)的人工調(diào)查更能準(zhǔn)確直觀地了解農(nóng)作物的分布特點，因此早已廣泛應(yīng)用于大尺度范圍的作物空間分布信息提取中[6].20世紀(jì)70年代，美國、日本、澳大利亞等國家逐漸開始對水稻、大小麥、大豆等農(nóng)作物空間分布進(jìn)行遙感監(jiān)測[7-8]，20世紀(jì)90年代，我國也逐步建立了專業(yè)的農(nóng)作物遙感監(jiān)測平臺.隨著MODIS(傳感器)到Landsat(平臺)、TM(傳感器)、SPOT(平臺)，再到國產(chǎn)高分系列衛(wèi)星的普及應(yīng)用，時間和空間分辨率不斷提升，為農(nóng)作物種植信息的精準(zhǔn)提取帶來了便利.

在農(nóng)作物的遙感識別提取方面，眾多學(xué)者進(jìn)行了深入研究.在方法上可以分為3類：1)基于單一影像識別法.此方法適用于作物種植結(jié)構(gòu)比較簡單的區(qū)域，通過單幅影像或兩到三幅影像對作物的種植信息進(jìn)行提取，其重點是找到待分作物的光譜信息特征與其他地物的區(qū)別.鄭長春等[9]基于SPOT-5數(shù)據(jù)，利用簡單決策樹結(jié)合歸一化植被指數(shù)NDVI、光譜波段和背景差異等特征量，對黑龍江省852農(nóng)場主要作物種植信息進(jìn)行提取.Mathur等[10-11]利用IRS-1D數(shù)據(jù)，采用支持向量機(jī)(SVM)的方法結(jié)合光譜波段分析，對印度旁遮普地區(qū)的棉花和水稻種植信息進(jìn)行提取.2)基于時間序列影像識別法.農(nóng)業(yè)土地系統(tǒng)是多種農(nóng)作物通過一種特定的規(guī)律種植并組合形成的空間表達(dá)，由于不同作物生長發(fā)育時間和階段不同，單一遙感影像難以覆蓋所有作物生長發(fā)育期和光譜識別的最佳時相[3].因此，利用時間序列遙感影像表現(xiàn)物候變化規(guī)律從而進(jìn)行農(nóng)作物遙感提取成為目前最為主流的技術(shù)方法[10-12].黃青等[13]基于MODIS-NDVI時間序列數(shù)據(jù)集，通過分析東北三省主要作物的物候歷及時序變化特征，構(gòu)建作物遙感提取模型并通過調(diào)整閾值的方法對東北三省主要作物種植信息進(jìn)行提取.Zhang等[14]根據(jù)時序MODIS影像，結(jié)合地表溫度產(chǎn)品LST(land surface temperature)，對中國和印度2000—2015年水稻種植面積進(jìn)行提取，并對兩國16年水稻面積的變化趨勢進(jìn)行分析.Atzberger等[15]基于AVHRR影像，通過神經(jīng)網(wǎng)絡(luò)算法自動計算NDVI時間序列曲線和端元豐度之間的非線性關(guān)系，對意大利托斯卡尼區(qū)域進(jìn)行作物種植信息自動提取.3)基于遙感影像數(shù)據(jù)與統(tǒng)計數(shù)據(jù)結(jié)合識別法.遙感技術(shù)為有效提取作物信息發(fā)揮了很大作用，但是由于遙感數(shù)據(jù)的獲取方式、大氣干擾、混合像元、分辨率等因素的影響，只采用遙感技術(shù)對大尺度區(qū)域進(jìn)行作物提取仍然有一定的缺陷[16-17].Leff等[18]采集已有全球范圍的遙感耕地數(shù)據(jù)，結(jié)合國家、省級和縣等多個行政單位的作物統(tǒng)計資料，通過非耕地區(qū)域掩膜、亞像素面積計算等方法，對20世紀(jì)90年代全球5分柵格尺度下的18種農(nóng)作物進(jìn)行空間分布信息提取.Fischer等[19]利用GAEZ模型，根據(jù)全球耕地分布圖、人口密度圖等信息，采用交叉信息熵理論和方法將作物統(tǒng)計信息分成5分柵格尺度像元，得到全球農(nóng)業(yè)生態(tài)區(qū)23類農(nóng)作物的空間分布特征.

針對目前農(nóng)作物空間分布信息遙感提取的研究現(xiàn)狀，對于大面積耕作地區(qū)多年份、多作物的研究較少，本文以整個黑龍江省為研究區(qū)，農(nóng)作物種植面積空間分布的快速提取為主要目標(biāo)，根據(jù)不同作物的波譜特性和物候特征，選取 MODIS-NDVI時間序列數(shù)據(jù)集的合適時相及特征參量，對黑龍江省主要農(nóng)作物進(jìn)行識別提取，并對3個年份的農(nóng)作物空間分布圖進(jìn)行分析，為研究區(qū)農(nóng)業(yè)遙感監(jiān)測和作物種植結(jié)構(gòu)調(diào)整等決策信息提供有效可行的參考和借鑒.

1 研究區(qū)概況

黑龍江省位于我國東北部，總面積約47.3萬km2，東北和西南部(三江平原、松嫩平原等地區(qū))地勢較低(圖1a)，有利于農(nóng)作物的種植.以大陸性季風(fēng)氣候為主，春冬季低溫少雨，夏季溫?zé)岫嘤?圖1b)，≥10 ℃有效積溫為1 600～3 300 ℃(圖1c).第一積溫帶分布在西南部地區(qū)，熱量充足，降水充沛，適合種植水稻等作物，而第六積溫帶緯度較高，水熱條件不足，作物呈零散種植分布.

黑龍江省作為我國重要的商品糧基地，農(nóng)業(yè)種植條件優(yōu)越，全省耕地近160 000 km2，居全國第一位.無霜期140～180 d，主要農(nóng)作物類型為水稻、玉米、大豆和大小麥等.2016年黑龍江省利用118 040 km2耕地面積生產(chǎn)了6. 057 8×107t糧食[20]，為全國提供了1/4的商品糧.

(a)數(shù)字高程圖

(b)月平均最高、最低氣溫和降水量

(c)積溫分布圖

2 數(shù)據(jù)與方法

2.1 數(shù)據(jù)來源

從美國地質(zhì)勘探局?jǐn)?shù)據(jù)中心(USGS EROS，https://ladsweb.modaps.eosdis.nasa.gov/)獲取的2001、2010和2017年46景8天合成MODIS地表反射率產(chǎn)品(MOD09Q1)，空間分辨率為250 m.利用MODIS產(chǎn)品批處理工具M(jìn)RT(MODIS reprojection tool)，將投影統(tǒng)一轉(zhuǎn)換成Albers等積投影，地理坐標(biāo)系為WGS-84.以研究區(qū)范圍為裁剪框?qū)D(zhuǎn)換后的數(shù)據(jù)進(jìn)行裁剪，并利用紅光波段數(shù)據(jù)和近紅外波段反射率數(shù)據(jù)計算歸一化植被指數(shù)NDVI，獲得3個年份的研究區(qū)MODIS-NDVI時間序列數(shù)據(jù)集.

輔助數(shù)據(jù)為黑龍江省2001、2010和2017年土地覆蓋類型產(chǎn)品MCD12Q1，精度驗證數(shù)據(jù)來源于黑龍江省歷年農(nóng)業(yè)統(tǒng)計年鑒.

2.2 研究方法

2.2.1 NDVI時間序列數(shù)據(jù)集

歸一化植被指數(shù)NDVI是目前應(yīng)用最廣泛的植被指數(shù)[21].完整的NDVI時間序列數(shù)據(jù)集能夠較好地反映農(nóng)作物從播種—生長—成熟—收割的全過程.由于農(nóng)作物的種類、生長特性和播種收割時間不同，表現(xiàn)出的光譜差異也不同.利用ENVI 5.2處理得到NDVI時間序列數(shù)據(jù)集.計算公式為

(1)

式中：ρnir為近紅外波段的反射率，ρred為紅波段的反射率.

2.2.2 MODIS-NDVI時間序列數(shù)據(jù)集重構(gòu)

由于云、氣溶膠以及地物的雙向性反射等因素會對地表覆蓋及陸地生態(tài)系統(tǒng)的監(jiān)測產(chǎn)生極大的影響[3]，因此需要對植被指數(shù)時間序列集進(jìn)行平滑處理，即重構(gòu)時序植被指數(shù).

目前常用的植被指數(shù)時間序列濾波去除噪聲分析方法主要有Savizky-Golay濾波(S-G)、非對稱高斯濾波(A-G)、雙邏輯調(diào)和函數(shù)濾波(D-L)和傅里葉諧波(Hants)等.J?nsson等[22]的研究表明Hants濾波在重現(xiàn)植被物候特征時具有較好的擬合度,平滑時間序列數(shù)據(jù)集的效果較好.其計算公式為

(2)

式中：諧波的余項A0等于序列的平均值；Aj為各諧波的振幅；ω=2jπ/N為各諧波的頻率；N為序列的長度；θj為各諧波的初相位；m為諧波的個數(shù)，m=N-1，這些正弦函數(shù)疊加構(gòu)成傅里葉序列.

2.2.3 農(nóng)作物遙感提取模型

在作物的整個生長周期內(nèi)，不同作物的光譜曲線特征都有一定程度的差異性，根據(jù)不同類型作物獨特的物候歷特征，抓住不同作物在不同時間段生長周期不同的物候差異，可以實現(xiàn)作物種植類型提取和空間分布提取[23].

根據(jù)對黑龍江省主要農(nóng)作物(玉米、水稻、大豆和大小麥)的物候歷(表1)及植被指數(shù)與物候期對應(yīng)關(guān)系的分析，利用MODIS-NDVI時序數(shù)據(jù)獲取不同農(nóng)作物在種植空間信息上的關(guān)鍵時間點，構(gòu)建研究區(qū)農(nóng)作物遙感提取模型(表2)，其中T1—T13為不同農(nóng)作物在不同生長發(fā)育期內(nèi)的NDVI閾值.

表1 黑龍江省主要農(nóng)作物(玉米、水稻、大豆和大小麥)物候特征Tab.1 Phenological characteristics of main crops (corn, rice, soybean and wheat) in Heilongjiang Province

表2 黑龍江省玉米、水稻、大豆和大小麥遙感提取模型Tab.2 Remote sensing extraction model for corn, rice, soybean and wheat in Heilongjiang Province

以大豆為例進(jìn)行說明.黑龍江省大豆種植時間大約在5月中上旬，5—6月出苗、開花，7月中下旬結(jié)莢，8—9月從鼓粒到成熟.隨著大豆的出苗、開花、結(jié)莢、成熟和收割，反映在遙感影像上的NDVI值隨著時間的變化先逐步增大，達(dá)到一個峰值后再逐漸減小.根據(jù)這一特性，大豆種植面積信息的識別要同時滿足表2中序號為3的所有條件，為了保證結(jié)果的準(zhǔn)確性，避免受到因各地種植習(xí)慣不同、氣候干擾、成熟和收割時間差異帶來的影響，模型中某一NDVI值的選擇包含了前后兩個時相的“或”.設(shè)置這些閾值參數(shù)并帶入上述模型后，能對農(nóng)作物的空間分布進(jìn)行提取.

3 結(jié)果與分析

3.1 作物種植結(jié)構(gòu)空間分布研究

根據(jù)表2作物空間分布遙感提取模型，得到黑龍江省主要農(nóng)作物2001、2010和2017年空間分布圖(圖2).

2001年

2010年

2017年

圖2 黑龍江省2001、2010和2017年主要農(nóng)作物種植結(jié)構(gòu)空間分布圖Fig.2 Spatial distribution map of main crop planting structures in Heilongjiang Province in 2001, 2010 and 2017

由圖2可見，自2001年開始，除了大小麥的種植面積不斷減少外，其他3種作物種植面積不斷擴(kuò)張，且耕作范圍有向三江平原逐步擴(kuò)張的趨勢.從行政區(qū)域上看，2001年作物種植主要分布在哈爾濱、齊齊哈爾和綏化，水稻種植區(qū)主要在南部，玉米和大豆主要分布于西南部和東北部；2010年，中東部的鶴崗、佳木斯和雙鴨山種植面積也逐步擴(kuò)大，多為玉米和大豆的混種模式；2017年，受到市場、政策等多方面因素的影響，種植結(jié)構(gòu)進(jìn)一步調(diào)整，大豆種植區(qū)主要分布于黑河和齊齊哈爾，佳木斯種植面積擴(kuò)大且結(jié)構(gòu)變化較為明顯，作物類型以玉米、水稻為主.

3.2 作物種植結(jié)構(gòu)時空變化分析

2001至2017年，黑龍江省主要農(nóng)作物種植面積顯著增長，種植結(jié)構(gòu)發(fā)生了明顯變化.為了探究種植結(jié)構(gòu)調(diào)整的空間分布格局，將整個研究區(qū)主要作物種植空間信息分布圖以市為行政單位進(jìn)行劃分，同時按經(jīng)度和緯度兩個方向進(jìn)行分析，以獲取其在時空變化上的規(guī)律.具體規(guī)則為：分市劃分時將4種作物的種植面積總和計為單位“1”，統(tǒng)計每個市范圍內(nèi)每種作物所占的百分比(圖3)；在經(jīng)緯度上，每隔0.5°統(tǒng)計各作物種植面積(圖4).

2001年

2010年

2017年

圖3 黑龍江省主要農(nóng)作物種植面積分市統(tǒng)計Fig.3 Statistical chart of Heilongjiang Province’s main crop planting surface integral city

由圖3可知，除了大小麥的種植面積明顯減少，玉米和水稻一直處于平穩(wěn)增長的狀態(tài)，而伊春和黑河兩市大豆種植占比顯著升高.玉米在近20年中的種植面積逐漸增大，并大幅向東部擴(kuò)張；水稻面積主要變化集中在黑龍江省南部和三江平原地區(qū)，原因是那里水熱條件好，積溫充足，更有利于水稻種植；大豆種植區(qū)主要向北緯48°—51°、東經(jīng)124°—126°遷移；大小麥的分布在經(jīng)緯上都呈現(xiàn)出少量分布到零星分布的趨勢(圖4).

3.3 精度評價分析

由于本次進(jìn)行作物提取選用的MODIS影像空間分辨率為250 m,不可避免地存在混合像元.但由于黑龍江省多以大面積成片種植為主，純凈像元較多，因此僅通過歷年農(nóng)業(yè)統(tǒng)計年鑒就可以對提取結(jié)果進(jìn)行精度評價.表3顯示，利用MODIS-NDVI時間序列數(shù)據(jù)集構(gòu)建的黑龍江省主要農(nóng)作物遙感提取模型最高精度為92.1%，最低精度77.0%，平均精度87.4%.

表3 黑龍江省主要農(nóng)作物遙感提取模型精度驗證Tab.3 Precision verification of remote sensing extraction model of main crops in Heilongjiang Province

主要從以下幾個方面對誤差精度進(jìn)行評價分析：1)在構(gòu)建黑龍江省主要農(nóng)作物種植面積遙感提取模型中，閾值T1—T13主要基于各農(nóng)場分布及選取的感興趣區(qū)域(ROI，region of interest)設(shè)定，且模型選取的時間節(jié)點并非完全準(zhǔn)確，有可能對提取結(jié)果產(chǎn)生影響.2)本次作物提取使用的MOD09Q1產(chǎn)品，時間分辨率為8 d，不能完全對應(yīng)研究區(qū)農(nóng)作物完整生長周期的關(guān)鍵時間節(jié)點，可能對作物面積信息提取的結(jié)果產(chǎn)生一定程度誤差.3)本研究中使用的MODIS影像空間分辨率為250 m，即一個像元內(nèi)存在像元“同物異譜”“同譜異物”的現(xiàn)象，對中小散戶混合種植信息提取不夠精確，從而影響作物的提取面積和精度.

4 結(jié)論

MODIS是具有中等空間分辨率、高時間分辨率且免費獲取的影像數(shù)據(jù)，在大尺度區(qū)域作物信息識別和提取等方面具有明顯優(yōu)勢.本文以2001、2010和2017年MODIS數(shù)據(jù)產(chǎn)品MOD09Q1為數(shù)據(jù)源，通過Hants濾波分析平滑MODIS-NDVI時間序列數(shù)據(jù)集，根據(jù)研究區(qū)主要農(nóng)作物的種植區(qū)劃、種植結(jié)構(gòu)和物候歷，分析玉米、水稻、大豆和大小麥的時間序列曲線特征，確定提取不同農(nóng)作物的特征點和參量，構(gòu)建農(nóng)作物遙感提取模型，設(shè)定閾值，利用物候關(guān)鍵特征參數(shù)值對黑龍江省主要作物進(jìn)行分類提取.結(jié)果顯示：1)黑龍江省在近20年中主要農(nóng)作物種植面積顯著提高，且擴(kuò)張范圍主要在東北三江平原和中北部地區(qū)；2)種植結(jié)構(gòu)方面，大小麥的種植面積持續(xù)減少，呈零散種植，玉米、水稻和大豆的種植面積持續(xù)增長，且玉米的增速最高；3)構(gòu)建的農(nóng)作物種植面積遙感提取模型的最高精度達(dá)到92.1%，最低精度77.0%，平均精度87.4%，表明通過MODIS影像可以有效對大范圍耕作區(qū)域進(jìn)行農(nóng)作物信息提取，能夠為區(qū)域農(nóng)情遙感監(jiān)測、作物種植結(jié)構(gòu)調(diào)整等提供可行的參考和借鑒.

本研究對于大范圍耕作區(qū)域農(nóng)作物種植面積提取效果較好，但在模型構(gòu)建、參數(shù)選擇、多源數(shù)據(jù)融合等方面還存在缺陷，未來的研究方向主要有：1)由于作物是經(jīng)過長期自然進(jìn)化和人工選擇后形成的[24]，而中國各地氣候類型和地域性差異較大，農(nóng)業(yè)生產(chǎn)利用面積廣且季節(jié)性較強(qiáng)，散戶、個體種植導(dǎo)致各地農(nóng)業(yè)種植習(xí)慣和種植結(jié)構(gòu)極其復(fù)雜，因此在作物面積提取前根據(jù)調(diào)查和統(tǒng)計信息及時更新不同地區(qū)的物候歷和種植結(jié)構(gòu)是農(nóng)作物種植信息提取本地化的基礎(chǔ)[14]；2)根據(jù)不同農(nóng)作物的物候歷特征及光譜特征，及時更新遙感提取模型中的關(guān)鍵數(shù)據(jù)閾值是保證模型精度的前提；3)由于MODIS影像為250 m中等空間分辨率，必然存在較多的混合像元，怎樣在作物種植信息提取中考慮混合像元分解問題、如何與其他中高分辨率的衛(wèi)星影像結(jié)合來提高提取精度都是未來研究的重點方向.