2000-2016年黃河源區植被NDVI變化趨勢及影響因素

劉啟興, 董國濤, 景海濤, 周俊利, 黨素珍, 賈培培

(1.河南理工大學 測繪與國土信息工程學院, 河南 焦作 454000；2.黃河水利委員會 黃河水利科學研究院 水利部 黃土高原水土保持流失過程與控制重點實驗室, 鄭州 450003)

植被作為土地覆被系統中的主要組分，是陸地生態系統存在的基礎條件，也是連接土壤、大氣、水分和人類土地利用的自然“紐帶”[1]。植被是陸地表面能量交換過程、生物地球化學循環過程和水文循環過程中重要的下墊層，在土地利用、覆被變化、全球變化研究中起著“指示器”的作用。歸一化植被指數(NDVI)被認為是反映植被生長狀態及植被覆蓋程度的最佳指示因子，是監測區域或全球植被和生態環境變化最有效指標[2-4]，監測植被動態變化以及分析其與氣候關系的響應情況己經成為全球變化研究中的重要應用方向之一[5-6]。黃河源區位于青藏高原腹地，是黃河流域上游最重作要的產流區和水源涵養區，也是我國生態屏障的重要保障。因此，研究黃河源區植被覆蓋的變化情況，對于整個黃河流域和青藏高原環境生態系統的演變具有重大意義[7-8]。

國內外學者利用NDVI數據來監測植被覆蓋變化以及探尋植被覆蓋與氣象因子之間存在的響應關系。對NDVI影響因素的研究有以下3個思路：將提取的研究區域植被指數面平均值與氣象要素進行簡單相關分析或回歸分析[9]；在像元尺度上進行簡單相關分析[10]；選取研究區氣象站點周圍的NDVI值與對應氣象站點氣象數據進行簡單相關分析和多元回歸分析[11]。研究發現，植被覆蓋與氣象因子之間存在的響應關系受地理因素的影響[12]，不同植被類型NDVI與氣溫、降水的響應程度存在一定差異，但所有植被類型與氣溫、降水都呈現為正相關關系[13-14]，并且這種響應關系存在一定的滯后性[15-16]。劉憲峰等基于MODIS NDVI數據分析了2000—2011年黃河源區植被覆蓋的時空變化特征及驅動因素，結果表明黃河源區植被呈增加趨勢，植被覆蓋的增加主要與氣候的暖濕化及生態保護工程的實施有關[17]。李輝霞等探討了三江源地區不同植被類型對氣候變化的響應，通過分離氣候變化與人類活動對植被NDVI的貢獻定量評估了生態工程對植被變化的影響，認為氣候要素對該區NDVI的貢獻大于人類活動[18]。目前尚未有研究對黃河源區NDVI進行像元尺度的多元回歸分析，缺乏定量化區分氣象因素與人類活動對黃河源區內植被指數的影響的手段。開展NDVI像元尺度的分析，可以準確分析研究區域內每個地點的時空變化及氣象因素對小范圍的影響，空間連續性和異質性能夠更好地體現和分析；利用多元回歸分析，同時考慮多個氣象因素對植被指數的影響，能準確分析氣象因素對植被指數的影響及貢獻，實現氣象因素與人類活動對植被指數影響的定量區分，比較準確地基于像元尺度預測未來植被覆蓋狀況[19-21]。

本文利用M-K檢驗及趨勢分析在像元尺度上對黃河源區2000—2016年植被NDVI進行時空變化分析，并利用相關分析法及復直線回歸分析方法，定量分析氣象因素和人為因素對植被動態變化的影響。

1 研究區概況

黃河源區(31.5°—36.5°N，95.5°—103.5°E)位于青藏高原東北部，面積約12.2萬km2,橫跨青海、四川、甘肅三省(圖1)，絕大部分在海拔3 000～5 000 m以上，區域內分布有高山、盆地、峽谷、草甸、湖泊和冰川、凍土等地貌。由于地處高海拔寒冷地區，黃河源是對氣候變化響應最敏感和生態環境最脆弱的地區。黃河源區自然環境類型多樣，高寒植被分布廣泛，有高寒草甸、高寒草原、高寒沼澤、高山稀疏植被、高寒灌叢、常綠針闊葉林等植被類型，其中主要以高寒草甸和高寒草原為主，其約占黃河源區總面積的70%以上。

