黃淮海地區(qū)耕地復(fù)種指數(shù)的時(shí)空格局演變

李 卓,劉淑亮,孫然好,*,劉維忠

1 中國(guó)科學(xué)院生態(tài)環(huán)境研究中心 城市與區(qū)域生態(tài)國(guó)家重點(diǎn)實(shí)驗(yàn)室,北京 100085 2 遼寧工程技術(shù)大學(xué)測(cè)繪與地理科學(xué)學(xué)院,阜新 123000 3 黑龍江第二測(cè)繪工程院,哈爾濱 150000

我國(guó)正處于國(guó)民經(jīng)濟(jì)發(fā)展的戰(zhàn)略機(jī)遇期,如何在新形勢(shì)下應(yīng)對(duì)城鎮(zhèn)化進(jìn)程加快和農(nóng)村經(jīng)濟(jì)結(jié)構(gòu)調(diào)整所造成的大量耕地持續(xù)向建設(shè)用地轉(zhuǎn)換的窘境,解決人口增長(zhǎng)與耕地資源緊缺的突出矛盾、確保耕作效率是當(dāng)務(wù)之急[1-2]。面對(duì)人多地少的基本國(guó)情,掌握作物物候規(guī)律,以“時(shí)間交換空間”的多熟種植方式是增加糧食產(chǎn)量的有效途徑之一[3]。復(fù)種指數(shù)作為耕種制度的重要內(nèi)容,反映了農(nóng)事活動(dòng)對(duì)農(nóng)業(yè)資源要素的利用方式,同時(shí)也是宏觀衡量耕地利用程度的關(guān)鍵指標(biāo)[4- 6]。以往的農(nóng)業(yè)決策多基于統(tǒng)計(jì)數(shù)據(jù)中的耕地面積和播種面積來計(jì)算復(fù)種指數(shù),這種計(jì)算方式時(shí)效性差并缺少空間維度信息,同時(shí)還摻雜了大量的人為誤差[7]。當(dāng)今,天地一體化的遙感技術(shù)日益成熟,為高效提取大尺度區(qū)域作物種植信息提供了優(yōu)異的條件,尤其是高時(shí)頻遙感數(shù)據(jù)的出現(xiàn)提高了植被指數(shù)時(shí)間序列的連續(xù)性和準(zhǔn)確性,進(jìn)一步促進(jìn)了復(fù)種指數(shù)遙感監(jiān)測(cè)的發(fā)展[8- 9]。Devendr等[10]研究發(fā)現(xiàn)亞洲地區(qū)糧食增產(chǎn)主要由耕地復(fù)種方式實(shí)現(xiàn);Verburg等[11]分析了1991年中國(guó)復(fù)種指數(shù)的區(qū)域性差異;丁明軍等[12]利用全國(guó)的SPOT-NDVI數(shù)據(jù)得到1999—2013年的耕地復(fù)種指數(shù)的時(shí)空變化曲線,并得出中國(guó)耕地復(fù)種指數(shù)從北至南從半濕潤(rùn)區(qū)到濕潤(rùn)區(qū)逐漸增大的結(jié)論。

