基于時變Budyko模型的滹沱河上游徑流變化歸因分析

苗正偉，路 梅，丁志宏

(1.河北水利電力學(xué)院 水利工程學(xué)院，河北 滄州 061001； 2.中水北方勘測設(shè)計研究有限責(zé)任公司，天津 300222)

1 研究背景

受氣候變化和人類活動等因素的影響，徑流的非平穩(wěn)性在全球范圍內(nèi)已普遍存在[1]，這對水資源的規(guī)劃與管理提出了嚴峻挑戰(zhàn)，識別徑流變化的影響因子是應(yīng)對這種挑戰(zhàn)的主要任務(wù)之一。近些年來，學(xué)者們通過水文模型法[2-4]、Budyko框架分析法[5-6]、經(jīng)驗統(tǒng)計法[7-8]等多種方法進行徑流影響量化分析，其中，基于Budyko模型方法相對簡單，且已被多項研究證實其結(jié)果是穩(wěn)定、可靠的[6，9]。敏感性分析法[10]和分解法[11]是基于Budyko框架分析流域年尺度上徑流變化原因的兩種主要方法。其中，分解法可以獨立地估計影響因子對徑流變化的貢獻率，而無需計算敏感性系數(shù)[12]，因此更為簡單、直觀。

在分解法中，影響徑流的因素主要是干旱指數(shù)(潛在蒸散發(fā)與降水的比值)和方程參數(shù)。很多研究[9，13-15]將這二者分別視為氣候變化和人類活動的表征。但是，干旱指數(shù)僅僅是潛在蒸散發(fā)與降水的比值，其他諸如降水的季節(jié)性、降水強度、氣溫等氣候特性也會對流域的產(chǎn)流過程產(chǎn)生影響，但它們并未包含于干旱指數(shù)，因此，干旱指數(shù)不能完全代表氣候變化。另一方面，單參數(shù)Budyko方程中的參數(shù)一般與地形、土壤、植被等流域特征有關(guān)[16]，它在某種程度上可表征人類活動對徑流的影響，但有研究表明，降水等氣候因素也可以通過改變方程參數(shù)而間接影響徑流變化，因此，該參數(shù)其實是人類活動和氣候變化對徑流影響的綜合表征[12，17]。鑒于此，基于單參數(shù)Budyko方程，本研究將徑流影響因素直接劃分為“干旱指數(shù)”和“流域特征”兩類，而不將之稱為氣候變化和人類活動。

以前基于Budyko假設(shè)的徑流變化歸因分析一般基于徑流序列的突變點劃分基準期和變化期，然后用各時段的水文氣候要素均值擬合Budyko方程，進而采用分解法或敏感性分析法進行徑流影響分離[13，18-19]。但該過程存在2點不足：①在基準期和變化期，Budyko方程參數(shù)均被視為常數(shù)，這意味著干旱指數(shù)和流域特征對徑流的影響是在突變點處突然發(fā)生的，這顯然不符合實際，兩者對徑流的影響是逐漸發(fā)生的，而不是一蹴而就的；②有些流域的徑流序列雖然呈現(xiàn)顯著趨勢，但沒有顯著突變點，此時，該方法適用性受限。

鑒于前人研究中存在的上述問題，本文構(gòu)建時變Budyko方程在年尺度上量化干旱指數(shù)和流域特征對徑流的影響。首先，在小波分析的基礎(chǔ)上對水文氣象序列進行滑動處理；然后以時間為協(xié)變量，描述Budyko方程參數(shù)及氣候因子的時間變化，進而構(gòu)建時變Budyko方程；最后基于分解法在年尺度及多年尺度上分離干旱指數(shù)和流域特征對滹沱河上游流域徑流變化的影響。

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

2.1 研究方法

2.1.1 時間序列的滑動平均處理

Budyko方程需要耦合水量平衡方程才能對徑流進行模擬，為了滿足流域蓄水變量可以忽略的前提條件，首先對原始水文氣象時間序列進行滑動平均處理。但對滑動窗口寬度的選擇，目前并無統(tǒng)一標準，主觀性較強[12，20-21]?？紤]到徑流變化歸因分析的本質(zhì)在于探究徑流影響因子對徑流長期變化趨勢的影響，這就需要消除時間序列的周期性波動和隨機性成分。因此，首先采用小波分析法對原時間序列進行周期性分析，然后進行滑動平均處理，這樣，既能弱化周期性成分的影響，又能滿足流域蓄水變化量近似為0的前提條件?；趶?fù)Morlet小波分析法[22]的結(jié)果表明：滹沱河上游流域的降水、潛在蒸散發(fā)、年徑流深系列的第一主周期尺度上的平均周期長度分別為8.3、7.8、10.2 a，在此基礎(chǔ)上，綜合考慮觀測數(shù)據(jù)長度以及便于表達的原則，最終確定對原序列進行11 a滑動平均處理，并將其記為