圖1 研究區位置及地形

2 數據源及分析方法

2.1 數據源

2.1.1 NDVI數據來源及處理 本文采用遙感數據MOD13Q1為NASA網站(https:∥ladsweb.modaps.eosdis.nasa.gov/search/)免費下載，時間范圍為2000—2016年。MOD13Q1遙感影像數據是由美國對地觀測計劃EOS/Terra衛星攜帶的中分辨率成像光譜儀MODIS(Moderate Resolution Imaging spectro radiometer)獲取的采用正弦曲線投影(SIN，sinusoidal Projection)方式的3級網格數據產品，具有1 km的空間分辨率和16 d的時間分辨率，單景影像覆蓋面積為1 200×1 200 km2，數據采用格式為HDF—EOS。黃河源區地理范圍覆蓋全球正弦曲線投影SIN系統中編號為h25v05和h26v05的兩個影像分帶。

原始影像通過MRT(MODIS Reprojection Tools)工具統一進行投影拼接轉換，轉換為UTM投影，坐標系為WGS84坐標，并利用黃河源區矢量邊界對處理后的數據進行裁剪。為了消除云層、大氣與太陽高度角等干擾影響，將每年NDVI數據采用最大值合成法(MVC，maximum Value Composites)進行合成，獲得從2000—2016年的17幅影像圖。

2.1.2 氣象數據來源及處理 氣溫和降水數據來源于中國氣象科學數據共享服務網(http:∥data.cma.cn/)下載的全國氣象站點數據。利用研究區內站點站號，通過IDL編程提取出黃河源區間內9個氣象站點的月降水量和氣溫數據。由于氣象因子存在空間不均勻性，同時為了逐像元地分析氣象因素和人類活動等不確定因素對植被覆蓋的影響，采用反距離權重插值法將研究區9個氣象站點實測降水、氣溫數據插值為與MODIS NDVI數據集相同的空間分辨率柵格影像，獲得每一個像元的降水和氣溫數據。采用反距離權重插值法將站點實測降水、氣溫數據插值為與MODIS NDVI數據集相同的空間分辨率柵格影像。

2.2 分析方法

2.2.1 Mann-Kendall檢驗法 Mann-Kendall檢驗法被廣泛應用在植被等的趨勢分析中。其優點在于檢測范圍寬、人為影響小、定量化程度高。該方法對于變化要素從一個相對穩定狀態變化到另一個狀態的變化檢驗非常有效[22]。利用Mann-Kendall檢驗法對黃河源區MODIS NDVI隨時間序列變化的特征進行分析，研究結果不僅從統計學意義上揭示植被隨時間的變化趨勢，還揭示植被發生突變的時間段。在Mann-Kendall突變檢驗中，對于具有n個樣本量的時間序列x，構造一個秩序列：

(1)

其中

(2)

可見，秩序列sk是第i時刻數值大于j時刻數值個數的累計數。

在時間序列隨機獨立的假定下，定義統計量：

(3)

式中：UF1=0；E(sk)，var(sk)是累計數sk的均值和方差，在x1,x2,…,xn相互獨立，且有相同連續分布時，它們可由下式算出：

(4)

(5)

給定顯著性水平α,將UFk和UBk兩個統計量曲線和顯著性水平線繪在同一個圖上，若UFk和UBk的值大于0，則表明序列呈上升趨勢，小于0則呈下降趨勢。當超過臨界直線時，表明上升或下降趨勢顯著，超過臨界線的范圍確定為突變的時間區域。如果UFk和UBk兩條曲線出現交點，那么交點對應的時刻便是突變開始的時間。