黃淮海地區(qū)作為我國(guó)糧食主產(chǎn)區(qū),近30年來為我國(guó)糧食增產(chǎn)貢獻(xiàn)45%[13],在《全國(guó)新增1000億斤糧食生產(chǎn)能力規(guī)劃(2009—2020年)》中,國(guó)務(wù)院分配164.5×108kg的新增糧食產(chǎn)能任務(wù),進(jìn)一步奠定了黃淮海糧食產(chǎn)區(qū)在保障我國(guó)糧食安全中的重要地位[14]。尤其是該區(qū)域北部京津冀地區(qū)作為我國(guó)經(jīng)濟(jì)的第三增長(zhǎng)極,加之特殊的地形地貌和人文氣候條件,傳統(tǒng)農(nóng)耕已逐漸向生鮮食品生產(chǎn)轉(zhuǎn)變,用來滿足特大城市群地區(qū)人口的生活所需。陳麗等[15]以黃淮海地區(qū)糧食均衡增產(chǎn)為研究目標(biāo),發(fā)現(xiàn)糧食增產(chǎn)潛力空間表現(xiàn)為南部高于西北部,并強(qiáng)調(diào)糧食生產(chǎn)布局合理優(yōu)化的重要性。閆慧敏等[16]利用遙感技術(shù)結(jié)合地面作物物候觀測(cè)數(shù)據(jù)分析了黃淮海地區(qū)二熟制作物的生長(zhǎng)過程、物候特征和作物植歷的空間差異;洪舒蔓等[14]以縣級(jí)區(qū)域?yàn)檠芯繂卧?討論了黃淮海地區(qū)耕地資源時(shí)空變化特征及其對(duì)糧食生產(chǎn)格局變化的影響。先前研究缺少基于長(zhǎng)時(shí)間序列復(fù)種指數(shù)的提取及其空間格局變化分析,不利于種植結(jié)構(gòu)的調(diào)整和種植方式的優(yōu)化。此外,隨著精耕細(xì)作的不斷推行,尤其是在山地、丘陵等地域,多熟制經(jīng)濟(jì)作物逐漸被廣泛耕種,以往只考慮一熟制、兩熟制顯然不夠充分。

鑒于此,本文以黃淮海地區(qū)為研究區(qū)域,使用2001—2015年1 km分辨率的MODIS NDVI 16 d影像,通過Savitzky-Golay濾波對(duì)時(shí)間序列曲線作平滑重構(gòu),采用二次差分法結(jié)合物候信息提取時(shí)間序曲線峰值個(gè)數(shù),進(jìn)而確定復(fù)種次數(shù),最后結(jié)合耕種區(qū)邊界利用R語(yǔ)言繪制復(fù)種指數(shù)空間分布圖,旨在為黃淮海地區(qū)糧食生產(chǎn)安全政策制定和土地利用規(guī)劃提供參考依據(jù)。

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

1.1 研究區(qū)概況及數(shù)據(jù)來源

為了充分辨識(shí)黃淮海地區(qū)復(fù)種指數(shù)的時(shí)空差異,同時(shí)兼顧研究區(qū)的完整性,本研究以行政區(qū)界限為基準(zhǔn)選取京、津、冀、魯、豫三省兩市為研究對(duì)象。該區(qū)域東臨黃海,西起太行山,南接大別山,北達(dá)燕山南麓,土地面積約為53.66×104km2,其中由黃河以南至淮河北岸的黃淮海平原占研究區(qū)總面積的53%。境內(nèi)屬溫帶季風(fēng)氣候,全年積溫約4500℃,無霜期最高達(dá)到230 d,年降水在600—800 mm之間,土壤以穩(wěn)產(chǎn)、高產(chǎn)的棕壤和褐土土壤為主,有機(jī)質(zhì)含量在15—25 g/kg[17]。良好的氣候條件加上優(yōu)質(zhì)的土壤資源,使其成為我國(guó)傳統(tǒng)農(nóng)業(yè)和農(nóng)耕文明的發(fā)源地之一。然而,該區(qū)域人口數(shù)量巨大,加之“京津冀一體化”、“山東藍(lán)色海洋經(jīng)濟(jì)區(qū)”等經(jīng)濟(jì)開發(fā)建設(shè)工程的推行,導(dǎo)致耕地面積縮減嚴(yán)重,未來如何在有限的耕地上提高復(fù)種指數(shù)成為當(dāng)今關(guān)注的焦點(diǎn)。

復(fù)種指數(shù)提取使用的遙感影像來源于NASA(美國(guó)國(guó)家航空和空間管理局)的MOD13A2數(shù)據(jù)集,該產(chǎn)品是空間分辨率為1 km的 MODIS(Terra星)16 d數(shù)據(jù)。本研究時(shí)間跨度為2001—2015年,在MRT(Modis Reprojection Tool)軟件中進(jìn)行鑲嵌、投影等預(yù)處理,最后通過研究區(qū)掩膜裁剪獲得345期影像。黃淮海耕種區(qū)邊界由中國(guó)西部環(huán)境與生態(tài)科學(xué)數(shù)據(jù)中心提供。

