基于NDVI時間序列的近12年中國植被覆蓋的單調趨勢分析

李 強，張 翀

（1.陜西師范大學 旅游與環境學院，陜西 西安710062；2.陜西學前師范學院 環境與資源管理系，陜西 西安710100）

植被作為陸地生態圈的重要組成，是氣候系統的重要元素［1］。眾所周知，全球氣候正處于一個持續變暖的階段，強烈影響陸地生物圈［2］，在這樣的背景下，掌握陸地植被覆蓋度年際間的變化規律，對評價陸地生態系統的環境質量、調節生態過程具有重要的理論和實際意義［3］。歸一化植被指數（normalized differ－ence vegetation index，NDVI）由紅波段與近紅外波段的反射率計算而來［4］，和植物的生產力密切相關［5］，并且NDVI趨勢可以用來衡量植被覆蓋的改善與退化［6］。NDVI趨勢被用在很多測算中，包括全球變暖的生態響應［7］、物候變化［8］、作物狀況［9］、土地覆蓋變化［10］以及沙漠化［11］。植被覆蓋的趨勢及年際變化會影響植被與大氣之間能量的交換［12］。一系列研究表明，在全球，特別是北半球，存在生長季的開始時刻提前以及生長季長度增大的趨勢［13］，生長季長度增大與氣溫的上升會加快地表水分蒸發，增大干旱脅迫與林火發生的可能［14］，并增加了固碳強度［15］。因此，生長季內植被覆蓋相對于全年更能作為反映諸如土壤退化等的指示器。利用NDVI序列進行趨勢分析時，由于數據集存在自相關會影響模型假設，線性模型應慎重使用。所以既要剔除序列的季節性，也需要進行非參數趨勢檢驗［16］。本文基于時間序列諧波分析法、線性趨勢和Kendall’sτ趨勢等方法對中國植被生長季起止時刻進行了確定，并計算得到中國生長季植被覆蓋的單調性趨勢以及生長季長度趨勢，并進一步解釋不同植被類型植被覆蓋趨勢的差異以及生長季變化規律，以期為中國當前的生態建設、修復提供有用的空間信息和理論支撐。

1 數據與方法

1.1 數據來源

本文所用的資料包括黃土高原1999—2010年SPOT VEGETATION旬值NDVI數據、131個臺站12a（1999—2010年）的年降水資料和年平均氣溫資料。NDVI數據來自于互聯網（http：∥free.vgt.vito.be／home.php），空間分辨率為1 000m。植被類型數據來自寒區旱區科學數據中心。本文對旬值NDVI數據經過時間序列諧波分析重構其旬值序列，并得到旬值異常值序列，進而對逐年異常值序列求均值得到年際異常值序列；利用NDVI旬值重構序列提取植被物候特征，確定植被生長季始末期，從而計算得到逐年生長季內NDVI的均值時間序列。

1.2 NDVI時間序列的諧波分析

時間序列諧波分析法（harmonic analysis of time series，HANTS）是平滑和濾波兩種方法的綜合，它能夠充分利用遙感圖像存在時間性和空間性的特點，將其空間上的分布規律和時間上的變化規律聯系起來。時間序列諧波分解法進行影像重構時充分考慮了植被生長周期性和數據本身的雙重特點，能夠用代表不同生長周期的植被頻率曲線重新構建時序NDVI影像，真實反映植被的周期性變化規律。剔除季節影響幾乎可以完全消除數據的季節性特征。剔除季節影響后，一期與另一期的對比將更有意義，而且可以幫助確定是否存在趨勢。HANTS分析中單獨年份的諧波分析與整個時間段的諧波分析之間的差異即為異常值序列，異常值序列剔除了數據的季節性特征［17］。

1.3 物候特征提取

