川西高原植被時空格局及其對氣候變化的響應

伍良旭, 王 晗, 邵懷勇, 周思彤

(成都理工大學 地球科學學院, 成都 610059)

植被是地球生態系統主要的組成部分，是連接氣候與土壤的紐帶，在全球物質能量交換和氣候變化中具有不可替代的作用[1]。植被動態變化在高原生態環境中具有敏感的指示作用，合理利用區域資源和有效治理生態環境提供依據，已成為生態環境研究的熱點[2]。NDVI是表征植被動態變化的一個重要指標，能較精確地反映地表植被的生長狀態及演變特征，在區域植被退化、濕地景觀格局等領域發揮著重要的作用[3]。川西高原位于我國青藏高原向川西南山地和四川盆地的過渡地區，是我國生態安全戰略格局中青藏高原生態屏障和黃土高原—川滇生態屏障的重要組成部分，地形復雜，且受高原背風坡和西南季風等系統的共同影響，是典型的生態環境脆弱區和生態氣候敏感區[2]，但目前針對該區域植被時空格局與氣候因子的相關研究較少，因此對該地區的植被時空格局及其對氣候變化響應的研究具有重要意義。

此外，現在的遙感數據通常不能同時具備較長的時間序列和較高的空間分辨率，如GIMMS NDVI3g是目前最長時間序列的NDVI數據集之一[4]，前人利用其進行了1982—2009年東北地區[5]、1982—2006年蒙古高原[6]和1982—2013年新疆地區[7]的植被變化研究，但8 km的空間分辨率應用在區域植被變化研究中過于粗糙。雖然如今高分辨率輻射計(AVHRR)、地球觀測系統(SPOT)-VEGETATION，MODIS和ENVISAT中分辨率成像光譜儀(MERIS)等衛星傳感器所采集的數據有較高的時間分辨率，通過合成一定時間內的圖像，能夠減少甚至消除云污染[4]，遙感數據質量不斷提高，使得區域植被趨勢分析在全球范圍內得到廣泛應用[8-10]，但如何提高數據的空間分辨率仍是亟待解決的問題。相對于GIMMS NDVI3g數據精度更高、更能表達植被時序變化的MODIS NDVI產品，卻因其時間覆蓋范圍(2001年7月至今)無法針對長時間序列進行研究。因此，如何提高早期遙感數據的分辨率對于長時間序列的植被變化監測研究顯得尤為重要。

針對遙感數據時空尺度的融合問題，前人做了大量的研究。如Gao等[11]提出的STARFM模型(spatial and temporal reflection fusion Model)可融合Landsat和MODIS地表反射率，生成具有高時空分辨率的合成反射率數據。Singh等[12]利用STARFM融合模型充分提取了農作物8 a生長季節的物候。Bhandari等[13]和Schmidt等[14]評估了STARFM合成長期時間序列數據的能力，并獲取不同植被的物候特征進行研究。Meng等[15]和Zhu等[16]對STARFM算法改進，分別提出了時空融合自適應植被指數模型(STAVFM)和增強型自適應時空融合模型ESTARFM，以更準確地表征地表的時空變化。由于在地形復雜地區，ESTARFM的融合效果比STARFM較好，使得ESTARFM算法被頻繁用于生成高時空分辨率的合成遙感圖像[17-18]。然而，目前較少有學者針對川西高原使用ESTARFM模型將GIMMS NDVI3g與MODIS NDVI進行時空融合并開展1995年至今植被時空格局的相關研究，為此本文基于ESTARFM模型開展川西高原GIMMS NDVI3g和MODIS NDVI數據的時空融合，這對開展其長期的植被變化監測具有重要意義。

本研究嘗試基于ESTARFM模型建立川西高原GIMMS NDVI3g與MODIS NDVI數據的函數關系，以提高GIMMS NDVI3g數據的空間分辨率。在此基礎上分析川西高原長時間的植被變化特征及其對氣候變化的響應，為當地的環境保護決策提供基礎數據和科學參考。

1 研究區概況

