2000-2015年太行山地區(qū)NDVI變化及驅(qū)動因素分析

李新元 王莉 景海濤 范晨雨 余汛

摘要：基于2000-2015年MODIS NDVI數(shù)據(jù)，結(jié)合同期氣象、土地利用等相關(guān)數(shù)據(jù)進(jìn)行研究。根據(jù)相關(guān)計(jì)算獲得NDVI變化數(shù)據(jù)和土地利用類型變化數(shù)據(jù)，對NDVI變化數(shù)據(jù)進(jìn)行反向研究，尋找影響NDVI變化的驅(qū)動因素。結(jié)果表明，2000-2015年研究區(qū)植被覆蓋情況整體得到改善，平均增速為0.0075/年;NDVI增加區(qū)域的驅(qū)動因素是氣象因素對植被生長的促進(jìn)作用;人類活動使耕地的農(nóng)作物長勢更加旺盛，直接導(dǎo)致了耕地范圍內(nèi)的NDVI提高。太行山綠化工程、退耕還林工程、農(nóng)村人口向城市轉(zhuǎn)移等生態(tài)保護(hù)措施對草地、林地、水域的植被生長也產(chǎn)生了積極影響，人類活動對離城市較遠(yuǎn)的自然區(qū)域植被起促進(jìn)作用;NDVI減少區(qū)域的驅(qū)動因素是人類活動使太行山城市周邊地區(qū)的耕地、草地面積減少，城市面積增加，對城市及其周邊地區(qū)的植被生長起減弱作用。

關(guān)鍵詞：NDVI;ArcGIS;植被覆蓋度;氣候因素;人類活動;太行山地區(qū)

中圖分類號：Q948.1 文獻(xiàn)標(biāo)識碼：A

文章編號：0439-8114（2020）09-0070-07

DOI：10.14088/j.cnki.issn0439-8114.2020.09.014

植被作為陸地生態(tài)系統(tǒng)的主體，它與地質(zhì)、地貌、氣候、水文、土壤、動物和微生物共同構(gòu)成了自然地理環(huán)境，是最能反映其他要素性質(zhì)的指示者[1]。同時，植被還是連結(jié)土壤、水分與環(huán)境的自然紐帶，在陸地碳平衡和氣候系統(tǒng)調(diào)節(jié)中發(fā)揮了重要作用。變化環(huán)境下的植被生態(tài)系統(tǒng)演變規(guī)律、驅(qū)動機(jī)制和預(yù)測調(diào)控研究是全球變化研究中的核心內(nèi)容之一[2-4]。監(jiān)測植被動態(tài)變化，評估人類活動和氣候變化對植被變化的影響程度，對確定合理的生態(tài)工程布局和適應(yīng)性管理對策具有重要的實(shí)用價(jià)值，已成為變化環(huán)境下植被生態(tài)系統(tǒng)可持續(xù)管理研究領(lǐng)域的前沿和熱點(diǎn)[5-7]。

歸一化植被指數(shù)（NDVI）是反映植被生長狀態(tài)的最佳指示因子，能夠表征植被覆蓋的整體狀況，且與植被覆蓋度之間存在著極為顯著的線性關(guān)系，因此利用NDVI計(jì)算植被覆蓋度的方法被廣泛使用。

太行山脈作為京津和華北大平原的天然屏障，既可抵御西北寒潮的襲擊，又可接納東南暖濕氣流，是海河水系的發(fā)源地及黃河中下游部分支流的上游。太行山區(qū)山高坡陡，巖石裸露，植被稀少，由于森林資源曾遭到嚴(yán)重破壞，水土流失嚴(yán)重，易暴發(fā)山洪。政府相關(guān)部門為改善太行山地區(qū)的植被覆蓋度，實(shí)施了太行山綠化工程等一系列植被恢復(fù)計(jì)劃。目前有關(guān)太行山地區(qū)的植被變化已經(jīng)有一部分研究成果，例如李薇等[8]關(guān)于太行山區(qū)不同坡度NDVI變化趨勢差異的分析研究，李曉榮等[9]對太行山區(qū)植被NPP時空變化特征及其驅(qū)動力的分析，胡實(shí)等[10]對中國典型山地植被垂直地帶性特征及其影響要素的研究，以及其他有關(guān)太行山地區(qū)的零散的、具有小區(qū)域特征的研究。通過研究分析發(fā)現(xiàn)已有的研究結(jié)果沒有深入分析太行山地區(qū)的植被覆蓋變化信息和植被變化的驅(qū)動因素。因此，本研究利用MODIS NDVI數(shù)據(jù)，結(jié)合氣溫、降水、DEM、居民地分布、土地利用類型等數(shù)據(jù)，對太行山地區(qū)2000-2015年16年間植被覆蓋度的時空變化情況進(jìn)行分析，并對影響植被變化的氣候因素和人類活動進(jìn)行相關(guān)分析。

