黃文潔，曾桐瑤，黃曉東

（蘭州大學草地農(nóng)業(yè)科技學院 / 草地農(nóng)業(yè)生態(tài)系統(tǒng)國家重點實驗室，甘肅 蘭州 730020）

物候?qū)W是研究生物圈動植物與氣候、水文、土壤等環(huán)境因子之間相互關(guān)系的科學[1]。植物物候是指植物由于受氣候和環(huán)境因素驅(qū)動而出現(xiàn)的以年為單位的周期現(xiàn)象[2-3]。物候變化直接影響生物生產(chǎn)，并在生態(tài)系統(tǒng)生產(chǎn)力和碳循環(huán)中發(fā)揮關(guān)鍵作用，許多學者形象地認為物候是季節(jié)性和內(nèi)生性的最敏感的綜合指標。近100年來，在變暖已成為事實的背景下[4]，植物物候因其對氣候變化響應的敏感性已經(jīng)成為生態(tài)系統(tǒng)對氣候變化反饋的重要感應器和全球變化中重點研究的對象[5]。

傳統(tǒng)的物候觀測采用野外人工目視觀察法，該方法耗費大量人力物力，并且不易開展大范圍長時間序列的監(jiān)測。遙感觀測具有覆蓋范圍廣、連續(xù)的空間、長時間序列的特點，能夠定量反映植物的季節(jié)生長和發(fā)育及其年際間的年際變化等特征，而成為傳統(tǒng)物候觀測手段的有效補充[6]。越來越多的學者以遙感技術(shù)為手段深入研究植被物候變化，使其逐步發(fā)展為一門獨立的學科并成為全球氣候變化研究的重要組成部分。宋春橋等[7]利用2001-2010年MODIS EVI數(shù)據(jù)對藏北高原典型樣區(qū)物候的時空變化進行了分析。田柳茜[8]選用GIMMS NDVI (1982-2006)、SPOT NDVI (1999-2013)和 MODIS NDVI(2001-2013)數(shù)據(jù)集探討了青藏高原返青期的時空特征。李鵬[9]選擇GIMMS NDVI(1982-2012)、 SPOT NDVI(1999-2012)和MODIS NDVI(2001-2012)3種植被指數(shù)數(shù)據(jù)來研究青藏高原植物枯黃期的時空轉(zhuǎn)變。近年來，關(guān)于草地植被物候時空變化特征在其他地區(qū)已有大量的研究，但針對青藏高原地區(qū)高寒草地植被物候的研究，尚存在地面驗證數(shù)據(jù)較少，時空尺度較小，以及不夠深入等問題。因此，本研究基于地面觀測數(shù)據(jù)以及遙感數(shù)據(jù)，分析青藏高原高寒草地的植被物候特征，以期為后期研究物候?qū)夂虻捻憫龊脺蕚洹?/p>

NDVI植被指數(shù)由于對葉綠素敏感，曲線變化幅度大，可消除部分太陽高度角和觀測條件帶來的噪聲等特點，使其成為物候研究中廣泛而有效的選擇[10]。利用遙感數(shù)據(jù)提取植物物候期主要有閾值法、最大斜率法、滑動平均法和經(jīng)驗公式法等方法[11]。這些數(shù)學方法各有優(yōu)劣，如丁明軍等[12]利用1999-2009年SPOT VGT-NDVI數(shù)據(jù)采用最大比率法對高寒草地物候的時空變化進行了分析。司文才和劉峻明[13]選用2000-2009年NDVI數(shù)據(jù)，采用動態(tài)閾值法研究了河北省冬小麥(Triticum aestivum)的返青期，得出冬小麥返青期由南到北逐漸推遲的結(jié)論。侯學會等[14]利用Logistic模型對NDVI累積曲線進行擬合，依據(jù)曲率極值方法提取植被物候，研究了近10年東北地區(qū)的森林物候。

青藏高原是全球氣候變化響應較為敏感的地區(qū)，其熱液條件處于生物生長極限，生態(tài)系統(tǒng)極其脆弱[15]。本研究針對青藏高原高寒草地，選用2001-2015年MODIS-NDVI數(shù)據(jù)，基于TIMESAT軟件提取植被返青期、枯黃期及生長季長度3個物候指標，分析其時空變化特征，探究其驅(qū)動因子，以期為青藏高原高寒草地物候監(jiān)測打下基礎。

