2000-2018年間青藏高原植被覆蓋變化及其與氣候因素的關系分析

韓炳宏， 周秉榮， 顏玉倩， 石明明， 蘇淑蘭， 趙恒和，牛得草， 傅 華

(1. 中國氣象局成都高原氣象研究所， 四川 成都， 610072； 2. 青海省防災減災重點實驗室， 青海 西寧， 810001；3. 青海省氣象科學研究所， 青海 西寧， 810001； 4. 青海省海南州氣象局， 青海 共和， 813099；5. 蘭州大學草地農業生態系統國家重點實驗室， 甘肅 蘭州， 730020)

植被作為全球陸地生態系統的重要組分，對全球物質能量循環、碳平衡調控及維持氣候穩定等方面起著重要作用[1]。光照、氣溫和降水等氣候要素是植被生長發育不可或缺的環境因子，對植被健康穩定的生長狀況具有顯著影響[2-3]。歸一化植被指數(Normalized Difference Vegetation Index,NDVI)作為評價植被生長狀況的重要因子，不僅可以反映全球或區域尺度上植被生長狀況和分配格局，也可以表觀長時間序列上的變化趨勢[4]。通常與植被生產力和葉面積指數等重要植被特征參數密切相關[5]。有些研究結果表明，1982—1999年北半球中高緯度和中國大多數地區的NDVI表現出增加趨勢[6-7]。而有些研究結果顯示，2000—2009年北半球多個地區植被NDVI呈減弱趨勢[8]。Piao等[9]對歐亞大陸溫帶和寒帶生長季植被NDVI變化趨勢進行了研究，認為1982—1997年顯著增加，而1998—2006年呈下降趨勢，尤其在春季和夏季表現最為明顯；Mohammat等[10]認為，亞洲內陸地區生長季植被由于受春季變冷和夏季干旱的影響，其NDVI增加趨勢于1990年停止；杜加強等[11]研究表明，1998—2012年新疆地區夏季植被NDVI由之前的極顯著增加轉變為顯著減少。

青藏高原是中國面積最大、世界上海拔最高的高原，堪稱地球“第三極”[12]，擁有獨特的生態系統類型，被認為是亞洲乃至北半球氣候變化的“感應器”[13]，不僅關系到高原周邊地區數億居民的供水安全，同時也對我國乃至東亞地區生態系統穩定起著屏障作用[12,14]。隨著全球氣候變暖和人類活動加劇，青藏高原正面臨生態環境壓力增加與生態安全屏障功能改變等風險[12]。由于受氣候變化和人類活動的雙重影響，青藏高原高寒植被變化受到越來越廣泛的關注。近30a來，尤其是21世紀以來，青藏高原植被覆蓋度呈不明顯的綠度增加或變綠趨勢[15-16]。趙紫薇[17]利用1982—2012年全球監測模型與制圖研究歸一化植被指數(Glaobal Inventory Modelling and Mapping Studies Normalized Difference Vegetation Index,GIMMS NDVI)和2001—2013年中尺度分辨率成像光譜歸一化植被指數(Moderate Resolution Imaging Spectroradiometer Normalized Difference Vegetation Index,MODIS NDVI)對近30 a青藏高原植被動態變化時序進行了對比分析，發現青藏高原地區植被呈整體改善趨勢；孟夢等[18]利用1982—2013年GIMMS NDVI數據研究了青藏高原NDVI變化趨勢及其對氣候變化的響應，指出近年來青藏高原植被長勢逐漸變良好，覆蓋度呈增加態勢。以上研究雖得到了較為相近的結果，但由于受遙感數據時序的限制，并未涉及近5a的草地植被變化狀況，加之研究方法和不同時段采用的遙感影像產品不同，獲取的結果僅能代表過去某一時段草地植被生長狀況的好壞，同時也缺少對不同時段草地植被與同期氣象因子的關系分析。為此，本文基于遙感監測分析，利用2000—2018年MODIS NDVI逐旬數據對青藏高原植被覆蓋時空變化特征及其驅動因素進行系統研究，以期為青藏高原植被在區域乃至全球氣候變化過程中的調控機理提供科學依據。