1.2 時(shí)序曲線的熟制分析

歸一化植被指數(shù)(Normalized Difference Vegetation Index,NDVI)作為當(dāng)前主流的植被生長(zhǎng)狀況反映指標(biāo),在尺度、時(shí)效以及表征植被的生物物理特征等方面具有明顯優(yōu)勢(shì)[18]。計(jì)算公式為:

(1)

式中,ρNIR為近紅外波段地表反射率,ρR表示紅波段地表反射率。選用的MODIS NVDI產(chǎn)品每年有23個(gè)時(shí)像,通過建立NDVI值和時(shí)間的二維坐標(biāo)系,可直觀反映農(nóng)作物的生長(zhǎng)周期。一年中農(nóng)作物從播種到收獲,NDVI值隨物候特征經(jīng)歷了升高-降低的動(dòng)態(tài)過程,從而形成了單峰或多峰時(shí)序曲線,每個(gè)波峰點(diǎn)代表了農(nóng)作物生長(zhǎng)周期的最旺盛點(diǎn)。因此,可以認(rèn)為峰值點(diǎn)是最佳的熟制判斷點(diǎn),進(jìn)而通過提取NDVI時(shí)間序列中的峰值頻數(shù)判定熟制[19]。

1.3 時(shí)序曲線平滑重構(gòu)

太陽(yáng)高度角、大氣氣溶膠以及積雪等外界因素影響近紅外波段在地表的反射強(qiáng)度,導(dǎo)致時(shí)序曲線斜率迅速下降[20]。突變形成的鋸齒狀“偽波峰”會(huì)削弱時(shí)序曲線的周期性變化趨勢(shì),影響復(fù)種指數(shù)提取的精度,所以對(duì)時(shí)序曲線進(jìn)行平滑重構(gòu)是提取復(fù)種指數(shù)的關(guān)鍵步驟。

Savitzky-Golay濾波理論簡(jiǎn)單,是一種利用局部多項(xiàng)式回歸模型平滑時(shí)序數(shù)據(jù)的時(shí)域低通濾波方法[21]。公式如下:

(2)

式中,Y*是函數(shù)擬合值,Y為原始數(shù)據(jù),N=2m+1為函數(shù)濾波窗口的值,Ci表示第i個(gè)點(diǎn)的權(quán)重。

該方法對(duì)時(shí)間序列曲線平滑去噪的保真性效果主要受活動(dòng)窗口大小和多項(xiàng)式階數(shù)的影響,本研究參照申健等[3]對(duì)關(guān)中地區(qū)復(fù)種指數(shù)的提取的方法,在TIMESAT軟件中,多次調(diào)試窗口大小和擬合次數(shù),最終確定(4,2)為最佳擬合參數(shù)值。此外,為了削弱S-G濾波的邊緣效應(yīng)(越冬作物的冬前峰現(xiàn)象),選取一年的完整數(shù)據(jù)、該年上一年的下半年和該年下一年的上半年組成一個(gè)研究時(shí)段,即通過2001—2003年、2004—2006年、2007—2009年、2010—2012年、2013—2015年5個(gè)時(shí)段,分別提取2002年、2005年、2008年、2011年、2014年5個(gè)年份復(fù)種指數(shù)。

1.4 二次差分法提取復(fù)種指數(shù)

對(duì)于大區(qū)域尺度時(shí)間序列NDVI峰值的提取,二次差分法應(yīng)用較為廣泛,并且取得了較好的探測(cè)效果[3- 4,22]。二次差分基本原理如下[23]:首先計(jì)算相鄰NDVI的差值,記為序列S1,公式如(3);其次根據(jù)公式(4),對(duì)S1序列判斷正負(fù),并重新賦值,記為序列S2;最后依次求S2序列前后元素差,記為序列S3。

S1i=NDVIi-NDVIi-1

(3)

(4)