1 研究區(qū)概況、數(shù)據(jù)來源與方法

1.1 研究區(qū)概況

太行山地區(qū)位于北緯34°34′-40°43′、東經(jīng)110°14′-114°33′，是中國東部地區(qū)的重要山脈。太行山脈從東北向西南延伸，北起拒馬河谷，南至?xí)x豫邊境黃河沿岸，西臨山西省汾河，東瀕華北平原，橫跨山西、河南、河北、北京4省市。太行山地理環(huán)境復(fù)雜，地勢西緩東陡，受河流切割，多橫谷，為東西交通要道，氣候?qū)儆谂瘻貛О霛駶櫞箨懶约撅L(fēng)氣候。

本研究中研究區(qū)的劃分與傳統(tǒng)意義上以行政區(qū)邊界的劃分方法不同，研究區(qū)的選擇更多考慮自然環(huán)境因素。研究區(qū)面積為12.7484萬km²，以山區(qū)為主，不包括太行山脈周邊的平原地區(qū)，在行政區(qū)域上無明顯界限（圖1）。

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

采用的NDVI數(shù)據(jù)來源于中國科學(xué)院計(jì)算機(jī)網(wǎng)絡(luò)信息中心國際科學(xué)數(shù)據(jù)鏡像網(wǎng)站（http：//www.gs-cloud.cn）MODNDIM中國500m NDVI月合成產(chǎn)品，數(shù)據(jù)是由MODNDID計(jì)算得到，計(jì)算方法為取月內(nèi)每天最大值。坐標(biāo)系EPSG：4326（WGS84），空間分辨率500m，時間分辨率每月。時間跨度為2000-2015年。

根據(jù)研究區(qū)矢量邊界數(shù)據(jù)對NDVI數(shù)據(jù)進(jìn)行裁剪。本研究所使用的MODIS NDVI數(shù)據(jù)已經(jīng)經(jīng)過水、云、氣溶膠等處理，減弱了部分因素對數(shù)據(jù)的干擾。為進(jìn)一步消除其他異常值對研究的影響，采用最大值合成法（Maximum value composites，MVC），將月最大NDVI合成年最大NDVI，以年最大NDVI來表示該地區(qū)的植被覆蓋情況，在一定程度上提高了數(shù)據(jù)的可靠性[11]。由此獲得研究區(qū)2000-2015年16幅影像圖。

氣象數(shù)據(jù)來源于中國氣象數(shù)據(jù)共享服務(wù)網(wǎng)，包括太行山區(qū)域內(nèi)18個氣象站2000-2015年日降水量數(shù)據(jù)。對氣象數(shù)據(jù)進(jìn)行剔除異常值等預(yù)處理，利用反距離權(quán)重插值法將每年的氣溫、降水?dāng)?shù)據(jù)插值為與MODIS NDVI數(shù)據(jù)集相同空間分辨率的柵格影像。研究區(qū)邊界矢量數(shù)據(jù)來源于中國科學(xué)院山地表生過程與生態(tài)調(diào)控重點(diǎn)實(shí)驗(yàn)室[12]。太行山DEM數(shù)據(jù)下載于地理空間數(shù)據(jù)云平臺。太行山綠化工程相關(guān)數(shù)據(jù)來源于中國林業(yè)網(wǎng)。土地利用類型數(shù)據(jù)來源于ESACCI-LC的2000年和2015年全球土地利用類型覆蓋圖，空間分辨率為300m。根據(jù)研究區(qū)范圍進(jìn)行裁剪后按照聯(lián)合國土地覆蓋分類系統(tǒng)（LCCS）進(jìn)行分類，分為耕地、草地、林地、城市、水域、裸地6類土地利用類型。

1.3 研究方法