1 材料與方法

1.1 研究區概況

青藏高原南起喜馬拉雅山脈南緣，北至昆侖山、阿爾金山和祁連山北緣，西部為帕米爾高原和喀喇昆侖山脈，東及東北部與秦嶺山脈西段和黃土高原相接，地處26°00′～39°47′ N，73°19′～104°47′ E，東西長約2.8×103km，南北寬約3.0×102～1.5×103km，總面積約2.5×106km2，地形上可分為藏北高原、藏南谷地、柴達木盆地、祁連山地、青海高原和川藏高山峽谷區6個部分[12-13]，主要包括中國西藏全部和青海、新疆、甘肅、四川、云南的部分地區(圖1)。其自然歷史發育極其年輕，受多種因素共同影響，東南部屬于暖濕性氣候，西北部屬干冷性氣候，具有太陽輻射強、氣溫低、氣溫日較差大、年變化小、干濕季節分明等特點；年降水量為486 mm，年均溫在—5.6～8.6℃之間，年日照時數為2 300～3 600 h，太陽總輻射為140～180 C·m-2。高原上凍土廣布，植被多為草地植被[19]。

圖1 研究區地勢及高程Fig.1 Terrin and elevation of study area

1.2 數據來源及處理

NDVI是通過紅外與近紅外波段的組合實現對植被信息狀態的表達，通常用來表征植被覆蓋程度、生長狀況、葉面積指數、生物量以及光合有效輻射等植被參數[19]，已被廣泛應用到全球及區域尺度上植被變化與氣候互作研究中。

