基于連續時間序列NDVI數據的中國生態狀況時空變化特征

宗加權, 白淑英, 馮朝陽, 蒲 陽

(1.南京信息工程大學, 南京 210044 ; 2.中國環境科學研究院,國家環境保護區域生態過程與功能評估重點實驗室, 北京 100012)

植被是陸地生態系統中的重要組成部分，在陸地表層能量交換、水分循環和生物地球化學循環過程中扮演著重要角色[1]，植被覆蓋度變化是生態環境變化的直接結果，很大程度上表示了生態環境總體狀況[2]。相對于傳統實地調查植被覆蓋特征費時費力，通過遙感手段獲取數據，優勢有空間范圍覆蓋面積大、時間序列長、方便獲取等。歸一化差值植被指數(NDVI)是目前常用的研究地表植被覆蓋特征的重要指標之一[3]，長時間序列監測植被覆蓋特征變化對更好地模擬陸地生態系統的動態變化特征具有重要意義[4]。AVHRR NDVI是目前時間序列最長的全球連續數據集，具有很高的科研價值和實際意義[5-6];MODIS NDVI是從2000年開始的數據集，提高了葉綠素敏感度和空間分辨率、排除了大氣水汽的干擾、改善了合成方法，但時間序列長度有限。探尋兩者間的關系，通過數據擬合構建長時間序列的NDVI數據集，是實現對植被長期觀測的關鍵[7-8]。一些學者基于不同NDVI數據集重構NDVI數據，并結合氣象數據等分析了植被狀況的時空變化特征，如:李凈等[9]基于重構后的數據產品和植被類型數據集，分析了不同時期植被NDVI變化特征，發現西北地區NDVI整體呈上升趨勢;李飛等[10]根據MODIS NDVI與AVHRR數據構建了較長時間序列NDVI數據集，采用相關性分析法分析了全國農業、森林、草地與稀疏植被區植被活動的變化特征，結果指出近30 a來中國陸地植被活動整體增強，這與方精云等[11]利用AVHRR NDVI研究中國地區植被活動的結果相同。趙玉萍等[12]利用1982—2003年GIMMS NDVI數據集和逐月氣象資料，用相關分析法研究了藏北高原草地生態系統NDVI與氣候因子的相關性及滯后性，得到NDVI與月平均氣溫等呈高度正相關(p<0.001)。劉可等[13]利用1982—2012年GIMMS NDVI3g和中國陸地生態系統類型數據，采用一元線性經驗模態分解和相關分析方法，研究了近30 a中國各生態系統NDVI的時空變化特征，結果表明近30 a中國植被活動顯著上升，農田、森林、草地和水體與濕地生態系統的NDVI總體為非穩定上升趨勢。中國陸地生態系統中的森林、草原、荒漠、水體與濕地和其他生態系統類型是對土地利用/覆蓋的分類判別匯總，對生態系統狀況有指示意義。水體與濕地在生態環境評價中起到重要作用，其水體與濕地狀況好壞和面積變化影響著生態環境狀況，水域面積的減少會導致生態環境惡化[14]。水體和濕地因為水源豐富、土壤肥沃，很容易被開墾為良田和進行水產的養殖，水體與濕地生態系統因此被嚴重破壞[15]，牛振國等[16]研究表明，近30 a，中國濕地類型由濕地轉化為非濕地主要是氣候變化和農業活動引起的。

以上研究存在一些問題，首先采用的AVHRR數據為一代數據產品，數據時間跨度不長，數據的空間分辨率不高(8 km×8 km);其次大多數研究多從行政區劃和地理單元探究植被活動特征，未考慮到不同生態系統間的鑲嵌性，因此難以詳細地探究不同生態系統間NDVI的差異。對此，本研究基于重構的1981—2018年長時間序列植被指數數據集(GIMMS NDVI3g和MODSI NDVI)，采用分區統計、NDVI趨勢分析法，分析中國大陸范圍(包含香港、澳門，未包含中國臺灣省)非耕地、非聚落和濕地生態系統NDVI的空間分布和時間動態變化特征，試圖為中國生態環境保護和可持續發展提供決策依據。

1 數據來源及預處理

1.1 數據來源

本研究利用2000年2月—2018年12月MODIS13 A3數據集植被指數產品的月合成數據(https:∥ladsweb.modaps.eosdis.nasa.gov/search/)，空間分辨率為1 km。還有美國國家航天航空局NASA Goddard Space Flight Center的全球監測與模型研究組(GIMMS)采用最大值合成法(MVC)合成的15 d、空間分辨率為8 km的最新GIMMS NDVI 3g.v1數據，時間范圍為1981年7月—2015年12月(https:∥ecocast.arc.nasa.gov/data/pub/gimms/3g.v1)。