1.3.1 像元二分模型 歸一化植被指數(shù)是一種由遙感傳感器所接收的地物光譜信息推算而得的反映地表植被狀況的定量值。根據(jù)像元二分模型，一個像元的ND VI可以表達(dá)為由綠色植被部分所貢獻(xiàn)的信息NDVI_veg與由無植被覆蓋（裸土）部分所貢獻(xiàn)的信息NDVI_soil這兩部分組成，像元植被覆蓋度（NDVI_FVC）表達(dá)式如下：

式中，NDVI_max為該像元上校正過的歸一化植被指數(shù);NDVI_soil為全部是裸土或無植被覆蓋區(qū)域的歸一化植被指數(shù);NDVI_veg為全部被植被所覆蓋的歸一化植被指數(shù)[13]。NDVI_veg和NDVI_soil取給定置信度的置信區(qū)間內(nèi)的最大值與最小值，分別取累積頻率為99.5%與0.5%時對應(yīng)的值作為NDVI_veg和NDVI_soil，小于NDVI_soil的像元賦值為0，大于NDVI_veg的像元賦值為1。

1.3.2 Mann-Kendall檢驗(yàn)Mann-Kendall（簡稱為M-K）非參數(shù)檢測法所用樣本不需要遵從一定的分布規(guī)則，也不受少數(shù)異常值的干擾，因此其被廣泛應(yīng)用于長時間序列的趨勢和突變檢驗(yàn)分析[14，15]。在M-K檢驗(yàn)中，將時間序列植被覆蓋率X作為獨(dú)立分布的樣本數(shù)據(jù)，X_i表示年際植被覆蓋率由區(qū)域內(nèi)所有像元對應(yīng)年份的NDVI求平均值計(jì)算得到，S表示秩序列。M-K檢驗(yàn)計(jì)算公式為：

趨勢變化由OF檢驗(yàn)值判斷：

式中，E（S）和var（S）分別為S的均值和方差。

若OF>0，表明序列呈上升趨勢UF<0，則表明序列呈下降趨勢。將植被覆蓋度的逆序列按上述步驟進(jìn)行計(jì)算，使|UB|=-UF，給定顯著性水平α=0.05，當(dāng)|UB|曲線超出臨界線時變化趨勢顯著，而OF和UB在臨界線之間出現(xiàn)交點(diǎn)，則交點(diǎn)對應(yīng)的時刻或區(qū)域就是突變的開始。

1.3.3 趨勢分析 一元線性回歸分析可以模擬每個柵格在研究時間內(nèi)的變化趨勢，對每個像元進(jìn)行基于空間位置的時間變化特征能夠反映整個研究對象的時空變化特征，綜合反映出區(qū)域時空格局的演變過程[16，17]。Slope為像元NDVI線性回歸方程的斜率（s），表示隨時間變化的植被指數(shù)變化趨勢[18]。其計(jì)算公式如下：

式中，s為像元NDVI回歸方程的斜率;n為年跨度，i為年序號，NDVI_i為第i年的NDVI。s>0表明NDVI的變化趨勢是增加的，反之則是減少的。采用F檢驗(yàn)進(jìn)行顯著性分析，公式為：實(shí)觀測值;y_i為未能擬合部分的y_i;y為覆蓋度的16年平均值，并將結(jié)果劃分為極顯著相關(guān)（P<0.01）、顯著相關(guān)（0.010.05）3個等級。

1.3.4 相關(guān)分析 利用研究區(qū)處理過的年NDVI數(shù)據(jù)分別與插值后的氣象數(shù)據(jù)進(jìn)行相關(guān)分析，計(jì)算相關(guān)系數(shù)，分析植被覆蓋變化與氣象因子變化之間的相關(guān)性[9]。

式中，i為年序號;NDVI_i為第i年的NDVI;NDVI是對應(yīng)像元16年的年平均NDVI;n為統(tǒng)計(jì)年數(shù)（16年）;x_i為氣象要素;x為氣象要素多年平均值;r為Pearson相關(guān)系數(shù)。r>0表示氣象要素對植被覆蓋度的影響呈正相關(guān);r<0表示氣象要素對植被覆蓋度的影響呈負(fù)相關(guān)。采用T檢驗(yàn)進(jìn)行顯著性分析，并將結(jié)果劃分為極顯著相關(guān)（P<0.01）、顯著相關(guān)（0.010.05）3個等級。

