雅魯藏布江流域植被覆蓋變化及其對(duì)氣候變化的響應(yīng)

韓先明，左德鵬，李佩君，徐宗學(xué)，高曉曦

(1. 北京師范大學(xué)水科學(xué)研究院城市水循環(huán)與海綿城市技術(shù)北京市重點(diǎn)實(shí)驗(yàn)室，北京 100875;2.水利部宣傳教育中心，北京 100053)

植被作為重要的自然資源，是人類生存環(huán)境中不可或缺的組成部分。植被覆蓋與氣象要素密切相關(guān)，通過水汽輸送的方式，植被適應(yīng)著氣候變化并在一定程度上表征氣候變化。同時(shí)，植被參與并影響全球水循環(huán)、生化循環(huán)、能量平衡等過程，對(duì)生態(tài)系統(tǒng)變化產(chǎn)生深遠(yuǎn)的影響。植被覆蓋是連接氣象和水文的載體，通過生態(tài)系統(tǒng)的反饋可對(duì)其他生態(tài)要素產(chǎn)生不同程度的作用。植被指數(shù)在評(píng)價(jià)植被覆蓋度、生長活力和生物量等植被狀態(tài)信息時(shí)發(fā)揮著重要作用，其中歸一化植被指數(shù)[1-3](normalized difference vegetation index,NDVI)是表征植被生長狀態(tài)和植被空間分布密度的指示因子，能夠很好地反映植被的動(dòng)態(tài)變化，在氣象、水文等領(lǐng)域的研究中被廣泛應(yīng)用。

雅魯藏布江作為我國西南地區(qū)最大的外流河，地處青藏高原腹地，區(qū)位海拔跨度極大。流域內(nèi)植被及生態(tài)環(huán)境十分脆弱，是我國重要的地理屏障，因此研究雅魯藏布江流域植被變化及水文過程，提高對(duì)該地區(qū)植被覆蓋變化的認(rèn)識(shí)，厘清植被覆蓋變化同氣象要素的關(guān)系對(duì)我國的水能開發(fā)、環(huán)境保護(hù)、地緣政治有著重大戰(zhàn)略意義[4-9]。許多學(xué)者[10-14]在不同地區(qū)對(duì)NDVI與氣象要素的關(guān)系進(jìn)行了探究，姜琳等[15]研究了降雨因子對(duì)雅魯藏布江流域植被生長的影響；姜紅濤等[16]利用NDVI、降水和氣溫?cái)?shù)據(jù)，分析了艾比湖流域NDVI變化趨勢(shì)及其對(duì)降水和氣溫的響應(yīng)特征，指出NDVI總體上與年降水量呈正相關(guān)，而與平均氣溫呈負(fù)相關(guān)。樸世龍等[17]發(fā)現(xiàn)近年來我國大部分地區(qū)NDVI呈增長趨勢(shì)，西北干旱區(qū)域和青藏高寒區(qū)域NDVI增長率在夏季達(dá)到最大，東部季風(fēng)區(qū)域NDVI增長率在春季達(dá)到最大。大多研究表明NDVI對(duì)降水的響應(yīng)存在一定滯后性[18-19]，同降水及氣溫存在正相關(guān)關(guān)系[20]；NDVI與海拔也存在一定關(guān)系，一般來說低海拔地區(qū)NDVI值相對(duì)較高[21]。就雅魯藏布江流域而言，Li等[22]通過遙感手段對(duì)1999—2013年雅魯藏布江流域NDVI年際變化及其沿海拔梯度變化特征進(jìn)行了分析，發(fā)現(xiàn)以闊葉林和針葉林為主的低海拔地區(qū)植被覆蓋增加明顯。陳斌[23]利用SPOT數(shù)據(jù)集對(duì)雅魯藏布江大峽谷自然保護(hù)區(qū)不同植被類型、海拔梯度下的NDVI時(shí)空變化規(guī)律進(jìn)行了研究。呂洋等[24]分析了雅魯藏布江流域植被覆蓋時(shí)空變化規(guī)律，發(fā)現(xiàn)NDVI與降水之間存在較為明顯的正相關(guān)關(guān)系。