S3i=S2i+1-S2i

(5)

式中,i代表序列中第i個(gè)元素。作物時(shí)序曲線上的波峰出現(xiàn)在序列S3中元素為-2且前后元素皆為0的位置。

盡管時(shí)序曲線經(jīng)過了S-G濾波的平滑重構(gòu),但仍有部分離散峰值點(diǎn)被提取出來,產(chǎn)生誤差。主要原因有兩個(gè):一是由于非耕種區(qū)NDVI值較低,平滑后時(shí)序曲線特征不明顯,微小的波動(dòng)都形成了“偽波峰”,造成了錯(cuò)誤提取;二是非生長(zhǎng)期的雜草以及其他植被的生長(zhǎng),會(huì)出現(xiàn)雙峰現(xiàn)象,影響提取精度。因此,還需要添加相應(yīng)的約束條件對(duì)探測(cè)到的波峰進(jìn)行取舍。通過分析黃淮海地區(qū)作物時(shí)序曲線和作物物候特征,本研究設(shè)定“橫縱二維約束”:(1)波峰的峰值要大于0.5;(2)兩相鄰的波峰之間至少包括4個(gè)時(shí)相。

獲取耕地像元峰值頻數(shù)后,使用公式(6)計(jì)算各行政單位復(fù)種指數(shù)。

(6)

式中,Mi表示從單個(gè)的NDVI時(shí)間序列曲線中提取的峰值頻數(shù),N表示研究區(qū)或者行政區(qū)劃內(nèi)的像元總個(gè)數(shù),C表示耕地復(fù)種指數(shù)。

2 結(jié)果與分析

2.1 NDVI時(shí)序曲線平滑處理效果

作物在不同的生長(zhǎng)階段(播種、出苗、拔節(jié)、抽穗、成熟),NDVI時(shí)序曲線隨之出現(xiàn)相應(yīng)的波動(dòng),這種鋸齒狀的波動(dòng)不利于熟制的判斷和提取,利用Savitzky-Golay濾波對(duì)2001—2015年MODIS NDVI時(shí)間序列曲線做平滑去噪重構(gòu)。對(duì)比原始曲線,如圖1所示,無論是一熟制,還是兩熟制,平滑處理在剔除干擾值的同時(shí),保持了原始曲線的基本特征,并突出刻畫了波峰與波谷相交替的狀態(tài),良好的表達(dá)了農(nóng)作物生長(zhǎng)的節(jié)律性,更適合于熟制提取。

圖1 一年一熟和一年兩熟NDVI時(shí)序曲線平滑前后效果對(duì)比Fig.1 Comparison between before and after smoothing time-series NDVI date of one crop per annual and two crops per annual

2.2 黃淮海地區(qū)復(fù)種指數(shù)及時(shí)空演變特征

復(fù)種指數(shù)的分布格局是受自然因素影響,經(jīng)過長(zhǎng)期耕種逐漸形成的。由于人為操控因素的主導(dǎo)作用,短時(shí)間內(nèi)復(fù)種指數(shù)格局不會(huì)發(fā)生較大改變,本研究在R環(huán)境下以3年數(shù)據(jù)為一期,分別繪制2002年、2005年、2008年、2011年、2014年共計(jì)5個(gè)年份復(fù)種指數(shù)空間分布圖(圖2)。從總體分布特征來看,一熟制和兩熟制提取結(jié)果基本符合我國(guó)復(fù)種指數(shù)隨緯度增加而減小的大趨勢(shì),三熟制由于受到人類需求影響,作物的選擇和利用方式差異性較大,分布特征不明顯。不同熟制空間分布特征主要有:一年一熟制主要分布于河南省西部、山東省中東部和河北省大部分地區(qū);一年兩熟制則集中于河北省南部、魯西南以及河南省東部,這里屬于華北平原農(nóng)耕區(qū)的中心地帶,受水分、熱量和經(jīng)濟(jì)等條件的影響,已形成歷史悠久的冬小麥-夏玉米的耕作制度;一年三熟制分布格局不集中,以河北省西北部為主要代表。

