2000-2015年太行山地區NDVI變化及驅動因素分析
2020-08-19 12:51:38李新元王莉景海濤范晨雨余汛
湖北農業科學 2020年9期
李新元 王莉 景海濤 范晨雨 余汛
摘要：基于2000-2015年MODIS NDVI數據，結合同期氣象、土地利用等相關數據進行研究。根據相關計算獲得NDVI變化數據和土地利用類型變化數據，對NDVI變化數據進行反向研究，尋找影響NDVI變化的驅動因素。結果表明，2000-2015年研究區植被覆蓋情況整體得到改善，平均增速為0.0075/年;NDVI增加區域的驅動因素是氣象因素對植被生長的促進作用;人類活動使耕地的農作物長勢更加旺盛，直接導致了耕地范圍內的NDVI提高。太行山綠化工程、退耕還林工程、農村人口向城市轉移等生態保護措施對草地、林地、水域的植被生長也產生了積極影響，人類活動對離城市較遠的自然區域植被起促進作用;NDVI減少區域的驅動因素是人類活動使太行山城市周邊地區的耕地、草地面積減少，城市面積增加，對城市及其周邊地區的植被生長起減弱作用。
關鍵詞：NDVI;ArcGIS;植被覆蓋度;氣候因素;人類活動;太行山地區
中圖分類號：Q948.1 文獻標識碼：A
文章編號：0439-8114（2020）09-0070-07
DOI：10.14088/j.cnki.issn0439-8114.2020.09.014
植被作為陸地生態系統的主體，它與地質、地貌、氣候、水文、土壤、動物和微生物共同構成了自然地理環境，是最能反映其他要素性質的指示者[1]。同時，植被還是連結土壤、水分與環境的自然紐帶，在陸地碳平衡和氣候系統調節中發揮了重要作用。變化環境下的植被生態系統演變規律、驅動機制和預測調控研究是全球變化研究中的核心內容之一[2-4]。監測植被動態變化，評估人類活動和氣候變化對植被變化的影響程度，對確定合理的生態工程布局和適應性管理對策具有重要的實用價值，已成為變化環境下植被生態系統可持續管理研究領域的前沿和熱點[5-7]。
歸一化植被指數（NDVI）是反映植被生長狀態的最佳指示因子，能夠表征植被覆蓋的整體狀況，且與植被覆蓋度之間存在著極為顯著的線性關系，因此利用NDVI計算植被覆蓋度的方法被廣泛使用。
太行山脈作為京津和華北大平原的天然屏障，既可抵御西北寒潮的襲擊，又可接納東南暖濕氣流，是海河水系的發源地及黃河中下游部分支流的上游。太行山區山高坡陡，巖石裸露，植被稀少，由于森林資源曾遭到嚴重破壞，水土流失嚴重，易暴發山洪。政府相關部門為改善太行山地區的植被覆蓋度，實施了太行山綠化工程等一系列植被恢復計劃。目前有關太行山地區的植被變化已經有一部分研究成果，例如李薇等[8]關于太行山區不同坡度NDVI變化趨勢差異的分析研究，李曉榮等[9]對太行山區植被NPP時空變化特征及其驅動力的分析，胡實等[10]對中國典型山地植被垂直地帶性特征及其影響要素的研究，以及其他有關太行山地區的零散的、具有小區域特征的研究。通過研究分析發現已有的研究結果沒有深入分析太行山地區的植被覆蓋變化信息和植被變化的驅動因素。因此，本研究利用MODIS NDVI數據，結合氣溫、降水、DEM、居民地分布、土地利用類型等數據，對太行山地區2000-2015年16年間植被覆蓋度的時空變化情況進行分析，并對影響植被變化的氣候因素和人類活動進行相關分析。
1 研究區概況、數據來源與方法
1.1 研究區概況
太行山地區位于北緯34°34′-40°43′、東經110°14′-114°33′，是中國東部地區的重要山脈。太行山脈從東北向西南延伸，北起拒馬河谷，南至晉豫邊境黃河沿岸，西臨山西省汾河，東瀕華北平原，橫跨山西、河南、河北、北京4省市。太行山地理環境復雜，地勢西緩東陡，受河流切割，多橫谷，為東西交通要道，氣候屬于暖溫帶半濕潤大陸性季風氣候。
本研究中研究區的劃分與傳統意義上以行政區邊界的劃分方法不同，研究區的選擇更多考慮自然環境因素。研究區面積為12.7484萬km2，以山區為主，不包括太行山脈周邊的平原地區，在行政區域上無明顯界限（圖1）。
1.2 數據來源與預處理
采用的NDVI數據來源于中國科學院計算機網絡信息中心國際科學數據鏡像網站（http：//www.gs-cloud.cn）MODNDIM中國500m NDVI月合成產品，數據是由MODNDID計算得到，計算方法為取月內每天最大值。坐標系EPSG：4326（WGS84），空間分辨率500m，時間分辨率每月。時間跨度為2000-2015年。
根據研究區矢量邊界數據對NDVI數據進行裁剪。本研究所使用的MODIS NDVI數據已經經過水、云、氣溶膠等處理，減弱了部分因素對數據的干擾。為進一步消除其他異常值對研究的影響，采用最大值合成法（Maximum value composites，MVC），將月最大NDVI合成年最大NDVI，以年最大NDVI來表示該地區的植被覆蓋情況，在一定程度上提高了數據的可靠性[11]。由此獲得研究區2000-2015年16幅影像圖。
氣象數據來源于中國氣象數據共享服務網，包括太行山區域內18個氣象站2000-2015年日降水量數據。對氣象數據進行剔除異常值等預處理，利用反距離權重插值法將每年的氣溫、降水數據插值為與MODIS NDVI數據集相同空間分辨率的柵格影像。研究區邊界矢量數據來源于中國科學院山地表生過程與生態調控重點實驗室[12]。太行山DEM數據下載于地理空間數據云平臺。太行山綠化工程相關數據來源于中國林業網。土地利用類型數據來源于ESACCI-LC的2000年和2015年全球土地利用類型覆蓋圖，空間分辨率為300m。根據研究區范圍進行裁剪后按照聯合國土地覆蓋分類系統（LCCS）進行分類，分為耕地、草地、林地、城市、水域、裸地6類土地利用類型。
1.3 研究方法
1.3.1 像元二分模型 歸一化植被指數是一種由遙感傳感器所接收的地物光譜信息推算而得的反映地表植被狀況的定量值。根據像元二分模型，一個像元的ND VI可以表達為由綠色植被部分所貢獻的信息NDVIveg與由無植被覆蓋（裸土）部分所貢獻的信息NDVIsoil這兩部分組成，像元植被覆蓋度（NDVIFVC）表達式如下：
式中，NDVImax為該像元上校正過的歸一化植被指數;NDVIsoil為全部是裸土或無植被覆蓋區域的歸一化植被指數;NDVIveg為全部被植被所覆蓋的歸一化植被指數[13]。NDVIveg和NDVIsoil取給定置信度的置信區間內的最大值與最小值，分別取累積頻率為99.5%與0.5%時對應的值作為NDVIveg和NDVIsoil，小于NDVIsoil的像元賦值為0，大于NDVIveg的像元賦值為1。
1.3.2 Mann-Kendall檢驗Mann-Kendall（簡稱為M-K）非參數檢測法所用樣本不需要遵從一定的分布規則，也不受少數異常值的干擾，因此其被廣泛應用于長時間序列的趨勢和突變檢驗分析[14，15]。在M-K檢驗中，將時間序列植被覆蓋率X作為獨立分布的樣本數據，Xi表示年際植被覆蓋率由區域內所有像元對應年份的NDVI求平均值計算得到，S表示秩序列。M-K檢驗計算公式為：
趨勢變化由OF檢驗值判斷：
式中，E（S）和var（S）分別為S的均值和方差。
若OF>0，表明序列呈上升趨勢UF<0，則表明序列呈下降趨勢。將植被覆蓋度的逆序列按上述步驟進行計算，使|UB|=-UF，給定顯著性水平α=0.05，當|UB|曲線超出臨界線時變化趨勢顯著，而OF和UB在臨界線之間出現交點，則交點對應的時刻或區域就是突變的開始。
1.3.3 趨勢分析 一元線性回歸分析可以模擬每個柵格在研究時間內的變化趨勢，對每個像元進行基于空間位置的時間變化特征能夠反映整個研究對象的時空變化特征，綜合反映出區域時空格局的演變過程[16，17]。Slope為像元NDVI線性回歸方程的斜率（s），表示隨時間變化的植被指數變化趨勢[18]。其計算公式如下：
式中，s為像元NDVI回歸方程的斜率;n為年跨度，i為年序號，NDVIi為第i年的NDVI。s>0表明NDVI的變化趨勢是增加的，反之則是減少的。采用F檢驗進行顯著性分析，公式為：實觀測值;yi為未能擬合部分的yi;y為覆蓋度的16年平均值，并將結果劃分為極顯著相關（P<0.01）、顯著相關（0.01
1.3.4 相關分析 利用研究區處理過的年NDVI數據分別與插值后的氣象數據進行相關分析，計算相關系數，分析植被覆蓋變化與氣象因子變化之間的相關性[9]。
式中，i為年序號;NDVIi為第i年的NDVI;NDVI是對應像元16年的年平均NDVI;n為統計年數（16年）;xi為氣象要素;x為氣象要素多年平均值;r為Pearson相關系數。r>0表示氣象要素對植被覆蓋度的影響呈正相關;r<0表示氣象要素對植被覆蓋度的影響呈負相關。采用T檢驗進行顯著性分析，并將結果劃分為極顯著相關（P<0.01）、顯著相關（0.01
1.3.5 殘差分析 利用復直線相關回歸模型擬合出受氣候因子影響的植被覆蓋度，進而計算其與真實植被覆蓋度之間的差值，忽略其他不確定因素，把差值作為人為因素來考慮，該方法即為殘差分析[19]。
對NDVI與溫度、降水之間進行復直線回歸分析，回歸方程為：
Z=a+bx+cy（8）
其中的參數公式為：
式中，x、y、z分別為氣溫數據x、降水數據y和NDVI數據z的均值;σx、σy、σz分別為各自數據系列的均方差;rxz、rxy、ryz分別為x與z，x與y、y與z的相關系數。
NDVI模擬殘差計算公式：
NDVIε=NDVI-NDVIP
對圖3進行統計分析，2000-2015年太行山地區Slope值大于0的區域面積為104728.1km2，占總面積的82.15%，植被變化表現為增加;Slope值小于0的區域面積為22755.9km2，占總面積的17.85%），植被變化表現為減少。采用F檢驗進行顯著性分析，結果為植被變化極顯著（P<0.01）區域占總面積的18.00%，植被變化顯著（0.01
2.2 太行山地區2015年氣象因素對NDVI貢獻值分析
利用2000-2015年NDVI的時間序列數據以及氣溫、降水數據求出參數a、b、c的值，通過確定NDVI與氣溫、降水之間的復直線相關回歸方程公式帶入相應的氣溫、降水數據得到2015年的NDVI預測值（圖4）。這個預測值就是氣象因素對NDVI的貢獻值。
將2015年氣溫和降水對植被NDVI的貢獻值空間分布圖（圖4）結合太行山地區的DEM圖像（圖1）進行分析發現，高海拔地區氣候因素對NDVI的貢獻值比低海拔地區更大，人類居住地及其周圍的區域氣候因子對植被NDVI的貢獻值較小。NDVI貢獻值在0.5以上的區域較多分布在地形較為復雜的山地，這源于地形復雜的區域人類活動較少，有比較完整的生態系統，所以在氣候條件適宜的情況下植被自然生長更加迅速。北部地區由于長期的礦物開采導致生態環境破壞嚴重，植被稀少，植被恢復比較困難，所以該地區自然狀態下的植被生長比其他地區恢復的更慢。
統計分析氣候因子對植被NDVI的貢獻值，結果（表2）顯示，2015年99%以上的區域氣候因子對NDVI的貢獻為正，氣候因子對NDVI貢獻值在0.5～1.0的區域占總面積的85.22%;而2015年NDVI實際值在0.5～1.0的區域占總面積的86.84%，二者在相同區間（0.5～1.0）面積相近。這說明氣象因素是太行山地區植被覆蓋變化的主要驅動因素之一，表現出對植被生長的積極促進作用。
2.3 太行山地區2015年人類活動對NDVI的影響分析
為定量化分析人類活動對太行山地區植被覆蓋的影響，將實際觀測值減去利用復相關回歸分析計算得到的預測NDVI值，即得到殘差值，認為該殘差值就是人類活動對NDVI貢獻值（圖5）。
統計分析人類活動對NDVI貢獻值，結果（表3）顯示，對植被NDVI起降低作用的區域面積為35237.68km2，占太行山地區總面積的27.64%：對太行山地區植被NDVI起促進作用的區域面積為92246.32km2，占太行山地區總面積的72.36%。
根據從歐空局獲取的預處理后獲得的2000年和2015年的不同土地利用類型數據進行疊加分析得到不同土地利用類型之間的流轉變化情況，進行統計分析得到2000-2015年土地利用類型變化的轉移矩陣（表4）。結果顯示，在6種土地利用類型中耕地的總面積減少了844.33km2，變化率為-2.31%;草地的總面積減少了1877.47km2，變化率為-3.76%;林地的總面積增加了206.99km2，變化率為0.57%;城市的總面積增加了2437.06km2，變化率為244.14%：水域的總面積增加了77.74km2，變化率為18.86%;裸地的總面積無變化。由表4還可以看出，耕地、草地、林地、水域都有轉入和轉出變化，但城市無轉出變化。耕地、草地面積減少顯著，林地面積略微增加，城市面積增加顯著，水域、裸地變化不明顯。耕地、草地、林地、水域都有部分區域轉入為城市用地，其中草地和耕地轉入面積最多。
對2000年和2015年各土地利用類型的矢量區域進行提取，結合研究區在研究時段內的NDVI和變化趨勢Slope值進行統計分析獲得不同土地類型的平均NDVI變化結果（表5）和Slope值小于。的區域在2015年不同土地類型分布圖（圖6）。其中裸地面積太小，平均NDVI值變化不具有代表性，不進行統計。由表5可知，2000-2015年耕地的平均NDVI增加了0.15;草地的平均1VDV1增加了0.10;城市的平均NDVI增加了0.10;水域的平均NDVI增加了0.06。由此可知，16年來為使農作物更加高產，人們通過品種改良、科學施肥、科學管理等手段使耕地種植的農作物長勢更加旺盛，直接導致了耕地范圍內的NDVI最大值提高，此結論與李曉榮等[9]對太行山地區植被NPP的研究結果相同。
由圖6可知Slope值小于。的區域在南部都是圍繞城市相對集中分布，在北部地區的分布比較零散。將城市的分布信息與各種植被變化信息結合分析可知，2000-2015年NDVI降低的區域大部分集中在城市及周圍區域，除北部自然環境較惡劣的區域外，氣象因素對NDVI貢獻值較低的區域都是城市及其周圍地區，人類活動的NDVI貢獻值為負的區域也大都分布在城市及其周圍地區。
3 小結與討論
3.1 小結
本研究基于2000-2015年MODIS NDVI數據，結合氣溫、降水、DEM、行政區域分布和土地類型數據，對太行山地區2000-2015年16年間的植被覆蓋度的時空變化特征進行研究，對影響其變化的氣候因素和人類活動對植被覆蓋度深層變化情況進行相關分析，再結合分析結果對太行山地區的植被變化情況及其原因做初步總結，結果如下。