1.3.5 殘差分析 利用復(fù)直線相關(guān)回歸模型擬合出受氣候因子影響的植被覆蓋度，進(jìn)而計(jì)算其與真實(shí)植被覆蓋度之間的差值，忽略其他不確定因素，把差值作為人為因素來考慮，該方法即為殘差分析[19]。

對NDVI與溫度、降水之間進(jìn)行復(fù)直線回歸分析，回歸方程為：

Z=a+bx+cy（8）

其中的參數(shù)公式為：

式中，x、y、z分別為氣溫?cái)?shù)據(jù)x、降水?dāng)?shù)據(jù)y和NDVI數(shù)據(jù)z的均值;σ_x、σ_y、σ_z分別為各自數(shù)據(jù)系列的均方差;r_xz、r_xy、r_yz分別為x與z，x與y、y與z的相關(guān)系數(shù)。

NDVI模擬殘差計(jì)算公式：

NDVI_ε=NDVI-NDVI_P

對圖3進(jìn)行統(tǒng)計(jì)分析，2000-2015年太行山地區(qū)Slope值大于0的區(qū)域面積為104728.1km²，占總面積的82.15%，植被變化表現(xiàn)為增加;Slope值小于0的區(qū)域面積為22755.9km²，占總面積的17.85%），植被變化表現(xiàn)為減少。采用F檢驗(yàn)進(jìn)行顯著性分析，結(jié)果為植被變化極顯著（P<0.01）區(qū)域占總面積的18.00%，植被變化顯著（0.010.05）區(qū)域占總面積的65.55%。NDVI呈降低趨勢的區(qū)域呈部分區(qū)域相對集中、整體相對零散的分布情況。NDVI呈增加趨勢的區(qū)域分布較為密集，其中植被NDVI變化趨勢0.001≤s<0.010的面積為54920.11km²，占總面積的43.08%;NDVI變化趨勢0.010≤s<0.030的面積為42809.13km²，占總面積的33.58%。結(jié)果表明，太行山地區(qū)在20000-2015年NDVI整體呈增加趨勢，且變化趨勢大都在0.001～0.030。

2.2 太行山地區(qū)2015年氣象因素對NDVI貢獻(xiàn)值分析

利用2000-2015年NDVI的時間序列數(shù)據(jù)以及氣溫、降水?dāng)?shù)據(jù)求出參數(shù)a、b、c的值，通過確定NDVI與氣溫、降水之間的復(fù)直線相關(guān)回歸方程公式帶入相應(yīng)的氣溫、降水?dāng)?shù)據(jù)得到2015年的NDVI預(yù)測值（圖4）。這個預(yù)測值就是氣象因素對NDVI的貢獻(xiàn)值。

將2015年氣溫和降水對植被NDVI的貢獻(xiàn)值空間分布圖（圖4）結(jié)合太行山地區(qū)的DEM圖像（圖1）進(jìn)行分析發(fā)現(xiàn)，高海拔地區(qū)氣候因素對NDVI的貢獻(xiàn)值比低海拔地區(qū)更大，人類居住地及其周圍的區(qū)域氣候因子對植被NDVI的貢獻(xiàn)值較小。NDVI貢獻(xiàn)值在0.5以上的區(qū)域較多分布在地形較為復(fù)雜的山地，這源于地形復(fù)雜的區(qū)域人類活動較少，有比較完整的生態(tài)系統(tǒng)，所以在氣候條件適宜的情況下植被自然生長更加迅速。北部地區(qū)由于長期的礦物開采導(dǎo)致生態(tài)環(huán)境破壞嚴(yán)重，植被稀少，植被恢復(fù)比較困難，所以該地區(qū)自然狀態(tài)下的植被生長比其他地區(qū)恢復(fù)的更慢。

統(tǒng)計(jì)分析氣候因子對植被NDVI的貢獻(xiàn)值，結(jié)果（表2）顯示，2015年99%以上的區(qū)域氣候因子對NDVI的貢獻(xiàn)為正，氣候因子對NDVI貢獻(xiàn)值在0.5～1.0的區(qū)域占總面積的85.22%;而2015年NDVI實(shí)際值在0.5～1.0的區(qū)域占總面積的86.84%，二者在相同區(qū)間（0.5～1.0）面積相近。這說明氣象因素是太行山地區(qū)植被覆蓋變化的主要驅(qū)動因素之一，表現(xiàn)出對植被生長的積極促進(jìn)作用。