圖2 黃淮海地區(qū)復(fù)種指數(shù)空間分布圖Fig.2 Spatial distribution maps of multiple cropping index in Huang-Huai-Hai regions

從行政區(qū)上講,河北省南北緯度跨度較大,北部積溫、降水等農(nóng)作物生長(zhǎng)的自然條件受到一定程度的限制,只能滿足一季作物的生長(zhǎng);兩熟制從南至北形成了一條逐漸變窄的“綠帶”;河北省是京津冀重要的生鮮食品供應(yīng)基地,三熟制分布在河北省西北和北京、天津城市周邊這里主要以生長(zhǎng)周期較短的經(jīng)濟(jì)作物為主,尤其是城市周邊地區(qū)耕地已經(jīng)由標(biāo)準(zhǔn)的北方玉米、小麥等簡(jiǎn)單的作物類型日益多樣化,耕地種植大量的草皮和花卉等“裝飾性的農(nóng)業(yè)”開始起步,并逐漸成熟[24]。山東省中東丘陵地帶,素有“水果之鄉(xiāng)”等稱譽(yù),以種植一年一熟的水果為主;西南部土壤肥沃,農(nóng)業(yè)生產(chǎn)人口比重較大,是一年兩熟耕種區(qū)。河南省東部大部分地區(qū)是國(guó)家高標(biāo)準(zhǔn)農(nóng)田,種植條件與山東西部相似,屬于一年兩熟制。北京、天津經(jīng)濟(jì)發(fā)達(dá),耕地資源有限,主要以一年一熟的水稻和經(jīng)濟(jì)作物為主。

圖3 黃淮海地區(qū)不同省市復(fù)種指數(shù)年際變化趨勢(shì)Fig.3 The inter-annual change of multiple cropping index in provence or city

黃淮海地區(qū)復(fù)種指數(shù)存在顯著的時(shí)空差異(圖3)。復(fù)種指數(shù)均值從大到小依次是河南省(165.2%)、山東省(140.5%)、河北省(128.8%)、北京市(104.5%)、天津市(103.4%)。各省市復(fù)種指數(shù)年際變化趨勢(shì)表現(xiàn)較為一致:2002—2005年升高、2005—2011年降低、2011—2014年升高。自然條件和農(nóng)村經(jīng)濟(jì)結(jié)構(gòu)是影響復(fù)種指數(shù)的主要因素,短時(shí)間內(nèi),自然條件通常不會(huì)發(fā)生劇烈變化。查閱資料驗(yàn)證,黃淮海地區(qū)在研究時(shí)段未出現(xiàn)大規(guī)模自然災(zāi)害,因此,復(fù)種指數(shù)的年際變化可能取決于耕地的收益狀況和輪作。2004年國(guó)家對(duì)從事糧食生產(chǎn)的農(nóng)民實(shí)行糧食直補(bǔ)政策,提高了農(nóng)民種糧的積極性,復(fù)種指數(shù)得到增加。但2005年后,糧食價(jià)格持續(xù)走低,農(nóng)民更傾向于經(jīng)濟(jì)價(jià)值更高的一季作物,導(dǎo)致復(fù)種指數(shù)降低。2011年后,國(guó)家宏觀調(diào)控糧食生產(chǎn)使復(fù)種指數(shù)得到提升。

2.3 與其他遙感監(jiān)測(cè)結(jié)果比較

由于參考數(shù)據(jù)有限,本文搜集整理了2008年之前的相關(guān)研究成果與本文進(jìn)行對(duì)比(表1)。從整體來看,相對(duì)誤差位于-9.3%—2.5%區(qū)間內(nèi),表現(xiàn)出較高的耦合性。其中,最大值出現(xiàn)在河南省(2008年,-9.3%),最小出現(xiàn)在北京市(2002年,-0.5%)。梁守真等[25]在研究2000—2009年環(huán)渤海地區(qū)復(fù)種指數(shù)時(shí),與本文使用了相同的數(shù)據(jù)源、不同的處理方法(SPLINE差值、領(lǐng)域比較法),對(duì)比結(jié)果,最大相對(duì)誤差僅6.6%,是對(duì)本文采用Savitzky-Golay濾波、二次差分技術(shù)反演黃淮海地區(qū)復(fù)種指數(shù)可行性的驗(yàn)證。另外,對(duì)比數(shù)據(jù)可以看出,本文提取的復(fù)種指數(shù)略低于前人的研究結(jié)果,是由于判斷熟制進(jìn)行波峰篩選設(shè)置的閾值過于嚴(yán)格引起的,但誤差均在可接受范圍內(nèi)。

