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

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

(1. 中國氣象局成都高原氣象研究所， 四川 成都， 610072； 2. 青海省防災(zāi)減災(zāi)重點(diǎn)實(shí)驗(yàn)室， 青海 西寧， 810001；3. 青海省氣象科學(xué)研究所， 青海 西寧， 810001； 4. 青海省海南州氣象局， 青海 共和， 813099；5. 蘭州大學(xué)草地農(nóng)業(yè)生態(tài)系統(tǒng)國家重點(diǎn)實(shí)驗(yàn)室， 甘肅 蘭州， 730020)

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

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

1 材料與方法

1.1 研究區(qū)概況

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

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

1.2 數(shù)據(jù)來源及處理

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

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

1.3 研究方法

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

NDVIi=Max(NDVIij)

(1)

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

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

(2)

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

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

ΔNDVI=MNDVI,i-MNDVI,j

(3)

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

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

(4)

式中變量i為1～19序號(hào)，MNDVI,i表示第i年的NDVI值。變化趨勢圖則反映了近19a青藏高原地區(qū)的NDVI的變化趨勢。某格點(diǎn)的趨勢線是這個(gè)格點(diǎn)19a的NDVI值用一元線性回歸模擬出來的一個(gè)總的變化趨勢，即θSlope為這條趨勢線的斜率。這個(gè)趨勢線并不是簡單的最后一年與第一年的連線。若θSlope>0，則說明NDVI在19a間的變化趨勢是增加的；若θSlope<0則是減少。另外，利用ArcGIS空間分析工具Reclass命令，將θSlope值進(jìn)行重分類，依次劃分為嚴(yán)重退化[—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個(gè)等級(jí)[19]；再通過Zonal Satistics命令，得出7個(gè)等級(jí)所對應(yīng)的像元數(shù)、面積、均值、標(biāo)準(zhǔn)差和面積百分比。

2 結(jié)果與分析

2.1 2000—2018年年際NDVI變化

通過2000—2018年的年平均NDVI進(jìn)行逐年差值分析，從變化幅度上看，逐年NDVI的差值結(jié)果顯示2001，2006和2012年分別在2000，2005和2011年的基礎(chǔ)上顯著增加，增加幅度分別為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年呈不同程度的增加態(tài)勢。整體來看，青藏高原地區(qū)近19a年平均NDVI呈逐漸增加趨勢(圖2)。

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

2.2 2000—2018年NDVI季節(jié)變化

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

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

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

2.3 趨勢分析

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

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

表1 2000—2018年MNDVI變化趨勢結(jié)果統(tǒng)計(jì)Table 1 Statistics of MNDVI change trend from 2000 to 2018

θSlope植被變化趨勢Vegetation variation trend像元數(shù)Pixelelements面積Areas/ km2平均值A(chǔ)verage標(biāo)準(zhǔn)差Standard deviation百分比Percentage %—0.030～—0.010嚴(yán)重退化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與氣候因子的關(guān)系

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

圖5 NDVI與氣象因子的關(guān)系Fig.5 Relationship between NDVI and meteorologic factors

3 討論

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

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

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

4 結(jié)論

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