(1)

式中：i為原始時間序列的時序；xi為第i年的原始水文氣象要素值；t為11 a滑動平均處理之后的時間序列的時序；yt為11 a滑動平均處理后的第t年的水文氣象要素值。

2.1.2 氣候變量模擬

滑動平均處理在某種程度上弱化了周期性波動的影響，而時間序列的隨機項也可能會使得分析結(jié)果具有較大的不確定性。為消除這種不確定性，首先以時間為協(xié)變量模擬降水和潛在蒸散發(fā)的長期變化趨勢，將二者表示為：

(2)

(3)

2.1.3 時變Budyko方程

單參數(shù)Budyko方程可表示為[20]

(4)

式中：E、P分別為實際蒸散發(fā)和降水量多年平均值(mm)；φ=Ep/P，為干旱指數(shù)，Ep為潛在蒸散發(fā)量多年平均值(mm)；w為流域特征參數(shù)；B(·)表示Budyko-type方程。

根據(jù)水量平衡原理，有

Q=P-E-ΔS。

(5)

式中：Q為多年平均徑流深(mm)；ΔS為流域蓄水變化量，在較長時間尺度上可忽略不計。將式(4)代入式(5)可得

Q=P[1-B(φ，w)] 。

(6)

為了捕捉各因素在年尺度上對徑流的影響，可對w進行協(xié)變量分析，前人的研究常選擇氣溫、降水、GDP、人口、灌溉面積等作為解釋變量[12，17，20]，但這種方法計算復(fù)雜，而且對數(shù)據(jù)要求很高，當研究區(qū)較小時(如滹沱河上游流域)，很難獲取高精度、長系列的GDP、灌溉面積等數(shù)據(jù)。因此，引入時序t作為方程參數(shù)w的解釋變量，并且考慮線性和二次多項式2種形式，從而得到時變方程為：

Rt=Pt[1-B(φt，wt)] ；

(7)

wt=β0+β1t；

(8)

wt=β0+β1t+β2t2。

(9)

式中：Rt為基于時變Budyko方程所估計的第t年的徑流深(mm)；時序t=1，2，…；β0、β1、β2為方程參數(shù)，采用極大似然法估計。通過Budyko方程參數(shù)的協(xié)變量分析，既可以描述流域特征的時間變化，也可以基于Budyko假設(shè)實現(xiàn)年徑流的模擬。單參數(shù)Budyko方程形式多樣，其中，Choudhury-Yang方程形式簡單、應(yīng)用廣泛[23-24]，本研究基于Choudhury-Yang方程開展后續(xù)研究，其表達式為

(10)

將式(10)代入式(7)可得

(11)

當式(11)中的wt關(guān)于時序t的函數(shù)形式為常數(shù)、線性、二次多項式時，分別對應(yīng)于常參數(shù)Budyko方程、線性Budyko方程、二次多項式Budyko方程。

2.1.4 基于分解法的徑流影響分離

分解法的基本原理參見文獻[11]、文獻[12]、文獻[21]。在年尺度上，干旱指數(shù)和流域特征對徑流的年際影響量分別為：

(12)

(13)

式中Δa，t、Δo，t分別為從第t-1到第t年干旱指數(shù)和流域特征對年徑流的年際影響量。由此可知，干旱指數(shù)、流域特征從開始年份t0到第t年對年徑流的累積影響量分別為：

(14)

(15)

式中Aa，t、Ao，t分別為從初始年份t0到第t年干旱指數(shù)和流域特征對年徑流的累積影響量。

2.1.5 多年尺度上干旱指數(shù)和流域特征對徑流的影響

前人的研究大多在多年尺度上(即基準期和變化期)分離徑流變化，為便于比較，本文在年尺度的基礎(chǔ)上推導(dǎo)多年尺度干旱指數(shù)和流域特征對徑流的影響。設(shè)研究時域的初、末年份分別為t0、t2，年徑流系列的突變點為t1。那么，第t年的年徑流Rt可表達為初始年份年徑流Rt0與累積影響量的和，即

