王衛(wèi)光，陸文君，邢萬秋，李進興，李長妮

(1. 河海大學(xué)水文水資源與水利工程科學(xué)國家重點實驗室，江蘇 南京 210098；2. 河海大學(xué)水文水資源學(xué)院，江蘇 南京 210098；3. 河海大學(xué)地球科學(xué)與工程學(xué)院，江蘇 南京 211100)

氣候變化和人類活動對水文循環(huán)產(chǎn)生的影響一直是水文學(xué)家關(guān)注的焦點[1-2]。蒸散發(fā)作為陸面水循環(huán)中最重要的水文過程之一，迄今仍然是研究中的薄弱點[3]。關(guān)于蒸散發(fā)研究的基礎(chǔ)理論較多，其中Budyko框架理論由于原理簡單、物理機制明確，在基于氣候-土壤-植被系統(tǒng)的氣象學(xué)、水文學(xué)、生態(tài)水文學(xué)等研究中得到了廣泛的應(yīng)用[4]。自1948年以來，眾多關(guān)于Budyko框架的經(jīng)驗公式被提出，其中含有參數(shù)的經(jīng)驗公式在區(qū)域上對 Budyko 曲線的擬合精度較高，運用較廣。1979年，崔啟武等[5]首次將反映流域下墊面特征的參數(shù)n引入到Budyko方程中。目前在國內(nèi)大多數(shù)的研究中，參數(shù)n都被當(dāng)作恒定值來處理[6-8]。但事實上，恒定不變的參數(shù)n并不能很好地模擬流域水文狀況，因為不僅氣候變化會影響參數(shù)n，人類對水資源的不斷開發(fā)利用、城市經(jīng)濟的發(fā)展等都會改變n值，影響降雨產(chǎn)匯流過程，進而改變水資源的時空分布格局。因此，研究參數(shù)n隨時間的演變規(guī)律及其驅(qū)動機制很有必要。

圖1 研究區(qū)的氣象站點和水文站點分布

黃河流域流經(jīng)不同氣候區(qū)，流域氣候和地貌差異顯著[9]，在氣候變化和人類活動的共同影響下，參數(shù)n變異的情況更復(fù)雜，但關(guān)于此流域參數(shù)n的演變規(guī)律和歸因分析卻未見報道。因此，本文基于花園口控制站以上黃河流域內(nèi)76個氣象站1956—2000年的逐日氣象觀測資料，和花園口水文站1956—2000年的逐月徑流深(圖1)，運用移動窗口法將Budyko方程的輸入變量進行平滑和濾波處理后推算參數(shù)n，并在總結(jié)文獻研究[10-12]的基礎(chǔ)上，篩選出7個代表氣候變化的因子：降水量(P)、溫度(T)、潛在蒸發(fā)量(E0)、逐年最大日降水量(Pmax)、降水季節(jié)性指標(biāo)(SI)、降雨集中度(CI)、氣候干旱指數(shù)(PDSI)和3個反映人類活動的因子：人口狀況(Pop)、生產(chǎn)總值(GGDP)和有效灌溉面積(Airr)，運用多元逐步回歸模型、敏感性系數(shù)法和貢獻評估法，定量評價黃河流域參數(shù)n的變化對氣候變化和人類活動的響應(yīng)。

1 研究方法

1.1 Budyko方程

Budyko方程描述了實際蒸發(fā)依賴于潛在蒸發(fā)量(熱能來源)和可利用水量(水量來源)的程度[13]，利用Budyko 框架可以有效評價氣候、植被和水文循環(huán)之間的相互影響[14]。其中經(jīng)量綱分析和數(shù)學(xué)推導(dǎo)并被廣泛應(yīng)用的傅抱璞公式[15]可較好地表達流域內(nèi)的水熱耦合狀態(tài)。Yang等[16]以傅抱璞公式為基礎(chǔ)，對Choudhury公式變換成另一種表達Budyko曲線的方程為

(1)

式中：E為流域?qū)嶋H蒸發(fā)量，mm。

本文采用式(1)結(jié)合流域水量平衡方程(E=P-R)來推算Budyko方程中的參數(shù)n。

1.2 參數(shù)n各影響因子及計算方法

1.2.1 參數(shù)n各影響因子

a. 潛在蒸發(fā)量(E0)。計算潛在蒸發(fā)的方法主要可分為綜合法、溫度法、輻射法和水面蒸發(fā)量法，其中Penman-Monteith(PM)法由于具有明確的物理機制而成為FAO唯一推薦的方法[17]。因此，本文選取綜合考慮了多種氣象要素并經(jīng)Shuttleworth于1993年修正的Penman-Monteith(PM)公式來計算E0。

