袁 雷，周刊社，張東東

（西藏自治區(qū)氣候中心，拉薩 850000）

0 引言

IPCC第六次評估報告第一組工作報告指出，2011—2020年全球地表溫度比工業(yè)革命時期上升了1.09℃。從未來20年的平均溫度變化預估來看，全球溫升預計將達到或超過1.5℃。在全球變暖的背景下，陸地生態(tài)系統(tǒng)正在受到氣候變化的影響，作為陸地生態(tài)系統(tǒng)的主要組分，植被對氣候變化十分敏感[1]。歸一化植被指數(shù)(normalized differential vergetation index，NDVI)是反映植被生長狀況的最佳指標之一[2]，廣泛用于植被生長狀況演變分析[3-4]，植被NDVI在不同時空尺度下的變化特征及與氣候變化的關(guān)系也有很多研究[5-7]。青藏高原屬于氣候變化的敏感區(qū)和生態(tài)脆弱帶，被認為是研究陸地生態(tài)系統(tǒng)對氣候變化響應機制的理想場所[8]。羌塘國家級自然保護區(qū)位于青藏高原腹地，由于其獨特的高寒荒漠生態(tài)系統(tǒng)，生態(tài)更為脆弱和敏感[9]，已經(jīng)有很多學者對羌塘的植被生態(tài)[9-11]、荒漠化[12-13]、土壤溫濕度[14-15]、降水[16]、持續(xù)干旱日數(shù)[17]、冰雹日數(shù)[18]、積雪[19]、土壤凍結(jié)天數(shù)[20]、地表濕潤度[21]等進行了研究。其中，吳曉萍等[9]利用2001—2010年的NDVI數(shù)據(jù)及氣象資料對羌塘國家級自然保護區(qū)東、中、西部各縣NDVI、氣溫、降水變化進行了分析，認為保護區(qū)東部NDVI呈顯著上升趨勢，中部和西部呈下降趨勢。東部年均NDVI與氣溫呈顯著正相關(guān)、與降水呈顯著負相關(guān)，中部反之；西部年均NDVI與氣溫、降水均呈顯著負相關(guān)。杜軍等[11]分析了羌塘氣溫、降水、水體以及NDVI變化趨勢，認為氣溫、降水和NDVI都呈增加趨勢。這些研究以區(qū)域平均NDVI分析植被變化，很難體現(xiàn)小區(qū)域的植被NDVI較詳細的變化特征，對于不同功能區(qū)植被NDVI時空變化特征及其差異沒有進行定量化研究，也少有羌塘植被NDVI對氣候變化響應的研究。鑒于此，筆者以MODIS傳感器的NDVI數(shù)據(jù)作為基礎(chǔ)，以羌塘國家級自然保護區(qū)各功能區(qū)為研究對象，逐像元定量計算其NDVI時空變化趨勢及其穩(wěn)定性，期望能夠為羌塘國家級自然保護區(qū)生態(tài)環(huán)境保護提供科學參考。

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

1.1 研究區(qū)概況

羌塘國家級自然保護區(qū)（以下簡稱保護區(qū)）是僅次于格陵蘭國家公園的世界第二大陸地自然保護區(qū)，位于西藏自治區(qū)西北部，平均海拔5000 m以上，被稱為“世界屋脊的屋脊”，保護區(qū)內(nèi)植被種類較少、群落結(jié)構(gòu)簡單，從東南向西北草地植被大體呈高寒草甸草原、高寒草原、高寒荒漠草原、高寒荒漠的分布[11]。高寒植被以紫花針茅(Stipa purpurea)、藏羊茅(Festuca wallichanica)和高山蒿草(Kobresia pygmaea)等為主[10]，有實驗區(qū)、緩沖區(qū)和核心區(qū)3個大的功能分區(qū)（圖1）。

圖1 羌塘國家級自然保護區(qū)地理位置

保護區(qū)年平均氣溫大都在0℃以下，年平均地溫為-5.0～-1.5℃，降水量為95.6～294.9 mm，光照條件充足，年日照時數(shù)為2800～3600 h。杜軍[11]研究表明1971—2017年自然保護區(qū)附近氣象站點年平均氣溫以0.46℃/10 a的速率顯著升高，明顯高于同期全球和亞洲地表溫度的升溫率。

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