中國陸地生態系統數據來源于中國科學院地理科學與資源研究所、中國科學院資源環境科學數據中心的“資源環境數據云平臺”。包括1980年、1990年、1995年、2000年、2005年、2010年和2015年，共7時期數據。該數據是在遙感解譯獲取的1∶10萬比例尺土地利用/土地覆蓋數據的基礎上，通過對各生態系統類型進行辨識和研究，經過分類處理形成多期中國陸地生態系統類型空間分布數據集，矢量數據柵格化生成的1 km柵格數據。生態系統類型依據土地利用/土地覆蓋遙感分類系統，歸分為7大生態系統類型:(1) 農田生態系統，主要包括水田和旱地;(2) 森林生態系統，主要包括密林地(有林地)、灌叢、疏林地、其他林地;(3) 草地生態系統，主要包括高覆蓋度草地、中覆蓋度草地、低覆蓋度草地;(4) 水體與濕地生態系統，主要包括沼澤地、河渠、湖泊、水庫、冰川與永久積雪、灘地;(5) 聚落生態系統，主要包括城鎮、農村居民地、工礦;(6) 荒漠生態系統，主要包括沙地、戈壁、鹽堿地;高寒荒漠;(7) 其他生態系統，主要包括裸土地和裸巖礫石地。

1.2 數據預處理

采用MRT軟件對MODIS NDVI數據進行拼接，投影轉換，數據格式轉換。采用Matlab軟件對GIMMS NDVI數據進行逐月最大值合成(MVC)，進一步減少云的影響，并降低月內物候循環的影響[17]，還原真值，格式轉換。由于GIMMS NDVI數據的空間分辨率約為8 km，將GIMMS NDVI由8 km重采樣為1 km的數據，不損失影像信息，以匹配MODIS NDVI。本文采用的地理坐標系是WGS-84,Albers投影坐標系，隨后統一采用IDL程序對兩套數據進行逐月裁切、掩膜、還原真值等操作。

2 研究方法

MODIS NDVI數據被稱為GIMMS NDVI數據的延續[18]，兩套數據集有一段時間的重疊部分，即2000年2月份-2015年12月份。已有學者[17,19]的研究表明可以利用兩套數據的重疊時間部分，以此建立逐月線性轉換方程，但未能從文章中找出適合本研究區域的轉換方程，本研究中建立的最佳轉換方程為y=0.970034x+0.016953，擬合方程的判定系數R2為0.88,p<0.01，擬合效果較好。利用轉換方程將1981年7月—2001年12月的GIMMS NDVI進行修正擬合，為了保留MODIS數據1 km的空間分辨率，將GIMMS NDVI數據也重采樣為1 km，并連接MODIS NDVI 2002年1月至2018年12月數據，最終獲得中國大陸地區的1981年7月至2018年12月長時間序列空間分辨率為1 km的NDVI數據集。

2.1 NDVI時空變化趨勢分析法

2.1.1 NDVI時間變化趨勢分析法 趨勢分析是對一組隨時間變化的變量進行回歸分析，得到其變化趨勢的方法。對多年植被指數數據逐像元NDVI值進行一元線性回歸，得到研究區內植被隨時間的變化趨勢，即NDVI的年際變化。計算公式為[20-21]:

(1)

式中:k為研究時間段內的年數,取1-n；YNDVIk是第k年某季相的NDVI平均值,如夏季NDVI平均值采用當年6月、7月、8月份NDVI值求得平均；slope>0表示植被覆蓋呈增加趨勢;slope小于0則相反;slope=0表示植被覆蓋沒有明顯的變化。

2.1.2 NDVI空間變化分析法 將大陸地區分為六大區域(東北、華北、華東、中南、西北、西南)，分區統計分析NDVI的總體變化特征，能整體上體現NDVI均值和趨勢變化情況。

2.2 主要生態系統類型的NDVI時空變化分析法

采用分區統計的方法來統計不同生態系統NDVI的變化情況，及其空間分異特征。由于生態系統數據為多年一幅的柵格數據，為充分利用生態系統數據，將連續時間序列的NDVI與最鄰近年的生態系統數據進行對應統計，表1是分區統計對應表。

表1 分區統計年份對應

3 結果與分析

3.1 中國植被NDVI時空變化分析

3.1.1 中國植被NDVI時間變化分析 圖1是整個中國大陸范圍內NDVI大于0值，2000年2月至2015年12月GIMMS NDVI平均值修正前后與MODIS NDVI平均值變化曲線，由圖1可知，原始GIMMS NDVI的周期性變化趨勢與MODIS NDVI相似，每年冬季原始GIMMS NDVI都小于MODIS NDVI，在植被生長季，兩者相差不大，修正后GIMMS NDVI與MODIS NDVI值相接近，具有一致性。值得注意的是MODIS NDVI數據在2000年初和2001年冬季，其值明顯低于修正前后GIMMS NDVI數據，可見MODIS數據在剛獲取的階段，數據質量可能不穩定，所以對于中國大陸地區，在使用MODIS NDVI數據時，可以從2001年7月份開始選用。

