黃河流域Budyko方程參數n演變規(guī)律及其歸因研究

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

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

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

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

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

1 研究方法

1.1 Budyko方程

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

(1)

式中：E為流域實際蒸發(fā)量，mm。

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

1.2 參數n各影響因子及計算方法

1.2.1 參數n各影響因子

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

(2)

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

b. 逐年最大日降水量(Pmax)。根據中國氣象數據網上的氣象數據統計出流域逐年最大日降水量(Pmax)。

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

(3)

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

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

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

Y=aXexp(bX)

(4)

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

(5)

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

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

1.2.2 影響關系的研究方法

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

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

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

(6)

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

(7)

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

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

(8)

令

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

(9)

其中

2 結果與分析

2.1 參數n的變化趨勢

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

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

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

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

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

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

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

2.2 參數n各影響因子的演變趨勢

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

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

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

(a) 徑流

(b) 降水

(c) 潛在蒸發(fā)

(d) 溫度

(e) 逐年最大日降水量

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

(g) 降雨集中度

(h) 帕爾默干旱指數

(a) 人口

(b) 國民生產總值

(c) 有效灌溉面積

2.3 參數n變化對各影響因素的響應關系

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

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

(10)

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

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

2.4 影響因子變化對參數n的影響評估

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

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

3 結 論

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

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

c. 在黃河流域，人類活動對Budyko方程中參數n變化的貢獻率比氣候變化的貢獻率更大，而在氣候因子中，降水和潛在蒸發(fā)對參數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]崔啟武，孫延俊.論水熱平衡聯系方程[J].地理學報，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假設的海河流域蒸散發(fā)量和徑流量估算研究[J].中國農村水利水電，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假設預測長江流域未來徑流量變化[J].水科學進展，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ū)中的應用研究[J].中國農村水利水電，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]韓艷利，婁廣艷，葛雷，等.黃河流域與水有關生態(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].大氣科學,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].水電能源科學，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))