(2)

式中：Δ為飽和水汽壓曲線斜率，kPa /℃；γ為干濕常數(shù)，kPa /℃；Rn為凈輻射，MJ/(m2·d)；G為土壤熱通量，MJ/(m2·d)；λ為潛熱，MJ/kg；U2為2 m高處的風(fēng)速，m/s；es為平均飽和水汽壓，kPa；RH為相對濕度，%。

b. 逐年最大日降水量(Pmax)。根據(jù)中國氣象數(shù)據(jù)網(wǎng)上的氣象數(shù)據(jù)統(tǒng)計出流域逐年最大日降水量(Pmax)。

c. 降水季節(jié)性指標(biāo)(SI)。選用Walsh等[18]的公式來計算降水季節(jié)性指標(biāo)：

(3)

式中：SIi是第i年的SI；Ri為第i年的降水量，mm；Xin代表第i年n月份的降水量，mm。

對公式(3)求得的逐年SI取多年平均值即可得到本文所需的降水季節(jié)性指標(biāo)。

d. 降雨集中度(CI)。采用Martin-vide[19]于2004年提出的計算降雨集中度方法，此方法用于研究區(qū)域的降水極值情況更直觀，且簡單易行[20]。根據(jù)降雨集中度的定義，累積降水量百分比(Y)和累積降水天數(shù)百分比(X)呈負的指數(shù)分布，即符合洛倫茲曲線分布：

Y=aXexp(bX)

(4)

式中的a、b是可用最小二乘法率定得到的系數(shù)。則

(5)

e. 氣候干旱指數(shù)(PDSI)。氣候干旱指數(shù)采用帕爾默干旱指數(shù)(Palmer drought severity index，PDSI)代表。PDSI是一個常用的干旱指標(biāo)，不僅考慮當(dāng)時的水分條件，而且考慮前期水分狀況、持續(xù)時間，是個定量描述旱情的較好的指標(biāo)，在水文、 氣象、農(nóng)業(yè)等領(lǐng)域廣泛應(yīng)用。

f. 人口狀況(Pop，萬人)、生產(chǎn)總值(GGDP，億元)和有效灌溉面積(Airr，103hm2)。數(shù)據(jù)整理時，首先根據(jù)統(tǒng)計局官網(wǎng)年報整理得黃河流域流經(jīng)各省的人口狀況Pop、生產(chǎn)總值GGDP和有效灌溉面積Airr。將這些數(shù)據(jù)分別除以各省的行政區(qū)域面積即得各省的Pop密度、GGDP密度和Airr密度，然后乘以黃河流域在各省份分布的面積，最后將所得結(jié)果相加就可以得到黃河流域的人口狀況(Pop)、生產(chǎn)總值(GGDP)和有效灌溉面積(Airr)。

1.2.2 影響關(guān)系的研究方法

a. 流域氣象要素空間插值。黃河流域的降水(P)、平均氣溫(T)、最高氣溫、最低氣溫、相對濕度、平均風(fēng)速和日照時數(shù)采用氣象站點相應(yīng)的氣象數(shù)據(jù)通過插值方法獲得。具體計算方法如下：①將研究區(qū)劃分成10 km×10 km大小的網(wǎng)格；②根據(jù)網(wǎng)格內(nèi)氣象站點的數(shù)據(jù)通過反距離權(quán)重法得到每個網(wǎng)格的各氣象要素值；③對覆蓋整個研究區(qū)的所有網(wǎng)格的各氣象要素值分別求平均，即可得到該研究區(qū)的氣象要素值。

b. 平滑與濾波處理。移動窗口法(moving window)是簡單且有效的平滑隨機變化方法，被廣泛應(yīng)用于處理非平穩(wěn)序列數(shù)據(jù)。在本文中，考慮到Budyko水熱耦合平衡方程適用于多年平均尺度(一般多于10年)，該尺度下的流域蓄水變化幾乎為0，并且受自然氣候變化的影響小，因此將移動窗口定為11年，以1961年代表1956—1966年段，以此類推，2000年代表1995年以后的時段，然后對Budyko水熱耦合平衡方程中的變量及各影響因子進行滑動平均。

c. 多元逐步回歸模型(MSR)。為確定各影響因子對Budyko水熱耦合平衡方程中參數(shù)n變化的影響，選用MSR模型來建立代表氣候變化(P、E0、T、Pmax、SI、CI、PDSI)和人類活動(Pop、GGDP、Airr)的各因子對參數(shù)n的響應(yīng)關(guān)系。