NDVI數(shù)據(jù)為2000—2020年的MODIS NDVI月數(shù)據(jù)，空間分辨率為0.01°，來自于國家氣象中心，已經(jīng)進行了大氣校正、輻射校正和幾何校正。在研究植被變化時，使用最大值合成法(maximum value composites，MVC)計算植被生長季（6—9月）的NDVI數(shù)據(jù)。氣象數(shù)據(jù)使用中國第一代大氣和陸面再分析產(chǎn)品(CRA/LAND)的2000—2020年月氣溫和月總降水量數(shù)據(jù)，該數(shù)據(jù)來自于國家氣象信息中心，使用該數(shù)據(jù)計算出植被生長季平均氣溫和植被生長季總降水量。由于使用的CRA/LAND數(shù)據(jù)空間分辨率為0.25°，因此在分析NDVI與氣象因子相關(guān)性時，將植被生長季NDVI數(shù)據(jù)重采樣為0.25°。植被類型數(shù)據(jù)來源于中科院1:100萬中國植被數(shù)據(jù)集。

1.3 分析方法

1.3.1 Sen趨勢分析 使用Sen趨勢分析來逐像元描述植被生長季NDVI及氣候因子的長期變化趨勢[22-23]。Sen趨勢分析又被稱為Sen斜率估計，計算時間序列斜率對的中值，抗噪聲能力強，是一種穩(wěn)健的非參數(shù)統(tǒng)計的趨勢計算方法。該方法計算效率高，對于測量誤差和離群數(shù)據(jù)不敏感，常被用于長時間序列數(shù)據(jù)的趨勢分析中。Sen斜率計算如式(1)所示。

式中，xj和xi為時間序列數(shù)據(jù)，使用趨勢度β來判斷時間序列趨勢的升降，當β＞0時時間序列呈上升的趨勢，當β＜0時時間序列呈下降的趨勢。

1.3.2 Mann-Kendall趨勢檢驗 使用Mann-Kendall(MK)檢驗法對趨勢分析結(jié)果進行檢驗[22-23]，趨勢檢驗法過程如下：對于時間序列Xt=( )x1,x2,...,xn，做如下假設(shè)：H0—假設(shè)序列中的數(shù)據(jù)為獨立同分布樣本，即無顯著趨勢，H1—假設(shè)序列存在上升或下降單調(diào)趨勢。檢驗統(tǒng)計量S由式(2)～(3)計算。

本文中時間序列長度為21（2000—2020年），采用Sen趨勢逐像元計算保護區(qū)植被生長季NDVI、植被生長季氣溫和植被生長季降水量變化趨勢；使用MK檢驗法逐像元計算保護區(qū)植被生長季NDVI、植被生長季氣溫和植被生長季降水量變化趨勢顯著性。最終結(jié)果分為5級[4]，即極顯著下降(slope＜0，P≤0.01)、顯著下降(slope＜0，0.01＜P≤0.05)、變化不顯著(P＞0.05)、顯著上升(slope＞0，0.01＜P≤0.05)、極顯著上升(slope＜P≤0，P≤0.01)。

1.3.3 變化趨勢率 使用最小二乘法分析研究區(qū)及各功能區(qū)平均NDVI的變化趨勢率，其計算如式(5)[24]。

式中，n為研究序列長度，在本研究中為21；i表示第i年，xi表示第i年的NDVI，θ表示變化趨勢率。

1.3.4 穩(wěn)定性 變異系數(shù)CV即標準差與平均值的比值[24]，反映變異程度，廣泛運用于波動水平的分析中，計算如式(6)。

式中，SDNDVI為逐年植被生長季NDVI的標準差，為逐年植被生長季NDVI的平均值。CV值消除了單位和平均值不同對2個或多個變量變異程度比較的影響。本研究采用CV分析保護區(qū)植被生長季NDVI逐個像元21年間的變異情況，以揭示植被生長季NDVI的穩(wěn)定性。并按幾何間隔法將穩(wěn)定性分為5類[4]，即高穩(wěn) 定 (0.017＜CV≤0.055)、較 高 穩(wěn) 定 (0.055＜CV≤0.058)、中 等 穩(wěn) 定 (0.058＜CV≤0.096)、較 低 穩(wěn) 定(0.096＜CV≤0.692)、低穩(wěn)定(0.692＜CV≤9.954)。

1.3.5 相關(guān)性檢驗 采用皮爾遜(Pearson)相關(guān)分析方法研究氣象要素與植被生長季NDVI之間的關(guān)系，使用T檢驗來檢驗相關(guān)系數(shù)顯著性。相關(guān)系數(shù)如式(7)所示[24]。

