黃河流域植被指數對氣候變化的響應及其與水沙變化的關系

郭 帥， 裴艷茜， 胡 勝,3， 楊冬冬， 邱海軍,3， 曹明明

(1.西北大學 城市與環境學院， 陜西 西安 710127； 2.西北大學 地表系統與災害研究院，陜西 西安 710127； 3.陜西省地表系統與環境承載力重點實驗室， 陜西 西安 710127)

植被是全球環境變化研究的重要環節，是鏈接大氣與土壤的重要物質載體[1-2]。植被分布與變化過程與全球氣候格局關系密切，植被既受全球大氣候背景的影響，也會對區域小氣候產生作用[3-4]。全球氣候變暖大背景下，探討不同區域植被對氣候變化響應過程，對于科學認識陸地植被變化特征，辨識不同區域植被對氣候變化的響應差異，以及指導生態環境政策科學制定具有重要意義。黃河是中華民族的母親河，黃河流域是中華文明的發源地，隨著黃河流域中部地區黃土高原退耕還林還草工程實施，黃河流域植被及河流徑流量、輸沙量發生了顯著變化，對整個黃河流域產生影響[5-8]。開展黃河流域植被對氣候變化響應及其與水沙關系研究，對于科學認識黃河流域生態環境變化過程及政策調控具有重要理論與現實意義。

黃河流域植被時空變化及其水熱響應、黃河水沙關系變化是多學科學者關注的熱點問題[9-10]。2019年9月18日，習近平總書記在“黃河流域生態保護和高質量發展”座談會上強調“共同抓好大保護、協同推進大治理、讓黃河成為造福人民的幸福河”[11]，為黃河流域生態文明建設提出了要求，指明了方向。黃河流域開展的相關研究較多，主要涉及黃河流域總體、黃河源區、黃河中游黃土高原等區域，其中，已有研究認為黃河流域總體植被呈改善趨勢，尤其是退耕還林還草工程實施以來，植被顯著改善[12-13]，但流域內部存在空間差異，如受氣候變化影響，黃河源區氣溫升高，蒸發增強，對植被生長產生不利影響[14-15]，而黃河中游流域黃土高原和下游植被變化更多受到人類活動影響，趨于改善[16-17]。黃河水沙時間序列變化及其與植被變化關系等方面的研究開展較多，其中黃河徑流量、輸沙量銳減尤為受到關注[18-19]，但水沙銳減與流域內植被變化間的物理機制尚未厘清[5,8-9]。

MODIS (moderate resolution imaging spectroradiometer) NDVI(normalized difference vegetation index)具有時間尺度短但空間分辨率高的特點，在三江源地區、黃土高原地區的研究中使用頻率較高[20-22]，而GIMMS (global inventory modelling and mapping studies) NDVI具有空間分辨率低但時間尺度長的特點，是開展大尺度、大范圍研究的首選數據來源，應用也較為廣泛[23-25]。本研究以1982—2015年黃河流域GIMMS NDVI，氣溫，降水和2000—2015年黃河上中下游徑流量、輸沙量數據為基礎，利用線性回歸方法、相關分析方法等，開展黃河流域NDVI時空變化和對氣溫、降水的響應特征，及與徑流量和輸沙量關系研究，以期為該流域生態環境政策調整提供參考。

1 研究方法

1.1 研究區概況

黃河流域位于95°50′29″—119°6′26″E和32°6′53″—41°48′18″N之間，涉及青海、四川、甘肅省、寧夏回族自治區、內蒙古自治區、陜西、山西、河南和山東省，面積總計8.09×105km2。海拔高度0～6 255 m，其中流域西部為青藏高原，中部主要為黃土高原，東部為華北平原，地形自西向東由高原轉為平原。黃河流域氣候類型多樣，黃河源區以高原氣候為主，上游地區主要為溫帶大陸性干旱—半干旱氣候，中下游地區為溫帶大陸性半干旱—半濕潤氣候。流域內部季節特征差異明顯，溫度東西向梯度特征高于南北向；降水集中在夏季，但分布不均勻，年際變化較大；總體上流域濕度小，蒸發大，無霜期短。

利用ArcGIS中水文分析模塊，基于90 m分辨率STRM (shuttle radar topography mission) DEM (digital elevation model)提取黃河流域各小流域范圍和河網，并以內蒙古自治區托克托縣、河南省孟津縣作為黃河流域上中下游的分界點，將各小流域進行歸并，獲得黃河流域上中下游范圍，其中上游范圍面積41.21×104km2，占流域總面積的50.94%；中游范圍面積3.29×105km2，占40.66%；下游范圍面積6.79×104km2，占8.40%。黃河流域上游地區尤其是青藏高原、高山地區，植被變化與氣候因素關系密切；中游黃土高原腹地除受氣候變化影響外，人類活動尤其是退耕還林還草工程的實施對植被恢復與重建效果明顯；下游地區河南、山東省境內植被變化主要受到人類活動的干擾[26-28]。