文中NDVI數據來自美國LP DAAC數據中心(http://lp-daac.usgs.gov/main.asp)MODIS儀器提供的2000年1月—2018年12月逐旬NDVI數據，分辨率為1km，數據格式為HDF-EOS，投影方式為Sinusoidal。數據的合成、再投影處理通過MRT和投影軟件ENVI實現。根據青藏高原范圍將下載的MODIS數據進行拼接，使其投影為標準的經緯度網格坐標，將軟件輸出的MODIS資料保存為GEOTIF格式；在ENVI里完成青藏高原矢量邊界裁剪，獲得研究區的NDVI數據。為更好地反映研究區不同時間尺度植被覆蓋狀況，將數據另存為ASCII格式，以便用于植被指數與氣象要素的相關分析。

1.3 研究方法

1.3.1最大值合成法 最大值合成法(Maximum Synthesis Method,MVC)是目前國際上通用的最大合成法，采用該方法取一個月每旬的數據最大值為月值，進一步消除云、大氣和太陽高度角的干擾[20-21]。此法假設一個月每旬中NDVI值最大的那一天天氣是晴朗的，不受云層的影響，就取這個最大值作為月NDVI值，計算公式如下：

NDVIi=Max(NDVIij)

(1)

式中NDVIi為第i月的NDVI值，NDVIij為第i月第j旬的NDVI值，可以認為NDVIi是一月內植被最豐盛時期NDVI值[21]。利用(1)方法依次獲得年和季節NDVI值，然后采用公式NDVI=0.004DN-0.1將其轉化為真實NDVI值。NDVI的取值范圍為—1.0～1.0，一般認為NDVI達到0.1以上表示有植被覆蓋，值越大植被覆蓋度越大；0.1以下則表示地表無植被覆蓋，如裸土、沙漠、戈壁、水體、冰雪和云[22-23]。

1.3.2均值法 均值處理是將某時間間隔內的NDVI數據求均值，以消除或降低由時間段端點年份氣候異常對植被生長狀況的影響；每年的NDVI值由各月最大NDVIi求均值獲得[24]，計算公式如下：

(2)

式中MNDVI,i表示第i年的NDVI，n表示月份，NDVIi是第i月的最大NDVI值。

1.3.3差值法 差值法是用于2個時期NDVI值的變化，即后一時期所有像元NDVI值與前一時期所有像元的差。宋怡等[25]利用該方法計算了年份間NDVI值的差；徐慧等[24]也同樣采用該方法計算了年份間NDVI值的差。不同時期NDVI差值的計算公式如下：

ΔNDVI=MNDVI,i-MNDVI,j

(3)

式中NDVI表示NDVI差值，MNDVI,i表示后一年NDVI值，MNDVI,j表示前一年NDVI值。

1.3.4NDVI變化趨勢(Slope)分析 一元線性回歸分析可以模擬每個柵格的變化趨勢，Stow等[20]就用該方法來模擬植被的綠度變化率(GreennessRateofChange,GRC)，GRC被定義為某時間段內的季節合成歸一化植被指數年際變化的最小次方線性回歸方程的斜率。此處同樣用該方法來模擬多年NDVI的變化趨勢，計算公式為：

(4)

式中變量i為1～19序號，MNDVI,i表示第i年的NDVI值。變化趨勢圖則反映了近19a青藏高原地區的NDVI的變化趨勢。某格點的趨勢線是這個格點19a的NDVI值用一元線性回歸模擬出來的一個總的變化趨勢，即θSlope為這條趨勢線的斜率。這個趨勢線并不是簡單的最后一年與第一年的連線。若θSlope>0，則說明NDVI在19a間的變化趨勢是增加的；若θSlope<0則是減少。另外，利用ArcGIS空間分析工具Reclass命令，將θSlope值進行重分類，依次劃分為嚴重退化[—0.030，—0.010]、中度退化[—0.010，—0.005]、輕微退化[—0.005，—0.003]、保持不變[—0.003，0.003]、輕微改善[0.003，0.005]、中度改善[0.005，0.010]和明顯改善[0.010，0.030]7個等級[19]；再通過Zonal Satistics命令，得出7個等級所對應的像元數、面積、均值、標準差和面積百分比。

2 結果與分析

2.1 2000—2018年年際NDVI變化

通過2000—2018年的年平均NDVI進行逐年差值分析，從變化幅度上看，逐年NDVI的差值結果顯示2001，2006和2012年分別在2000，2005和2011年的基礎上顯著增加，增加幅度分別為0.126，0.106和0.11；2005，2007，2008和2009年青藏高原年均NDVI較2004，2006，2007和2008年分別減少了0.02，0.032，0.052和0.046。另外，除2001和2006年外，年平均NDVI在2000—2009年呈減少趨勢；2010—2018年呈不同程度的增加態勢。整體來看，青藏高原地區近19a年平均NDVI呈逐漸增加趨勢(圖2)。

圖2 2000—2018年年平均NDVI差值Fig.2 Difference of annual average NDVI from 2000 to 2018

2.2 2000—2018年NDVI季節變化

本文采用氣候學上公認的季節劃分方法，3—5月為春季、6—8月為夏季、9—11為秋季、12—翌年2月為冬季，分別將其NDVI平均獲得季節NDVI，再將19個年份同一季度NDVI求均值，從而獲得4張季節NDVI影像，以分析青藏高原2000—2018年NDVI平均生長狀況在季節上的空間分布及差異(圖3)。

整體來看，2000—2018年青藏高原NDVI自西北向東南依次增加，且四季NDVI空間分布差異較大(圖3)。從地理位置來看，低植被覆蓋區(NDVI<0.2)主要分布于西藏大部、新疆和甘肅局部以及青海西北部地區；中植被覆蓋區(0.2≤NDVI≤0.3)主要位于青海與甘肅、西藏、四川和云南的交界區；高植被覆蓋區(NDVI>0.3)主要分布在四川和云南大部、青海和甘肅以及西藏東南局部地區。從草地類型分布來看，低植被覆蓋區主要以荒漠、戈壁和裸地為主，中植被覆蓋區以高寒草原和溫性草原為主，高植被覆蓋區以高寒草甸為主。

圖3 2000～2018年NDVI季節變化Fig.3 NDVI seasonal changes from 2000 to 2018

2.3 趨勢分析

從NDVI趨勢分析結果(圖4)可以看出，在四川和云南局部、青海中東部及南部、西藏東南部及北部部分地區，NDVI的變化趨勢呈明顯的上升趨勢，說明近19a中這些區域的植被情況得到比較好的改善。而在新疆和甘肅局地、青海東南部及東北部、西藏東南部及其西部的邊遠地區，NDVI呈下降趨勢。尤其在和田南部局地、林芝大部、祁連山部分地區及果洛藏族自治州大部，降低的趨勢比較明顯，即在這些地區植被退化情況較為嚴重。另外，由表1可知19a年來，草地植被退化總面積為60.03 km2，恢復改善總面積達69.13 km2。故整個青藏高原地區植被生長狀況得到改善的區域面積大于植被退化的面積。其中，植被改善區域面積占整個青藏高原總面積的27.35%，退化區域占23.75%。因此，青藏高原植被雖局部惡化，但整體仍處于恢復狀態。這也進一步驗證了前面的分析結果。

圖4 2000—2018年MNDVI變化趨勢Fig.4 Change trends of MNDVI from 2000 to 2018

表1 2000—2018年MNDVI變化趨勢結果統計Table 1 Statistics of MNDVI change trend from 2000 to 2018

θSlope植被變化趨勢Vegetation variation trend像元數Pixelelements面積Areas/ km2平均值Average標準差Standard deviation百分比Percentage %—0.030～—0.010嚴重退化Serious degradation48 8073.89×1040.170.141.54—0.010～—0.005中度退化Moderate degradation145 01711.00×1040.150.154.35—0.005～—0.003輕微退化Slight degradation566 20945.14×1040.240.1617.86—0.003～0.003保持不變Remain unchanged1 551 107123.65×1040.140.1248.910.003～0.005輕微改善Slight improvement616 03649.11×1040.180.1419.430.005～0.010中度改善Moderate improvement193 92016.03×1040.530.136.340.010～0.030明顯改善Serious improvement50 1423.99×1040.670.061.58

2.4 2000—2018年NDVI與氣候因子的關系

全球氣候變暖有可能加速土壤有機質的分解而促進較為寒冷的地區植被凈初級生產力的提高已基本形成共識[1,12-13]。在大部分相對干旱的區域，植被生長與降水呈正相關；在濕潤地區，則為負相關。本研究發現19a中青藏高原地區NDVI與氣溫和降水的關系較好(圖5a，b)。其中，NDVI與氣溫的相關性系數達0.909；NDVI與降水的相關性雖不及氣溫，但相關系數也能達到0.793。水在植物的營養物質輸送、結構維持和各種生理生化過程中起著十分重要的作用，因此植被的生長與降水量的多少具有明顯的正相關關系；植物葉片的凈光合速率隨溫度的升高而增加，氣溫過高或過低都會影響植物的生理生化過程。從兩者的線性關系可以看出，當氣溫每升高1℃，NDVI增加0.128；降水每增加100 mm，NDVI相應增加0.172。

圖5 NDVI與氣象因子的關系Fig.5 Relationship between NDVI and meteorologic factors

3 討論

青藏高原作為亞洲生態安全屏障和全球氣候變化敏感區，其植被生長狀況對區域或全球尺度上氣候變化的響應具有重要的調控作用。為此，本文利用2000—2018年MODIS NDVI數據從季節、年際水平對青藏高原不同地區、不同時間尺度上植被覆蓋變化及其驅動因素進行了分析。結果表明，近19a來，青藏高原地區NDVI雖局部有所惡化，但整體仍呈穩定恢復趨勢；與卓嘎等[19]的研究結論一致。然而，Shen等[26]利用增強型植被指數(Enhanced Vegetation Index,EVI)研究了2000—2012年生長季青藏高原植被覆蓋變化，認為青藏高原植被覆蓋變化呈不顯著減少趨勢，主要表現為植被減少速率大于植被增加速率；且不同植被類型對氣候變化的敏感性不同。其研究結論與本文不一致，這主要是因為Shen等[26]的研究時段以及研究方法不同。另外，我們發現自21世紀以來，特別是2000—2009年期間，青藏高原植被覆蓋呈現不顯著的減少趨勢，即植被退化面積和速率大于植被恢復，2010—2018年間植被覆蓋不同程度的增加，這可能與近10年青藏高原地區氣溫和降水的明顯增多有關[19]。宋怡等[25]基于1998—2004年SPOT_VGT NDVI數據對中國西北地區植被變化特征進行了系統研究，研究表明：西北地區植被覆蓋在1998—2004年間出現退化趨勢，尤其以2000—2002年間的變化幅度較大；青藏高原地區植被雖局部區域有改善，但總的改善幅度小于退化幅度。由此可見，青藏高原植被覆蓋在2010年之后出現了較大幅度的改善，且逐年恢復的速率和面積大于退化速率和面積。

植物對氣溫和降水的響應機制主要表現在其對植物生理生化過程中酶活性的影響方面，水是植物代謝途徑中的主要生理生化過程(糖酵解、三羧酸循環、磷酸戊糖、乙醛酸循環和乙醇酸氧化等)的參與者，而氣溫則對植物代謝途徑中的一系列酶(氧化酶、還原酶和脫氫酶等)的活性產生重要影響。以往研究表明，青藏高原氣候條件與NDVI的相關存在明顯的區域性差異。青藏高原東北部氣溫和降水與NDVI具有較好的正相關關系，西南部呈負相關關系[19]。孟夢等[18]研究發現青藏高原植被覆蓋與空氣中的水汽含量顯著相關，對最低氣溫的響應最為敏感。我們發現青藏高原地區NDVI與氣溫和降水的相關性較好。其中，NDVI與氣溫間的相關性大于其與降水間的相關性，這與以往研究[18-19,26-28]結論一致；這也說明在全球氣候變暖背景下，青藏高原地區氣候在暖干化向暖濕化演變的過程中，青藏高原地區植被覆蓋對氣溫響應的敏感性高于降水。

由于青藏高原下墊面類型分布多樣，地貌特征復雜，降水格局不均等特點，青藏高原地區植被覆蓋和類型表現出明顯的區域性差異。另外，青藏高原地區氣象站點主要集中于高原中東部，高原西北部地區站點相對較少，甚至個別區域無觀測站點，故利用本文所得結論對整個青藏高原地區植被覆蓋進行描述仍存在一定差異。因此，今后需采用氣候格點數據進一步探討青藏高原不同地區植被覆蓋差異及其對氣候變化的響應及互饋機制。

4 結論

近19a來，青藏高原地區植被NDVI主要經歷了2個顯著階段，即前10a呈減少趨勢，后9a呈不同程度的增加態勢。NDVI整體呈逐漸增加趨勢。從地理位置來看，低植被覆蓋區主要分布于西藏大部、新疆南部和甘肅局部以及青海西北部地區；中植被覆蓋區主要位于青海與甘肅、西藏、四川和云南的交匯區域；高植被覆蓋區主要分布在四川和云南西北部大部分地區、青海和甘南以及西藏東南局部地區。NDVI趨勢分析結果顯示，除高原局部地區植被有所退化，但大部地區植被生長狀況明顯改善，且草地植被改善的面積大于退化。因此，青藏高原植被雖局部惡化，但整體仍處于恢復狀態。青藏高原地區NDVI與同期氣溫和降水間均有較好的正相關關系，但氣溫要好于降水。且氣溫每升高1℃，NDVI增加0.128；降水每增加100 mm，NDVI相應增加0.172。