多種處理方法被用來描述NDVI時間序列的生長季的起始時刻（start of growing season，SOS），如半極大值［18］、10% 閾 值 法［19］、拐 點 法［20］、最 大 曲 率法［21］、滑動 平 均 法［22］。 本 文 按 照 Whitetffu 等［8］的方法，采用HANTS平滑NDVI的一階導數，一階導數的最大值對應時刻為SOS（最大NDVI增長），生長季的結束時刻（end of growing season，EOS）為SOS過后首次NDVI下降到與SOS相等的時刻。該方法與其他方法比較相對可靠。本文主要采用各年NDVI年內變化的三階導數，但是這種方法只適合于一年一熟作物區［17］，所以對于農耕區的一年多熟區域采用10%閾值法進行計算［19］。求得中國區域范圍內每個柵格的每年的SOS（一階導數的極大值）與EOS（一階導數的極小值），SOS至EOS即為生長季，兩者之間的時間長度即為生長季長度（length of growing season，LOS）。

1.4 趨勢分析

1.4.1 線性趨勢 運用線性趨勢線分析植被覆蓋的變化趨勢，即以時間為自變量，NDVI為因變最，利用最小二乘法計算斜率值K。K值的符號反映上升或下降的變化趨勢，K＜0表示在計算時段內呈下降趨勢，K＞0表示呈上升趨勢。K值絕對值的大小可以度量其演變趨勢上升、下降的程度［23］。

1.4.2 Kendall’sτ趨勢 Mann等［24］首次建議用Mann—Kendall法來檢驗時間趨勢Kendall’sτ的顯著性，該方法已經被應用在季節數據上，主要是水文分析方面，近年應用在 NDVI數據上［6，25］。為了更好顯示出Kendall’sτ的空間分布，本文將其分為8個等級：極顯著持續退化（－1～－0.75）、較顯著持續退化（－0.75～－0.5）、顯著持續退化（－0.5～－0.25）、不顯著持續退化（－0.25～0）、不顯著持續改善（0～0.25）、顯著持續改善（0.25～0.5）、較顯著持續改善（0.5～0.75）、極顯著持續改善（0.75～1）。

2 結果與分析

2.1 剔除季節性的NDVI趨勢

圖1為中國范圍內內剔除季節性NDVI數據年際變化曲線，植被覆蓋呈明顯增加趨勢（0.03／10a），低于中國 1998—2007年 植 被 覆 蓋 速 度 0.048／10a［26］。（1）1999—2000年中國植被覆蓋減小，2000—2004年植被覆蓋增加，2000年植被覆蓋最低，2001年開始上升，主要是由于2000和2001年屬特大干旱年，且2000年干旱程度較2001年嚴重，該時間中國大范圍降水 偏 少，發 生 了 建 國 以 來 最 嚴 重 的 旱 災［27－28］；（2）2004年植被覆蓋達到峰值，隨之在2005年下降，2005年以后植被覆蓋呈穩定上升趨勢。

圖1 中國1999－2010年剔除季節性NDVI數據的年際變化

1999—2010年剔除季節性旬NDVI數據小波變化系數如圖2所示，根據其規則交替，顯示出在80～110頻率尺度上存在5～7a（180～252旬）的周期，在1～30頻率尺度上存在1～2a（36～72旬）的周期。圖1也顯示出中國植被覆蓋變化呈現出一定的周期性特點：4a的穩定增長期（2000—2004與2005—2009年）、1a的突然下降期（1999—2000，2004—2005與2009—2010年），即存在5a的周期。

圖2 中國1999－2010年剔除季節性NDVI旬數據的小波變換系數