1.2 數據來源與處理

研究使用數據包括：①90 m分辨率STRM DEM數據，下載自地理空間數據云(http:∥www.gscloud.cn/)，用于黃河流域小流域提取及參與上、中、下游范圍的劃分。②1982—2015年15d 8 km分辨率GIMMS NDVI數據，下載自美國國家航空航天局(http:∥neo.sci.gsfc.nasa.gov/)。對每年24張影像，采用最大值合成(maximum value composition， MVC)方法獲得年值序列數據。③1980—2015年1 km分辨率中國逐年平均氣溫、降水量空間插值數據集，該數據集是利用全國2 400多個氣相站點日監測數據，應用澳大利亞ANUSPLIN插值軟件生成。下載自資源環境數據云平臺(http:∥www.resdc.cn/)。利用黃河流域邊界裁切得到黃河流域年平均氣溫和年降水量空間插值數據，并重采樣為8 km分辨率進行相關分析。④2000—2015年黃河流域《水資源公報》和《泥沙公報》，下載自水利部黃河水利委員會網站(http:∥www.yrcc.gov.cn/)。用于統計黃河及上中下游(頭道拐、三門峽、利津水文站)年徑流量、年輸沙量數據。

1.3 研究方法

1.3.1 線性回歸方法 線性回歸方法用于分析1982—2015年NDVI時間序列和空間序列(像元尺度)變化過程及2000—2015年NDVI與徑流量、輸沙量線性擬合關系，計算公式如下[29]：

Y=aX+b

(1)

在分析NDVI變化趨勢時，式中Y為NDVI，X為年份；在分析NDVI與徑流量和輸沙量關系時，Y為徑流量或輸沙量，X為NDVI。a為線性回歸方程系數，b為常數。a值計算公式如下[29]：

(2)

在分析NDVI變化趨勢時；x為年份(1982—2015年)；y為1982—2015年NDVI影像；i為樣本數(i=1,2,…,34)。其中x1=1982，…，x34=2015；y1為1982年NDVI影像，…，y34為2015年NDVI影像。在分析NDVI與年徑流量和輸沙量關系時，x為年徑流量或年輸沙量(2000—2015年)；y為2000—2015年對應年份NDVI影像平均值；i為樣本數(i=1,2,…,16)。其中x1為2000年NDVI影像平均值，…，x16為2015年NDVI影像平均值；y1為2000年徑流量或輸沙量，…，y16為2015年徑流量或輸沙量。a值同前。

1982—2015年NDVI變化趨勢顯著性采用F檢驗，計算公式如下[29]：

(3)

(4)

(5)

式中：x為年份(1982—2015年);y為1982—2015年NDVI影像;i為樣本數(i=1,2,…,34)，n=34;y′為y擬合值;x′為x的平均值;xi和yi含義同前。

1.3.2 相關系數 用于分析1982—2015年黃河流域NDVI與年平均氣溫、年降水量的相關關系，計算公式如下[29]：

(6)

式中：r為相關系數；x為1982—2015年NDVI影像；y為年平均氣溫或年降水量，i為樣本數(i=1,2,…,34)，其中x1為1982年NDVI影像，…，x34為2015年NDVI影像；y1為1982年平均氣溫或降水量插值數據，…，y34為2015年平均氣溫或降水量插值數據。

2 結果與分析

2.1 NDVI時空變化

2.1.1 NDVI隨時間變化特征 黃河流域東西長南北窄，自西向東自然地理條件與社會經濟條件差異極大，也造成流域內NDVI具有一定的差異。黃河流域1982—2015年NDVI平均值為0.520 1，上、中、下游流域NDVI平均值分別為0.470 1,0.538 6和0.725 0，即NDVI平均值：下游流域>中游>流域總體>上游流域，表現出與降水自東向西逐漸減少的一致性，其變化趨勢見圖1。黃河流域總體、上游、中游和下游流域在1982—2015年NDVI均表現為增加趨勢，線性遞增速率分別為0.001 6/a(p<0.001)，0.001 0/a(p<0.001)，0.002 6/a(p<0.001)和0.000 7/a(p>0.05)，顯著性檢驗結果表明，僅有下游流域增加趨勢不顯著。