川西高原(27.59°—34.31°N，97.36°—104.62°E)地屬青藏高原與四川過渡地帶，覆蓋阿壩州和甘孜州等高原藏區，區域面積約為25.38萬km2。研究區內地勢高，地形起伏大，平均海拔約4 000 m。境內屬亞熱帶高原季風氣候，地域差異明顯，多年降水量約435.6～1 584.1 mm。近年來，由于全球氣候變化以及人類活動的增加，當地植被覆蓋逐漸變化，植被的退化對境內生態系統安全會造成極大威脅[19]，因此，研究分析川西高原植被覆蓋時空分布特征，并掌握主要氣候因素對植被變化的影響非常必要[20]。

2 數據與方法

2.1 數據來源與處理

(1) 遙感數據：MODIS 13Q1植被指數數據集產品，空間分辨率500 m×500 m，時間分辨率30 d。首先進行Savitzky-Golay濾波處理降低異常值等噪聲影響，隨后利用最大值合成法合成2005年和2015年MODIS NDVI數據。該方法能減少云霧等大氣效應對數據的影響，提高數據精度[21]。

1995年、2005年和2015年的GIMMS NDVI3g數據空間分辨率為8 km×8 km，時間分辨率為15 d。采用最大值合成法合成3個時間段的年GIMMS NDVI3g數據，并將其空間分辨率重采樣為500 m×500 m，為后續使用ESTARFM模型進行數據融合做準備。

(2) 氣象數據：年均氣溫、年降水量空間插值數據集產品。薄板光順樣條法采用平滑樣條函數對多變量數據進行分析和插值，即使用函數逼近曲面，是一種較好的氣候要素插值的方法[22]。因此，該數據集采用基于薄板光順樣條法開發的ANUSPLIN插值軟件(Hutchinson，2001)對氣候要素進行插值，其中將氣象站點的經度、緯度和海拔高程作為年降水量空間插值的自變量，而年均溫度空間插值則考慮了溫度隨垂直高度的變化[22]。

MODIS NDVI數據來源于中國科學院計算機網絡信息中心地理空間數據云平臺(http:∥www.gscloud.cn/)；GIMMS NDVI3g數據來源于ECOCAST網站(https:∥ecocast.arc.nasa.gov/data/pub/gimms/3 g.v1/)；氣象數據和高程數據來源于中國科學院資源環境科學數據中心(http:∥www.resdc.cn)；為提高試驗結果精度，轉換所有數據到統一的投影坐標系WGS_1984_Zone_48N。

2.2 研究方法

2.2.1 ESTARFM算法 為生成高時空分辨率的遙感圖像，使用ESTARFM模型融合MODIS NDVI和GIMMS NDVI3g數據。ESTARFM模型是在STARFM模型[11]基礎上根據像元的異質性調整賦權方法，通過設置轉換系數來改進預測結果，異質性越大，預測精度越高，可以保存更多的空間特征細節[23]。

本文根據第1時期m時刻MODIS NDVI與tp時刻的GIMMS NDVI3g得到融合后tp時刻的GIMMS NDVI3g數據，記為mm=(Xw/2,yw/2,tp,B)，或利用第2時期n時刻的觀測數據得到融合后的tp時刻GIMMS NDVI3g數據，記為mn=(Xw/2,yw/2,tp,B)。通過tm和tn與時刻tp間檢測到的GIMMS NDVI3g變化幅度來計算權重，如公式(1)所示[16]：

(k=m,n)

由公式(2)計算中心像元的值：

M(xw/2,yw/2,tp,B)=Tm·Mm(xw/2,yw/2,tp,B)+

Tn·Mn(xw/2,yw/2,tp,B)

(2)

由于篇幅限制，具體融合原理請見參考文獻[16]。統計2005年和2015年融合后的GIMMS NDVI3g和MODIS NDVI對應的像元值繪制散點圖(圖1)。融合的GIMMS NDVI3g數據與MODIS NDVI數據具有高度相似性，各像元值大致相同，兩者基本分布在y=x的趨勢線上(p<0.01)，預測效果較好。

圖1 真實MODIS NDVI與融合后的GIMMS NDVI3g數據散點圖