1 研究區(qū)概況

圖1 青藏高原草地類型及牧試臺站空間分布Figure 1 Distributions of the grasslands and climate stations on the Tibetan Plateau

青藏高原位于我國西南部 73°19′-104°47′ E，26°00′-39°47′ N(圖 1)，是世界上海拔最高的高原，有世界“第三極”的美譽[16]。青藏高原長約 2 800 km，寬 300～1 500 km，總面積達 257.24 萬 km2，占我國陸地面積的26.8%[17]。青藏高原區(qū)域地形復雜多樣，地勢落差大，海拔由西北向東南逐漸降低。全年氣溫低、區(qū)域差異大，降水少且分布不均勻，東南多、西北部少，空氣稀薄，具有很強的太陽輻射。由于海拔的原因使得青藏高原地區(qū)自然條件與其他中緯度暖溫帶和亞熱帶地區(qū)大有差異，青藏高原的草地類型主要有高寒草甸、高寒草原、山地草甸等17種。青藏高原的氣候變化(水熱、干濕)間接地影響了青藏高原植被的生長過程，而植被的生長反過來也會影響碳收支、氣候變化、土地覆蓋、農(nóng)作物估產(chǎn)及放牧等經(jīng)濟活動[18]。研究青藏高原物候時空變化特征能夠揭示全球變暖背景下高寒草地生態(tài)碳循環(huán)及其對氣候變化的響應[19]。

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

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

