基于NDVI時序數據的華北地區耕地物候參數時空變化特征*

查 燕，宋 茜，衛 煒，陳仲新，楊 鵬※

(1.中國農業科學院農業資源與農業區劃研究所/農業部農業遙感重點實驗室，北京 100081；2.農業部規劃設計研究院，農業部耕地利用遙感重點實驗室，北京 100125)

0 引言

植被物候的動態變化是生物圈對全球氣候變化以及陸地水文循環機制變化的動態響應，與陸地生態系統初級生產力和碳循環息息相關，因而成為全球變化的研究熱點之一，在全球變化、生態學和農業生產應用等領域發揮了重要作用[1-3]。農作物物候期數據是重要的農業信息，是農業生產、田間管理、計劃決策等的重要依據，也是作物模擬模型的重要參數[4-5]。農業物候期信息可以通過田間觀測法、積溫法和遙感法等多種途徑獲得[6]。由于衛星遙感具有宏觀、高效和便捷的特點，能夠提供大范圍的多時相重復觀測數據，為群落尺度和區域尺度的物候研究提供了有利條件，逐漸成為物候信息獲取的重要手段[7]。目前，歸一化植被指數(normalized difference vegetation index，NDVI)是應用最為廣泛的植被指數，它對植被生長狀況、植被覆蓋類型等信息非常敏感，能以最直觀的形式反映作物從播種、出苗、抽穗、成熟和收割等物理過程，可用于跟蹤作物的季節性動態變化[8-9]。

對于大區域而言，SPOT/VGT NDVI 長時間序列數據對植被生長期特征監測有較好地反應[10]，其逐旬的時間分辨率可以彌補空間分辨率的不足，可有效地對農作物長勢進行動態監測[11]，因此，SPOT/VGT NDVI長時間序列數據已經在作物生長過程監測、作物面積提取、作物類型識別和作物種植制度判定等方面得到了廣泛應用[12-14]。利用SPOT/VGT 的旬合成NDVI數據對區域尺度的耕地復種指數研究較多，如范錦龍等[15]利用諧函數分析法(HANTS)提取了1999～2002年全國耕地復種指數； 唐鵬欽等[16]利用小波變換提取了2007年華北平原的耕地復種指數； 丁明軍等[17]利用Savitzky-Golay濾波法研究了1999～2013年中國耕地復種指數的時空變化。對于耕地典型物候期的識別，李正國等[18-19]采用非對稱高斯函數法(Asymmetric Gaussian，AG)對華北地區(1998～2007年)和東北地區(1998～2009年)耕地物候期進行了識別； 衛煒等[20]基于2005年SPOT/VGT數據，提取了我國北方耕地物候信息。作為地表物候的重要組成部分，不同于草原、灌叢和森林等自然植被在年內完成一個生長過程，耕地由于受人類活動影響大和多熟種植制度的存在，耕地物候比自然植被物候更為復雜[9]。然而，目前耕地物候的研究都只是停留在某一年或過去一個時期上，很少進行長時間的動態監測。因此開展基于時序植被指數提取耕地生長季特征研究，能更好地反映長時間區域尺度上作物物候期特征的時空差異。

文章以我國華北地區5省(市)為研究區域，基于1km×1km 分辨率的SPOT/VGT逐旬NDVI數據，利用非對稱高斯函數擬合法對NDVI數據進行平滑去噪，重建作物生長植被指數曲線，通過比例閾值法提取了1999年和2013年華北地區耕地物候參數(生長季開始期和結束期)，并用波峰閾值法對耕地種植制度進行了識別。研究1999～2013年華北地區物候參數和種植制度的時空變化過程，為中國農業政策制定及管理提供科學依據。

1 研究區域與數據

1.1 研究區域概況