利用線性趨勢原理，得到中國剔除季節性NDVI的線性趨勢空間分布（如附圖8所示）。植被明顯改善區集中在黃土高原—壩上高原以南，橫斷山脈—黃土高原西緣以東的大面積區域，另外新疆綠洲區也呈改善趨勢；明顯退化區，在中國東部及東南部，主要分布于黃河三角洲平原、長江三角洲平原、珠江三角洲平原與騰格里沙漠以南—隴中高原以北的區域，其中黃三角植被結構簡單，生態系統年輕化的特點，并且是國內外少有的資源富集區，開發潛力大，脆弱的生態環境在近12a來經濟快速發展中呈退化趨勢，而長江三角與珠江三角以及騰格里沙漠—隴中高原之間的區域，主要是由于快速經濟發展、開發利用以及城市化等要素導致植被覆蓋的下降，樸世龍等［29］發現長三角與珠三角在過去18a間也呈下降趨勢。另外在烏蘭察布高原—渾善達克沙地南緣—錫林郭勒高原東南部一線的條帶地區退化較為明顯，與穆少杰等［3］研究結果一致。

圖3反映出，中國各種植被類型均呈改善趨勢。（1）常綠針葉林、常綠闊葉林、稀疏灌叢、低山草原、耕地為明顯改善區，趨勢值在0.004～0.007，年際波動性最強。常綠針葉林與常綠闊葉林主要分布在中國華東南部、華南和西南地區的山地丘陵區，尤以橫斷山脈、大巴山改善明顯；稀疏灌叢分布在兩廣南部的地勢平坦區；低山草原主要分布在黃土高原及其周邊低山地區；耕地分布較為廣泛，但陜北、河南—安徽的黃淮平原的植被改善趨勢明顯高于其他地區，其中陜北自1999年退耕還林還草開始以來，林草替換原有耕地，植被覆蓋明顯增加［30］，黃淮平原近50a來氣溫變化不明顯，降水呈增加趨勢［31］，降水量的增加促使該地區植被覆蓋的增加。（2）落葉針葉林、落葉闊葉林、密集灌叢、高山亞高山牧場草地、典型草原、牧場草地和沼澤濕地的改善趨勢次之，趨勢值在0.001～0.004，波動性相對較弱。落葉針葉林集中分布在大小興安嶺以及新疆的阿爾泰山與天山山脈；落葉闊葉林分布于秦嶺北麓、小興安嶺以及長白山脈與壩上高原；密集灌叢分布于青藏高原東部、云貴高原、羅霄山等山地區。（3）海岸濕地、荒漠草原與溫帶高山—亞高山高寒草原改善趨勢不明顯，趨勢值在0～0.001，波動性最小，為生態脆弱區。

綜上所述，對于趨勢及波動性：常綠林＞落葉林、稀疏灌叢＞密集灌叢、低山草原＞高山亞高山牧場草地、典型草原與牧場草地＞荒漠草原與溫帶高山—亞高山高寒草原、耕地＞沼澤濕地＞海岸濕地。反映出中國植被覆蓋受緯度地帶性影響較為明顯。緯度較低，生物多樣性程度較高，則植被覆蓋改善趨勢明顯，相反則生物多樣性程度下降，改善趨勢次之。輕微改善的地區多為生態環境脆弱區，其生物多樣性很低，易受人類活動等因素的影響，趨勢值處于0值附近。

2.2 生長季NDVI的單調趨勢

附圖8b為中國1999—2010年生長季植被覆蓋的Kendall’sτ趨勢空間分布。Kendall’sτ趨勢平均值為0.2873，持續改善的面積比重為85.85%，持續退化的比重占14.15%，說明中國生長季植被覆蓋總體上呈持續改善趨勢。

空間分布上（附圖8b）：（1）中國生長季植被覆蓋顯著、較顯著與極顯著持續退化區比重為3.91%，主要分布于長江三角、珠江三角、騰格里沙漠—隴中高原之間的區域、柴達木盆地東南部與烏蘭察布高原—渾善達克沙地南緣—錫林郭勒高原東南部一線的條帶地區［3］，植被覆蓋退化區還零散分布在青藏高原整個地區，尤以南部為主的區域。（2）不顯著持續性變化區為易變化區域，比重為38.25%，集中分布與新疆—內蒙古的沙漠區，該地區NDVI小于0.1，變化可能與到達傳感器的能量的微小差異引起，這種微小差異與傳感器本身以及大氣的影響有關；四川盆地、東南丘陵以及退化區周邊地區為不顯著持續性變化區。（3）持續改善區分布范圍遍及全國，持續性較高的區域（極顯著與較顯著區域）主要分布在105°E以東，30°—40°N 的區域，比重為23.52%，該區域北界大致為中溫帶與暖溫帶過渡處，南界大致為中亞熱帶與北亞熱帶的過渡處，西界至高原溫帶邊緣；另外在河西走廊、橫斷山區、新疆綠洲區以及青藏高原東北部也有分布。

