內蒙古大青山國家級自然保護區植被歸一化指數時空變化及其與環境因子的關系

胡爾查,王曉江,鐵 牛,洪光宇,蘇 和,張艷楠

1 內蒙古自治區林業科學研究院,呼和浩特 010010 2 內蒙古大青山森林生態系統國家定位觀測研究站,呼和浩特 010010 3 內蒙古農業大學,呼和浩特 010010 4 內蒙古自治區林業和草原監測規劃院,呼和浩特 010010

植被作為陸地生態系統的主體,在全球氣候變化加劇,陸地生態系統發生顯著變化的大背景下[1—3],對環境變化的響應十分敏感[4—6],常被作為表征生態環境變化的綜合指示器,很大程度上代表著區域環境的總體狀況[7—9]。植被時空變化受地形地貌、氣候、人類活動等多重因素的影響[10—13],植被變化及其驅動因素之間的響應機制是研究全球變化與陸地生態系統關系中的重要內容, 成為諸多學者所關注的熱點研究問題[14—17]。

歸一化植被指數(Normalized Difference Vegetation Index,NDVI)能較好地反映區域植被覆蓋與植被生長狀況,且在時效、尺度方面都具有明顯優勢,通常被用來進行區域尺度的植被評價和植被變化研究,并在全球范圍內得到廣泛應用[18—22]。MODIS NDVI空間分辨率較高,同時具備高時間分辨率,適用于較大尺度上研究植被的動態變化。近年來,借助遙感NDVI來研究較長時間序列的區域植被覆蓋變化的成果較為豐富,主要集中于植被覆蓋的整體時空變化特征、植被對地形地貌[13,23]、氣候變化[6,24]、人類活動[20,25]的響應等方面的內容,研究方法也不斷改進,常用的方法有趨勢分析法、Sen+Manna-Kendall非參數檢驗方法、Hurst指數法、線性相關分析、殘差分析法、TSS-RESTREND(Time Series Segmentation and Residual Trend)算法等。利用NDVI時序數據研究植被覆蓋變化,對評價區域經濟建設與生態環境保護具有重要的現實意義,可以為全球變化提供重要的理論依據[26—28]。

自然保護區以其特殊的功能和代表性的自然生態系統對促進國民經濟持續發展和科技文化事業發展中發揮著重要作用。內蒙古大青山國家級自然保護區處于半干旱典型草原區,是陰山山地生態系統最完整、生物多樣性最集中的區域,在維護和保持區域生態系統穩定性,涵養水源、保持水土,屏護山前河套平原乃至華北平原具有重要的意義。由于自然保護區的植被變化反映了保護區的總體保護效果,決定了自然保護區保護功能的發揮程度,自然保護區植被變化及其驅動因素研究已成為保護區研究的重點內容[29]。通過對保護區生長季NDVI時空動態變化及其主要影響因素的研究,不僅可以直觀反映保護區環境變化趨勢,也為進一步采取相關措施提供重要依據。2000年內蒙古大青山自然保護區管理局成立,2008年1月晉升為國家級自然保護區,但迄今為止,鮮有學者對大青山自然保護區成立以來的植被變化及其影響因素進行系統研究。大青山植被垂直帶普明顯,地形復雜多樣,因此,本研究從海拔、坡度、坡向等地形因素角度對大青山自然保護區建立以來保護區生長季NDVI時空變化規律進行研究,并探討了NDVI與氣溫和降水量的相關關系,以期在氣候變化大背景下為如何更好地開展保護區植被保育、保護區管理等方面提供科學依據和信息參考。

1 研究區概況