圖1 1982-2015年黃河流域及上、中、下游流域歸一化植被指數(NDVI)隨時間變化過程特征

2.1.2 NDVI隨空間變化特征 在像元尺度上，將計算出的1982—2015年黃河流域NDVI變化速率劃分為增加和減少2類，并將計算出的F檢驗結果根據閾值(4.15)劃分為顯著和不顯著2類，最后將變化速率和F檢驗結果進行疊加，則可將NDVI變化趨勢劃分為4類，分別為顯著增加、(不顯著)增加、(不顯著)減少、顯著減少(見圖2)。

圖2 1982-2015年黃河流域NDVI隨空間變化趨勢

由圖2可知，1982—2015年黃河流域NDVI以增加趨勢為主，占流域總面積的91.66%，其中顯著增加面積占75.18%，集中分布在中游及其附近區域；顯著減少面積僅占1.87%，集中分布在中游關中平原、上游青海省境內，下游河南省、山東省黃河河灘、上游內蒙古自治區境內也有分布。上游流域NDVI呈增加趨勢面積總流占域面積的90.05%，其中顯著增加趨勢占70.37%；中游NDVI增加趨勢占95.04%，其中顯著增加趨勢占85.33%；下游NDVI增加趨勢占84.97%，其中顯著增加趨勢占55.09%。

2.2 NDVI與氣溫及降水的關系

氣溫、降水是植被生長發育的必須條件，充足的熱量和水分條件有利于植被的生長。黃河流域自南而北熱量遞減，水分自東向西遞減。黃河源區地處青藏高原區東北部，相對流域其他區域較濕潤，黃河流域中北部和西北部(陜西北部、內蒙古南部、甘肅中北部、寧夏回族自治區等)由于降水相對較少，屬于半干旱和干旱地區。以相關系數0值為界，將1982—2015年黃河流域NDVI與年平均氣溫、年降水量的相關關系劃分為負相關和正相關，再以相關系數臨界值(0.349 4)將負相關和正相關劃分為顯著負相關、不顯著負相關、不顯著正相關、顯著正相關4類。1982—2015年黃河流域NDVI與年平均氣溫相關性分布(見圖3a)表明，流域總體上NDVI與年平均氣溫呈顯著正相關的面積占22.39%，集中分布在陜西省北部、內蒙古自治區南部和寧夏回族自治區東部，其他地區有零散分布；顯著負相關比例較低(占3.41%)。黃河流域上游、中游和下游地區NDVI與年平均氣溫呈顯著正相關的面積分別為占21.23%,25.97%,11.88%。黃河流域NDVI與年降水量相關性分布(見圖3b)表明，總體上NDVI與年降水量呈顯著正相關面積占21.99%，集中分布在流域中北部(陜西省北部、內蒙古自治區境內)；顯著負相關比例不足1%，集中分布在流域的西部(青海省境內)和中南部(關中平原等)。黃河流域上游、中游和下游地區NDVI與年降水量呈顯著正相關面積分別占22.09%,25.72%,3.25%。

圖3 1982-2015年黃河流域NDVI與氣溫及降水相關性的空間分布特征

2.3 NDVI變化與年徑流量及年輸沙量的關系

黃河流域水沙變化是多種因素共同作用的結果。本文僅從統計學角度探討NDVI變化與年徑流量、年輸沙量關系，是對黃河流域植被變化與水沙關系的粗淺認識，無法揭示其物理過程。采用2000—2015年利津水文站數據表征黃河流域總體情況，頭道拐水文站數據表征黃河流域上游情況，三門峽水文站與頭道拐水文站差值表征黃河流域中游情況，利津水文站與三門峽水文站差值表征黃河流域下游情況。

利用黃河流域及上中下游NDVI年值和年徑流量、年輸沙量制作散點圖，并進行線性擬合，結果見圖4。黃河流域及上游、中游和下游流域NDVI與年徑流量之間均正相關關系，即隨著NDVI增加，年徑流量呈增加趨勢。線性擬合表明，黃河流域NDVI增加0.001，則年徑流量增加1.39×108m3(p<0.1)(圖4a)，上游、中游和下游流域年徑流量分別增加1.85×108m3(p<0.001)、2.79×107m3和5.28×107m3(圖4b，4c，4d)，其中僅黃河流域總體和上游流域分別通過0.1和0.001水平顯著性檢驗。