由圖3b發現，各種植被類型的Kendall’sτ趨勢與線性趨勢規律具有很強的相似性，即緯度地帶性明顯，低緯度地區植被多樣性程度高，持續性強，反之持續性弱，從而也可以表明趨勢的持續性與趨勢的大小存在很強的相關性，即趨勢高則持續性大，反之持續性弱。

圖3 中國不同植被類型的統計值注：1.落葉針葉林；2.常綠針葉林；3.常綠闊葉林；4.落葉闊葉林；5.密集灌叢；6.稀疏灌叢；7.海岸濕地；8.高山—亞高山牧場草地；9.低山草原；10.典型草原；11.荒漠草原；12.牧場草地；13：沼澤濕地；14.耕地；15.溫帶高山—亞高山高寒草原

2.3 生長季長度趨勢

受物候變化影響的生長季Kendall’sτ趨勢由SOS與LOS的變化決定。如果植物生長季穩定，則年際SOS與EOS基本不變。附圖8c表明LOS的線性趨勢與其變異系數的空間變化，中國植被生長季長度平均趨勢為0.763 6，總體上呈增加趨勢。

附圖8c—d顯示出，LOS趨勢值主要集中在－1～3.5，減小與增加趨勢的比重相當，生長季長度減小的區域比重占49.88%。（1）連片減小區域主要分布在青海東部及南部與天山山脈，波動性較?。▓D3d）；陜北高原—呂梁山—太行山以及壩上高原大興安嶺兩側、小興安嶺、長白山周圍的低山區，波動性很小（圖3d）；內蒙古的呼倫貝爾高原波動性較?。▓D3d），以及山東丘陵，主要分布著耕地，波動性較高（圖3d），其生長季長度也呈減小趨勢。（2）生長季長度增加的比重為50.12%，連片增加區域主要分布在高度較高的山地區，包括阿爾泰山脈、秦巴山區、大小興安嶺與長白山脈，該地區除常綠林以外的其他植被類型的生態環境較為脆弱，波動性很高，值大于1.5（圖3d）。

SOS與EOS的年際變化以及線性趨勢顯示，SOS呈減小趨勢，而EOS成增加趨勢，并且EOS趨勢的增大程度大于SOS減小程度（0.589 6＞0.021 1），SOS與EOS的年際變化特征基本相反，峰值對應谷值，即SOS提前—EOS推遲或SOS推遲—EOS提前，并且后者情況只有在2007年表現明顯，所以中國LOS總體上的增長趨勢由SOS提前和EOS推遲造成，主要是因為在氣候增溫［30］的趨勢下，春季氣溫會提前達到植被生長所需的適宜溫度，而在秋季，氣溫的下降趨勢得到延緩。

3 結論

（1）中國植被覆蓋呈明顯增加趨勢，各種植被類型年際變化規律基本相同。中國植被覆蓋變化呈現出一定的周期性特點：4a的穩定增長期，1a的突然下降期。

（2）中國生長季植被覆蓋Kendall’sτ趨勢平均值為0.287 3，總體上呈持續改善趨勢，持續改善的面積比重為85.85%，持續退化的比重占14.15%，持續性較高的改善區主要分布在105°E以東，30°—40°N之間的區域，比重為23.52%，該區域北界大致為中溫帶與暖溫帶過渡處，南界大致為中亞熱帶與北亞熱帶的過渡處，西界至高原溫帶邊緣。