已有研究對(duì)雅魯藏布江流域植被覆蓋變化進(jìn)行了相關(guān)分析，然而不同植被類型下的植被生長狀況及其持續(xù)性分析還未做充分探討。考慮到研究區(qū)測(cè)站稀疏，屬于缺資料地區(qū)，本文基于MODIS高精度遙感影像，采用趨勢(shì)分析、Hurst指數(shù)、偏相關(guān)性分析對(duì)研究區(qū)2000—2016年生長季、年內(nèi)和年際植被覆蓋空間異質(zhì)性做詳細(xì)分析，探討植被覆蓋對(duì)氣候變化的響應(yīng)，以期為研究區(qū)生態(tài)環(huán)境保護(hù)和氣候變化研究提供參考。

圖1 研究區(qū)概況及氣象站點(diǎn)分布

1 研究區(qū)概況和數(shù)據(jù)

1.1 研究區(qū)概況

雅魯藏布江流域[25]位于北緯27°8′～31°2′，東經(jīng)82°0′～97°1′之間，地處我國西藏自治區(qū)境內(nèi)，是世界上平均海拔最高的流域。流域發(fā)源于西藏西南部喜馬拉雅山脈北麓的杰馬央宗冰川，上游稱為馬泉河[26]。流域海拔高度絕大部分在3 000 m以上，從源頭杰馬央宗冰川到下游河谷,由干旱、半干旱半濕潤、濕潤地區(qū)組成。流域東西狹長、南北窄短，南有喜馬拉雅山脈，北有岡底斯山及念青唐古拉山脈，南北之間為藏南谷地，在米林以下的下游地區(qū)分布著著名的大拐彎地帶[27]。流域降水年際變化不大，多年平均降水量300～500 mm，但空間分布差異極大且年內(nèi)分配不均，部分地區(qū)夏季降水可達(dá)全年降水量80%。

通過咨詢西藏水利專家并查閱相關(guān)文獻(xiàn)[28]，將雅魯藏布江流域分為上、中、下游，上游為日喀則市拉孜縣以上寬谷區(qū)，植被覆蓋多以草甸和草原為主；中游為河谷區(qū)，植被覆蓋多以高山植被和草甸為主；米林縣派鎮(zhèn)以下為下游峽谷區(qū)，植被覆蓋多以灌叢、針葉林和闊葉林為主。研究區(qū)水文氣象站點(diǎn)主要位于流域中、下游且分布不均，屬于缺資料流域，氣象站分布、高程及流域概況如圖1所示。

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

1.2.1植被數(shù)據(jù)