式中，rxy相關(guān)系數(shù)，n為年份，x為自變量（氣溫、降水等氣象要素），y為因變量（植被生長季NDVI）。

我小學五年級因抗戰(zhàn)爆發(fā)、家鄉(xiāng)淪陷、學校停辦而失學，一直在老家農(nóng)村種地。實際上即使我在上學的時候，也是一直跟著大人下地勞動的，農(nóng)村的孩子，一般十來歲早就下地勞動了。

2 結(jié)果與分析

2.1 保護區(qū)植被生長季NDVI空間分布特征

從整個保護區(qū)來看，植被生長季NDVI空間分布差異較大（圖2），表現(xiàn)為由東南向西北逐漸減小的空間分布特點。2000—2020年植被生長季NDVI的平均值為0.138，根據(jù)植被生長季NDVI值大小將保護區(qū)格點像元劃分為5個等級[25]，分別為≤0.1、0.1～0.2、0.2～0.3、0.3～0.4、＞0.4。植被生長季NDVI≤0.1為無植被覆蓋區(qū)，植被生長季NDVI處于0.1～0.2區(qū)間為低植被覆蓋區(qū)，植被生長季NDVI處于0.2～0.4為中等植被覆蓋區(qū)，植被生長季NDVI大于0.4為高等植被覆蓋區(qū)[25]，無植被覆蓋區(qū)占保護區(qū)面積的22.97%，主要分布于北緩沖區(qū)昆侖山南麓，植被群系為墊狀駝絨藜荒漠；低植被覆蓋區(qū)占保護區(qū)面積的64.39%，植被群系為青藏薹草草原、紫花針茅草原和部分沙生針茅草原；中等植被覆蓋區(qū)占保護區(qū)面積的12.56%，主要分布于實驗區(qū)東部和南緩沖區(qū)中東部，植被群系為紫花針茅草原、三指雪蓮花和西藏扁芒菊稀疏植被；高等植被覆蓋區(qū)只占保護區(qū)面積的0.08%，零星分布于保護區(qū)東部。

圖2 生長季NDVI空間分布圖

在保護區(qū)的3個功能區(qū)中，低植被覆蓋區(qū)面積占比均達到58.50%以上，高植被覆蓋區(qū)面積占比均在0.12%以下，表明3個功能區(qū)整體植被條件均較差；實驗區(qū)、緩沖區(qū)和核心區(qū)的中等植被覆蓋區(qū)域面積占比分別為21.98%、11.21%和8.30%，表明實驗區(qū)的植被條件在保護區(qū)中最好。3個功能區(qū)植被生長季NDVI均是東部好于西部。

2.2 趨勢分析

2.2.1 趨勢變化空間特征 從近21年保護區(qū)植被生長季NDVI變化趨勢空間分布來看（圖3a），保護區(qū)86.39%區(qū)域植被生長季NDVI呈緩慢的上升趨勢，變化趨勢率為(0～0.01)/10 a的區(qū)域占保護區(qū)面積的57.51%，在保護區(qū)廣泛分布；變化趨勢率為(0.01～0.03)/10 a的區(qū)域占保護區(qū)面積的27.06%，主要分布于保護區(qū)西北部。從變化趨勢顯著性檢驗來看（圖3b），植被生長季NDVI上升達到顯著性水平(P＜0.05)的區(qū)域占保護區(qū)面積的45.15%，主要分布在保護區(qū)中部和北部，增加趨勢多在(0～0.03)/10 a。植被生長季NDVI下降達到顯著性水平(P＜0.05)的區(qū)域占保護區(qū)面積的1.11%，零星分布在緩沖區(qū)南部?？傮w來看，保護區(qū)植被狀況明顯改善。結(jié)合圖2～3可以看出，植被生長季NDVI變化趨勢率基本上表現(xiàn)為：植被生長季NDVI值小的區(qū)域，變化上升趨勢率?。恢脖簧L季NDVI值大的區(qū)域，變化上升趨勢率大；植被生長季NDVI值處于0.1～0.2范圍時保護區(qū)西北部變化趨勢率高于東南部。實驗區(qū)、緩沖區(qū)和核心區(qū)植被生長季NDVI上升達到顯著性水平（顯著上升和極顯著上升）的面積占比分別是24.47%、47.43%、57.27%，表明核心區(qū)植被生長季NDVI改善程度優(yōu)于實驗區(qū)和緩沖區(qū)。

