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

1 內(nèi)蒙古自治區(qū)林業(yè)科學研究院,呼和浩特 010010 2 內(nèi)蒙古大青山森林生態(tài)系統(tǒng)國家定位觀測研究站,呼和浩特 010010 3 內(nèi)蒙古農(nóng)業(yè)大學,呼和浩特 010010 4 內(nèi)蒙古自治區(qū)林業(yè)和草原監(jiān)測規(guī)劃院,呼和浩特 010010

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

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

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

1 研究區(qū)概況

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

2 數(shù)據(jù)來源與研究方法

2.1 數(shù)據(jù)來源與預(yù)處理

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

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

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

2.2 研究方法

2.2.1趨勢分析

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

(1)

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

2.2.2變異系數(shù)

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

(2)

2.2.3相關(guān)分析

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

(3)

3 結(jié)果與分析

3.1 NDVI隨時間變化特征

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

圖2 大青山自然保護區(qū)NDVI變化趨勢

3.2 NDVI空間分布特征

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

圖3 大青山自然保護區(qū)NDVI空間分布

表1 NDVI分布面積比例

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

3.2.1各功能區(qū)NDVI空間分布

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

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

3.2.2不同地形NDVI空間分布

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

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

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

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

3.3 NDVI變化趨勢分析

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

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

3.3.1各功能區(qū)NDVI變化趨勢

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

表2 各功能區(qū)NDVI變化趨勢

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

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

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

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

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

3.4 植被覆蓋穩(wěn)定性分析

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

表3 大青山自然保護區(qū)NDVI變異系數(shù)(CV)統(tǒng)計

圖9 大青山自然保護區(qū)NDVI變異系數(shù)分布

3.5 氣象因子對NDVI的影響

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

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

4 討論與結(jié)論

4.1 討論

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

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

4.2 結(jié)論

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

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

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

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

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

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