Rt=Rt0+Aa，t+Ao，t。

(16)

基準期(即t0～t1)的平均年徑流Rpre、變化期(即t1+1～t2)的平均年徑流Rpost表達如下：

(17)

同理，

(18)

基準期到變化期的徑流變化量為

Rdiff=Rpost-Rpre=

(19)

式(19)右側(cè)的第1項、第2項分別為從基準期到變化期干旱指數(shù)和流域特征對徑流的影響量。

2.1.6 其他方法

采用Sen斜率法[25-26]及Mann-Kendall趨勢檢驗法[27]估計水文氣象系列的趨勢及其顯著性；同時采用滑動t檢驗、Mann-Kendall突變檢驗[28]對水文氣象系列進行突變檢測以提高結(jié)果的可靠性，通過2種方法檢驗(顯著性水平均為0.05)的突變點才可判定為顯著的突變點[29]。

2.2 研究區(qū)概況與數(shù)據(jù)來源

滹沱河全長587 km，是海河水系子牙河的主要支流之一，其上游的小覺水文站控制面積超過14 000 km2，是其主要控制站。滹沱河上游有云中河、牧馬河、陽武河等眾多支流，是崗南水庫的主要產(chǎn)水區(qū)。崗南水庫總庫容為17.04億m3，以防洪為主，兼有供水、灌溉等功能。滹沱河上游流域多年平均降水量為471.3 mm，多年平均氣溫為8.2 ℃，具有超過70%的植被覆蓋度，屬溫帶季風(fēng)氣候[30]。

本文所用主要數(shù)據(jù)包括：①小覺水文站1961—2016年逐月徑流數(shù)據(jù)，摘錄于水文年鑒；②五臺山、原平、五寨、太原這4個國家基本氣象站1961—2016年逐日氣象數(shù)據(jù)，包括氣溫、降水、風(fēng)速、相對濕度等，數(shù)據(jù)來自中國氣象數(shù)據(jù)網(wǎng)，氣象站分布如圖1所示。

圖1 滹沱河上游流域概況Fig.1 Topography，river networks，and meteorological stations in the upper reaches of Hutuo River basin

3 結(jié)果與分析

3.1 水文氣象要素的時間變化趨勢及擬合

經(jīng)11 a滑動平均處理之后的水文氣象要素序列的Sen趨勢檢驗結(jié)果如圖2所示：降水、年徑流深都呈顯著(P<0.001)下降趨勢，潛在蒸散發(fā)與干旱指數(shù)則呈顯著(P<0.001)上升趨勢。經(jīng)M-K與滑動t突變檢測，四者均未發(fā)生突變。但原始年徑流序列在1981年發(fā)生突變，限于篇幅，原始水文氣象序列的趨勢及突變檢測結(jié)果省略。

圖2 滹沱河上游流域水文氣象要素變化趨勢Fig.2 Trends of hydrology and meteorology elements in the upper reaches of Hutuo River basin

考慮線性和二次多項式2種函數(shù)形式對降水和潛在蒸散發(fā)的長期變化趨勢進行模擬，并計算了貝葉斯信息準則BIC、均方根誤差RMSE和Nash-Sutcliffe模型效率系數(shù)NSE[31]3個指標，結(jié)果如表1所示。由此可見，降水、潛在蒸散發(fā)的二次多項式模型的BIC、RMSE、NSE 3個指標均明顯優(yōu)于線性模型，因此，本研究將采用二次多項式模擬降水和潛在蒸散發(fā)的長期變化趨勢。

表1 氣候因子的協(xié)變量分析結(jié)果Table 1 Results of covariate analysis of precipitation and potential evapotranspiration

3.2 時變Budyko模型的有效性

常參數(shù)(包括進行、不進行滑動平均處理2種情況)、線性、二次多項式4種Budyko方程的評價結(jié)果如表2所示。若未對原始水文氣象數(shù)據(jù)進行滑動處理，常參數(shù)Budyko模型(即表2中的“常參數(shù)A”)的RMSE高達22.039 6 mm，NSE僅為0.416 83，說明不能進行有效的年徑流模擬，主要原因可能是在年尺度上流域蓄水變量不為0。而在進行11 a滑動平均處理之后，常參數(shù)Budyko模型(即表2中的“常參數(shù)B”)的RMSE降低為6.708 1 mm，NSE提升到將近0.748 46，具有較好的年徑流模擬能力，說明將時間序列進行滑動平均處理是必要的。時變Budyko模型較之常參數(shù)Budyko模型表現(xiàn)更優(yōu)異，BIC、RMSE大幅減少，NSE大幅提升，高達0.88以上，表明時變Budyko模型擁有很好的年徑流模擬能力。2個時變Budyko模型相比，二次多項式模型略占優(yōu)勢，因此，本文以二次多項式Budyko方程進行后續(xù)計算。