（3）中國植被覆蓋受緯度地帶性影響較為明顯，緯度低，生物多樣性程度較高，則植被覆蓋改善趨勢明顯，持續性強，反之改善趨勢及持續性較低；黃三角和風沙區等生態脆弱區，長江三角、珠江三角與蘭州市周邊地區等經濟快速發展區，植被覆蓋退化嚴重。

（4）中國植被生長季長度平均趨勢為0.763 6，總體上呈增加趨勢，主要由EOS的增大趨勢引起。本文之所以采用2000年土地覆蓋數據，就是因為退耕還林從1999年開始，而1999—2010年耕地的植被覆蓋增大趨勢明顯，而且呈持續性相對較高，這意味著在此期間有很大面積的耕地轉化為了林地［30］；從1978年“三北防護林”建設以來，西北、華北北部、東北西部植被覆蓋的改善趨勢及持續性較弱，且退化區面積仍然較大，但是總體上是改善的，與王強等［32］的結論非常接近，并且沙漠以及荒漠草原均呈持續改善趨勢，表明1999—2010年沙漠與荒漠草原均有轉化為植被覆蓋度較高的林草地的區域，植樹種草、禁牧輪牧和防沙治沙生態恢復措施的廣泛實施起到了重要作用［3］。

［1］ Foley J A，Levis S，Costa M H，et al.Incorporating dynamic vegetation cover within global climate models［J］.Ecological Applications，2000，10（6）：1620－1632.

［2］ IPCC.4th Assessment report of the intergovernmental panel on climate change［M］.Synthesis Report，2007：52.

［3］ 穆少杰，李建龍，陳奕兆.2001—2010年內蒙古植被覆蓋度時空變化特征［J］.地理學報，2012，67（9）：1255－1268.

［4］ Tucker C J.Red and photographic infrared linear combinations for monitoring vegetation［J］.Remote Sensing of Environment，1979（8）：27－150.

［5］ Prince S D，Tucker C J.Satellite remote sensing of rangelands in Botswana（Ⅱ）：NOAA AVHRR and herbaceous vegetation［J］.International Journal of Remote Sensing，1986，7（11）：1555－1570.

［6］ Alcaraz－Segura D，Chuvieco E，Epstein H E，et al.Debating the greening vs.browning of the North American boreal forest：Differences between satellite datasets［J］.Global Change Biology，2009，16（2）：760－770.

［7］ Pettorelli N，Vik J O，Mysterud A，et al.Using the satellite－derived NDVI to assess ecological responses to environmental change［J］.Trends in Ecology ＆ Evolution，2005，20（9）：503－510.

［8］ White M A，de Beurs K M，Didan K，et al.Intercomparison，interpretation，and assessment of spring phenology in North America estimated from remote sensing for 1982—2006［J］.Global Change Biology，2009，15（10）：2335－2359.

［9］ Tottrup C，Rasmussen M S.Mapping long－term changes in savannah crop productivity in Senegal through trend analysis of time series of remote sensing data［J］.Agriculture，Ecosystems ＆ Environments，2004，103（3）：545－560.

［10］ Hüttich C，Herold M，Schmullius C，et al.Indicators of Northern Eurasia’s land－cover change trends from SPOT－VEG ETATION time－series analysis 1998—2005［J］.International Journal of Remote Sensing，2007，248（18）：4199－4206.

［11］ Symeonakis E，Drake N.Monitoring desertification and land degradation over sub－Saharan Africa［J］.International Journal of Remote Sensing，2004，25（3）：573－592.

［12］ Baldocchi D，Falge E，Gu L，et al.FLUXNET：A new tool to study the temporal and spatial variability of ecosystem－scale carbon dioxide，water vapor，and energy flux densities［J］.Bulletin of the American Meteorological Society，2001，82（11）：2415－2434.