(6)

d. 敏感性系數(shù)法。為識別各因子對Budyko方程參數(shù)n變化的貢獻，選用敏感性系數(shù)法來估算參數(shù)n的變化對各因子的敏感程度，計算公式為

(7)

式中：xi代表第i個自變量；S(xi)為參數(shù)n的變化對變量xi的敏感性系數(shù)；Δxi為第i個自變量的變化值；Δn為參數(shù)n的變化值。

e. 貢獻評估法。數(shù)學(xué)上，對于函數(shù)y=f(x1,x2,…)，t時間序列的因變量y的變化可用微分方程表示：

(8)

令

分別代表y和xi的長期序列的坡度值，式(8)可轉(zhuǎn)換為

(9)

其中

2 結(jié)果與分析

2.1 參數(shù)n的變化趨勢

以P、R和E0的11年滑動平均值作為Budyko水熱耦合平衡方程(式(1))的輸入數(shù)據(jù)推算參數(shù)n。圖2(a)為變化的參數(shù)n值序列與恒定n值(用實測P、R和E0的1956—2000年均值通過Budyko計算得到恒定n值)的比較。從圖2(a)可以看出，從20世紀(jì)70年代初期到80年代中期滑動平均的參數(shù)n值在恒定n值上下波動，而20世紀(jì)70年代初期以前和80年代中期以后滑動平均的參數(shù)n值分別明顯小于和大于恒定n值?；瑒悠骄膮?shù)n值在1961—2000年間呈現(xiàn)顯著的上升趨勢。

整個研究時段中變參數(shù)n的最大值為2.746，最小值為1.783，分別出現(xiàn)在1964年和1999年。將它們分別代入式(1)，計算得蒸發(fā)率(E/P)相對于干旱指數(shù)(E0/P)的分布(圖2(b))。從圖2(b)可以看出，在相同的干旱指數(shù)下，隨著水熱耦合狀態(tài)改變(n的提升)，流域的蒸發(fā)率越大，即在降水量不變的條件下，黃河流域的n越大，蒸發(fā)越大，根據(jù)流域水量平衡方程(E=P-R)可知，徑流越小，意味著整個流域越來越干旱。

(a) 恒定n值與滑動平均n值的比較

(b) 蒸發(fā)率(E/P)相對于干旱指數(shù)(E0/P)的分布

表1 黃河流域上反映氣候變化和人類活動的10個因子的變化趨勢及其顯著性特征

注：“*”代表顯著性通過0.05顯著性檢驗。

此外，將滑動平均的參數(shù)n值與恒定n值分別代入到Budyko方程中計算的E與根據(jù)水量平衡(E=P-R)用實測P和R的滑動平均值計算的E進行比較，并選用納西系數(shù)(NSEC)和平均相對誤差(RE)來評價模型的模擬效果(圖3)。由滑動平均的參數(shù)n計算的E與實測滑動平均值計算的E基本吻合，RE為0，NSEC達到1。而由恒定參數(shù)n計算的滑動平均窗口的E與實測滑動平均值計算的E有較大差距，RE為-0.031，NSEC為0.683。這說明恒定不變的Budyko水熱耦合平衡方程參數(shù)n在模擬流域水文狀況時效果較差，證實了研究變參數(shù)n的必要性。

2.2 參數(shù)n各影響因子的演變趨勢

表1為黃河流域代表氣候變化和人類活動特征的因子的變化趨勢及其顯著性檢驗。由表1可知，黃河流域1956—2000年的R、E0和PDSI呈現(xiàn)顯著下降趨勢，并且MK檢驗和T檢驗的結(jié)果一致。對于T，T檢驗結(jié)果顯示其呈顯著上升趨勢，而MK檢驗結(jié)果并不顯著。在其他代表氣候變化的因子中，除了P呈非顯著下降趨勢外，Pmax、SI和CI均呈現(xiàn)非顯著上升趨勢。表明雖然黃河流域降水量在過去幾十年有微弱的下降，但降水極值(Pmax、CI等)變化卻愈發(fā)劇烈。

圖4展示了各水文氣象變量的時間序列以及11年滑動平均值。從圖4可以看出，花園口水文站R從1985年開始持續(xù)下降，流域E0在1985年之后由下降趨勢轉(zhuǎn)為輕微的上升趨勢，流域T則在1987年之后上升趨勢更為明顯，SI在1990年之后呈現(xiàn)明顯的增加趨勢,PDSI呈現(xiàn)持續(xù)下降的趨勢，其他因子變化較為平緩。