內蒙古大青山國家級自然保護區位于陰山山脈中段,西起包頭市昆都侖河,北與包頭市固陽縣、呼和浩特市武川縣相接,東至烏蘭察布市卓資縣頭道北山山脊,南連土默川平原,地理坐標為109°47′—112°17′E, 40°34′—41°14′N。保護區東西長度217km,南北平均寬度18km,保護區總面積為38.89萬hm2,主要保護對象為邊緣物種群落為代表的山地森林和瀕危珍稀物種等(圖1)。大青山山地為塊狀斷裂的中等高度山地,山體呈東西走向,山地南北的地貌形態非常不對稱,北部比較平緩,與蒙古高原之間沒有明顯的分界線,南坡則以巨大的正斷面與黃河平原即土默川平原截然分開。保護區處于半干旱的典型草原地帶,屬典型的大陸性半干旱季風氣候,年降水量約360mm左右,年均溫約5.9℃左右。保護區植被具有明顯的垂直分布特征,自下而上由三個垂直帶組成: 山麓草原帶、山地森林灌叢帶、山地草甸帶[30]。

2 數據來源與研究方法

2.1 數據來源與預處理

NDVI數據來源于美國航空航天局LAADS(http://ladsweb.modaps.eosdis.nasa.gov/search/)的MOD13Q1植被指數數據產品,空間和時間分辨率分別為250m和16d,時間范圍為2000—2020年。利用MRT、ArcGIS等軟件進行格式和投影轉換、剪裁、重采樣等處理,獲取研究區NDVI數據集,采用最大值合成法(Maximum Value Composite,MVC)合成月NDVI數據,消除云層、大氣和太陽高度角的干擾,最后利用均值法獲取逐年平均的NDVI數據,消除極端年份氣候異常對植被生長產生的影響[31]。根據研究區氣候和植被生長特征,本文選取每年生長季內4—10月份的NDVI數據來合成年均NDVI數據,來表征該年植被覆蓋狀況。

氣象數據來源于研究區周邊26個氣象站點的月降水量和月平均溫度數據,采用反距離權重內插法(Inverse Distance Weighted,IDW),對年降水量和年均溫度數據進行插值并生成與NDVI像元大小一致的2000—2020年時間序列的氣象柵格數據,本研究采用年降水量和年均溫度作為影響植被生長狀況的氣候因子,分析降水量和氣溫對植被變化的影響。

研究區數字高程模型(DEM)數據來源于中國科學院資源環境科學數據中心(http://www.resdc.cn/DOI)的SRTM(Shuttle Radar Topography Mission)數據,空間分辨率為90m,利用ArcGIS軟件進行裁剪、投影變換、重采樣等處理,重采樣后DEM數據空間分辨率為250m×250m,與NDVI空間分辨率一致。運用ArcGIS表面分析功能,提取大青山自然保護區海拔、坡度和坡向數據。

2.2 研究方法

2.2.1趨勢分析

采用一元線性回歸方法擬合每個像元NDVI隨時間的變化趨勢,計算公式如下:

(1)

式中,θslope為各像元NDVI變化趨勢斜率；i為年序號；n為監測時段累計年數,即21年；NDVIi為第i年的NDVI值。趨勢的顯著性進行雙尾t檢驗,根據檢驗結果將變化趨勢分為5個等級：極顯著增加(P<0.01,θslope>0),極顯著減少(P<0.01,θslope<0),顯著增加(0.010),顯著減少(0.010.05)。

2.2.2變異系數

變異系數的大小反映數據系列的波動程度,通過計算基于像元的NDVI變異系數來描述植被覆蓋度的穩定性,其計算公式如下：

(2)

2.2.3相關分析

本文基于像元尺度,分析NDVI與降水量和平均氣溫的相關性,采用Pearson相關系數表示兩因子之間的相關性。相關系數r的計算公式為：

(3)

3 結果與分析

3.1 NDVI隨時間變化特征

2000—2020年,內蒙古大青山國家級自然保護區生長季NDVI均值整體上表現出波動上升趨勢,增長速率為0.058/10a,生長季 NDVI均值最大值出現在 2019年、2020 年,為 0.63,最小值出現在2001年,為0.39。2004年、2008年、2013年和2016分別出現了峰值,與當年降水量數據對比發現,這幾年保護區年均降水量均達到了400mm以上,表明NDVI受降水量影響較明顯。2017年后,NDVI呈相對穩定增長趨勢(圖2)。

圖2 大青山自然保護區NDVI變化趨勢

3.2 NDVI空間分布特征

利用 2000—2020 年21a的年均 NDVI 數據,計算得到大青山自然保護區平均 NDVI 空間分布圖(圖3),可以看出,大青山自然保護區 NDVI 整體較高,NDVI值小于0.30的面積占總面積的0.50%,NDVI值0.40—0.50的面積占比最高,占總面積的23.56%(表1)。

圖3 大青山自然保護區NDVI空間分布

表1 NDVI分布面積比例

空間分布呈現出西部低、東部高的從西到東逐漸增加的特征。從行政區域看,包頭市昆都侖區、青山區、石拐區、固陽縣、土默特右旗北部,呼和浩特市武川縣、土默特左旗中部和南部,烏蘭察布市卓資縣南部等區域NDVI值較低。包頭市土默特右旗東北部,呼和浩特市土默特左旗北部、回民區、新城區,烏蘭察布市卓資縣西部和西北部區域NDVI值較高。

3.2.1各功能區NDVI空間分布

大青山自然保護區劃分為6個核心區、6個緩沖區和1個實驗區(圖1),空間上,各功能區NDVI由西向東逐漸增加。核心區、緩沖區和實驗區平均NDVI值分別為0.59、0.57和0.54,核心區、緩沖區位于大青山中高海拔區域,也是植被長勢較好的區域,NDVI值也相對較高(圖4)。

圖4 大青山自然保護區各功能區平均NDVI

3.2.2不同地形NDVI空間分布

大青山自然保護區生長季NDVI受地形因素影響較明顯(圖5),保護區生長季NDVI均值隨著海拔的升高而明顯增加,這與大青山地形起伏和山地植被垂直帶譜密切相關。大青山處于半干旱典型草原帶,山地植被自下而上由三個垂直帶組成,即山麓典型草原帶、山地森林灌叢帶和山地草甸帶,NDVI值較高的山地灌叢、喬木林和草甸主要分布在1400m以上的相對較高海拔區域,山麓典型草原NDVI值相對較低,分布在相對低海拔區域。

圖5 大青山自然保護區不同海拔、坡度和坡向生長季NDVI均值的變化

隨著坡度的增加,保護區生長季NDVI均值由0.49增加到0.63,喬木林、山地灌叢主要分布在坡度較大的中高海拔區域,NDVI值較高；坡度平緩的南北麓和山脊以地帶性草原植被為主,NDVI值也相對較低。

不同坡向上,生長季NDVI均值大小依次為陰坡>半陰坡>半陽坡>陽坡,陰坡、半陰坡立地條件相對好,土壤水分和養分條件較好,發育有較大面積的白樺(Betulaplatyphylla)、蒙古櫟(Quercusmongolica)、落葉松(Larixgmelinii)等喬木林及虎榛子(Ostryopsisdavidiana)灌叢,植被覆蓋度較高,NDVI值也相對較高。

3.3 NDVI變化趨勢分析

利用線性趨勢分析法對 2000—2020年各年平均 NDVI 值與其相對應的時間序列進行逐像元回歸分析,得到該時間段內植被覆蓋變化趨勢(圖6),結果表明,21年間大青山自然保護區植被覆蓋改善的區域面積占總面積的61.75%,其中,顯著改善、極顯著改善面積占總面積的41.13%、20.62%；基本保持不變的面積占總面積的38.17%；退化區域面積只占0.07%。保護區植被覆蓋整體呈明顯的增加趨勢,增長率為0.058/10a。在空間分布上,植被覆蓋顯著及極顯著增加的區域主要分布在保護區西部和東部,中部以變化不顯著為主。

圖6 2000—2020年大青山自然保護區植被覆蓋度變化趨勢

3.3.1各功能區NDVI變化趨勢

從各功能區NDVI均值變化趨勢看(表2),NDVI值極顯著增加、顯著增加、變化不顯著面積百分比從大到小的順序為實驗區>核心區>緩沖區。實驗區相對于較容易受到干擾或植被覆蓋變化較大的區域,長期以來,因在大青山南北麓低海拔地區放牧等人為干擾頻繁,水土流失和植被退化較嚴重,但隨著保護區的建立和國家林業重點工程項目的實施,保護區的實驗區植被覆蓋度變化明顯高于核心區和緩沖區。核心區植被覆蓋顯著增加,表明建立保護區以來,大青山典型植被類型得到有效保護和恢復。

表2 各功能區NDVI變化趨勢

3.3.2地形因素對NDVI變化的影響

海拔是影響山地植被分布的主要因素之一,大青山自然保護區生長季 NDVI變化面積隨海拔升高表現出單峰曲線的變化趨勢(圖7),NDVI顯著增加和極顯著增加面積比例在海拔1600—1800m處最高(圖8),NDVI變化不顯著面積在海拔1400—1600m處最高。海拔低于1400m時,NDVI變化不顯著面積比例最高。

圖7 2000—2020年不同海拔、坡向、坡度NDVI變化面積占比

圖8 2000—2020年不同海拔、坡向和坡度NDVI變化趨勢

大青山自然保護區各坡向面積比例從大到小順序為陽坡>半陰坡>半陽坡>陰坡。陰坡NDVI變化趨勢不顯著面積比例最大,主要分布在保護區中部偏北區域,由于陰坡植被以喬木林和高覆蓋草甸為主,NDVI變化較小；其他坡向上,NDVI變化趨勢面積比例大小順序為極顯著增加>不顯著>顯著增加,顯著增加和極顯著增加區域主要分布在保護區西部和東部,變化不顯著區域主要分布在保護區中部。

3.4 植被覆蓋穩定性分析

從像元尺度計算分析了大青山自然保護區生長季NDVI的變異系數,模擬其穩定性。2000—2020年,穩定性整體表現為“不同波動變化并存,較低、中等和較高波動變化居多,低波動和高波動變化面積占比相對較小”的特點(表3)。空間格局表現為：低波動變化區域和較低波動變化區域主要分布在中部和東部高海拔區域的保護區核心區和緩沖區(圖9),對應于高NDVI值區域的主要植被類型為白樺林、虎榛子灌叢等天然林和高山草甸,植被覆蓋度較高,生長穩定,NDVI波動較小,其中,較低波動變化區環簇于低波動變化區邊緣并向外延伸分布；中等波動變化面積最大,占總面積的27.26%,分散分布在低波動變化和高波動變化區域之外的整個保護區范圍內；較高波動變化分布區與中等波動變化分布類似,與中等波動變化相間分布；高波動變化區域主要分布在保護區西部、中部和實驗區,對應于NDVI值較低的區域,植被以山麓典型草原和低覆蓋度灌叢,如柄扁桃(Prunuspedunculata)灌叢等為主,處于保護區邊緣區域水土流失、植被退化較嚴重區域,對人為干擾和降水量的影響更敏感,NDVI波動變化較強烈。

表3 大青山自然保護區NDVI變異系數(CV)統計

圖9 大青山自然保護區NDVI變異系數分布

3.5 氣象因子對NDVI的影響

干旱、半干旱地區,降水量是影響植被生長狀況的最主要的氣象因子。從2000—2020年時間序列逐像元的保護區年降水量、年均溫度與NDVI的相關性分析可以看出,年降水量與NDVI的變化主要表現為正相關,正相關面積比例達到99.79%,其中P<0.05水平顯著正相關面積比例占59.25%、P<0.10水平顯著正相關面積比例占19.10%、不顯著正相關面積比例占21.44%；不顯著負相關面積占總面積的0.21%。年均溫度與NDVI相關性以不顯著為主,不顯著面積占總面積的97.91%,溫度對大青山自然保護區植被覆蓋無顯著影響。大青山自然保護區處于半干旱的典型草原地帶,保護區山地植被以天然植被為主,降水量是影響總體植被覆蓋狀況和變化的最主要因素(圖10)。

圖10 大青山自然保護區降水量、溫度空間分布及其與NDVI的相關性檢驗

4 討論與結論

4.1 討論

通過區域與像元尺度的NDVI時空變化分析,發現自2000年內蒙古大青山自然保護區管理局成立以來,大青山植被總體覆蓋狀況良好且呈上升趨勢,植被恢復效果顯著,這與保護區大力實施國家林業重點工程和保護區建設有關。NDVI變化受氣候、地形地貌、人為干擾和歷史等多種因素共同影響,大青山地處半干旱典型草原區,自然生態環境較脆弱,植被覆蓋更容易受到降水量波動和人為干擾的影響,尤其是人為因素突出了變化趨勢的不確定性[32]。本研究的起始時間與保護區建立時間一致,因此,研究時間段內人為影響可視為正向干擾,故主要分析了氣象因子和地形因子對NDVI時空變化的影響,人為因素對大青山自然保護區植被變化時空動態過程中相對貢獻的量化評定還有待于進一步研究。

植被覆蓋狀況與氣候因子有著密不可分的聯系,大青山自然保護區生長季NDVI受降水量空間分布影響明顯,從西到東,隨著降水量的增加,NDVI也逐漸增加。在研究時間段內,NDVI與降水量呈顯著正相關,而與溫度之間不存在顯著相關關系,降水量是大青山植被空間分布和變化的主要影響因素,這與干旱、半干旱區相關研究結果一致[33—35]。而從相關性空間分布看(圖10),保護區范圍內,降水量越高,相關性越不顯著,相反,降水量相對低的區域,NDVI與降水量相關性越顯著。這是由于保護區降水量低的區域,植被以稀疏的灌叢或草本群落為主,植被生長受有效降水量的影響更明顯,而降水量較高的區域,植被以覆蓋度較高的天然喬木林和灌木林為主,如白樺林和虎榛子灌叢等,受降水波動的影響較小所導致。隨著植被覆蓋變化及其對氣候響應的深入研究,不同區域進行研究得到的結論也有所差異[12,29,36]。相對而言,在干旱和半干旱地區的植被覆蓋對降水變化更為敏感,但在濕潤區,降水量高反而影響植物的生長[29],氣溫對植被覆蓋的影響突顯出來,植被覆蓋與氣候因子之間的響應作用是一個復雜的過程,受研究尺度(空間和時間尺度)、地形、區位等多種因素的影響。下一步應綜合分析植被恢復對氣候變化和人類活動的響應,量化其相對貢獻值,揭示多個因素相互作用機制,以此為保護區生態恢復和管理提供理論依據。

4.2 結論

本文選取2000—2020年MODIS NDVI數據及同期的氣象數據,采用趨勢分析、變異系數、相關分析等方法,結合DEM數據,從海拔、坡度和坡向的角度探討了內蒙古大青山國家級自然保護區生長季NDVI的時空變化特征及其與氣候因子的關系,得出以下結論:

(1)2000—2020年,內蒙古大青山國家級自然保護區植被總體覆蓋狀況良好且呈上升趨勢,NDVI增長速率為0.058/10a。

(2)從NDVI空間分布上看,保護區 NDVI空間分布呈現出西部低、東部高、中部偏南和偏北區域較低的分布特征；各功能區NDVI平均值大小順序為核心區>緩沖區>實驗區。

(3)從NDVI變化趨勢分析看,2000—2020 年大青山自然保護區植被覆蓋改善的區域面積占總面積的61.75%,其中,顯著改善、極顯著改善面積分別占總面積的41.13%、20.62%；基本保持不變的面積占總面積的38.17%；退化面積占總面積的0.07%；在空間分布上,植被覆蓋顯著、極顯著增加的區域主要分布在保護區西部和東部,中部以變化不顯著為主。保護區各功能區NDVI均值變化趨勢看,NDVI極顯著增加、顯著增加面積百分比大小順序為實驗區>核心區>緩沖區。

(4)植被覆蓋度穩定性分析看,2000—2020年,穩定性整體表現為“不同波動變化并存,較低、中等和較高波動變化居多,低波動和高波動變化面積占比相對較小”的特點,植被以持續恢復為主。

(5)從保護區生長季NDVI與年降水量和年均氣溫逐像元的相關性分析看,2000—2020年降水量與NDVI的變化主要表現為正相關,正相關面積比例達到99.79%；溫度與NDVI相關性以不顯著為主,不顯著面積占總面積的97.91%,溫度對大青山自然保護區植被覆蓋變化的影響較小。