［13］ Sparks T H，Aasa A，Huber K，et al.Changes and patterns in biologically relevant temperatures in Europe 1941—2000［J］.Climate Research，2009，39（3）：191－207.

［14］ Zhang Ke，Kimball J S，Mu Qiaozhen，et al.Satellite based analysis of northern ET trends and associated changes in the regional water balance from 1983to 2005［J］.Journal of Hydrology，2009，379（1／2）：92－110.

［15］ White M A，Running S W，Thornton P E.The impact of growing－season length variability on carbon assimila－tion and evapotranspiration over 88years in the eastern US deciduous forest［J］.International Journal of Biometeorology，1999，42（3）：139－145.

［16］ De Beurs K M，Henebry G M.Trend analysis of the Path finder AVHRR Land（PAL）NDVI data for the deserts of central Asia［J］.Geoscience and Remote Sensing Letters，IEEE，2004，1（4）：282－286.

［17］ Rogier de Jong，Sytze de Bruin，Allard de Wit.Analysis of monotonic greening and browning trends from global NDVI time－series［J］.Remote Sensing of Environment，2011，115（2）：692－702.

［18］ White M A，Thornton P E，Running S W.A continental phenology model for monitoring vegetation responses to interannual climatic variability［J］.Global Biogeochemical Cycles，1997，11（2）：217－234.

［19］ Jones M O，Kimball J S，Jones L A，Satellite passive microwave detection of North America start of season［J］.Remote Sensing of Environment，2012，123：324－333.

［20］ Moulin S，Kergoat L，Viovy N，et al.Global－scale assessment of vegetation phenology using NOAA／AV HRR satellite measurements［J］.Journal of Climate，1997，10（6）：1154－1170.

［21］ Zhang Xiaoyang，Friedl M A，Schaaf C B，et al.Monitoring vegetation phenology using MODIS［J］.Remote Sensing of Environment，2003，84（3）：471－475.

［22］ Reed B C，White M，Brown J F.Remote sensing phenology［M］∥Schwartz M D.Phenology：An Integrative Environmental Science.Dordrecht，The Nether－lands：Kluwer Academic Publishing，2003.

［23］ 徐建華.現代地理學中的數學方法［M］.北京：高等教育出版社，2002.

［24］ Hirsch R M，Slack J R，Smith R A.Techniques of trend analysis for monthly water quality data［J］.Water Resources Research，1982，18（1）：107－121.

［25］ Alcaraz－Segura D，Liras E，Tabik S，et al.Evaluating the consistency of the 1982—1999NDVI trends in the Iberian Peninsula across four time－series derived from the AVHRR Sensor：LTDR，GIMM S，FASIR，and PAL（Ⅱ）［J］.Sensors，2010，10（2）：1291－1314.

［26］ 邱海軍，曹明明.基于SPOT VEGETATION數據的中國植被覆蓋時空變化分析［J］.資源科學，2011，33（2）：335－340.

［27］ 程殿龍.2000年旱災與抗旱工作［J］.防汛與抗旱，2001（1）：30－35.

［28］ 民政部救災處.2001年全國自然災害和救災工作情況［J］.中國減災，2002（1）：30－33.

［29］ 樸世龍，方精云.最近18年來中國植被覆蓋的動態變化［J］.第四紀研究，2001，21（4）：294－302.

［30］ 李雙雙，延軍平，萬佳.近10年陜甘寧黃土高原區植被覆蓋時空變化特征［J］.地理學報，2012，67（7）：960－970.

［31］ 何太蓉，莊紅娟，劉存東.秦嶺—黃淮平原交界帶中東部近50年氣候變化特征與趨勢［J］.安徽農業科學，2009，37（14）：6532－6534.

［32］ 王強，張勃，戴聲佩.基于GIMMS AVHRR NDVI數據的三北防護林工程區植被覆蓋動態變化［J］.資源科學，2011，33（8）：1613－1620.