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

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

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

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

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

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

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

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

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

1.2 時序曲線的熟制分析

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

(1)

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

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

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

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

(2)

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

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

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

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

S1i=NDVIi-NDVIi-1

(3)

(4)

S3i=S2i+1-S2i

(5)

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

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

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

(6)

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

2 結(jié)果與分析

2.1 NDVI時序曲線平滑處理效果

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

圖1 一年一熟和一年兩熟NDVI時序曲線平滑前后效果對比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ù)及時空演變特征

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

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

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

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

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

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

由于參考數(shù)據(jù)有限,本文搜集整理了2008年之前的相關(guān)研究成果與本文進行對比(表1)。從整體來看,相對誤差位于-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ù)據(jù)源、不同的處理方法(SPLINE差值、領(lǐng)域比較法),對比結(jié)果,最大相對誤差僅6.6%,是對本文采用Savitzky-Golay濾波、二次差分技術(shù)反演黃淮海地區(qū)復(fù)種指數(shù)可行性的驗證。另外,對比數(shù)據(jù)可以看出,本文提取的復(fù)種指數(shù)略低于前人的研究結(jié)果,是由于判斷熟制進行波峰篩選設(shè)置的閾值過于嚴格引起的,但誤差均在可接受范圍內(nèi)。

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

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

3 結(jié)論與討論

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

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

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

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