對比2005年和2015年GIMMS NDVI3g、融合后的GIMMS NDVI3g和真實MODIS NDVI數據細節(圖2)可以看出，融合后的GIMMS NDVI3g與真實MODIS NDVI數據像元分辨率較好，地物紋理清晰，部分特征地區的細節能夠凸現出來，能較好地再現區域植被覆蓋細節特征。

注：A為2005年GIMMS NDVI3g；B為2005年融合后的GIMMS NDVI3g；C為2005年MODIS NDVI；D為2015年GIMMS NDVI3g；E為2015年融合后的GIMMS NDVI3g；F為2015年MODIS NDVI。圖2 融合數據對比圖

2.2.2 一元線性回歸模型 基于融合處理后的GIMMS NDVI3g數據(值域為[-1,1])，使用像元二分模型計算植被覆蓋度，由于篇幅限制，具體的植被覆蓋度計算公式和原理請見參考文獻[24]。

植被覆蓋度變化主要用于分析植被覆蓋度隨時間變化的關系，他們的關系滿足一元線性回歸模型的條件，以局部像元時間變化特征表示整個空間的變化規律[24]，其計算方式見公式(3)：

(3)

式中:n為總監測年數;fi為第i年植被覆蓋度;Slope為植被覆蓋度線性擬合斜率。利用兩者的相關關系來判斷植被覆蓋度時間變化的顯著性,斜率為負表示植被覆蓋度下降,反之則表示植被覆蓋度上升,擬合斜率絕對值越大,變化越明顯[24-25]。

2.2.3 偏相關性分析法 地理系統是一個由多要素組成的復雜系統，偏相關分析可剔除其他變量的影響，只研究一個因素與另一個因素的相關程度[25]。本文基于像元對川西高原1995年、2000年、2015年植被覆蓋度與年降水量、年均氣溫兩個氣候因素進行偏相關性研究，見公式(4)[25]：

(4)

式中:rab,c表示控制變量c的值不變,變量a，b間的偏相關系數。rab，rbc，rbc分別表示變量a和b，b和c，a和c的簡單相關系數。其中a，b和c分別為植被覆蓋度、年均溫度和年降水量。rab,c>0時表示兩變量呈正相關,rab,c<0時表示兩變量呈負相關,rab,c的絕對值越大,兩變量的相關性越大。

簡單相關系數計算見公式(5)[25]：

(5)

3 結果與分析

3.1 植被覆蓋特征分析

3.1.1 植被覆蓋空間分布特征 參照前人研究[26]將川西高原植被覆蓋度分為極低覆蓋度(0～0.2)、低覆蓋度(0.2～0.4)、中覆蓋度(0.4～0.6)、中高覆蓋度(0.6～0.8)和高覆蓋度(0.8～1.0)5個等級，得到川西高原1995年、2005年和2015年植被覆蓋度空間分布圖(圖3)。

圖3 不同等級植被覆蓋度空間分布

分析可得，1995年、2005年和2015年研究區極低、低和中覆蓋度區域主要分布在西部高海拔的白玉、巴塘、理塘和東部低海拔的瀘定、康定、九龍、理縣等縣區，南部的稻城、和西北部的瑪多縣主要為中高覆蓋度，其余地區多屬于高覆蓋度區域。西部地區海拔較高，多為高山植被，西部德格、甘孜一帶和東部瀘定縣區有荒漠和草甸區，其余地區水熱條件良好，主要為針葉林和闊葉林。2005年川西高原高覆蓋度區域面積相對1995年有所減少，中植被覆蓋和中高植被覆蓋區域面積增加。2015年高植被覆蓋區域覆蓋面積繼續降低，但趨勢減緩，中、中高植被覆蓋區域面積增加，且后者相對前者增加明顯。

川西高原整體的植被覆蓋較高，呈現較為明顯的區域差異。高植被覆蓋區主要分布在北部和南部地區，低值區分布在中部地區，此外沿海拔梯度，植被覆蓋呈先升高后降低的垂直分布特點。