NDVI數(shù)據(jù)由美國航空航天局NASA提供的MODIS(moderate resolution image spectroradiometer)遙感產(chǎn)品(MOD13A3)(https://modis.gsfc.nasa.gov/)處理所得，已經(jīng)過質(zhì)量控制，該產(chǎn)品時(shí)長為2000年2月至2016年12月。MOD13A3是MODIS的三級(jí)產(chǎn)品，是在二級(jí)產(chǎn)品的基礎(chǔ)上，對(duì)由傳感器成像過程產(chǎn)生的邊緣畸變(Bowtie效應(yīng))進(jìn)行校正所得。該產(chǎn)品空間分辨率為1 km×1 km，時(shí)間分辨率為月，原始數(shù)據(jù)為EOS-HDF格式，有效值域范圍為-0.2～1.0，數(shù)值越大說明植被覆蓋狀況越好。數(shù)據(jù)通過MRT軟件批量處理成GEO TIFF格式，并通過ArcGIS投影拼接裁剪處理，其中合成方法選用最大合成法。植被類型數(shù)據(jù)來源于國家自然科學(xué)基金委員會(huì)中國西部環(huán)境與生態(tài)科學(xué)數(shù)據(jù)中心(http://westdc.westgis.ac.cn)，空間分辨率為1∶100萬。

1.2.2氣象數(shù)據(jù)

氣象數(shù)據(jù)包括逐月降水和平均氣溫，均來自中國氣象科學(xué)數(shù)據(jù)共享服務(wù)網(wǎng)(http://data.cma.cn/)的中國地面氣候資料月值數(shù)據(jù)集。該數(shù)據(jù)集已通過質(zhì)量控制，月年值各要素?zé)o錯(cuò)誤數(shù)據(jù)。主要選取雅魯藏布江流域內(nèi)15個(gè)氣象站點(diǎn)(當(dāng)雄、拉孜、南木林、日喀則、尼木、拉薩、墨竹貢卡、瓊結(jié)、澤當(dāng)、江孜、嘉黎、波密、加查、林芝、米林)自2000年2月至2016年12月的逐月降水和平均氣溫?cái)?shù)據(jù)，經(jīng)去趨勢(shì)處理后用于和NDVI進(jìn)行偏相關(guān)分析。研究區(qū)氣象站點(diǎn)分布如圖1所示。

2 研究方法

2.1 偏相關(guān)系數(shù)法

偏相關(guān)系數(shù)[29]是考察多變量相關(guān)程度時(shí)采取的研究方法，是在控制其他變量的條件下，衡量其中某兩個(gè)變量直接線性相關(guān)關(guān)系的指標(biāo)。其中，一階偏相關(guān)系數(shù)的計(jì)算公式為

(1)

式中：r12-3為控制第3個(gè)變量的情況下第1個(gè)變量和第2個(gè)變量的相關(guān)系數(shù)；r12為第1個(gè)變量和第2個(gè)變量的相關(guān)系數(shù)；r13為第1個(gè)變量和第3個(gè)變量的相關(guān)系數(shù)；r23為第2個(gè)變量和第3個(gè)變量的相關(guān)系數(shù)。r12-3越接近1，表明偏相關(guān)性越高。

圖2 2000—2016年生長季平均NDVI空間分布

2.2 Hurst指數(shù)法

基于重標(biāo)極差(R/S)分析方法的Hurst指數(shù)[30]是定量描述時(shí)間序列信息長期依賴性的有效方法，該方法是由英國水文學(xué)家Hurst對(duì)尼羅河進(jìn)行的長期水文觀測(cè)歸納發(fā)明，其基本原理是：給定一時(shí)間序列{ξt}，對(duì)于任意正整數(shù)τ，定義均值序列：

(2)

(3)

(4)

式中：〈ξ〉τ為給定時(shí)間序列的均值序列；Xt,τ為時(shí)間序列與其均值序列的累計(jì)離差；Rτ為累計(jì)離差的極差；Sτ為時(shí)間序列的標(biāo)準(zhǔn)差。若存在R/S∝τH，則說明時(shí)間序列存在Hurst現(xiàn)象，H值稱為Hurst指數(shù)，其值可在雙對(duì)數(shù)坐標(biāo)系(lnτ, ln(R/S))中用最小二乘法擬合得到。Hurst指數(shù)主要有3種形式：①當(dāng)0.5

2.3 一元線性回歸法

一元線性回歸法是根據(jù)自變量和因變量的相關(guān)關(guān)系，建立兩者的線性回歸方程對(duì)因變量進(jìn)行分析的方法[31-32]。一元線性回歸法的模型為

Yt=a+bΧt

(6)

(7)

(8)

式中：Χt為自變量時(shí)間序列；Yt為因變量時(shí)間序列；a和b為一元線性回歸方程的參數(shù)；n為系列長度。當(dāng)b>0時(shí)，自變量與因變量呈正線性相關(guān)；當(dāng)b<0時(shí)，自變量與因變量呈負(fù)線性相關(guān)；當(dāng)b=0時(shí)，自變量與因變量無線性相關(guān)關(guān)系。

3 雅魯藏布江流域NDVI特征分析及其對(duì)氣候變化的響應(yīng)

3.1 植被覆蓋空間分布特征

3.1.1NDVI空間分布特征

從基于最大值合成法計(jì)算的2000—2016年雅魯藏布江流域生長季NDVI空間分布(圖2)可以看出，NDVI高值區(qū)主要分布于流域中游和下游地區(qū)，尤其是下游地區(qū)，米林、波密等站附近NDVI值大多在0.8以上；而拉孜以上上游地區(qū)由于海拔較高，不少地區(qū)冰雪終年覆蓋，罕有植被生長，故多為NDVI低值區(qū)，基本小于0.4。和圖1對(duì)比可以發(fā)現(xiàn)高海拔地區(qū)NDVI值普遍較低，低海拔地區(qū)NDVI值保持穩(wěn)定高值，說明植被覆蓋分布特征與海拔存在有一定相關(guān)關(guān)系。

3.1.2NDVI沿高程變化特征

以250 m為分度值對(duì)雅魯藏布江流域NDVI值沿高程變化特征進(jìn)行分析(圖3)?？梢钥闯?，NDVI值呈現(xiàn)一定程度的海拔分異性，隨海拔升高，NDVI值顯著下降。中低海拔地區(qū)為NDVI高值區(qū)，海拔高度為1 000～2 500 m范圍內(nèi)NDVI值基本保持同一水平，均值為0.91。2 500 m以上地區(qū)NDVI值逐步下降，在3 000 m附近呈微弱上升趨勢(shì)之后急劇下降。6 000 m以上NDVI值較低，均值為0.03，且隨海拔高度增加NDVI值基本保持不變。

圖3 NDVI沿高程分布變化

圖4 2000—2016年NDVI 10 a變化趨勢(shì)空間分布

3.2 植被覆蓋時(shí)間趨勢(shì)特征

3.2.1NDVI年際變化分析

雅魯藏布江流域狹長且空間差異性較大，不同地區(qū)NDVI年際變化差別顯著。對(duì)研究時(shí)段內(nèi)流域NDVI年際變化趨勢(shì)進(jìn)行空間分析(圖4)，可以看出，流域內(nèi)NDVI年際變化空間異質(zhì)性顯著，流域上游和下游地區(qū)年際NDVI變化趨勢(shì)中增加比例分別占56.82%和58.60%，而流域中游地區(qū)NDVI變化趨勢(shì)中增加比例僅占42.07%。流域上游和下游地區(qū)NDVI 10 a變化趨勢(shì)在-1%～1%之間的比例分別為73.17%和58.70%，而流域中游地區(qū)NDVI 10 a變化趨勢(shì)在-1%～1%之間的比例為51.72%。結(jié)果表明流域上游和下游地區(qū)植被覆蓋以改善為主且絕大部分區(qū)域變化輕微，流域中游地區(qū)植被覆蓋以惡化為主且部分區(qū)域惡化程度較重。雅魯藏布江流域人類活動(dòng)對(duì)下墊面植被變化的影響不容小覷，下一步可以聯(lián)系土地利用變化等因素進(jìn)一步探討其驅(qū)動(dòng)機(jī)制，在更長時(shí)間尺度上進(jìn)行植被覆蓋對(duì)氣溫和降水的時(shí)間滯后性分析，探究下墊面變化對(duì)濕度和蒸散量等更多氣象要素的響應(yīng)關(guān)系，并分析研究結(jié)果的不確定性，并繼續(xù)補(bǔ)充數(shù)據(jù)源，例如采用SPOT VGT和AVHRR遙感數(shù)據(jù)，應(yīng)用多源遙感數(shù)據(jù)時(shí)空融合技術(shù)研究下墊面長期時(shí)空演變特征。

分析2000—2016年雅魯藏布江流域上游、中游、下游地區(qū)NDVI年均值變化(圖5)，可以看出，NDVI總體呈上升趨勢(shì)，增長率為0.002%/10a(顯著性水平p<0.001)，表明近17 a流域植被覆蓋總體呈緩慢好轉(zhuǎn)趨勢(shì)，尤其是中游和下游地區(qū)NDVI自2015年后出現(xiàn)顯著增長。NDVI自上游、中游至下游呈依次增大趨勢(shì)，NDVI均值分別為0.21、0.36和0.55，全流域NDVI均值為0.41，這與上中下游地區(qū)海拔差異有一定關(guān)系。上游大部分地區(qū)位于雪線以上或山谷邊緣，坡度變化劇烈，植被生長稀疏；而下游地區(qū)海拔較低，氣溫較高，降水豐沛，在適宜的雨熱條件下植被生長旺盛。因此，可以看出雅魯藏布江流域植被分布因地形和海拔分布不均，且具有很大空間差異性。

圖5 2000—2016年NDVI年際變化

3.2.2NDVI年內(nèi)變化分析

對(duì)2000—2016年各月NDVI均值進(jìn)行計(jì)算后的空間分布如圖6所示，為了更加清楚地展現(xiàn)2000—2016年各月NDVI均值的最大值和最小值空間分布情況，制作最值統(tǒng)計(jì)圖如圖7所示。從圖6和圖7可以看出：2月NDVI低值區(qū)面積最大，中上游NDVI值均小于0.4，僅有流域下游植被覆蓋面積較大，NDVI值維持在0.9左右；8月NDVI高值區(qū)面積最大，中上游NDVI值均在0.4至0.5之間，流域中游植被覆蓋面積達(dá)年內(nèi)最大值，維持在0.7左右。

圖6 2000—2016年NDVI年內(nèi)變化空間分布

圖7 2000—2016年NDVI年內(nèi)最值空間分布

分別統(tǒng)計(jì)研究區(qū)上游、中游、下游和全流域NDVI年內(nèi)變化如圖8所示。可以看出NDVI值隨月份變化呈先增大、后減小趨勢(shì)，全流域1—3月NDVI值維持在0.22左右，在3月達(dá)最小值0.21，之后穩(wěn)步上升在7月達(dá)最大值0.48，9月之后下降較為明顯，在12月下降為0.25。流域內(nèi)各分區(qū)和全流域年內(nèi)變化趨勢(shì)基本一致，但下游地區(qū)NDVI值在1—3月呈明顯下降趨勢(shì)。

圖8 2000—2016年NDVI年內(nèi)變化

3.2.3NDVI季節(jié)變化分析

為了準(zhǔn)確反映各季節(jié)植被生長情況，按照春季(3—5月)、夏季(6—8月)、秋季(9—11月)和冬季(12—2月)對(duì)雅魯藏布江流域各季NDVI采用最大值合成法進(jìn)行統(tǒng)計(jì)分析，結(jié)果見圖9。從圖9可以看出，2000—2016年間流域夏季NDVI值最大，均值達(dá)0.40，其中2013年的0.42和2010年的0.37分別為最大和最小值；其次為秋季，NDVI均值為0.36，在2000—2005年間呈緩慢上升趨勢(shì)，并在2005年達(dá)最大值0.39，此后波動(dòng)下降，在2008年達(dá)最小值0.33，2009年后基本維持在0.35左右；春季和冬季NDVI值相對(duì)較小，均值均為0.21。夏季和秋季NDVI呈下降趨勢(shì)，10 a變幅均為-0.000 3%(夏季p=0.876，秋季p=0.897)，而春季和冬季NDVI呈上升趨勢(shì)，10 a變幅分別為0.003%(p=0.131)和0.002%(p=0.097)。

圖9 2000—2016年各季NDVI變化過程

3.2.4NDVI變化持續(xù)性分析

由上文分析可以發(fā)現(xiàn)雅魯藏布江流域NDVI變化具有明顯的周期性，因此有必要對(duì)該流域NDVI變化的持續(xù)性進(jìn)行研究。通過對(duì)NDVI月值逐格點(diǎn)進(jìn)行Hurst指數(shù)計(jì)算，結(jié)果空間分布如圖10所示，相關(guān)統(tǒng)計(jì)信息如表1所示。

圖10 2000—2016年逐月NDVI的Hurst指數(shù)空間分布

表1 雅魯藏布江流域逐月NDVI持續(xù)性統(tǒng)計(jì)

從圖10可以發(fā)現(xiàn)，NDVI的Hurst指數(shù)分布具有一定的空間差異性，在流域中游墨竹貢卡附近拉薩河流域及流域上游地區(qū)空間持續(xù)性較強(qiáng)，表明未來與過去變化趨勢(shì)趨同概率較高。全流域Hurst指數(shù)均值為0.51，說明雅魯藏布江流域NDVI呈弱持續(xù)性，未來與過去的變化趨勢(shì)趨同概率較低。如表1所示，研究區(qū)呈弱持續(xù)性比例最大，達(dá)43.93%，流域面積為105 633 m2；其次是弱反持續(xù)性，面積為85 826 m2，占35.69%；強(qiáng)持續(xù)性與強(qiáng)反持續(xù)性所占比例最小，分別為12.69%和7.69%，所占面積分別為30 526 m2和18 498 m2。

3.3 不同植被類型NDVI時(shí)空變化特征

雅魯藏布江流域植被類型數(shù)據(jù)來源于寒區(qū)旱區(qū)數(shù)據(jù)中心提供的1∶100萬中國植被圖集，通過對(duì)不同植被類型年際及年內(nèi)NDVI進(jìn)行疊加并采用最大值合成法進(jìn)行統(tǒng)計(jì)分析，結(jié)果如圖11所示。由圖11可知，不同植被類型年際和年內(nèi)分布從大到小依次為闊葉林、草叢、針葉林、栽培植被、灌叢、草甸、草原、高山植被和其他，均值分別為0.88、0.86、0.76、0.56、0.53、0.37、0.33、0.27和0.13。其中，2009年栽培植被和草原NDVI值較其他年份小，灌叢和高山植被NDVI值較其他年份大。植被非生長期栽培植被和草原NDVI值趨同，生長期開始后，4—7月草原和闊葉林NDVI值趨同，7—10月栽培植被和灌叢NDVI值趨同，這體現(xiàn)生長期草原和闊葉林、栽培植被和灌木分布差距較小的特點(diǎn)，草叢、灌叢相較于闊葉林、栽培植被等植被類型受植物節(jié)律影響，年內(nèi)變化顯著的特點(diǎn)。研究表明[20]高海拔地區(qū)不同高程植被覆蓋與植被類型有關(guān)，通常在3 600 m左右植被覆蓋開始出現(xiàn)急劇下降。對(duì)流層平均氣溫直減率為0.65℃/100 m，3 600 m以上絕大部分地區(qū)平均溫度在0℃以下，植被類型從闊葉林針葉林向草原荒漠過渡可能是不同海拔高度植被覆蓋度差異的原因。

圖11 2000—2016年不同植被類型NDVI年際和年內(nèi)變化過程

圖12 NDVI與氣象要素相關(guān)性分布

3.4 植被覆蓋對(duì)氣候變化的響應(yīng)關(guān)系

在雅魯藏布江流域植被與氣象要素的相關(guān)性分析中，為了探討區(qū)域NDVI與站點(diǎn)氣象要素的關(guān)系，根據(jù)研究區(qū)氣象站點(diǎn)位置，參考已有研究[33-37]，利用ArcGIS提取以各站點(diǎn)為中心，半徑為1.5 km的區(qū)域NDVI平均值，并對(duì)2000—2016年NDVI、氣溫和降水?dāng)?shù)據(jù)進(jìn)行去除趨勢(shì)處理，然后分別計(jì)算NDVI與氣溫和降水的偏相關(guān)系數(shù)，結(jié)果如圖12所示。NDVI與氣溫的偏相關(guān)系數(shù)在0.27～0.85之間，均值為0.47；而NDVI與降水的偏相關(guān)系數(shù)在-0.28～0.57之間，均值為0.28，且通過了顯著水平為0.05的顯著性檢驗(yàn)。下游同中上游相比，NDVI和氣溫的偏相關(guān)系數(shù)較降水大，說明流域下游地區(qū)植被覆蓋受氣溫影響更大。全流域NDVI與氣溫的平均偏相關(guān)系數(shù)大于降水，說明整體來說雅魯藏布江流域植被生長受氣溫影響要強(qiáng)于降水。