2.2.2 趨勢分析 趨勢分析是通過對一組隨時間變化的變量進行線性回歸分析，從而預測其變化趨勢的方法。該方法可計算研究區域監測時段內每個柵格點NDVI的變化趨勢，反映該區域植被NDVI變化的方向和速率。線性回歸方法是研究植被變化趨勢的常用方法，其優點在于通過利用各個年份的數據值進行擬合，可以消除研究時段內偶發性異常因素對植被長勢的影響，更真實地反映長時間序列植被變化趨勢。

(6)

(7)

(8)

式中：xi表示研究時段內第i個年份；b為變化趨勢線的斜率；a為截距；NDVIg為各個像元的NDVI；R為單個像元NDVIg變化率；x1和xn分別為研究時段的起始和終止年份。上述公式被廣泛應用在NDVIg的時間序列分析，具有較好的穩定性[20-22]。

2.2.3 相關性分析 采用person相關性分析法與t檢驗對NDVI與同期氣象數據進行相關分析和顯著性檢驗，計算公式如下：

(9)

(10)

2.2.4 復直線回歸分析 本文對NDVI與溫度、降水之間進行復直線回歸分析，回歸方程為：

z=a+bx+cy

(11)

參數計算公式為：

(12)

(13)

(14)

NDVI模擬殘差計算公式：

residual=NDVI-NDVIp

(15)

式中：residual為NDVI復直線回歸方程的殘差；NDVI為NDVI時間數據集；NDVIp是根據復直線回歸模型預測的NDVI值。

3 結果與分析

3.1 黃河源區植被時空變化特征

3.1.1 黃河源區植被NDVI在時間上的動態變化分析 圖2為利用matlab對黃河源區2000—2016年的MODIS NDVI進行MK趨勢檢驗圖。

圖2 2000-2016年黃河源區NDVI非參數檢驗

黃河源區2000—2016年的植被NDVI趨勢分析表明：UF統計量在2004—2016年大于零，說明黃河源區NDVI整體呈增加趨勢；雖然在2002—2004年UF統計量小于零，并且UF線與UB線在2004年出現交點，說明在這期間NDVI的變化趨勢發生較大變化，但交點統計量并未超過顯著性α<0.05的臨界值1.96，大體可以判斷出黃河源區NDVI在2000—2016年期間整體比較平穩，沒有特別明顯的突變年份；從2010年起UF統計量超過了顯著性α<0.05的臨界值1.96，說明黃河源區NDVI在2010年后增長趨勢顯著。

3.1.2 黃河源區植被NDVI在空間上的動態變化分析 圖3為源區植被NDVI年際變化趨勢空間分布圖，由圖可以看到大部分區域呈淺色，說明2000—2016年黃河源區大部分區域NDVI是增加的，大部分區域的變化趨勢為0～0.004/a，變化趨勢最大達到0.213/a；NDVI在2000—2016年在減少的區域主要分布在扎陵湖和鄂陵湖西部與南部、黃河源區的中南部以及若爾蓋濕地的北部。變化趨勢最小值為-0.031/a。源區整體植被NDVI變化趨勢均值為0.001 2/a。

圖3 2000-2016年黃河源區NDVI年際變化趨勢空間

表1 黃河源區NDVI一元線性回歸分析變化

表1為黃河源區NDVI不同級別變化趨勢面積統計表。植被NDVI年變化趨勢<0的面積為3.54萬km2，占黃河源區總面積的29.54%，其中變化趨勢在-0.001～0之間的面積為1.68萬km2，占黃河源區總面積的14.02%；植被NDVI年變化趨勢>0的面積為10.31萬km2，占黃河源區總面積的70.46%，其中變化趨勢在0～0.004/a之間的面積為70 189.20 km2，占黃河源區總面積的58.55%。結果表明黃河源區在2000—2016年期間NDVI整體呈增加趨勢，且變化趨勢都在0～0.004/a，黃河源區部分地區植被覆蓋情況正逐年好轉。

3.1.3 2000年和2016年植被NDVI變化比較 黃河源區2016年植被NDVI年最大值影像減去2000年植被NDVI年最大值影像，得到兩個年份NDVI差值空間分布圖(圖4)，并且做出不同差值范圍面積比例統計表(表2)。2016年植被NDVI同2000年相比，大部分區域植被NDVI差值在-0.2～0，0～0.2這兩個范圍。說明近17 a來，黃河源區植被覆蓋并未出現特別顯著的變化。NDVI差值分布圖中，NDVI減少地區域分布和NDVI年際變化趨勢空間分布中變化趨勢<0的區域基本相同。