黃河流域及中游流域NDVI與年輸沙量之間呈負相關關系，而上游和下游流域呈正相關關系。線性擬合表明，NDVI每增加0.001，黃河流域及中游流域年輸沙量下降3.80×105t和1.94×106t(p<0.1)(圖4e，4g)，上游和下游流域分別增加4.90×105t(p<0.1)和1.65×106t(圖4f，4h)，其中僅有上游和中游流域通過0.1水平顯著性檢驗。

圖4 2000-2015年黃河流域NDVI與徑流量及輸沙量的關系

3 討 論

植被變化主要受制于氣候變化和人類活動。盡管總體上黃河流域及上、中、下游流域NDVI均呈增加趨勢，但在空間上存在明顯差異。NDVI表現為下降的區域主要集中在黃河源區、流域北部內蒙古境內、關中平原和下游河灘區，而中部大面積區域呈增加過程。黃河源區1982—2015年NDVI下降明顯的區域(圖2)，其年平均氣溫多呈升高趨勢，而年降水量多呈下降趨勢，呈暖干化特征。氣候的暖干化抑制了黃河源區植被的生長發育。近些年黃河源區暖干化趨勢有所下降[15,21]，抑制作用有所減弱并表現出海拔梯度差異[30]。氣候變化對上游流域影響高于中游和下游流域。人類活動在中游和下游流域產生了截然不同的結果，其中黃河流域中部大面積區域的NDVI增加，與人類活動密不可分，退耕還林還草工程實施，促使大面積坡耕地轉為林草地，極大的改善了下墊面狀況[16-17,22]，人類活動對植被變化起到了積極作用，但隨之引起的其他自然要素連鎖反應需要深入開展相關研究，如下墊面改善引起的蒸散發變化、黃河水沙變化等問題；而NDVI下降的區域集中在流域北部內蒙古境內的呼和浩特市附近、關中平原的西安市附近、及下游河灘區的鄭州市和濟南市附近區域，這些區域NDVI的下降與城市的快速擴張具有密切關系[20,31]。

黃河流域水沙變化尤其是水沙銳減是學者們一直關注的問題，但其物理機制尚不清晰[9]。本研究涉及的NDVI與徑流量、輸沙量關系僅從統計學角度開展了分析，結果反映出黃河流域在總體植被改善情況下，黃河流域和上游流域年徑流量呈增加趨勢，而輸沙量在上游流域表現為增加趨勢，在中游流域表現為下降趨勢。黃河流域內部植被變化與徑流量、輸沙量變化關系的空間異質性反映出黃河流域水沙變化驅動機制的復雜性。在植被與水沙變化方面開展的相關研究表明植被改善對產沙影響高于產流[32]，退耕還林還草工程引起的植被覆蓋變化和水土流失治理工程實施，削弱了流域水沙動力關系[33]，黃河中游植被增加引起徑流系數下降[17]、河龍區間輸沙銳減半數由植被恢復所引起[34]，而黃河流域植被覆蓋達到25%時具有最高的水土保持效益[19]，這些研究反映出植被變化與徑流量、輸沙量變化之間具有一定的關系，但僅從統計學角度無法認識這種復雜關系，需要利用水文模型開展黃河流域植被變化與水沙變化關系研究[35-36]。

4 結 論

(1) 1982—2015年黃河流域NDVI平均值表現為：下游>中游>流域總體>上游流域。黃河流域及上、中游流域NDVI均呈顯著線性增加趨勢，其中中游流域增速最快，高于流域總體增速，達到0.002 6/a。黃河流域NDVI顯著增加區域面積占75.18%，集中分布在中游及其附近區域；上游、中游和下游流域NDVI顯著增加區域 面積分別占流域總面積的70.37%，85.33%和55.09%。

(2) 1982—2015年黃河流域NDVI與年平均氣溫、年降水量均以正相關為主，其中顯著正相關分別占22.39%和21.99%，集中分布在中北部區域(陜西省北部、內蒙古自治區南部和寧夏回族自治區東部)。上、中和下游流域中，以中游流域NDVI與年平均氣溫和年降水量的顯著正相關比例最高，分別達到25.97%和25.72%，均高于黃河流域總體。

(3) 隨著黃河流域植被總體改善，2000—2015年黃河流域NDVI與年徑流量、年輸沙量關系表現出明顯區域差異。黃河流域及上游流域年徑流量表現出隨NDVI增加而增加的趨勢，而上游和中游流域年輸沙量隨NDVI增加的表現趨勢不一致，上游流域表現為輸沙量增加趨勢，而中游流域表現為下降趨勢。黃河流域水沙變化是多種因素共同作用的結果，單純進行植被指數與水沙關系的統計學分析，不能揭示其物理過程機制，需要深入開展相關研究工作。