表1 與其他遙感監(jiān)測(cè)研究結(jié)果對(duì)比

相對(duì)誤差計(jì)算公式:(監(jiān)測(cè)值-文獻(xiàn)值)/監(jiān)測(cè)值。a代表數(shù)據(jù)來自梁守真[25],b代表數(shù)據(jù)來自朱孝林[23],c代表數(shù)據(jù)來自左麗君[26],d代表數(shù)據(jù)來自唐鵬欽[27]

3 結(jié)論與討論

本研究以黃淮海地區(qū)為例,基于2005—2015年MODIS NDVI時(shí)序數(shù)據(jù),采用Savitzky-Golay濾波和二次差分技術(shù)提取耕地復(fù)種指數(shù)及其空間分布特征,以期為辨識(shí)人類活動(dòng)方式和土地利用強(qiáng)度提供參考,同時(shí)也可以為黃淮海地區(qū)農(nóng)業(yè)種植結(jié)構(gòu)調(diào)整和種植方式優(yōu)化提供決策依據(jù)。研究結(jié)果主要有:

(1)黃淮海地區(qū)三省兩市中,河南省復(fù)種指數(shù)最高,在2014年達(dá)到了169.3%,山東省次之,天津市最小。各省市在研究時(shí)段內(nèi),復(fù)種指數(shù)年際變化趨勢(shì)大體一致,經(jīng)歷了升高-降低-升高的過程。自然條件和農(nóng)村經(jīng)濟(jì)結(jié)構(gòu)是影響復(fù)種指數(shù)的主要因素,短時(shí)間內(nèi),自然條件通常不會(huì)發(fā)生劇烈變化。因此,復(fù)種指數(shù)的年際變化主要取決于耕地的收益狀況和輪作。從空間分布特征來看,黃淮海地區(qū)耕地熟制具有明顯的地域性差異,二熟制主要集中于南部,東部和北部受地形和緯度影響主要以一熟制為主。

(2)與前人監(jiān)測(cè)結(jié)果比較,最大相對(duì)誤差(-9.3%)驗(yàn)證了本研究設(shè)置含有閾值的二次差分算法在提取黃淮海地區(qū)復(fù)種指數(shù)的可行性。同時(shí),根據(jù)三熟制的零星散布格局可以看出,將三熟制納入研究?jī)?nèi)容,并未影響到黃淮海地區(qū)南部傳統(tǒng)農(nóng)耕區(qū)的熟制識(shí)別,佐證了研究成果的可靠性。

中國(guó)農(nóng)業(yè)耕作具有典型的小農(nóng)特點(diǎn),農(nóng)作物復(fù)雜多樣且在地表不規(guī)則分布,不同作物之間或作物與裸地組成的混合像元是影響復(fù)種指數(shù)提取精度的主要限制因素。針對(duì)這些問題,未來利用多平臺(tái)、多分辨率遙感數(shù)據(jù)以及地面觀測(cè)數(shù)據(jù),確定不同區(qū)域的最優(yōu)提取算法及其本地化參數(shù),是實(shí)現(xiàn)復(fù)種指數(shù)的高精度監(jiān)測(cè)的有效途徑之一。此外,復(fù)種指數(shù)受經(jīng)濟(jì)形勢(shì)和農(nóng)民意愿等因素影響較大,因此,政策和管理措施的量化指標(biāo)有助于判斷復(fù)耕指數(shù)的閾值和范圍。