2016年NDVI值較2000年減少的區域面積為57 362 km2，占黃河源區總面積的47.84%，其中NDVI差值在-0.2～0的區域面積為57 073 km2，占總面積的47.6%。NDVI差值大于0的區域面積為62 517 km2，占黃河源區總面積的52.15%；其中NDVI差值在0～0.2的區域面積為61 931 km2，占總面積的51.66%。結果表明2016年黃河源區植被覆蓋較2000年基本保持不變有輕微改善。

圖4 2016年與2000年NDVI差值分布

表2 2016年與2000年NDVI差值分布

3.2 黃河源區NDVI動態變化驅動因素定量分析

已有研究表明，植被NDVI對降水、氣溫的響應關系存在一定的滯后性。為了較為準確地定量分析植被NDVI對降水和氣溫的響應關系，本文在獲得每個像元的降水和氣溫數據之后，分別求出黃河源區植被NDVI與歷年降水、氣溫數據的6月、7月、8月、6月和7月均值、7月和8月均值及6—8月3個月均值的相關系數，并進行t分布檢驗。結果表明黃河源區NDVI與7月份的降水總量和6月份的平均氣溫相關系數最大，且均通過了α>0.05的置信度檢驗。所以本文選擇了2000—2016年9個站點每年7月份的降水總量和6月份的平均氣溫作為降水和氣溫的研究序列。

圖5 黃河源區2016年植被NDVI氣候因素貢獻分布

在復直線回歸分析中，利用NDVI的2000—2016年時間序列數據和降水氣溫數據求出a，b，c參數，接著利用確定的NDVI與降水、氣溫之間的回歸方程以及相應的降水和氣溫數據得出NDVI預測值。這個值即為氣候因素對NDVI的貢獻值，用實際NDVI減去NDVI預測值得到的殘差即為人為因素和其他不確定因素對NDVI的影響。

圖5為通過復直線回歸分析得到的2016年黃河源區NDVI預測值即降雨和氣溫對NDVI貢獻值空間分布，黃河源區東部和東南部區域NDVI預測值多在0.5以上，而在北部和西北部區域NDVI預測值較小。這是由于黃河源區屬于大陸性高原氣候，氣溫降水總體上呈現為由東南向西北逐漸遞減的空間分布特征且北部和西北部地區靠近柴達木盆地，土壤水分不足，植被自我更新和自我生長能力較差，植被生長受降水、氣溫等氣候條件限制更大。

表3是氣溫、降水對黃河源區2016年植被NDVI貢獻值。99%以上的區域NDVI氣候因素貢獻值為正值，其中貢獻值在0.5～0.6的面積為17 054 km2，占總面積的14.23%；貢獻值在0.6～0.7的面積為21 651 km2，占總面積的18.06%；貢獻值在0.7～0.8的面積為34 565 km2，占總面積的28.83%，三者總計面積為73 271 km2，占總面積的61.12%。這說明僅在氣候因素的作用下，2016年黃河源區NDVI大部分在0.5～0.8的水平，而2016年NDVI實際觀測值在0.5～0.6的面積為16 733 km2，占黃河源區總面積的13.96%；觀測值在0.6～0.7的面積為21 449，占總面積的17.89%；觀測值在0.7～0.8面積為34 443 km2，占總面積的28.73%，觀測值和實際值的差別不大，可以看出目前氣候因素仍是黃河源區植被覆蓋的主要驅動因素。

為了定量化分析人類活動對黃河源區植被覆蓋的影響，將實際觀測的NDVI減去復直線回歸分析預測的NDVI，即計算出殘差，該殘差即可認為是人類活動對植被NDVI影響的部分(圖6和表4)。人類活動對植被NDVI起削減作用的區域面積為53 265 km2，占黃河源區總面積的44.43%；主要分布在黃河源區西北部的扎陵湖和鄂陵湖周邊以及東南部的若爾蓋盆地南部等地區。人類活動對黃河源區植被NDVI起促進作用的區域面積為66 613 km2，占黃河源區總面積的55.57%，其中貢獻值在0～0.02的區域面積為64 190 km2，占總面積的53.54%。說明這些區域人類活動等因素促進了植被覆蓋度增加。

