基于Budyko假設(shè)的黃河源區(qū)徑流演變與量化歸因分析

馬明衛(wèi), 王召航, 王文川, 羅黨

(1.華北水利水電大學(xué) 水資源學(xué)院,河南 鄭州 450046; 2.華北水利水電大學(xué) 水利學(xué)院,河南 鄭州 450046;3.華北水利水電大學(xué) 數(shù)學(xué)與統(tǒng)計學(xué)院,河南 鄭州 450046)

在全球氣候變化和人類活動影響加劇的背景下,黃河源區(qū)徑流量日趨減少,威脅著流域的生態(tài)安全,可能損害區(qū)域經(jīng)濟的高質(zhì)量發(fā)展,這一問題引起了政府決策部門和社會公眾的廣泛關(guān)注[1]。影響徑流形成和變化的因素主要有三大類,即氣候因素(降水、氣溫和蒸散發(fā)等)、流域下墊面條件(地形、土壤和植被狀況等)以及人類活動(淤地壩建設(shè)、過度農(nóng)業(yè)灌溉和直接取用水等)。其中,植被作為陸地生態(tài)系統(tǒng)的重要組成部分,在過去幾十年經(jīng)歷了巨大的變化,北極地區(qū)、歐亞北部地區(qū)、北美和發(fā)展中國家(特別是中國和印度),出現(xiàn)了植被綠化的現(xiàn)象[2]。植被變化不僅影響氣候,而且可以通過改變截留蒸發(fā)、植被蒸騰、降水入滲、地下水補給等水文過程,進而影響徑流的形成與演變[3]。因此,在徑流變化的成因分析中,準確把握植被變化對徑流的影響十分重要。

目前,關(guān)于徑流演變與量化歸因分析的研究大多采用以下兩類方法:第一類是基于過程的水文模型法,雖然水文模型法分析精度高,物理機制強,可以對降水、蒸散發(fā)、徑流等要素進行水文過程的模擬,但模型參數(shù)率定困難且受人的主觀因素影響較大,增加了結(jié)果的不確定性[4];第二類方法是基于Budyko假設(shè)的彈性系數(shù)法,與基于過程的水文模型法相比,結(jié)構(gòu)簡單、使用方便,而且對參數(shù)的選擇和模型構(gòu)建的要求相對較低,它在降低不確定性方面具有一定優(yōu)勢[5]。以往人們通常采用基于Budyko假設(shè)的彈性系數(shù)法來定量區(qū)分氣候變化和人類活動對徑流的影響,其中人類活動對徑流變化的貢獻是通過徑流變化總量與氣候?qū)搅髯兓绊懛至康牟钪祦砉烙嫷?即流域下墊面變化(包括植被狀況變化)也包含在人類活動影響中,這樣就不能精細地定量分析植被變化是如何對徑流產(chǎn)生影響的[6]。在Budyko方程中,下墊面參數(shù)n反映的是流域地形、土壤、土地利用和植被覆蓋等綜合效應(yīng)[7],有學(xué)者通過計算下墊面參數(shù)對徑流變化的貢獻率,來粗略地分析植被變化對徑流演變的影響[8-10]。但是,關(guān)于植被變化對徑流影響的相關(guān)研究還比較少且不夠深入,不能精細地計算出植被變化對徑流變化的貢獻量。

歸一化植被指數(shù)(Normalized Difference Vegetation Index,NDVI)是國際上公認的反映大尺度地表植被覆蓋和生長狀況的有效指標(biāo),通常用來表征植被活動的強弱[11]。本文選取近年來植被恢復(fù)較為明顯的黃河源區(qū),根據(jù)彈性系數(shù)法的定義,在Budyko框架中引入植被指數(shù)NDVI,將影響徑流的因素進一步細化,擴展Budyko框架下的具體應(yīng)用公式,量化分析氣候因素、流域下墊面條件(包括植被狀況)和人類活動對徑流變化的貢獻,以期為黃河源區(qū)水資源保護和利用提供參考。