4 結(jié) 論

a. 流域植被覆蓋自上游、中游至下游呈遞增趨勢(shì)，海拔2 500 m以上植被覆蓋開始下降，至6 000 m以上后保持在極低水平。流域上游和下游植被覆蓋以改善為主且絕大部分區(qū)域變化輕微，流域中游植被覆蓋以惡化為主且部分區(qū)域惡化程度較重。

b. 植被覆蓋在2月低值區(qū)面積最大，8月高值區(qū)面積最大。NDVI夏季最大，秋季次之，春季和冬季相對(duì)較小。夏季和秋季NDVI呈下降趨勢(shì)，而春季和冬季呈上升趨勢(shì)。

c. 流域Hurst指數(shù)空間差異性顯著，均值為0.51,說明流域植被覆蓋整體呈弱可持續(xù)性，墨竹貢卡附近拉薩河流域及流域上游空間持續(xù)性較強(qiáng)。整體來說流域植被覆蓋受氣溫影響要強(qiáng)于降水，下游最為明顯。

d. 不同植被類型NDVI年際和年內(nèi)分布從大到小依次為闊葉林、草叢、針葉林、栽培植被、灌叢、草甸、草原、高山植被和其他。生長期草原和闊葉林、栽培植被和灌木分布差距較小，非生長期栽培植被和草原分布較為相近。