2.3 太行山地區(qū)2015年人類活動對NDVI的影響分析

為定量化分析人類活動對太行山地區(qū)植被覆蓋的影響，將實(shí)際觀測值減去利用復(fù)相關(guān)回歸分析計(jì)算得到的預(yù)測NDVI值，即得到殘差值，認(rèn)為該殘差值就是人類活動對NDVI貢獻(xiàn)值（圖5）。

統(tǒng)計(jì)分析人類活動對NDVI貢獻(xiàn)值，結(jié)果（表3）顯示，對植被NDVI起降低作用的區(qū)域面積為35237.68km²，占太行山地區(qū)總面積的27.64%：對太行山地區(qū)植被NDVI起促進(jìn)作用的區(qū)域面積為92246.32km²，占太行山地區(qū)總面積的72.36%。

根據(jù)從歐空局獲取的預(yù)處理后獲得的2000年和2015年的不同土地利用類型數(shù)據(jù)進(jìn)行疊加分析得到不同土地利用類型之間的流轉(zhuǎn)變化情況，進(jìn)行統(tǒng)計(jì)分析得到2000-2015年土地利用類型變化的轉(zhuǎn)移矩陣（表4）。結(jié)果顯示，在6種土地利用類型中耕地的總面積減少了844.33km²，變化率為-2.31%;草地的總面積減少了1877.47km²，變化率為-3.76%;林地的總面積增加了206.99km²，變化率為0.57%;城市的總面積增加了2437.06km²，變化率為244.14%：水域的總面積增加了77.74km²，變化率為18.86%;裸地的總面積無變化。由表4還可以看出，耕地、草地、林地、水域都有轉(zhuǎn)入和轉(zhuǎn)出變化，但城市無轉(zhuǎn)出變化。耕地、草地面積減少顯著，林地面積略微增加，城市面積增加顯著，水域、裸地變化不明顯。耕地、草地、林地、水域都有部分區(qū)域轉(zhuǎn)入為城市用地，其中草地和耕地轉(zhuǎn)入面積最多。

對2000年和2015年各土地利用類型的矢量區(qū)域進(jìn)行提取，結(jié)合研究區(qū)在研究時段內(nèi)的NDVI和變化趨勢Slope值進(jìn)行統(tǒng)計(jì)分析獲得不同土地類型的平均NDVI變化結(jié)果（表5）和Slope值小于。的區(qū)域在2015年不同土地類型分布圖（圖6）。其中裸地面積太小，平均NDVI值變化不具有代表性，不進(jìn)行統(tǒng)計(jì)。由表5可知，2000-2015年耕地的平均NDVI增加了0.15;草地的平均1VDV1增加了0.10;城市的平均NDVI增加了0.10;水域的平均NDVI增加了0.06。由此可知，16年來為使農(nóng)作物更加高產(chǎn)，人們通過品種改良、科學(xué)施肥、科學(xué)管理等手段使耕地種植的農(nóng)作物長勢更加旺盛，直接導(dǎo)致了耕地范圍內(nèi)的NDVI最大值提高，此結(jié)論與李曉榮等[9]對太行山地區(qū)植被NPP的研究結(jié)果相同。

由圖6可知Slope值小于。的區(qū)域在南部都是圍繞城市相對集中分布，在北部地區(qū)的分布比較零散。將城市的分布信息與各種植被變化信息結(jié)合分析可知，2000-2015年NDVI降低的區(qū)域大部分集中在城市及周圍區(qū)域，除北部自然環(huán)境較惡劣的區(qū)域外，氣象因素對NDVI貢獻(xiàn)值較低的區(qū)域都是城市及其周圍地區(qū)，人類活動的NDVI貢獻(xiàn)值為負(fù)的區(qū)域也大都分布在城市及其周圍地區(qū)。

3 小結(jié)與討論

3.1 小結(jié)

本研究基于2000-2015年MODIS NDVI數(shù)據(jù)，結(jié)合氣溫、降水、DEM、行政區(qū)域分布和土地類型數(shù)據(jù)，對太行山地區(qū)2000-2015年16年間的植被覆蓋度的時空變化特征進(jìn)行研究，對影響其變化的氣候因素和人類活動對植被覆蓋度深層變化情況進(jìn)行相關(guān)分析，再結(jié)合分析結(jié)果對太行山地區(qū)的植被變化情況及其原因做初步總結(jié)，結(jié)果如下。