3.1.2 植被覆蓋變化趨勢特征 為研究1995—2015年川西高原植被覆蓋的時間變化規律，對植被覆蓋與時間進行一元線性回歸分析(圖4)，根據趨勢變化的方向和程度，劃分為5個等級(表1)[24]。植被覆蓋呈增長趨勢的面積占17.63%，其中明顯增加面積占15.53%，輕度增加面積占2.10%，主要分布在高海拔地區的瑪多、石渠、白玉等縣區，以瑪多縣分布尤為密集。植被覆蓋呈減少趨勢的面積占27.36%，其中明顯減少面積占18.29%，輕度減少面積占9.07%，減少區域廣泛分布在各個區域，以理塘、康定縣區分布較廣，結合前人研究分析認為[3]，理塘縣自20世紀80年代開始，因過度挖采、放牧等人為因素、鼠蟲害、氣候變化、暴雨沖刷等多種因素綜合并長期作用使理塘縣草原退化、沙源裸露導致植被覆蓋減少；康定縣由于亂砍濫伐，開荒種地，旅游開發等人為干擾因素和自然災害等多因素引起大面積土地退化，山坡地水土流失十分嚴重，植被減少趨勢明顯。植被覆蓋基本穩定的區域面積占55.01%，分布廣泛，主要在研究區的東北部的阿壩、若爾蓋縣區，該地區主要為沼澤、草甸等，植被生長狀況穩定，植被變化不明顯。總體而言，1995—2015年川西高原植被覆蓋呈增長趨勢面積占17.63%，呈增長趨勢面積占27.36%，基本穩定面積占55.01%，植被整體變化趨勢穩定。

圖4 川西高原1995-2015年植被覆蓋變化趨勢

表1 1995-2015年川西高原植被覆蓋變化面積及比例

3.2 植被對氣候變化的響應

3.2.1 年均氣溫與年降水量空間分布特征 植被的變化受到多方面因素的影響，其中年降水量和年均氣溫是對植被變化影響較大的兩個氣候因素[25]。計算1995年、2005年和2015年的年均氣溫和年降水量(圖5)。

圖5 1995-2015年年均氣溫和年降水量分布

分析可得，1995年、2005年和2015年川西高原年均氣溫為3.3℃，最高溫度為17.9℃，位于最南部的鹽源縣，最低溫度為-13.1℃，位于西北部的瑪多縣；年均降水量為848.6 mm，最大降水量為1 584.1 mm，位于西南部的瀘定縣；最小降水量為435.7 mm，位于西北部的瑪多縣。研究區內年均氣溫和年降水量的分布具有空間異質性，瑪多縣所在的西北部和瀘定縣所在的東南部區域氣溫整體低于中部地區，瀘定縣所在的東南地區，位于青藏高原和四川盆地過渡區，受東南季風和西南季風影響，降水充沛，氣溫不高；瑪多縣地區靠近青藏高原，從研究區中部到西北部地勢逐漸增高，氣溫逐漸下降。

3.2.2 植被覆蓋對氣候變化的響應 計算1995—2015年年均氣溫和年降水量與植被覆蓋的偏相關系數(圖6A—B)。年植被覆蓋與年均氣溫的偏相關系數變化范圍為[-0.71，0.71]，年植被覆蓋與年降水量的偏相關系數變化范圍為[-0.71，0.71]。統計分析可得年植被覆蓋與年均氣溫呈正相關的區域面積占14.9%，年植被覆蓋與年降水量呈正相關的區域面積占35.4%。根據年均氣溫和年均降水量與年植被覆蓋的偏相關系數方向，對研究區域進行分區討論[27]，年植被覆蓋與年均氣溫負相關，與年降水量負相關，歸于分區一；年植被覆蓋與年均氣溫負相關，與年降水量正相關，歸于分區二；年植被覆蓋與年均氣溫正相關，與年降水量負相關，歸于分區三；年植被覆蓋與年均氣溫正相關，與年降水量正相關，歸于分區四；年植被覆蓋與年均氣溫和年降水量的偏相關系數較低(都為0)，歸于分區五(圖6C)。