本研究采用2001-2015年青藏高原歸一化植被指 數(shù) (normalized difference vegetation index, NDVI)分析青藏高原高寒草地物候變化趨勢。NDVI能夠很好地反映植被的生物量、初級生產(chǎn)力以及植被覆蓋狀況，在植被生長和趨勢變化的研究中得到廣泛應用[20-21]。本研究所采用的MODIS(moderate resolution imaging spectroradiometer)V006 MOD13A1 NDVI數(shù)據(jù)集，時間為2001-2015年，空間分辨率為500 m，時間分辨率為16 d，每年23期。DEM數(shù)據(jù)從地理空間數(shù)據(jù)云(http://www.gscloud.cn/)中下載得到，空間分辨率為500 m。

研究采用了青海省氣象科學研究所提供的共計20個牧試站點的地面牧草物候觀測資料，用于驗證遙感提取的物候參數(shù)(表1)。站點觀測記錄表明，2003-2015年，青藏高原牧草返青期主要集中在每年4月初至6月中旬(4月10日-6月12日)，返青期在空間上表現(xiàn)出顯著的差異，不同區(qū)域不同草地類型的返青期相差最大達到2個月?？蔹S期主要集中在8月末至10月初(8月29日-10月14)，相比返青期，青藏高原牧草進入枯黃期的時段相對比較集中。

2.2 時間序列重建

16 d 最大合成法 (maximum value composite, MVC)雖然消除了部分云、大氣等對NDVI的影響，但NDVI時間序列仍然存在噪聲，氣溶膠、云層等噪聲很可能影響提取植被物候的準確性[22]。因此，在提取物候參數(shù)之前需要對NDVI數(shù)據(jù)進行平滑處理。TIMESAT 3.2平臺提供了3種時間序列重建的方法：非對稱高斯函數(shù)擬合法(AG)[23]、雙Logistic函數(shù)擬合方法(DL)[24]和Savizky-Glolay過濾算法(S-G)[25]。S-G算法最符合原始NDVI曲線的細節(jié)，在降低噪聲的同時保持曲線的范圍與形狀不發(fā)生偏移[26]，能夠真實地反映出青藏高原NDVI數(shù)據(jù)的周期性變化[27]，因此，本研究選用S-G過濾算法擬合NDVI時間序列數(shù)據(jù)，S-G濾波窗口的值設為4。

表1 牧試臺站經(jīng)緯度與海拔Table 1 Latitude, longitude, and altitude of the pasture test stations

2.3 物候參數(shù)提取

植被的 3 個參數(shù)，生長季開始 (the start of growth season, SOG)、生長季結(jié)束 (the end of growth season,EOG)和生長季長度 (length of growth season, LOG)被用來反映植物的物候變化。 采用儒略日 (Julian days)換算方法，將物候出現(xiàn)日期轉(zhuǎn)化為距1月1日的實際日數(shù)，即年序列累積日數(shù)，得到各物候期的時間序列?；谶b感技術(shù)的SOG是指大多數(shù)植物開始擴大葉片并正常生長的日期，此后植被進入快速生長階段。EOG指大多數(shù)植物不能正常生長并且葉子開始變黃的日期，之后植物進入休眠或死亡階段。LOG等于EOG減去SOG，這與傳統(tǒng)概念上生長季節(jié)是指“一年中植物生長的天數(shù)”不同。目前，針對時間序列的遙感數(shù)據(jù)發(fā)展了不同的物候期提取方法，包括閾值法[28-29]、曲線擬合法[30]、最大變化速率法[31]和經(jīng)驗回歸法[32]等，但未對最佳提取方法形成普遍共識[33]。結(jié)合該區(qū)域的前期研究，本研究采用閾值法提取青藏高原高寒草地植被物候期。其中，當NDVI增加至擬合函數(shù)左半部分振幅20%時定義為植物返青期，當NDVI降低至擬合函數(shù)右半部分振幅60%時定義為植物枯黃期[9]，生長季的長度則為返青期與枯黃期時間之差，據(jù)此利用TIMESAT軟件逐年提取青藏高原高寒草地植被物候參數(shù)。

2.4 趨勢分析方法

為了研究青藏高原高寒草地近15年植被物候時空變化特征，本研究采用趨勢分析法來擬合2001-2015年物候的年際變化[34-35]。年變化率(Slope)大于零表示返青期推遲、枯黃期推遲或生長季延長的趨勢，年變化率小于零表示返青期提前、枯黃期提前或生長季縮短的趨勢。公式如下：

式中：y為第i年的物候；i為年變量，從1年到n年。

3 結(jié)果與分析

3.1 精度驗證

利用20個牧試站提供的青藏高原牧草物候期觀測數(shù)據(jù)作為真值，采用回歸系數(shù)(R2)、平均誤差(Bias)和均方根誤差(RMSE)對MOD13A1產(chǎn)品提取的物候期進行精度驗證。2003-2015年20個驗證站點MOD13A1識別的返青期和枯黃期與地面觀測返青期和枯黃期的精度評價結(jié)果 (圖2)表明，2003-2015年，青藏高原高寒草地牧試站觀測返青期主要集中在第110-160天，而MOD13A1遙感數(shù)據(jù)提取的返青期主要集中在第100-161天。遙感識別的返青期與地面觀測值的R2= 0.61，Bias=-1.23，RMSE= 8.31，誤差超過 10 d 的占 19.6%。2003-2015年，青藏高原高寒草地牧試站觀測枯黃期主要集中在第250-290天，而 MOD13A1 遙感數(shù)據(jù) 提取的枯黃期主要集中在第240-288天。遙感識別的枯黃期與牧試站觀測值的R2= 0.64，Bias=5.64，RMSE=8.39，誤差超過 10 d 的占 22.2%。返青期和枯黃期的驗證結(jié)果表明，基于MOD13A1產(chǎn)品提取的物候參數(shù)可以較為準確地反映青藏高原高寒草地返青期與枯黃期的時空變化特征。

圖2 MOD13A1產(chǎn)品反演與地面觀測的返青期和枯黃期對比Figure 2 Comparison between the phenology observed from the ground stations and phenology based on the MOD13A1 retrieval

3.2 物候期的空間分布

通過柵格處理計算得出青藏高原高寒草地2001-2015年的物候空間分布特征圖(圖3)。青藏高原高寒草地返青期存在明顯的東早西遲的趨勢，從東南到西北，返青期開始時間逐漸推遲。多年平均返青期主要集中在儒略年第110-170天，這與牧試站地面觀測數(shù)據(jù)基本一致。青藏高原高寒草地平均返青期的時間為第137天，東部返青期主要在第130天左右，而西部則主要在第150天左右開始。返青期在高原東南部地勢較低的河谷地區(qū)最早，在緯度或海拔較高的區(qū)域，牧草地返青期較晚。

青藏高原高寒草地枯黃期在空間分布上也存在明顯的差異，由中部向西南和東南方向逐漸推遲，高原中部枯黃期在第260-280天，西部邊緣枯黃期遲于東部，約在10月中下旬枯黃。從枯黃期空間分布特征圖來看，整個青藏高原高寒草地的平均枯黃期時間為第286天左右，草地集中在9月初到10月底進入枯黃期。在返青期和枯黃期的共同作用下，牧草生長季長度自東南至西北逐漸縮短，這一趨勢符合研究區(qū)地勢和氣候的特征。生長季長度集中在第100-170天，平均生長季長度為 149 d。

3.3 物候期的年際變化

通過線性擬合青藏高原2001-2015年高寒草地物候期的年際變化趨勢(圖4)。返青期的趨勢并不是固定的，在不同的時間段內(nèi)趨勢不盡相同。2001-2009年間，返青期有明顯的提前趨勢，返青期提前 18 d；2009-2015 年間，返青期推遲 4 d?？傮w而言，青藏高原高寒草地的返青期呈現(xiàn)出提前的趨勢，2009年草地的返青期最早。青藏高原高寒草地植被枯黃期呈現(xiàn)波動提前的趨勢，這一趨勢較之返青期而言并不明顯，2002、2008和2013年牧草進入枯黃期較早，2001-2015年間青藏高原高寒草地植被枯黃期整體提前5 d。將2001-2015共15年分為兩個時間段，在2001-2008年期間枯黃期以每年0.857 d的速度提前，在隨后的2008-2015年又出現(xiàn)了推遲的趨勢，每年推遲0.286 d(圖4)。從整個研究期來看，2001-2015年枯黃期呈現(xiàn)出微弱的提前趨勢。在青藏高原高寒草地返青期和枯黃期年際變化的共同影響下，青藏高原高寒草地植被生長期呈現(xiàn)波動延長的趨勢，平均生長季長度延長了9 d，其中2007年的平均生長季長度最長，達到158 d。

圖4 2001-2015年青藏高原高寒草地植被物候年際變化Figure 4 Inter-annual variation of the vegetation phenology in the alpine grasslands of the Tibetan Plateau from 2001 to 2015

選取青藏高原山地草甸類、高寒草原類、高寒草甸類3種典型的高寒草地類型，分析其物候期差異。青藏高原不同草地類型的物候差異顯著(圖5)，結(jié)果表明，山地草甸類植物返青期較早，返青期集中在第124-134天；高寒草甸類的植物返青期次之，集中在第132-141天；高寒草原類的植物返青期最晚，在第142-152天。不同草地類型的枯黃期相比返青期時間相對比較集中，其中高寒草原類與高寒草甸類在第270-282天進入枯黃期，山地草甸類在第278-288天進入枯黃期。在返青期與枯黃期的共同影響下，山地草甸類的生長季長度最長，范圍在145～162 d，高寒草原類的生長季長度最短，范圍在123～140 d，高寒草甸類草地的生長季長度在133～148 d。2001-2015年間，不同草地類型的返青期都有提前的趨勢，但趨勢并不顯著，而枯黃期高寒草原類與高寒草甸類表現(xiàn)出一定的提前趨勢，山地草甸類基本沒有表現(xiàn)明顯的提前或推遲趨勢。

圖5 青藏高原不同草地類型物候年際變化特征Figure 5 Inter-annual variation of the vegetation phenology in the different grassland types of the Tibetan Plateau

2001-2015年，青藏高原高寒草地物候年際變化在空間分布上也呈現(xiàn)出顯著的差異(圖6)。牧草返青期中呈提前趨勢的占71.71%，呈推遲趨勢的占28.29%，在青藏高原高寒草地東部和西北地區(qū)返青期呈現(xiàn)出推遲的趨勢。其中，返青期顯著提前的占5.12%，顯著推遲的占2.11%。牧草枯黃期中呈提前趨勢的占65.18%，呈推遲趨勢的占34.82%，枯黃期提前的區(qū)域集中在青藏高原高寒草地中部?？蔹S期顯著提前的占4.05%，顯著推遲的占2.81%。從牧草生長季長度變化來看，具有增長趨勢的占57.81%，具有縮短趨勢的占42.19%，整體呈現(xiàn)延長趨勢。生長季長度中約2.67%的區(qū)域顯著提前，約6.15%的區(qū)域顯著推遲。

3.4 物候期的海拔差異

以100 m間隔為海拔梯度提取不同草地類型的物候參數(shù)，對物候與海拔建立線性擬合關(guān)系，反映出2001-2015年青藏高原高寒草地物候受海拔因素影響的趨勢。結(jié)果表明，青藏高原2001-2015年高寒草地植被物候受海拔因素影響(圖7)。在海拔2 000-3 500 m，返青期呈現(xiàn)不規(guī)則波動，海拔對不同草地類型的返青基本沒有任何影響；在3 500-5 100 m，返青期隨海拔的升高而推遲，在 5 100 m以上，隨著海拔的升高返青期反而提前?？蔹S期以5 200 m為界，出現(xiàn)兩種不同的趨勢，在2 000-5 200 m范圍內(nèi)，隨著海拔升高枯黃期呈現(xiàn)波動提前的趨勢；在5 200 m以上，枯黃隨著海拔的升高而推遲。對于生長季長度而言，以3 500和5 000 m為節(jié)點，小于3 500 m生長季長度變化不規(guī)則，海拔在 3 500-5 000 m，生長季長度縮短。海拔大于5 000 m 時，生長季長度反而延長。在海拔 5 000 m以上出現(xiàn)了植物物候特征與低海拔相反的變化趨勢，可能是由于受積雪覆蓋地表的影響，造成草地植被物候提取出現(xiàn)偏差。

綜上所述，海拔因素很大程度上影響了青藏高原高寒草地物候期。當海拔在3 500-5 000 m時，海拔每升高500 m，牧草返青期推遲4.3 d，枯黃期提前 1.3 d，生長季長度縮短 5.6 d。在海拔 3 500 m以下，牧草返青期和生長季長度與海拔并無明顯的相關(guān)關(guān)系(圖7)，這可能是受到人類活動的影響[13]。在海拔5 000 m以上，植被物候隨海拔變化的關(guān)系呈現(xiàn)出相反的趨勢，因為在5 000 m以上的區(qū)域，主要以寒漠或苔原為主，該區(qū)域植被生命活動短暫，對物候參數(shù)進行提取基本是不可行的，造成分析結(jié)果可信度較差。

青藏高原不同草地類型對海拔變化的響應有所不同，本研究選取了高寒草甸類、高寒草原類、山地草甸類3種草地類型來分析其在不同海拔高度植物的返青期、枯黃期和生長季長度的變化。分析結(jié)果表明，在 3 500-5 000 m 海拔區(qū)間，3 種草地類型的返青期都隨著海拔的升高而推遲。其中高寒草原類的推遲幅度最大，達到9.41 d。枯黃期隨海拔變化的趨勢不明顯，波動較大，山地草甸類與高寒草甸類與海拔呈負相關(guān)的關(guān)系，說明隨著海拔的升高，枯黃期逐漸提前，但高寒草原類隨著海拔的升高卻呈現(xiàn)枯黃期推遲的趨勢，R2達到 0.513 6。生長季長度在海拔 3 000-5 000 m 區(qū)間，都表現(xiàn)出隨著海拔的升高而縮短的趨勢，且趨勢顯著(圖8)。

4 討論與結(jié)論

圖7 2001-2015年青藏高原高寒草地物候期與海拔的關(guān)系Figure 7 Relationship between the phenology and altitude of the alpine grasslands of the Tibetan Plateau from 2001 to 2015

圖8 青藏高原不同草地類型物候期與海拔的關(guān)系Figure 8 Relationship between the phenology and altitude of the different grassland types of the Tibetan Plateau

青藏高原高寒草地植被物候的變化趨勢不是單純的線性變化，在不同的時間段，有不同的趨勢，呈現(xiàn)出的變化非常復雜。因此，推測這是由于物候提取方法不同，即物候期是由多種因素共同決定導致的，如水熱因素和人為因素等。丁明軍等[12]認為，1999-2009年間，青藏高原高寒草地返青期整體上呈提前趨勢，枯黃期呈推遲趨勢，生長期長度延長，在重合的時間段本研究也得出相同的結(jié)果。田柳茜[8]選用GIMMS (1982-2006)、SPOT (1999-2013)、MODIS (2001-2013)3種NDVI數(shù)據(jù)提取青藏高原植被返青期，盡管在絕對數(shù)值上不相同，均得出返青期提前的趨勢。馬曉芳等[36]選用1982-2005年GIMMS NDVI數(shù)據(jù)提取青藏高原高寒草地物候，得出返青期和枯黃期逐漸提前，生長季長度縮短的結(jié)論。呂燦賓[37]認為1980-2010年間生長季結(jié)束時間呈提前趨勢。受到時間尺度和不同遙感數(shù)據(jù)源精度的影響各研究得出的物候變化存在一定的差異，整體上，本研究結(jié)論與前人的研究基本吻合[7-9]，但在分析物候期與海拔之間關(guān)系方面，本研究得出了一些更有意義的結(jié)論。物候期并非隨著海拔的增加而呈線性變化，以5 000 m為分界點，在海拔3 500-5 000 m 間，草地物候與海拔之間有明顯的相關(guān)性，隨著海拔的升高牧草返青期推遲，枯黃期提前，生長期長度縮短。當海拔大于5 000 m之后，寒漠或苔原受極低溫度的影響，不具有明顯的物候期[14]。

不同的研究方法將影響草地關(guān)鍵物候參數(shù)的提取，本研究采用閾值法，提取青藏高原高寒草地植被物候期。由于遙感數(shù)據(jù)分辨率較低和存在云層、大氣的噪聲等問題，并且數(shù)據(jù)預處理方法與物候參數(shù)提取方法都將影響研究結(jié)果，所以很有必要結(jié)合地面實測數(shù)據(jù)對遙感數(shù)據(jù)提取的物候期進行精度評價。例如，Chen等[38]利用中國科學院的地面實測數(shù)據(jù)對1982-1993年中國東部溫帶植被物候期進行了驗證。本研究采用2003-2015年的牧試站物候觀測值對遙感數(shù)據(jù)進行均方根誤差與偏差評價。驗證結(jié)果表明，地面觀測數(shù)據(jù)與遙感監(jiān)測結(jié)果之間的誤差在可接受范圍內(nèi)，本研究提取的青藏高原高寒草地植被物候期參數(shù)的精度較高，可以反映青藏高原的草地植被物候的時空變化。

本研究利用MOD13A1產(chǎn)品分析了青藏高原高寒草地2001-2015年植被物候時空變化特征并探究了海拔對物候期影響。結(jié)果表明：1)青藏高原高寒草地物候在空間上存在較大差異性。由于受地形抬升和水熱條件惡化的影響，呈現(xiàn)出由東南向西北返青期逐漸推遲，枯黃期逐漸提前，生長季長度逐漸縮短的趨勢。2)不同草地類型物候期存在明顯差異，山地草甸類植物返青期較早，高寒草原類植物返青期較晚。高寒草甸草原類的枯黃期最早，山地草甸類的枯黃期最遲。3)青藏高原高寒草地植被物候年際變化的趨勢表現(xiàn)為返青期和枯黃期提前，整體上，生長季長度增長。具體分析結(jié)果表明，2001-2015年間，植被返青期提前了 14 d，枯黃期提前 5 d，生長季長度延長了 9 d。4)青藏高原高寒草地物候空間異質(zhì)性與海拔關(guān)系密切。存在著 3 500 和 5 000 m 的分界線，3 500 m以下物候隨海拔變化規(guī)律不明顯，在3 500-5 000 m之間，二者關(guān)系顯著，隨著海拔的升高，返青期推遲，枯黃期提前，生長季長度縮短。