表2 Budyko方程參數(shù)w協(xié)變量分析結(jié)果Table 2 Results of covariate analysis of parameter w of the Budyko-type equation

3.3 與傳統(tǒng)分段分解法的比較

因為原始觀測徑流在1981年發(fā)生突變，所以，以1981年為分割點劃分基準期和變化期，基于常參數(shù)和時變Budyko模型進行2個時段間的徑流變化分析，結(jié)果如表3所示。由表3可知：在對數(shù)據(jù)進行滑動平均處理的情況下，無論是徑流影響量還是各因素的貢獻率，二次多項式Budyko模型與常參數(shù)B模型的結(jié)果都很接近；若不進行滑動平均處理，常參數(shù)A模型所得的徑流減少量大于二次多項式Budyko模型，原因可能是二次多項式Budyko模型在某種程度上消除了徑流系列的周期性波動和隨機性成分，使得徑流的年際變化趨于平緩，盡管2種模型的貢獻率也有一定差異，但二者均顯示徑流減少的主要因素是干旱指數(shù)。以上對比分析表明，本文構(gòu)建的時變Budyko模型是可靠的。

表3 從基準期到變化期的徑流影響分離Table 3 Separation of annual runoff variation from prechange period to postchange period

3.4 滹沱河上游流域徑流變化歸因分析

基于二次多項式Budyko模型的年徑流變化歸因分析的結(jié)果如圖3所示。由圖3(a)可知，干旱指數(shù)和流域特征均導(dǎo)致年徑流持續(xù)遞減，至2011年，干旱指數(shù)和流域特征分別使得年徑流累積減少29.74 mm和12.14 mm。由圖3(a)曲線斜率變化可知：干旱指數(shù)對徑流的影響強度(即年際影響量的絕對值，mm/a，下同)逐年降低，而流域特征在20世紀90年代之前對徑流影響強度有較明顯的增大趨勢，之后基本保持穩(wěn)定。從貢獻率(圖3(b))的角度來說，干旱指數(shù)和流域特征對滹沱河上游流域年徑流減少的相對影響是一個此消彼長的過程，其中，干旱指數(shù)的貢獻率由86.79%持續(xù)降低至71.01%，而流域特征的貢獻率由13.21% 增加至28.09%，但是，滹沱河上游流域徑流減少的主導(dǎo)因素始終是干旱指數(shù)。

圖3 基于二次多項式Budyko模型的徑流影響分離Fig.3 Separation of the cumulative impacts of aridity index and catchment properties on runoff by the quadratic time-varying Budyko-type equation

但值得注意的是：圖3(b)顯示的流域特征對徑流變化的貢獻率越來越高，并非是反映了流域特征對徑流影響的強度持續(xù)增加，更沒有反映人類活動對滹沱河上游流域水循環(huán)的干預(yù)日益增強。實際上，從圖4可以看出，流域特征對年徑流的影響強度的變化相對較小，只在20世紀90年代之前有過一定幅度的增加，在20世紀90年代之后基本保持穩(wěn)定，2000年之后甚至呈現(xiàn)減弱趨勢，反觀干旱指數(shù)，其對年徑流的影響強度在整個研究時域持續(xù)減小，且減幅非常明顯，換言之，圖3(b)所顯示的流域特征對徑流變化貢獻率的增加，主要是因為干旱指數(shù)對年徑流的影響強度持續(xù)、明顯減弱所致，尤其是在20世紀90年代之后。這與現(xiàn)實情況相符，在20世紀90年代之前，經(jīng)濟的快速增長所導(dǎo)致的下墊面變化(比如土地利用變化)、水工程建設(shè)等均對水循環(huán)的干預(yù)有所增強，表現(xiàn)為圖4所示的20世紀90年代之前流域特征對年徑流的影響強度增加；而在20世紀90年代之后，尤其是進入21世紀之后，在可持續(xù)發(fā)展、人與自然和諧相處等理念的影響之下，人們越來越重視水資源的可持續(xù)利用，相繼采取了退耕還林、退耕還湖等一系列措施，因此，盡管城鎮(zhèn)化建設(shè)、灌區(qū)建設(shè)等一系列因素還在影響水循環(huán)，但是年徑流的減小程度應(yīng)是受到了一定程度的遏制，正如圖4所示；2000年之后，流域特征對年徑流的影響強度在基本保持穩(wěn)定的同時還有所減弱(計算結(jié)果顯示，流域特征對年徑流的影響強度在1999年達到最大，值為0.302 mm，此后呈小幅減少趨勢，至2011年減小至0.291 mm)。這種年際尺度的變化環(huán)境下的水文效應(yīng)，是基于常參數(shù)Budyko方程的傳統(tǒng)分解法或敏感性方法所不能揭示的。