圖6 1995-2015年年均氣溫和年降水量與植被覆蓋度偏相關分析及偏相關性分區

分區一占總面積的6.4%，主要分布在石渠縣、甘孜縣和色達縣的北部，植被類型主要為高寒草甸，分析認為該區域鼠災泛濫和人為過度放牧導致植被覆蓋稀疏，氣溫和降水反而加速了水土流失，抑制植被生長。分區二占總面積的30.8%，主要分布在西部的白玉縣、理塘縣和東北部的黑水縣等地，植被類型主要為草甸、常綠闊葉林、針葉林和灌叢，地形復雜，氣溫較高而降水量低，水土蒸發量大，水土流失較嚴重[28]，氣溫升高減弱其光合作用，抑制植被的生長，適當的降水可增加土壤含水量，促進植被生長。分區三占總面積的10.3%，主要分布在馬爾康縣、道孚縣、九龍縣一帶，該區域氣溫相對較低而降水量高，降水持續增加引發泥石流、滑坡等地質災害，并使植被光照條件下降，減弱其光合作用，抑制植被生長，從而使植被覆蓋降低[29-30]。分區四占總面積的4.6%，主要分布在西北部的瑪多縣及東南部的汶川縣和康定縣地區。汶川地區自2008年汶川地震后，該區域的地質活動較為頻繁，生態系統也在逐漸恢復，氣溫和降水作為植被生長的兩個重要因素，促進該區域生態重建[31]。瑪多縣位于高海拔地區，氣溫低、降水量少、人類活動影響程度小，氣溫和降水成為植被覆蓋變化的主要誘因[32]。分區五占總面積的43.8%，主要分布在若爾蓋縣和紅原縣等東北部高寒高緯度地區，土地沙化嚴重，植被變化的主導驅動因素是人類活動，因此，該區域的年植被覆蓋變化與年均氣溫和年降水量的相關性不大[33]。

在研究區不同的區域，年植被覆蓋與年均氣溫的偏相關性和年植被覆蓋與年均降水量的偏相關性的方向與程度不同，整體而言，植被覆蓋對降水的變化相較于氣溫更加敏感，此外研究區內年植被覆蓋與年降水量呈負相關的區域，往往與年均氣溫呈正相關，這與仙巍等[20]研究川西高山高原過渡帶植被對氣候變化響應的結果相符。

4 結 論

(1) 運用ESTARFM模型融合的2005年和2015年的GIMMS NDVI3g數據與對應的真實MODIS NDVI數據吻合度較高，能凸顯地區特征，融合結果可用于地形復雜的川西高原植被變化的研究，為早期的大尺度植被變化研究提供了一種新的思路。

(2) 1995年、2005年和2015年川西高原植被覆蓋結構整體較好，空間上呈現沿海拔先升高后降低的垂直分布特點，時間上植被變化趨勢穩定。植被覆蓋度較低的區域主要分布在低海拔區，近20 a來研究區植被覆蓋度整體變化趨勢比較穩定，其中以高覆蓋度區域面積的逐年減少和較低覆蓋度區域面積的逐年增加趨勢較為明顯。

(3) 1995年、2005年和2015年川西高原年均氣溫和年降水量存在空間異質性，氣溫空間差異表現為研究區北部氣溫低于南部、西部氣溫低于東部，降水空間差異表現為研究區西部降水量低于東部。

(4) 年植被覆蓋與年均氣溫呈正相關的區域面積占14.9%，年植被覆蓋與年均降水量呈正相關的區域面積占35.4%，此外研究區內年植被覆蓋與年降水量呈負相關的區域，往往年植被覆蓋與年均氣溫呈正相關，植被覆蓋度對降水的變化相較于氣溫的變化更加敏感。

本研究中使用的數據時間間隔較長，可能會忽略一些細節。此外，由于數據獲取等原因，僅選取了年均氣溫和年降水量兩個主要的氣候要素與植被進行響應分析，未考慮地表因素和人為因素等其他因素對植被的響應。未來的研究將會考慮更短的時間間隔和更多因素對植被變化的影響研究。