研究區域位于我國糧食主產區之一的華北平原地區，主要包括北京市、天津市、河北省、河南省和山東省，耕地面積約1 700萬hm2。大部分區域屬于暖溫帶大陸性季風氣候，氣候條件良好，年平均溫度13.0℃，年均總輻射量5 100～5 300MJ/m2，>0℃年積溫4 200～5 500℃，無霜期為170～220d，熱量條件可以滿足一年一熟和一年兩熟農作需要。年降水量為500～900mm，主要集中在夏季的7～9月份，季節分配不均。主要農作物包括冬小麥、玉米、棉花和大豆。

1.2 數據來源

該文所采用的遙感數據是比利時佛萊芒技術研究所(Flemish Institute for Technological Research Vito)免費提供的SPOT VGT 逐旬 NDVI 最大合成數據(http：//www.vito-eodata.be)，空間分辨率為1km，時間范圍為1999年和2013年全年，共72景影像。該數據已經完成了一系列幾何校正、輻射校正、地圖投影、狀態標識以及大氣校正等處理。該數據集包含了各波段的反射率、NDVI和太陽方位角等多個數據層，利用VITO提供的處理工具VGTExtract提取其中的NDVI數據層，并將影像由經緯度轉換為雙標準緯線Albers投影，并通過與區域行政邊界掩膜處理，得到研究區域內各時相的遙感影像。數據產品中NDVI的DN值范圍為0～250，通過公式NDVI=DN×0.004-0.1轉換成標準植被指數值。

耕地數據來自中國科學院遙感應用研究所等7個單位2005年共同完成的中國1： 25萬土地覆蓋遙感調查與監測數據庫。該數據庫內容包括森林、草地、農田、聚落、濕地與水體、荒漠等6個一級類型和25個二級類型，其中農田數據可分為水田、水澆地和旱地，空間分辨率為100m。研究將華北地區的農田數據提取出來作為耕地圖層，轉換其地圖投影為阿爾伯斯等面積圓錐投影并將空間分辨率重采樣到1 km使之與VGT-S10 NDVI數據相匹配，從而得到研究區域的耕地分布數據(圖1)。

圖1 研究區域耕地分布

2 研究方法

2.1 NDVI時序數據預處理

盡管VGT-S10 NDVI產品經過最大值合成處理后能夠消除部分云、氣溶膠、太陽高度角等干擾因素的影響，但仍然殘留了許多噪聲，不能清晰地反映出耕地像元作物的生長過程，因此必須對NDVI時間序列數據進行去噪和平滑等處理。

該文采用非對稱高斯函數來分段模擬植被的生長過程，最后通過平滑連接各段高斯擬合曲線實現時間序列的重構[21]。該方法是一種由局部最優化擬合到全局擬合的方法，具有一定的靈活性，避免了整體數據對局部擬合的干擾，使得擬合后的NDVI曲線可以較好描述NDVI時序數據中復雜的和微小的變化，更加接近真實情況[22]。具體步驟如下：

首先，提取原始時序數據曲線中的谷值和峰值，采用高斯函數分別擬合曲線的左右部分。局部擬合函數為：

f(t)=f(t；c1，c2，a1，…，a5)=c1+c2g(t；a1，…，a5)

(1)

式(1)中，g(t；a1，…，a5)為高斯函數；c1和c2為控制曲線的基準和幅度；a1是對應時間變量t的峰值或谷值的位置參數；a2、a3和a4、a5分別控制曲線右半部分和左半部分的寬度和平度。

局部擬合函數可以很好地描述最大值與最小值間隔內的NDVI數據變化曲線。但是在曲線突出部，擬合效果不是很好。在此將變化曲線劃分成左邊谷值間隔區、中部峰值間隔區與右邊谷值間隔區，分別以不同的局部擬合函數來描述，最后再利用各局部擬合函數構建整體擬合函數，從而較好地描述整個植被生長期的NDVI變化過程。整體擬合函數為：

(2)

式(2),[tL，tR]區間是NDVI數據中待擬合部分的變化區間；fL(t)、fC(t)和fR(t)分別代表[tL，tR]區間內左邊谷值、中間峰值及右邊谷值對應的局部擬合函數；α(t)和β(t)為介于0～1的剪切系數。

2.2 作物物候特征參數提取

經過去噪和平滑處理的NDVI時間序列曲線可以很好地反映作物生長的年內動態變化特征，因此根據NDVI時間序列曲線能提取出作物主要的物候信息，如作物的出苗期或返青期、抽穗期和收獲期等特征參數。動態閾值法是目前提取作物物候特征參數常用的方法，其考慮了NDVI季節變化幅度，在像元尺度上對閾值動態設定，從而一定程度上消除了土壤背景值和植被類型的影響。

在我國北方地區的耕地生長季開始前和結束后，可能存在一些影響物候參數準確提取的因素，如冬小麥種植區小麥冬前峰的影響。為避免這些因素的干擾，該文根據先驗知識對NDVI時間序列數據進行“掐頭去尾”處理，即將1月、11月和12月的數據從參與分析計算的數據中加以剔除。另外由于耕地物候不同于自然植被物候，其受人類耕作活動的影響較大，參考我國北方地區的作物物候歷將生長季開始和結束的閾值分別設置為10%和50%[20]。因此，當NDVI時間序列曲線達到的上升階段振幅的10%時即認為是生長季開始。類似地，定義NDVI時間序列曲線達到下降階段振幅的50%時為生長季結束。

2.3 作物種植制度提取

每種作物生長周期都對應著NDVI時間序列曲線上的一個波峰，多熟種植制度下會有多個生長周期，波峰的個數即代表了生長周期的個數。其中一年一季作物的NDVI在一年內形成明顯的單峰曲線，一年兩季作物的NDVI形成雙峰曲線，一年三季作物的NDVI形成三峰曲線。因此，可以通過NDVI時間序列曲線上波峰的數目來確定作物種植制度[23]。然而在實際的NDVI時間序列曲線上除了表征作物生長周期的有效波峰外，還存在著一部分由噪聲引起的干擾波峰。這些干擾波峰使得擬合曲線上的波峰個數比實際偏多，導致對生長周期個數的誤判。前人研究表明，對于農作物而言其生長季的NDVI峰值通常會達到0.4以上[4， 7]。因此，該文選擇0.4作為波峰閾值對探測到的峰值進行判定取舍，只有峰值大于0.4的波峰才被認為是能夠真實反映作物生長周期的有效波峰，每個耕地像元的種植制度由該像元NDVI時間序列曲線上有效波峰的個數來確定。

3 結果與分析

3.1 農作物典型物候期的空間分布特征

3.1.1 作物返青/出苗期

圖2 華北地區作物返青/出苗期空間分布

圖2a和2b分別描述了1999年華北地區第一季作物返青/出苗期和第二季作物出苗期空間分布。從分布特征看，第一季、第二季作物生長開始期都存在明顯的空間差異，具有從南向北逐漸推遲的空間特征。第一季作物返青/出苗期最早發生在2月上旬，主要是河南南部少數地區； 大部分一年兩熟區，如河南中部、北部，河北中部、南部的低平原區、太行山和燕山山前平原區，以及山東西部、西北部平原區，作物出苗期或返青期多開始于3月上旬至4月上旬。河北省北部和東部一年一熟區的作物生長季開始約在4月下旬至6月上旬，而山東省沂蒙山地丘陵區和魯東丘陵區，作物返青/出苗期多開始于5月上旬至6月下旬，也有少部分區域作物返青/出苗期開始較早，大約在3月上旬至4月下旬。北京和天津地區作物返青/出苗期多發生在4月上旬至6月上旬。第二季作物的出苗期多在6 月下旬至7月上旬，河南南部少數地區發生5 月下旬至6 月上旬。

圖3 夏玉米生長季開始期估算值與地面站點觀測數據對比散點

2013年華北地區第一季、第二季作物返青/出苗期空間分布與1999年相比，河南大部分地區、山東魯西、魯西北部平原區，河北南部地區，第一季作物返青/出苗期明顯提前，多發生在2月上旬至3月上旬。河北北部、北京、天津等地區，第一季作物返青/出苗期變化不大。而第二季作物出苗期明顯提前的是河南南部和中部地區，由1999年6月下旬7月上旬，提前至6月上旬。有研究表明，華北平原近50年來，年平均氣溫整體呈現顯著上升趨勢，尤其春季和冬季增溫對年均溫增溫的貢獻較大[24]。華北平原夏玉米播種至出苗期的平均溫度在近30 年的上升幅度為0.51℃/10年(P<0.05)，其中河北中部的石家莊、保定、廊坊和唐山該期的平均溫度增加顯著，平均增幅達1.04℃/10年(P<0.05)[25]。因此，華北地區作物生長開始期提前可能與氣候變暖關系密切。

圖4 華北地區作物成熟期空間分布

圖5 華北地區農作物種植制度的空間分布

圖6 華北地區農作物種植制度變化的空間分布 圖7 1999～2013年華北地區各省冬小麥播種面積與一年兩熟制耕地像元數對比

冬小麥夏玉米是華北地區一年兩熟制的典型作物，而夏玉米是最主要的第二季作物。利用2013年SPOT/VGT數據提取的第二季作物生長開始期與25個地面站點夏玉米觀測數據對比如圖3所示，二者相關系數為0.40，平均絕對誤差(Mean Absolute Error)和均方根誤差(Root Mean Square Error)分別為5.2和7.09，說明遙感提取結果與實際情況比較吻合，可以用來分析典型作物物候期的空間變化。

3.1.2 作物成熟期

華北地區作物成熟期空間分布也存在明顯的地域差異，河南除豫西和豫南的山地丘陵區以外，大部分地區作物成熟比較早。河北中南部的太行山和燕山山前平原區，山東的魯西和魯西北平原區等，由于第一季作物生長期開始早，因此其成熟期結束也較早，主要在5月下旬至6月下旬。第一季作物成熟早為第二季作物完成整個生長周期提供了足夠的熱量資源保障。而河北的中部和北部、山東中部和東部，以及北京市和天津市等地區，作物返青/出苗期較晚，其生長結束期也就較晚，多發生在9月上旬至10月中上旬(圖4a，圖4c)。華北地區第二季作物成熟期比較集中，基本都集中在9月上旬至10月上旬，南部地區總體上要早于北部地區(圖4b，圖4d)。華北地區作物成熟期空間分布差異與該區域耕地第一季作物的生長開始期差異大有關外，還與不同區域的作物種類和作物生長周期不一樣有關[4]。

3.2 農作物種植制度的空間分布特征

圖5是利用NDVI數據分別提取的1999年和2013年華北地區耕地種植制度空間分布圖。從圖5可以看出，華北地區耕地種植制度以一年兩熟為主，一年一熟區主要分布在河北省北部和中部地區，山東省中部的沂蒙山地丘陵區和東部的魯東丘陵區，以及河南省西部和南部的山地丘陵區。一年兩熟區主要分布在河北省南部、山東省西部、西北部的平原區，以及河南省北部地區。北京市和天津市耕地主要為一年一熟，兩熟制區域面積很小。華北地區作物種植制度空間分布格局和該區域自然地理分布特征具有很好的一致性，因為耕地種植制度在很大程度上受水、熱、光和土等自然資源共同作用和影響[18]。

由于受農業技術、社會經濟發展，以及全球氣候變化等因素的綜合影響，經過15年的耕種(1999～2013年)，華北地區農作物種植制度的空間格局發生了很大的變化(圖6)，河北省中部和南部，以及河南省南部有部分區域由一年兩熟變成了一年一熟； 而河南西部、山東中部和東部地區，有零星的區域由一熟制變成了兩熟制； 華北地區北部受熱量資源制約，仍舊保持一年一熟制不變。與1999年相比， 2013年兩熟制種植面積下降了21.1%，而一熟制種植面積增加了38.7%。楊婷等[12]利用遙感數據研究中國近些年熟制變化結果表明，中國兩熟制區中2 819個像元變回了一熟制，這種現象在河北和山東較明顯。

華北地區絕大多數一年兩熟區，冬小麥是最主要的第一季作物。因此，為了驗證上述種植制度空間提取結果的可靠性，該文利用《中國農業統計資料》華北地區各省15年(1999～2013年)的冬小麥播種面積與各省一年兩熟制耕地像元的數量進行對比。從圖7可以看出，冬小麥播種面積與一年兩熟區具有極顯著的線性相關關系，R2為0.97，說明遙感提取結果符合實際情況，可以用來分析區域種植制度的空間變化。

4 討論

基于NOAA/AVHRR、SPOT/VGT和EOS/MODIS等中高分辨率NDVI時間序列數據已經在作物生長過程監測、作物類型識別、作物種植制度判定等方面得到了廣泛應用。NOAA/AVHRR的NDVI數據產品其空間分辨率通常為8km，由于我國華北地區耕地種植制度復雜多樣， 8km 的空間分辨率在一定程度上影響了耕地種植制度和物候信息結果的準確性[4]。MODIS提供了500m和1km兩種相對較高空間分辨率的NDVI數據產品，但是其時間分辨率為16d。SPOT/VGT的VGT-S10 NDVI產品具有適中的1km的空間分辨率和10d的時間分辨率，相對上述兩種數據而言更加適合我國北方地區的耕地物候研究[20]。

由于小麥冬前峰的存在可能使得重構的NDVI擬合曲線形狀發生變化，從而導致提取的生長季開始期比實際提前，而生長季結束期比實際推遲。在該研究中根據先驗知識將耕地生長季之外的1月、11月和12月的NDVI時間序列數據加以剔除，這種“掐頭去尾”處理，能夠避免小麥冬前峰的影響，進而準確獲取耕地物候信息。

耕地種植制度與物候分布格局和自然地理環境緊密相關，尤其是不同區域的溫度、降水和光照等氣候資源的稟賦和匹配程度對耕地種植結構和作物布局有很大影響[4， 26]。該文僅對物候期和種植制度時空變化特征進行了分析，未對這些影響因素進行深入探討。因此，未來計劃通過分析區域內主要農業氣候資源，如≥0℃積溫、平均氣溫、溫度生長期天數、降水等的時空變化特征，進一步完善耕地物候期對農業氣候資源特征的響應機制研究。

5 結論

該文基于SPOT/VGT數據提取了我國華北地區耕地典型物候期的空間格局，并對該區域種植制度進行了識別。研究結果表明，華北地區作物生長開始期和結束期都存在明顯的空間差異， 15年來河北北部、北京、天津等地區，第一季作物返青/出苗期變化不大。而河南南部和中部地區， 2013年第二季作物出苗期明顯提前，由1999年6月下旬7月上旬，提前至6月上旬。我國華北地區種植制度仍以一年二熟制為主，華北地區北部受熱量資源制約，仍舊保持一年一熟制不變。與1999年相比， 2013年兩熟制種植面積下降了21.1%，而一熟制種植面積增加了38.7%。總之，華北地區耕地物候的時空變化與種植制度密切相關，同時也受到自然資源和人類活動的共同影響。

[1] 侯學會， 牛錚，高帥，等.基于SPOT-VGT NDVI 時間序列的農牧交錯帶植被物候監測.農業工程學報， 2013, 29(1): 142～150

[2] Cleland E E，Chuine I，Menzel A，et al.Shifting plant phenology in response to global change.Trends in Ecology and Evolution， 2007， 22： 357～365

[3] Weiss J L，Gutzler D S，Coonrod J E A，et al.Seasonal and inter-annual relationships between vegetation and climate in central New Mexico，USA.Journal of Arid Environments， 2004, 57(4): 507～534

[4] 吳文斌， 楊鵬，唐華俊，等．基于NDVI 數據的華北地區耕地物候空間格局研究．中國農業科學， 2009, 42(2): 552～560

[5] 辛景峰， 宇振榮，Driessen P M.利用NOAA NDVI數據集監測冬小麥生育期的研究.遙感學報， 2001, 6(5): 442～447

[6] 張峰， 吳炳方，劉成林，等.利用時序植被指數監測作物物候的方法研究.農業工程學報， 2004, 20(1): 155～159

[7] Wu Wenbin，Yang Peng，Tang Huajun，et al.Characterizing spatial patterns of phenology in cropland of China based on remotely sensed data.Agricultural Sciences in China， 2010, 9(1): 101～112

[8] 張峰， 吳炳方，劉成林，等.區域作物生長過程的遙感提取方法，遙感學報， 2004, 8(6): 515～528

[9] 李正國， 唐華俊，楊鵬，等．植被物候特征的遙感提取與農業應用綜述.中國農業資源與區劃， 2012, 33(5): 20～28

[10]Chuai X W，Huang X J，Wang W J，et al.NDVI，temperature and precipitation changes and their relationships with different vegetation types during 1998～2007 in Inner Mongolia，China.International Journal of Climatology， 2013， 33： 1696～1706

[11]李楊， 江南，侍昊，等.基于SPOT/VGT NDVI 的大區域農作物空間分布.農業工程學報， 2010, 26(12): 242～247

[12]楊婷， 趙文利，王哲怡，等.基于遙感影像 NDVI 數據的中國種植制度分布變化.中國農業科學， 2015, 48(10): 1915～1925

[13]歐立業， 羅烈琴，易明華.基于 NDVI時間序列數據的江西省水稻種植制度變化研究.測繪與空間地理信息， 2016, 39(4): 63～66

[14]Li Zhengguo，Tang Huajun，Yang Peng，et al.Spatio-temporal responses of cropland phenophases to climate change in Northeast China.Journal of Geographical Sciences， 2012, 22(1): 29～45

[15]范錦龍， 吳炳方.復種指數遙感監測方法.遙感學報， 2004, 8(6): 628～636

[16]唐鵬欽， 吳文斌，姚艷敏，等.基于小波變換的華北平原耕地復種指數提取.農業工程學報， 2011, 27(7): 220～225

[17]丁明軍， 陳倩，辛良杰，等.1999～2013年中國耕地復種指數的時空演變格局.地理學報， 2015, 70(7): 1080～1090

[18]李正國， 楊鵬，周清波，等.基于時序植被指數的華北地區作物物候期/種植制度的時空格局特征.生態學報， 2009, 29(11): 6216～6225

[19]李正國， 唐華俊，楊鵬，等．基于時序植被指數的東北地區耕地生長季特征識別與應用研究．北京大學學報：自然科學版， 2011, 47(5): 882～892

[20]衛煒， 吳文斌，李正國，等.基于SPOT/VGT數據的中國北方耕地物候提取研究.中國農業資源與區劃， 2016, 37(4): 77～86

[21]J?nsson P，Eklundh L.Seasonality extraction and noise removal by function fitting to time-series of satellite sensor data，IEEE Transactions of Geoscience and Remote Sensing.2002, 40(8): 1824～1832

[22]J?nsson P，Eklundh L.TIMESAT a program for analyzing time-series of satellite sensor data.Computers & Geosciences， 2004， 30： 833～845

[23]閆慧敏， 曹明奎，劉紀遠，等.基于多時相遙感信息的中國農業種植制度空間格局研究.農業工程學報， 2005, 21(4): 85～90

[24]阿多， 熊凱，趙文吉，等.1960～2013年華北平原氣候變化時空特征及其對太陽活動和大氣環境變化的響應.地理科學， 2016, 36(10): 1555～1564

[25]胡洵， 王靖.華北平原夏玉米各生育階段農業氣候要素變化特征.干旱地區農業研究， 2015, 33(4): 251～258

[26]李正國， 唐華俊，楊鵬，等.東北三省耕地物候期對熱量資源變化的響應.地理學報， 2011, 66(7): 928～939