圖4 基于二次多項式的Budyko方程的年徑流年際影響分離Fig.4 Separation of the annual impacts of aridity index and catchment properties on runoff by the quadratic time-varying Budyko-type equation

4 討論與結(jié)論

4.1 討 論

4.1.1 時變Budyko模型的優(yōu)點

本研究基于Wang和Hejazi[11]于2011年提出的分解法構(gòu)建了一種年徑流變化歸因分析的新方法，與傳統(tǒng)的基于單參數(shù)Budyko方程的分解法或者敏感性系數(shù)法相比，新方法的優(yōu)點在于：

(1)對Budyko方程參數(shù)w進行了協(xié)變量分析，從而可以描述流域特征的時間變化，而傳統(tǒng)方法假定在某固定時段內(nèi)(如變化期或者基準期)流域特征保持不變，這與現(xiàn)實不符。

(2)基于w的協(xié)變量分析所構(gòu)建的時變Budyko方程，可以在年尺度上分離干旱指數(shù)和流域特征對徑流的影響，這可以幫助我們更深入地理解變化環(huán)境下的水文響應(yīng)。

(3)時變Budyko方程較之常參數(shù)方程擁有更強的年徑流模擬能力。

4.1.2 模型的改進

由于資料限制，本研究僅以時間為協(xié)變量構(gòu)建了時變Budyko模型，雖然這種做法在不少研究中均有采用，但其缺乏物理基礎(chǔ)，其穩(wěn)定性和可靠性有待進一步驗證，因此，在資料允許的情況下，建議選擇具有實際物理意義的解釋變量。此外，僅以時間為協(xié)變量，仍然無法在Budyko方程參數(shù)中分離人類活動和氣候變化對徑流的影響，這使得我們無法精確分離氣候變化和人類活動的水文效應(yīng)。因此限制了結(jié)果的適用性，如果分別選擇代表氣候變化和人類活動的兩類變量為解釋變量，如用北大西洋濤動NAO、太平洋年代際振蕩IPO、大西洋多年代際振蕩AMO等大尺度氣候因子反映氣候變化，而用GDP、人口、耕地面積、土地利用等因子反映人類活動，由此對w進行協(xié)變量分析，即可從中分離氣候變化和人類活動對徑流的影響，這不僅使得結(jié)果的可解釋性更強，而且其實際意義也會更加明顯。

4.2 結(jié) 論

本研究基于Budyko方程參數(shù)w的協(xié)變量分析，構(gòu)建了時變Budyko模型，進一步改進了分解法，實現(xiàn)了在年尺度上分離干旱指數(shù)和流域特征對徑流的影響。將該方法應(yīng)用于滹沱河上游流域，所得的主要結(jié)論如下：

(1)時變Budyko模型的NSE遠高于傳統(tǒng)常參數(shù)模型，因此時變Budyko模型大大提高了年徑流的模擬能力?；诰€性方程的時變Budyko模型具有和二次多項式Budyko模型相似的結(jié)果。

(2)基于時變Budyko模型的分段研究的結(jié)果與傳統(tǒng)常參數(shù)Budyko模型的結(jié)果相近，一定程度上說明了本文所提出的方法是可靠的。

(3)干旱指數(shù)和流域特征均導(dǎo)致滹沱河上游流域年徑流減少，其中，干旱指數(shù)的貢獻率由86.79%持續(xù)遞減到71.01%，而流域特征的貢獻率由13.21%持續(xù)增加至28.09%，因此，滹沱河上游流域徑流減少的主導(dǎo)因素是干旱指數(shù)。

(4)干旱指數(shù)對年徑流的影響強度逐年降低，而流域特征的影響強度相對穩(wěn)定，因此，流域特征對年徑流變化貢獻率的增加，主要是干旱指數(shù)對年徑流的影響強度明顯減弱所致。