圖3 生長季NDVI變化趨勢(a)及顯著性檢驗(b)

2.2.2 NDVI年際趨勢變化特征 從保護區(qū)植被生長季各等級NDVI面積比例逐年變化分析（圖4），植被生長季NDVI值處在0.1～0.2范圍的面積比例在每年均最高，變化區(qū)間在57.7%～67.8%，總體表現(xiàn)出明顯的增加趨勢，平均每年增加0.33個百分點(P＜0.001)；植被生長季NDVI≤0.1的面積比例呈現(xiàn)明顯的減少趨勢，平均每年減少0.56個百分點(P＜0.001)，到2020年只占保護區(qū)面積的21.0%，在2018年甚至只占到保護區(qū)面積比例的17.7%；對于植被生長季NDVI值在0.2～0.3范圍的面積比例呈現(xiàn)震蕩上升趨勢，在2017年和2018年達到保護區(qū)植被面積的16.7%和16.3%，分析表明保護區(qū)無植被覆蓋區(qū)減少，低植被覆蓋區(qū)和中等植被覆蓋區(qū)增加，植被總體呈現(xiàn)變好趨勢。

圖4 保護區(qū)生長季NDVI各等級面積比例

2.3 時空穩(wěn)定性分析

使用CV來分析2000—2020年植被生長季NDVI的年際間穩(wěn)定性程度，保護區(qū)植被生長季NDVI時空穩(wěn)定性見圖5。較低穩(wěn)定度分布面積最大，占整個保護區(qū)植被面積的47.86%；中等穩(wěn)定度分布面積也較大，占整個保護區(qū)植被面積的38.8%；較低穩(wěn)定度和中等穩(wěn)定度總共占整個保護區(qū)面積的86.66%，表明2000年以來保護區(qū)的NDVI穩(wěn)定性較差，這是可能由于保護區(qū)植被以草地植被和灌叢為主，受氣象條件影響明顯，氣象條件的較小波動就能引起植被狀態(tài)的較大變動。較高穩(wěn)定度和高穩(wěn)定度主要分布于保護區(qū)中東部及西部邊沿。保護區(qū)中東部和南部氣象條件較好，利于植被生長，因此植被生長季NDVI穩(wěn)定性較高；保護區(qū)西部植被生長季NDVI穩(wěn)定性較高的原因可能是該區(qū)域雪山融水較豐富，利于植被生長。

結(jié)合圖2和圖5可以看出，植被生長季NDVI≤0.1的大部分區(qū)域處于較低穩(wěn)定度，表明該區(qū)域易受環(huán)境影響，對外界干擾反應敏感，這與王青霞[25]、徐維新[26]的研究結(jié)果較為一致。

圖5 生長季NDVI穩(wěn)定性空間分布

在各功能區(qū)中，實驗區(qū)、緩沖區(qū)和核心區(qū)高穩(wěn)定度區(qū)域面積占比分別是11.47%、9.48%和9.49%，中等穩(wěn)定度區(qū)域面積占比分別是40.44%、37.65%和39.47%，較低穩(wěn)定度區(qū)域面積占比分別是44.82%、49.39%和47.61%，實驗區(qū)高穩(wěn)定度和中等穩(wěn)定度區(qū)域面積占比均略高于緩沖區(qū)和核心區(qū)，而較低穩(wěn)定度面積占比略低于緩沖區(qū)和核心區(qū)，說明實驗區(qū)植被穩(wěn)定度好于緩沖區(qū)和核心區(qū)。原因可能是實驗區(qū)在保護區(qū)南部和東南部，氣象條件較好，植被穩(wěn)定度高。

2.4 相關(guān)性分析

2000—2020年，保護區(qū)植被生長季平均氣溫變化率為 0.01～0.87℃/10 a（圖 6a），平均值是 0.43℃/10 a，增溫趨勢明顯，東南角至西北角的條帶增溫最明顯，每10年增加0.5～0.6℃，該結(jié)論與杜軍[11]的研究結(jié)果基本一致；植被生長季累積降水量變化率為-3.55～27.52 mm/10 a（圖6b），平均值為9.92 mm/10 a，保護區(qū)中部降水量呈減少趨勢、東部與西部降水呈增加趨勢，其中，東南部達到每10年增加17～27 mm，西部大部分區(qū)域也達到了每10年增加10～17 mm?？傮w來看，保護區(qū)植被生長季氣候趨向于“暖濕化”。