圖1 2000年2月至2015年12月MODIS NDVI、原始與修正后GIMMS NDVI時間序列曲線對比

由表2看出，修正后GIMMS NDVI在3種數據中的平均值和最小值都為最大，平均值最大為0.373，最小值最大為0.273;除了2000年和2001年，修正后GIMMS NDVI與MODIS NDVI的年平均值環比增長率和最大值有相似的變化規律，且有相同的NDVI線性趨勢(0.000 8/a)，最大值最小僅相差0.002。這說明，除了2000年、2001年，修正后GIMMS NDVI在年尺度上與MODIS NDVI相似的變化規律。

表2 不同數據年平均值、最大值、最小值和年平均值環比增長率

從表3中看出，MODIS NDVI和修正后GIMMS NDVI在2月至5月NDVI平均值環比增長率逐漸增加，而從5月至6月NDVI值月平均環比增長率有所下降，這是因為冬小麥成熟收割，新播種的作物未生長成熟造成的NDVI值下降現象，原始GIMMS NDVI未能捕捉到這個現象;修正后GIMMS NDVI平均值和最大值與MODIS NDVI有相似的變化規律，GIMMS NDVI月均值環比增長率最大為0.266%,MODIS NDVI為0.273%;修正后GIMMS NDVI最小值大于MODIS NDVI,2月份最大差值為0.096，這可能是由于2000年2月份MODIS NDVI數據剛剛獲取，數據質量不穩定導致NDVI值偏低，10月份相差值最小為0.008。

表3 不同數據月平均值、最大值、最小值和月平均值環比增長率

3.1.2 中國植被NDVI空間變化分析 從圖2中可以看出中國大陸夏季的NDVI的空間分布整體呈現東南和東北高、西北低的格局。在沿海的華北平原、長江三角洲和珠江三角洲等經濟發達地區呈現低NDVI值的現象，這是由于人類自身發展的需要，NDVI值隨著城鎮化的發展而降低;在新疆的阿爾泰山南坡、天山山脈北坡的NDVI值卻反常高，最高可達到0.896，是因為在準格爾盆地西部有缺口，大西洋和北冰洋的水汽能夠影響到兩側的山脈，形成迎風坡，對于植被生長有利;另外，在西北新疆塔里木盆地、準格爾盆地、甘肅西北部、青海西北部、內蒙古西部有沙漠、荒漠戈壁等，西藏自治區內的高原有高原荒漠和裸巖，氣候惡劣，不利于植物生長的環境，所以NDVI值普遍偏低，平均值為0.275。華北平原主要有農田，冬小麥收割之后，夏季播種農作物還未完全長成，摻雜土壤背景信息，因此NDVI值較低，平均值為0.471。位于華東地區的浙江、福建有丘陵，貴州、重慶、湖南和湖北西部有山地適宜森林生長，NDVI值較高，平均值為0.668。內蒙古東北部和東北三省內有大興安嶺、小興安嶺和長白山山脈，夏季森林和農作物生長旺盛NDVI值較高，平均值為0.775。

圖2 2018年夏季中國大陸NDVI平均值

表4 2018年夏季中國大陸六大分區NDVI均值

由圖3和表5，6可以看出過去38 a中，有49.923%的NDVI穩定不變，49.846%的NDVI輕微改善，中國大陸地區整體植被沒有發生劇烈變化。西北地區NDVI趨勢整體無明顯變化(0.000 098/a)，但新疆天山山脈北坡和塔里木河流域處發現NDVI的變化較為劇烈，最大為0.030 272/a，呈上升趨勢;內蒙古陰山南、甘肅省內祁連山、寧夏賀蘭山等NDVI變化為上升趨勢;在經濟發達的長江三角洲、珠江三角洲、京津翼地區NDVI趨勢呈現下降現象。在西藏藏南地區、內蒙古大興安嶺、東北三省內的小興安嶺和長白山，夏季NDVI值較高，而年際變化趨勢中呈現下降趨勢，說明植被覆蓋濃郁，多年以來的植被覆蓋有所下降。

表5 1981-2018年六大分區NDVI趨勢均值

表6 1981-2018年NDVI變化趨勢分類級別

圖3 1981-2018年中國大陸夏季NDVI變化趨勢

3.2 主要生態系統類型的NDVI變化分析

中國陸地生態系統具體分為農業、森林、草地、水體與濕地、聚落、荒漠和其他共7種生態系統類型。農田和聚落受人類活動影響較大，近30 a來農作物總產和單產都有所提高[22]，對于全國范圍內植被NDVI增加是正向作用，但想要更詳細地分析中國大陸的生態狀況則需要中國陸地生態系統分類數據的支持。自然生態系統有森林、草地、水體與濕地、荒漠或者其他生態系統，這些生態系統類型受人類影響較小，對生態系統狀況有指示意義，按表1對上述生態系統分區統計其NDVI長時間序列變化，有助于中國生物多樣性保護和開展生態環境狀況評價。