1 研究區(qū)概況

黃河源區(qū)位于青藏高原東北部,地理位置位于32°10′N～36°05′N和95°50′E～103°30′E范圍內(nèi),如圖1所示。其中唐乃亥水文站以上流域面積約12.2萬km2,占黃河流域總面積的16.2%;源區(qū)年徑流量為227.5億m3(2018年),占黃河流域年徑流量的35%。源區(qū)內(nèi)平均海拔在3 000 m以上,以山地地貌為主,山脈綿延,地勢高聳,地形復(fù)雜。區(qū)域氣候?qū)儆诟咴箨懶詺夂?表現(xiàn)為冷熱兩季交替,干濕兩季分明,年溫差小,日溫差大,日照時間長,輻射強烈,無四季區(qū)分的氣候特征,多年平均氣溫為-0.35 ℃,年均降水量為361.60 mm[12]。

圖1 黃河源區(qū)概況圖

2 數(shù)據(jù)與研究方法

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

黃河源區(qū)徑流資料為唐乃亥水文站1956—2018年的逐日實測流量數(shù)據(jù),氣象數(shù)據(jù)來自中國氣象局國家氣象中心(http://www.data.cma.cn/)。綜合考慮黃河源區(qū)所有氣象站的氣象數(shù)據(jù)序列長度、數(shù)據(jù)完整性和空間位置等特性,最終選用了河南、久治、達日、瑪曲、若爾蓋、瑪多和清水河7個氣象站1956—2018年的降水和氣溫觀測數(shù)據(jù)。

歸一化植被指數(shù)NDVI能夠表征陸面植被覆蓋和植被生長情況,是植被宏觀檢測的重要指標(biāo),故采用NDVI來表征植被變化。本研究使用的GIMMS NDVI數(shù)據(jù)來自美國國家航空航天局(http://ecocast.arc.nasa.gov/data/pub/gimms/),時間范圍為1982—2015年,空間分辨率為8 km,時間分辨率為15 d。根據(jù)研究區(qū)的空間范圍,提取相應(yīng)網(wǎng)格單元逐月NDVI數(shù)據(jù),再對每年12個月的月尺度NDVI數(shù)據(jù)取平均值得到黃河源區(qū)1982—2015年的年均NDVI序列。

潛在蒸散發(fā)數(shù)據(jù)采用Hamon公式計算求得。Hamon方法認為氣溫是驅(qū)動蒸散發(fā)的主要動力,同時也考慮了日照時長與飽和蒸氣壓對蒸散發(fā)的影響,該方法的優(yōu)點在于只需要氣溫數(shù)據(jù)便可以估算潛在蒸散發(fā)量[13]。

2.2 研究方法

2.2.1 趨勢及突變檢驗

針對長時間序列數(shù)據(jù),采用線性趨勢法和Mann-Kendall(M-K)非參數(shù)統(tǒng)計檢驗法綜合分析黃河源區(qū)氣溫、降水、潛在蒸散發(fā)、徑流和NDVI的變化趨勢。利用Mann-Kendall突變檢驗法識別相應(yīng)氣象、水文、植被要素演變中的突變現(xiàn)象,結(jié)合滑動T突變檢驗法[14]來去除M-K突變檢驗法中的干擾點,從而進一步確定合理的突變點。根據(jù)年徑流序列的突變點,劃分基準期和變化期。

2.2.2 植被變化對徑流影響的研究方法

流域水量和能量平衡是Budyko理論的基礎(chǔ),其中水量平衡表達式為:

R=P-Ea-ΔS。

(1)

式中:R為年徑流深,mm;P為年降水量,mm;Ea為年實際蒸散發(fā)量,mm;ΔS為年蓄水變化量,mm。對于閉合流域,長時間尺度的蓄水變化量很小(近似為0),為了減少黃河源區(qū)流域ΔS的影響,取各變量時間尺度為10年的滑動平均值。

Budyko理論認為實際蒸散發(fā)量主要取決于降水量和潛在蒸散發(fā)量之間的平衡,其兩個邊界條件為:在極端干燥的環(huán)境下,所有降水量都將蒸發(fā);在極端濕潤的環(huán)境下,蒸散發(fā)消耗的能量都將轉(zhuǎn)化為潛熱[15]。

YANG H B等[16]在Budyko理論的基礎(chǔ)上,推導(dǎo)出了實際蒸散發(fā)量的經(jīng)驗方程:

(2)

式中:E0為年潛在蒸散發(fā)量,mm;n是綜合反映流域下墊面特征的無量綱參數(shù)(簡稱下墊面參數(shù)),與地形、土壤、植被覆蓋和土地利用等一系列因素有關(guān)。

公式(2)中n反映的是流域下墊面變化情況,所以n值的大小在年或者多年尺度上應(yīng)該是變化的。對于一個特定流域,土壤和地形等狀況較少會發(fā)生變化,故下墊面參數(shù)n主要受流域植被狀況變化的影響,根據(jù)LI D等[17]的研究結(jié)果,下墊面參數(shù)n與植被指數(shù)NDVI有如下線性關(guān)系:

n=a·NDVI+b。

(3)

其中,對于特定流域,a和b是常數(shù),a的大小表征了n對NDVI變化的敏感性程度,b則反映除了NDVI以外其他因素對n的影響。

2.2.3 彈性系數(shù)法與貢獻率計算

根據(jù)彈性系數(shù)法的定義[18]推導(dǎo)出徑流量R對降水量P、潛在蒸散發(fā)量E0和下墊面參數(shù)n的彈性系數(shù)計算公式:

(4)

(5)

(6)

式中:εP、εE0、εn分別為徑流量R對降水量P、對潛在蒸散發(fā)量E0、對流域下墊面參數(shù)n的彈性系數(shù)。

彈性系數(shù)為正值,表示徑流量會隨著該影響因子的增加而增加;為負值,則說明徑流量會隨著該影響因子的增加而減小。彈性系數(shù)的絕對值越大,該影響因子對徑流量的影響程度越大。

根據(jù)全微分的概念,可以用微分方程表示各要素對徑流量R變化的貢獻量,計算公式如下:

(7)

將公式(7)兩邊同時除以R,結(jié)合公式(4)(5)(6)得到每個因子對徑流量R相對變化量的表達式:

(8)

式中:等號右邊3項分別是P、E0、n變化引起徑流量R的相對變化量。

本文根據(jù)彈性系數(shù)的定義,并結(jié)合公式(3)將植被指數(shù)NDVI引入到基于Budyko假設(shè)的彈性系數(shù)法中,以Budyko方程中的下墊面參數(shù)n作為中間變量,來定量分析植被變化對徑流變化的影響程度,其中徑流量R對NDVI的彈性系數(shù)為:

(9)

(10)

各要素對徑流變化的貢獻率ηxi的計算公式如下:

(11)

式中:ΔRxi為要素xi變化引起的徑流量變化;ΔR為所有要素變化引起徑流量變化的總和。

3 結(jié)果分析

3.1 水文氣象要素演變趨勢與徑流突變點檢測

黃河源區(qū)1956—2018年年平均氣溫和年平均潛在蒸散發(fā)量的年際變化趨勢如圖2所示。由圖2可知:氣溫變化傾向率為0.032 ℃/年,潛在蒸散發(fā)量變化傾向率為1.026 mm/年;通過氣溫和潛在蒸散發(fā)量的線性擬合結(jié)果可以看出,黃河源區(qū)1956—2018年氣溫和潛在蒸散發(fā)量上升趨勢明顯。M-K趨勢檢驗結(jié)果見表1,該表顯示,黃河源區(qū)1956—2018年氣溫和潛在蒸散發(fā)量呈顯著上升趨勢(顯著水平p>0.05),與線性趨勢法檢驗得出的結(jié)果相一致。

表1 黃河源區(qū)各要素M-K趨勢檢驗結(jié)果

黃河源區(qū)1956—2018年降水量和徑流深的年際變化趨勢如圖3所示。黃河源區(qū)在過去60年間降水相對豐沛,多年平均降水量達到361.65 mm,且有逐年上升的趨勢,但是黃河源區(qū)1956—2018年的徑流量卻并沒有隨著降水量的增加而增加,其變化傾向率為-0.189 mm/年。這說明降水量的變化不是導(dǎo)致黃河源區(qū)徑流量減少的因素,潛在蒸散發(fā)、其他環(huán)境因素和人類活動可能是造成徑流量減少的原因。

圖3 黃河源區(qū)1956—2018年降水量和徑流深變化趨勢

為了進一步研究黃河源區(qū)徑流的變化特性,采用M-K和滑動T兩種突變檢驗方法綜合識別徑流變化的突變點。M-K突變檢驗結(jié)果如圖4所示,UF和UB兩條曲線相交于1991年和1993年兩點,均在0.05顯著水平范圍內(nèi)。結(jié)合滑動T突變檢驗結(jié)果(圖5,子序列長度大于10年),最終確定1991年為黃河源區(qū)徑流變化的突變點。從圖5中可以看出,徑流量在1991年前后發(fā)生突變后開始持續(xù)下降。為了進一步厘清其中的機理,下文將研究時段劃分為1956—1991年和1992—2018年兩個時段進行詳細討論。

圖4 黃河源區(qū)1956—2018年徑流M-K突變檢驗結(jié)果

圖5 黃河源區(qū)1956—2018年徑流滑動T突變檢驗結(jié)果

3.2 下墊面參數(shù)n與植被變化的關(guān)系

為了定量分析黃河源區(qū)植被變化對下墊面參數(shù)n的影響,首先對黃河源區(qū)植被變化趨勢進行分析。圖6給出了黃河源區(qū)1982—2015年年均NDVI變化情況。由圖6可知:年均NDVI呈波動增加趨勢;在研究時段內(nèi),黃河源區(qū)年均NDVI的變化范圍是0.337 2～0.383 5,最小值和最大值分別出現(xiàn)在1983年和2010年。總體上,黃河源區(qū)1982—2015年年均NDVI線性變化斜率為0.000 4/年(p>0.05),植被恢復(fù)較為明顯,黃河源區(qū)整體上呈變綠趨勢。

圖6 1982—2015年黃河源區(qū)年均NDVI年際變化

通過公式(1)和公式(2)聯(lián)立,并由已知的徑流量R、降水量P和潛在蒸散發(fā)量E0,即可推求出黃河源區(qū)在不同年份的下墊面參數(shù)n。然后根據(jù)公式(3),利用植被指數(shù)NDVI和下墊面參數(shù)n擬合出方程中的系數(shù)a和b,結(jié)果如圖7所示。經(jīng)擬合,得到a=25.292,b=-6.989 4,即下墊面參數(shù)n與植被指數(shù)NDVI的關(guān)系為n=25.292NDVI-6.989 4(p>0.05)。p=0.017 4>0.05,達到95%的置信度,判定系數(shù)R2=0.836 6,說明模型有效且擬合的效果較好。由參數(shù)n與NDVI的關(guān)系式可以看出,植被變化是流域下墊面變化的主要控制因子。根據(jù)系數(shù)a(a=25.292)可知,NDVI每增加1%,下墊面參數(shù)n就增加25%以上,這意味著黃河源區(qū)流域產(chǎn)流性能對植被變化有著極強的敏感性。

圖7 黃河源區(qū)年均NDVI與下墊面參數(shù)n的關(guān)系

3.3 徑流演變量化歸因

根據(jù)上文研究結(jié)果,將研究時段劃分為基準期(1956—1991年)和變化期(1992—2018年)。表2詳細列出了黃河源區(qū)基準期和變化期的多年平均徑流深R、降水量P、潛在蒸散發(fā)量E0和下墊面參數(shù)n的具體數(shù)值以及相應(yīng)的變化量。其中,變化期的多年平均徑流深和降水量較基準期的減少,多年平均潛在蒸散發(fā)量和下墊面參數(shù)較基準期的增加。

表2 徑流突變點前后各要素的變化量

為了進一步對徑流變化進行歸因分析,根據(jù)公式(4)(5)(6)和公式(9)計算徑流量對降水量、潛在蒸散發(fā)量、下墊面參數(shù)和植被指數(shù)的彈性系數(shù),結(jié)果見表3。由表3可以看出,徑流量變化與降水量變化呈正相關(guān),與潛在蒸散發(fā)量、下墊面參數(shù)和植被指數(shù)變化呈負相關(guān)。具體表現(xiàn)為,當(dāng)降水量、潛在蒸散發(fā)量、下墊面參數(shù)和植被指數(shù)NDVI均增加10%時,在變化期,徑流量將分別增加15.9%、-6.2%、-13.0%和-9.8%。

表3 徑流突變點前后各要素的彈性系數(shù)

根據(jù)公式(11)求得氣候因素、下墊面條件和人類活動對徑流變化的貢獻率見表4。結(jié)果表明,降水量對徑流量的減少為負貢獻,潛在蒸散發(fā)量、植被指數(shù)和人類活動對徑流量的減少為正貢獻。其中,1992—2018年降水量對徑流量減少的貢獻率為-35.45%,潛在蒸散發(fā)量對徑流量減少的貢獻率為9.24%,下墊面條件對徑流量減少的貢獻率為55.61%(植被指數(shù)NDVI占25.30%),人類活動對徑流量減少的貢獻率為70.60%?？傮w上看,人類活動是黃河源區(qū)徑流量減少的首要因素。氣候因素對徑流量減少的貢獻率變化不大,對徑流量減少為負貢獻,抑制了徑流量減少的趨勢。流域下墊面變化是造成徑流量減少的另一個重要因素。由表4可知,植被恢復(fù)導(dǎo)致徑流量的減少抵消了大部分降水量增加引起的徑流量增加,并且由于潛在蒸散發(fā)量的增加,導(dǎo)致過去幾十年徑流量的總體小幅下降。這一結(jié)果強調(diào)了植被變化對控制黃河源區(qū)長期徑流演變的重要性,需要引起一定的重視。

表4 各影響因子對徑流量變化的貢獻率 %

4 討論

近幾十年來,黃河源區(qū)徑流量總體上呈下降趨勢[19]。張成鳳等[20]采用基于Budyko假設(shè)的彈性系數(shù)法定量識別黃河源區(qū)徑流變化的原因,其研究結(jié)果表明,氣候因素使徑流量先減少后增加,其他環(huán)境驅(qū)動因素是黃河源區(qū)徑流量減少的主要因素,但其研究只能將冰川凍土和人類活動歸類為其他環(huán)境驅(qū)動因素,不能更精細地量化分離各要素是如何對徑流產(chǎn)生影響的。商瀅等[21]認為,降水不是引起黃河源區(qū)徑流變化的原因,降水對徑流量減少的貢獻率為6.1%,人類活動和其他氣候因子對黃河源區(qū)徑流產(chǎn)生了巨大影響。葉培龍等[22]的研究顯示,從年際波動來看,降水是黃河源區(qū)年徑流量變化的主要影響因子,而生態(tài)植被、凍土退化、水儲量變化和人類活動對徑流量的影響也不容忽視。由此可見,以往的研究在徑流演變驅(qū)動機制方面存在一定的分歧。

自1998年我國實施退耕還林還草工程以來,黃河源區(qū)整體上呈變綠趨勢,目前還較少有學(xué)者研究植被變化是如何對徑流演變產(chǎn)生影響的。本文根據(jù)彈性系數(shù)法的定義,在Budyko框架中引入植被指數(shù)NDVI,將影響徑流的因素進一步細化。其中,氣候因素、下墊面條件和人類活動是影響徑流變化的3個主要因素,將氣候因素劃分為降水和潛在蒸散發(fā)兩個因素,將下墊面條件劃分為植被指數(shù)NDVI和其他會對下墊面變化造成影響的因素。通過量化分離造成徑流變化的影響因素,發(fā)現(xiàn)植被指數(shù)NDVI對黃河源區(qū)徑流的演變存在不可忽視的作用。具體表現(xiàn)為,植被恢復(fù)導(dǎo)致徑流的減少量可以抵消大部分降水增加引起徑流的增加量。

需要強調(diào)的是,本研究仍存在一些不足之處。在本研究中,流域下墊面條件對黃河源區(qū)徑流減少的貢獻率為55.61%,其中植被指數(shù)NDVI占25.30%,其他占30.31%。這里的“其他”可能包括流域內(nèi)的土壤地形變化、湖泊水庫變化和土地利用變化等。一般情況下對于特定的流域,其土壤、地形和湖泊水庫較難發(fā)生改變,但近幾十年來,黃河源區(qū)的土地利用卻發(fā)生了較為明顯的變化(圖8)。

圖8 1980年、2015年黃河源區(qū)土地利用情況

從圖8中可以看出:黃河源區(qū)主要土地利用類型是耕地、林地、草地、水域、建設(shè)用地和未利用土地,其中又以草地、林地和未利用土地為主;2015年相較于1980年,耕地、林地和建設(shè)用地面積分別增加了46、461、23 km2,相對增幅分別為10.27%、5.59%、62.16%;草地和水域面積面積分別減少了669、162 km2,相對減幅分別為0.72%、5.85%。本研究中的植被指數(shù)NDVI對徑流變化的影響可以反映出耕地、林地增加和草地減少對徑流變化的影響,卻還不能得出建設(shè)用地增加和水域減少是如何對黃河源區(qū)徑流變化產(chǎn)生影響的。在未來的研究中,應(yīng)進一步細化“其他”因素對徑流變化產(chǎn)生的影響,使研究結(jié)果更加準確。

5 結(jié)論

本文利用黃河源區(qū)水文氣象觀測數(shù)據(jù)和衛(wèi)星遙感NDVI數(shù)據(jù),采用基于Budyko假設(shè)的彈性系數(shù)法,定量分析了黃河源區(qū)氣候變化、流域下墊面條件變化和人類活動對徑流演變的相對貢獻,主要結(jié)論如下:

1)1956—2018年黃河源區(qū)降水量顯著增加(p>0.05),但徑流量呈現(xiàn)不顯著下降趨勢,下降速率為-0.189 mm/年,徑流量在1991年前后發(fā)生突變。

2)在多年時間尺度上優(yōu)化了Budyko假設(shè)下的具體應(yīng)用公式,擬合得到黃河源區(qū)下墊面參數(shù)n與植被指數(shù)NDVI的函數(shù)關(guān)系為n=25.292NDVI-6.989 4(p>0.05)。根據(jù)系數(shù)a為25.292可知NDVI每增加1%,下墊面參數(shù)n則增加25%以上,這意味著黃河源區(qū)產(chǎn)流性能對植被變化有著極強的敏感性。

3)黃河源區(qū)的徑流變化主要是由人類活動引起的,在徑流量的總變化中占70.60%;氣候變化對徑流量減少的貢獻率為-26.21%,抑制了徑流量減少的趨勢;流域下墊面條件變化是造成徑流量減少的另一個重要因素,植被恢復(fù)可能會使黃河源區(qū)產(chǎn)流量和徑流量進一步減少。