圖6 保護區(qū)2000-2020年生長季氣溫和降水量變化趨勢

逐像元分別計算生長季NDVI與生長季平均氣溫、生長季總降水量的相關(guān)系數(shù)。結(jié)果顯示空間分布差異較為明顯，保護區(qū)大部分區(qū)域生長季NDVI與氣溫呈正相關(guān)性（圖7a），占保護區(qū)總面積的80.9%。相關(guān)系數(shù)集中分布在0～0.5范圍，占保護區(qū)總面積的72.9%。與氣溫呈負相關(guān)性的區(qū)域主要分布在保護區(qū)西南部。保護區(qū)所有像元生長季NDVI均與氣溫達到顯著相關(guān)(P＜0.05)，99.4%的區(qū)域達到極顯著相關(guān)(P＜0.001)。保護區(qū)大部分區(qū)域生長季NDVI與降水量呈正相關(guān)性（圖7b），占保護區(qū)面積的80.7%。相關(guān)系數(shù)集中分布在0～0.5范圍，占保護區(qū)總面積的77.4%。與降水量呈負相關(guān)性的區(qū)域空間分布特征不明顯。保護區(qū)所有像元生長季NDVI均與降水量均達到極顯著相關(guān)(P＜0.001)。表明氣溫與降水量對保護區(qū)生長季NDVI都有重要的正向影響，在保護區(qū)西南部的部分區(qū)域，可能是氣溫升高導致蒸發(fā)加大使生長季NDVI與氣溫呈現(xiàn)出負相關(guān)性。

圖7 保護區(qū)2000—2020年逐年生長季NDVI與氣溫(a)和降水量(b)的相關(guān)性

3 結(jié)論與討論

（1）保護區(qū)多年植被生長季NDVI均值為0.02～0.55，平均值為0.138，呈現(xiàn)由東南向西北逐漸減小的空間分布特征。

（2）保護區(qū)86.39%區(qū)域植被生長季NDVI呈緩慢的上升趨勢，并且大部分達到顯著水平，雖然有13.61%的植被生長季NDVI處于下降趨勢，但只有很少的區(qū)域達到了顯著性水平；植被總體呈現(xiàn)上升變好趨勢。

（3）2000年以來保護區(qū)的生長季NDVI穩(wěn)定性較差，較低穩(wěn)定度和中等穩(wěn)定度占整個研究面積的86.66%，植被生長季NDVI≤0.1的大部分區(qū)域處于較低穩(wěn)定度。

（4）在植被生長季，保護區(qū)氣候趨向于“暖濕化”；氣溫與降水量對保護區(qū)生長季NDVI都有重要的正向影響。

青藏高原植被生長季NDVI在0.1～0.9，多年均值為0.49[1]，青藏高原草地生長季NDVI變化范圍在0.431～0.471[27]，羌塘保護區(qū)植被生長季多年NDVI的平均值為0.138，明顯低于青藏高原整體植被NDVI和青藏高原草地NDVI，說明羌塘保護區(qū)植被狀況較差。

從全球來看，有56.3%的陸地區(qū)域植被NDVI呈增加趨勢，且具有顯著季節(jié)變化趨勢[28]。在中國，NDVI呈增加趨勢的面積大約占53.8%，達到顯著水平的面積占29.3%[29]。在青藏高原，NDVI升高區(qū)域面積約占57.1%[1]。羌塘保護區(qū)植被生長季NDVI呈增加趨勢的面積占植被總面積的86.4%，其中達到顯著水平的面積占保護區(qū)植被面積的45.15%，明顯高于青藏高原整體水平。

綜合NDVI變化趨勢和氣象要素變化趨勢來看，NDVI下降區(qū)域與氣溫升高、降水量減少區(qū)域重疊較大，而NDVI上升區(qū)域主要與氣溫升高、降水量增加區(qū)域重疊較大，徐增讓[30]的研究表明，保護區(qū)居民點和人口在實驗區(qū)，緩沖區(qū)和核心區(qū)基本上無分布，說明保護區(qū)影響植被生長的人為干擾小，水熱條件是保護區(qū)植被生長最重要影響因素，氣候“暖濕化”有利于保護區(qū)植被改善。