森林和農田生態系統在過去38 a中NDVI呈現波動式上升，上升趨勢分別為0.001 4/a，0.001 1/a，森林生態系統NDVI最低值為1982年0.692，最高值為2017年0.766。1987年5月6日至6月2日大興安嶺發生特大森林火災;2000年6月17，18日黑龍江大興安嶺呼中區呼源、松嶺區發生特大森林火災，同年8月6日新疆裕民縣巴旦杏林保護區發生重大森林火災;2002年7月28日和2004年6月22日內蒙古大興安嶺北部原始林區分別發生特大、重大森林大火;2018年6月1日和2日內蒙古大興安嶺汗馬國家級自然保護區和北部原始林區發生森林火災等，這些森林大火使NDVI值在時間序列上呈現出偏低或者谷。在過去38 a中草地生態系統的NDVI在0.358至0.402間變化，草地生態系統NDVI變化較為劇烈，說明NDVI對于不同覆蓋度草地的變化敏感，1995年、1997年錫林浩特市草原分別發生3次和1次特大火災，該年份NDVI值有所下降，草地生態環境狀況變差。荒漠生態系統的NDVI在過去38 a中呈下降趨勢(-0.000 6/a)，最低為0.094，最高為0.126，荒漠生態系統主要包含沙地、戈壁、鹽堿地、高寒荒漠，這些土地類型環境惡劣，十分脆弱，NDVI能夠反映出其變化狀況，在1981—2001年NDVI幾乎沒有大的變化，處于相對穩定的狀態，而在2002—2018年NDVI值的上升趨勢為0.000 9/a，說明近17 a的生態環境呈現向好的態勢發展。其他生態系統的NDVI在過去38 a中呈現下降趨勢(-0.000 7/a)，最低為0.138，最高為0.184，與荒漠生態系統類似，在1981—2001年，NDVI也變化不大，2002—2018年，其NDVI的上升趨勢與荒漠生態系統一樣，生態系統呈現向好態勢。水體與濕地生態系統NDVI在過去38 a中呈現下降趨勢(-0.001 7/a)，最低為0.311，最高為0.401，在1981—2001年NDVI呈現較穩定狀態，在2002—2018年中呈現上升趨勢(0.001 6/a)。

圖4 1981-2018年森林和農田生態系統NDVI變化

4 討 論

本研究利用NDVI探究生態系統狀況，但可能有所欠缺，另如氣候、地形等多種因素可以加以綜合考慮;趨勢線分析方法能很好地體現NDVI的變化趨勢，但是如何尋找合適的NDVI變化趨勢分級區間，需要進一步研究。劉可等[13]僅使用2010年一時期的中國陸地生態系統類型數據統計近30 a GIMMS NDVI3g數據，指出森林、農田、草地和水體與濕地總體NDVI非穩定上升，本文利用生態系統類型數據統計重構后的NDVI數據，指出森林和農田生態系統NDVI整體呈增強趨勢，草地生態系統NDVI呈非穩定上升趨勢，這可能是因為本文利用的是多期生態系統類型數據。利用生態系統分區統計不同年份的NDVI數據時，雖然盡可能利用生態系統數據，但最終統計的結果可能與實際的有差別，荒漠、水體與濕地和其他生態系統的NDVI曲線都在2000—2001年處有下降的陡坡，這可能是因為對GIMMS NDVI數據進行重采樣，仍無法解決其分辨率粗的問題，但若將MODIS NDVI數據進行降分辨率，則會損失大量信息。

5 結 論

(1) 構建數據集時，修正后GIMMS NDVI和MODIS NDVI相似的年份和月份變化規律;中國大陸地區MODIS NDVI在2000年、2001年的冬季數據質量不佳，可以從2001年7月開始選用數據。中國大陸地區夏季NDVI整體分布特征呈現西北低，東南和東北高，在經濟發達區域，由于人類經濟社會活動因素出現低NDVI值現象。

(2) 在過去38 a中，對受人類活動影響較小的生態系統，森林生態系統植被呈上升趨勢，森林火災對森林的影響在NDVI序列變化能夠反映出來;草地生態系統呈現穩定趨勢，NDVI變化較為劇烈，NDVI序列可以作為監測草地生態狀況變化的指示;荒漠、水體與濕地、其他生態系統的都呈下降趨勢，三者在1981—2001年NDVI變化都較為穩定，且2002年以后NDVI呈現向好態勢發展，這與MODIS傳感器分辨率和波段寬度設置優于AVHRR傳感器有關。