圖5為黃河流域反映人類活動特征的因子的變化趨勢。結(jié)合表1可知流域Pop、GGDP和Airr均呈顯著上升趨勢(均通過MK檢驗和T檢驗)。流域Pop基本呈現(xiàn)線性增加趨勢，GGDP在1985年之后迅速增加，而Airr則從20世紀(jì)80年代初增長變緩。

(a) 徑流

(b) 降水

(c) 潛在蒸發(fā)

(d) 溫度

(e) 逐年最大日降水量

(f) 降水季節(jié)性指標(biāo)

(g) 降雨集中度

(h) 帕爾默干旱指數(shù)

(a) 人口

(b) 國民生產(chǎn)總值

(c) 有效灌溉面積

2.3 參數(shù)n變化對各影響因素的響應(yīng)關(guān)系

為進一步研究各變量對參數(shù)n變化的影響，運用MSR模型將代表氣候變化的因子和人類活動的因子分別作為模型的自變量，參數(shù)n為因變量進行回歸模擬。結(jié)果如下：

n=-1.886+0.00 416P+0.000 248Pop+
0.000 088 6GGDP-0.00 013Airr

(10)

由式(10)可知，與參數(shù)n變化最為相關(guān)的變量分別為P、Pop、GGDP和Airr。用實測P、Pop、GGDP和Airr代入公式(10)模擬的參數(shù)n和滑動參數(shù)n的結(jié)果見圖6。從圖6可以看出MSR模型對n值模擬效果較好，趨勢線y=1.02x-0.06，R2為0.98，同時說明參數(shù)n確實受P、Pop、GGDP和Airr的影響較大。

圖6 MSR模型模擬n值和滑動n值的比較

2.4 影響因子變化對參數(shù)n的影響評估

依據(jù)研究區(qū)MSR模型結(jié)果(式(10))進行敏感性系數(shù)和貢獻率的計算，結(jié)果見表2。由表2可知，Budyko方程參數(shù)n的變化對P最為敏感，敏感性系數(shù)值達到1.169；而對Pop、GGDP和Airr的敏感性系數(shù)值分別為0.795、0.038和-0.152。表2同時列出了這4個因子對參數(shù)n變化的貢獻率，Pop和GGDP對參數(shù)n的增加起到促進作用，并且Pop的貢獻最大，貢獻率為128.7%。P和Airr對參數(shù)n值的增加起到抑制作用，貢獻率分別為-20.2%和-44.3%。總之，氣候變化和人類活動對參數(shù)n值變化的貢獻分別為-20.2%和120.2%，氣候變化中P抑制了Budyko方程參數(shù)n值的增加，而人類活動的加劇促進了參數(shù)n的增加。

表2 P、Pop、GGDP、Airr對參數(shù)n變化的敏感性系數(shù)和貢獻率

3 結(jié) 論

a. 黃河流域的Budyko方程參數(shù)n在1956—2000年間呈現(xiàn)增加的趨勢，相應(yīng)地，蒸散發(fā)也增加，從而出現(xiàn)了黃河流域的徑流在1956—2000年間減小的現(xiàn)象。

b. 黃河流域Budyko方程參數(shù)n的變化與氣候變化和人類活動均相關(guān)。P的變化對參數(shù)n的影響最大，對參數(shù)n的升高起到反作用。人類活動(Pop、GGDP和Airr)對參數(shù)n的變化也較為敏感，人類活動的加劇對參數(shù)n的增加起到促進作用。

c. 在黃河流域，人類活動對Budyko方程中參數(shù)n變化的貢獻率比氣候變化的貢獻率更大，而在氣候因子中，降水和潛在蒸發(fā)對參數(shù)n變化的貢獻率最大。

[1]李洋洋，白潔芳，周維博，等.灞河流域氣候因子對水沙變化的影響[J].水資源保護，2017，33(5):98-105.(LI Yangyang,BAI Jiefang,ZHOU Weibo,et al.Impact of climate factors on variations of water and sediment in Bahe River Basin [J].Water Resources Protection,2017,33(5):98-105.(in Chinese))

[2]郭巧玲，韓振英，丁斌，等.窟野河流域徑流變化及其影響因素研究[J].水資源保護，2017，33(5):75-80.(GUO Qiaoling,HAN Zhenying,DING Bin,et al.Study of runoff variation characteristics and influence factors in Kuye River [J].Water Resources Protection,2017,33(5):75-80.(in Chinese))

[3]ZHANG Yongqiang,LIU Changming,TANG Yanhong,et al.Trends in pan evaporation and reference and actual evapotranspiration across the Tibetan Plateau [J].Journal of Geophysical Research:Atmospheres,2007,112(D12):1103-1118.