表3 氣溫降水對黃河源區2016年NDVI貢獻值

表4 人類活動對黃河源區NDVI貢獻值

圖6 人類活動對黃河源區2016年NDVI的貢獻值空間分布

黃河源區地區政府自2003年起開展退牧還草和減畜工程，2015年平均牲畜存欄量減少至1 959.8萬頭，減幅達到8.60%[23-24]。同時2000年國家批準成立三江源國家級自然保護區和一期工程，對黃河源區等地的生態環境產生了積極影響[25]。但仍有接近50%的區域殘差小于零，說明在全球氣候變化和人類活動的影響下，對全球變化高度敏感的高寒生態系統已經開始加速退化，導致高寒草原覆蓋度不斷降低，景觀破碎化加劇，草地水源涵養能力降低，水土流失加劇，這也與眾多研究結果相吻合[26-28]。

4 討論與結論

黃河源區超過70%的區域NDVI是增加趨勢，增加速率在0～0.004/a。2016年同2000年相比，黃河源區較大部分地區的植被NDVI是增加的，但在扎陵湖、鄂陵湖西部和南部、黃河源區的中南部以及若爾蓋濕地的北部等相當大的一部分地區，植被NDVI減少了0～0.2。

氣候因素對植被生長變化起到了主要的促進作用，西部和西北地區貢獻較小，西部和西北部地區由于臨近柴達木盆地，土壤水分不足，植被自我更新和自我生長能力較差。植被覆蓋與氣候因子之間的響應作用是一個復雜的交互系統，本文目前只考慮了氣候因素中的降水和氣溫，忽略了日照、風速、濕度等氣候因子對植被覆蓋的影響，深刻揭示其他各種相關因素對植被生長作用機理的研究還有待進行。

從復直線回歸分析結果可以看出：人類活動對黃河源區整體植被覆蓋產生一定顯著的影響，源區高寒草地的退化趨勢未得到有效遏制。黃河源區草場資源豐富，但可利用耕地面積少，是青藏高原地區重要的畜牧業基地之一。黃河源區在歷史上很長的一段時間內一直處于嚴重超載過牧的狀態。過牧是引起該區草地生態系統退化的主要因素。公路建設期間大量土石方的搬運工作造成的沿線50 m以內植被的破壞，仍然沒有恢復，而且隨著歷次公路的改擴建活動，破壞的面積和范圍在進一步擴大，從公路開始修建到歷次的改擴建活動，由于其間隔時間很短，加之這一地區降雨、溫度較低和植被的生長期短，自我恢復能力差，使得其沿線植被得不到很好的恢復，進而造成了水土流失、泥石流和滑坡等地質災害的頻發。因此要加強黃河源區生態保護力度，國家已經成立三江源自然保護區，并于2005年頒布《青海三江源自然保護區生態保護和建設總體規劃》，采取了牧區產業結構調整、轉變經濟發展方式、牧民進城以及禁牧補貼等一系列有效措施，提高了群眾的環保意識，大幅降低人為因素對長江黃河源區高寒草地破壞，延緩草地的退化。

農作物由于人為管理存在，導致NDVI值較大，存在對黃河源區NDVI分析結果產生影響的可能。通過分析2000年、2005年以及2010年三期黃河源區土地覆被數據發現，黃河源區耕地面積僅占黃河源區總面積的0.4%，因此，農作物對分析黃河源區植被NDVI的影響可以忽略不計。

通過植被NDVI動態分析可以有效地監測黃河源區封禁、退牧還草等政策落實和實施效果；且在現階段人類活動變化幅度不大的情況下，可以根據各種氣候模式預測的未來氣溫、降水數據預測黃河源區未來植被覆蓋狀況的理論值，以此為黃河源區生態保護和生態修復提供參考。