黃河流域植被指數(shù)對(duì)氣候變化的響應(yīng)及其與水沙變化的關(guān)系

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

(1.西北大學(xué) 城市與環(huán)境學(xué)院， 陜西 西安 710127； 2.西北大學(xué) 地表系統(tǒng)與災(zāi)害研究院，陜西 西安 710127； 3.陜西省地表系統(tǒng)與環(huán)境承載力重點(diǎn)實(shí)驗(yàn)室， 陜西 西安 710127)

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

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

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

1 研究方法

1.1 研究區(qū)概況

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

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

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

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

1.3 研究方法

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

Y=aX+b

(1)

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

(2)

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

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

(3)

(4)

(5)

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

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

(6)

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

2 結(jié)果與分析

2.1 NDVI時(shí)空變化

2.1.1 NDVI隨時(shí)間變化特征 黃河流域東西長南北窄，自西向東自然地理?xiàng)l件與社會(huì)經(jīng)濟(jì)條件差異極大，也造成流域內(nèi)NDVI具有一定的差異。黃河流域1982—2015年NDVI平均值為0.520 1，上、中、下游流域NDVI平均值分別為0.470 1,0.538 6和0.725 0，即NDVI平均值：下游流域>中游>流域總體>上游流域，表現(xiàn)出與降水自東向西逐漸減少的一致性，其變化趨勢見圖1。黃河流域總體、上游、中游和下游流域在1982—2015年NDVI均表現(xiàn)為增加趨勢，線性遞增速率分別為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)，顯著性檢驗(yàn)結(jié)果表明，僅有下游流域增加趨勢不顯著。

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

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

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

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

2.2 NDVI與氣溫及降水的關(guān)系

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

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

2.3 NDVI變化與年徑流量及年輸沙量的關(guān)系

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

利用黃河流域及上中下游NDVI年值和年徑流量、年輸沙量制作散點(diǎn)圖，并進(jìn)行線性擬合，結(jié)果見圖4。黃河流域及上游、中游和下游流域NDVI與年徑流量之間均正相關(guān)關(guān)系，即隨著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水平顯著性檢驗(yàn)。

黃河流域及中游流域NDVI與年輸沙量之間呈負(fù)相關(guān)關(guān)系，而上游和下游流域呈正相關(guān)關(guān)系。線性擬合表明，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水平顯著性檢驗(yàn)。

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

3 討 論

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

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

4 結(jié) 論

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

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

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