[4]ZHANG Lu,POTTER N,HICKEL K,et al.Water balance modeling over variable time scales based on the Budyko framework-Model development and testing [J].Journal of Hydrology,2008,360(1/2/3/4):117-131.

[5]崔啟武，孫延俊.論水熱平衡聯(lián)系方程[J].地理學(xué)報，1979，34(2):9-17.(CUI Qiwu,SUN Yanjun.On the correlative equation of the heat and water balance [J].Acta Geographica Sinica,1979,34(2):9-17.(in Chinese))

[6]栗銘，陳喜.基于Budyko假設(shè)的海河流域蒸散發(fā)量和徑流量估算研究[J].中國農(nóng)村水利水電，2014(6):107-111.(LI Ming,CHEN Xi.Research on the estimation of the evaporation and the runoff in the Haihe River Based on Budyko hypothesis [J].China Rural Water and Hydropower,2014(6):107-111.(in Chinese))

[7]郭生練，郭家力，侯雨坤，等.基于Budyko假設(shè)預(yù)測長江流域未來徑流量變化[J].水科學(xué)進展，2015，26(2):151-160.(GUO Shenglian,GUO Jiali,HOU Yukun,et al.Prediction of future runoff change based on Budyko hypothesis in Yangtze River basin [J].Advances in Water Science,2015,26(2):151-160.(in Chinese))

[8]許珊珊，金亦，劉金濤，等.Budyko理論在水文氣候分區(qū)中的應(yīng)用研究[J].中國農(nóng)村水利水電，2017(3):17-20.(XU Shanshan,JIN Yi,LIU Jintao,et al.The application of Budyko theory to hydro-climatic classification [J].China Rural Water and Hydropower,2017(3):17-20.(in Chinese))

[9]韓艷利，婁廣艷，葛雷，等.黃河流域與水有關(guān)生態(tài)補償框架的探討[J].水資源保護，2016，32(6):142-150.(HAN Yanli,LOU Guangyan,GE Lei,et al.Discussion of water-related ecological compensation framework for Yellow River Basin [J].Water Resources Protection,2016,32(6):142-150.(in Chinese))

[10]YANG Dawen,SUN Fubao,LIU Zhiyu,et al.Analyzing spatial and temporal variability of annual water-energy balance in nonhumid regions of China using the Budyko hypothesis [J].Water Resources Research,2007,43(4):436-451.

[11]SHAO Quanxi,TRAYLEN A,ZHANG Lu.Nonparametric method for estimating the effects of climatic and catchment characteristics on mean annual evapotranspiration [J].Water Resources Research,2012,48(3):949-960.

[12]WILLIAMS C A,REICHSTEIN M,BUCHMANN N,et al.Climate and vegetation controls on the surface water balance:synthesis of evapotranspiration measured across a global network of flux towers [J].Water Resources Research,2012,48(6):327-40.

[13]BUDYKO M I.Climate and life [M].New York：Academic Press,1974.

[14]DONOHUE R J,RODERICK M L,MCVICAR T R.On the importance of including vegetation dynamics in Budyko's hydrological model [J].Hydrology and Earth System Sciences,2006,11(2):983-995.

[15]傅抱璞.論陸面蒸發(fā)的計算[J].大氣科學(xué),1981,5(1):23-31.(FU Baopu.Calculation on land surface evaporation [J].Scientia Atmospherica Sinica,1981,5(1):23-31.(in Chinese))

[16]YANG Hanbo,YANG Dawen,LEI Zhidong,et al.New analytical derivation of the mean annual water-energy balance equation [J].Water Resources Research,2008,44(3):893-897.

[17]ALLEN R G,PEREIRA L S,RAES D,et al.Crop evapotranspiration guidelines for computing crop water requirements [M].Rome:Fao Irrigation and Drainage Paper,1998:56.

[18]WALSH R P D,LAWLER D M.Rainfall seasonality:description,spatial patterns and change through time [J].Weather,1981,36(7):201-208.

[19]MARTIN-VIDE J.Spatial distribution of a daily precipitation concentration index in peninsular Spain [J].International Journal of Climatology,2004,2(8):2709-2729.

[20]邢萬秋，王衛(wèi)光，吳楊青，等.淮河流域降雨集中度的時空演變規(guī)律分析[J].水電能源科學(xué)，2011，29(5):1-5.(XING Wanqiu,WANG Weiguang,WU Yangqing,et al.Change properties of precipitation concentration in Huaihe River Basin [J].Water Resource and Power,2011,29(5):1-5.(in Chinese))