基于Budyko假設的黃河下游非一致性徑流模擬

叢培月， 牟獻友， 冀鴻蘭， 張寶森

(1.內蒙古農業大學 水利與土木建筑工程學院， 內蒙古 呼和浩特 010018； 2.黃河水利委員會 黃河水利科學研究院， 河南 鄭州 450003)

1 研究背景

黃河是我國第二大河流，多年平均徑流量為580×108m3，約占全國總徑流量的2%[1]。 2014年，IPCC全球氣候變化第五次評估報告顯示：全球氣候變暖已被廣泛認可，變暖體現在地表溫度上升、極端氣候事件頻率增加等方面，1880-2012年全球平均溫度已升高0.85℃，氣候變化的嚴重性遠超于以往的認知[2]。氣候變化會改變降水及蒸散發，造成水文氣象序列出現非一致性[3]。水文氣象的非一致性會影響陸地水循環過程，驅動徑流量的變化[4]。由于歷史資料的缺失，黃河徑流量演化規律的分析和判別將受到一定的限制[5]，因此，流域徑流關系演變分析、徑流模擬應運而生，成為當前水文水資源領域研究的前沿與熱點之一，也是國際上一個具有挑戰性的問題[6]。

近50年來，黃河徑流量較以往有著顯著差異，國內外許多學者針對中國各江河流域的徑流模擬，進行了定量分析研究[7-8]。顧朝軍等[9]分析了黃河陜縣站1919-1997年天然徑流量變化，指出黃河天然徑流量年際變化存在明顯的階段性及周期性特征；王雁等[10]通過對徑流和氣象觀測數據的統計分析，發現黃河流域由于氣候變化引起了徑流量的減少。張建云等[11]提出，自20世紀70年代以來，受氣候變化的影響，黃河中游河川徑流呈現明顯的下降趨勢。孫福寶等[12]基于Budyko假設的傅報璞公式檢驗了干旱區水熱耦合理論。熊立華[13]分析了Budyko假設在中國南部濕潤地區的適用性，并提出Budyko公式改進方法，提高了徑流模擬精度。張麗梅等[14]利用Budyko假設分析渭河年徑流量變化特征。對徑流的模擬大多是基于水文氣象序列概率分布及其參數隨時間不變的假定[15]，張波等[16]指出，在環境變化的大背景下，一致性假設的水文理論無法精準地揭示水資源演變規律，在此背景下提出了非一致性年徑流模擬方法，使結果更加真實、可靠；唐雄朋等[17]探討了雅魯藏布江拉孜以上流域非一致性徑流過程模擬與預測，徑流模擬精度較高；鄧曉宇等[18]對氣候變化和人類活動綜合影響下的撫河流域進行徑流模擬研究，發現納西系數達到0.82以上。Milly等[19]指出，在氣候變化和人類活動的影響下，年徑流序列通常是非一致性的，導致許多現有方法徑流模擬的精度較低；Merz等[20]根據奧地利273個盆地41 a的歷史數據，通過校準模型參數來進行非一致性徑流模擬，提高了徑流模擬精度；Pathiraja等[21]通過使用集合卡爾曼濾波器同時估計模型參數，對西澳大利亞實驗區進行徑流預測，得出預測結果可被接受。

目前，對基于Budyko公式改進的Zhang模型在黃河下游流域的適用性研究尚有空缺。本文利用灰色關聯法[22]與滑動樣本熵等方法對黃河下游水文氣象序列進行非一致性識別，將改進的Zhang模型應用于黃河下游，研究其在非一致性水文氣象序列中模擬徑流的適用性。

2 資料與方法

2.1 研究區概況

河南鄭州桃花峪以下的黃河河段為黃河下游，地處中國東部，流經河南、山東，于山東墾利縣注入渤海。黃河下游段位于東經110°～120°、北緯30°～40°之間，河長約786 km，為典型游蕩型河段，河道寬淺，流域面積約2.3×104km2，占全流域面積的3%左右，總落差93.6 m，平均比降0.12%。黃河下游屬于溫帶季風氣候，夏季高溫多雨(夏季季平均降水量為371.76 mm，夏季季平均氣溫為26.39℃)，冬季寒冷干燥(冬季季平均降水量為56.71 mm，冬季季平均氣溫為1.34℃)，多年平均降雨量約617.52 mm，年內降水分布不均，多年平均氣溫約14.49℃，多年平均徑流量約253.45×108m3，多年平均日照時數約2 548.77 h。根據降雨徑流分帶，大部分為半濕潤地帶。

2.2 數據來源

收集了黃河下游 7 個水文站(花園口、夾河灘、高村、孫口、艾山、濼口、利津)1979-2015年共37 a的逐日降雨、逐日氣溫、逐日流量等水文資料以及周邊4個氣象站(鄭州、開封、濟南、惠民)的1979-2015年逐日日照時數等氣象資料，數據均來源于黃河流域歷年水文年鑒、中國氣象數據網(http://data.cma.cn)。

2.3 研究方法

本文利用灰色關聯法分析徑流復雜現象與降水、氣溫發展態勢的關聯、接近的程度，從而確定徑流形成的主控因素。由于水文氣象序列的非一致性會直接影響降水和潛在蒸散量，造成氣候-徑流關系的方式不同于過去的觀察[23]，因此在非一致性水文氣象序列下，選用基于歷史經驗的水文模型將影響模擬精度。本文首先進行水文氣象序列的一致性驗證，利用累積距平法和滑動平均法進行初步檢驗，并基于滑動樣本熵法[24]進行突變分析，聯合多種分析方法對水文氣象一致性進行診斷。由于實際蒸發量難以通過觀測直接獲得，通常利用潛在蒸散發量來計算實際蒸發量[25]，本文采用在濕潤與半濕潤區優先選擇的Hargreaves[26]法估算潛在蒸發量，通過潛在蒸發量利用Budyko假設計算蒸發量，最后利用水量平衡法反算徑流量。

2.3.1 Budyko假設 Budyko(1974年)假設基本原理是：假定流域內降水和潛在蒸散發量決定了多年平均蒸散發，蒸散發量為降水與潛在蒸散發量之間的函數，即可以對降水、潛在蒸散發和徑流之間定量關系進行較好的擬合[27]。但是由于植被類別、土壤水力特性及地形的區別，使得Budyko假設的廣泛性受到限制，Zhang等在Budyko基礎上引進模型參數ω，改進了Budyko假設公式。

Budyko假設可簡單表達為：

(1)

(2)

Zhang公式：

(3)

式中：E為年平均蒸發量，mm；P為年平均降水量，mm；ω為植被水分利用系數；F為干燥指數；E0為年平均潛在蒸散發量，mm。

水量平衡法：

R=P-E

(4)

式中：R為徑流深，mm；P為降水量，mm；E為蒸發量，mm。

2.3.2 模擬精度分析 本文采用納西系數(NSE)對年徑流量的模擬精度進行評價。NSE∈(-∞,1)，其中當NSE=1時，模擬精度最優；NSE∈(0.01,1)時，模擬精度為可接受。

(5)

3 結果與分析

3.1 降水、氣溫與徑流灰色關聯度

以黃河下游7個水文站為代表站，設定降雨徑流灰色關聯度為γ1，氣溫徑流灰色關聯度為γ2。通過計算得出1979-2015年7個水文站的γ1、γ2，如表1所示。7個站的γ1、γ2均大于0.5，表示徑流量與降水和氣溫關聯度較高，即氣溫和降水均為黃河下游徑流量的主控因素。此結果與張國宏等[28]的研究結論相似，黃河流域徑流量與氣溫和降水均有較好的相關性，其中各站均為γ2>γ1，表明在黃河下游徑流形成過程中，氣溫與徑流量的關聯度較大，對徑流的影響所占的權重大于降水。

表1 1979-2015年降水、氣溫與徑流的灰色關聯度

3.2 降水、氣溫變化趨勢分析

表2為1979-2015年黃河下游7個水文站的降水與氣溫突變年份。由表2可知，1979-2015年黃河下游流域降水和氣溫突變點不一致，如夾河灘站降水時間序列發生突變的時間為1999年，氣溫時間序列發生突變的時間為1993年；突變時間均在20世紀90年代前后，水文氣象時間序列突變時間具有較高的同步性。本文以花園口站降水、氣溫變化趨勢分析為例，其余6個站不再贅述。

表2 1979-2015年黃河下游7個站降水與氣溫突變年份

黃河下游1979-2015年間多年平均降水量為617.52 mm，降水量總體偏少。由線性擬合結果(圖1(a))，降水量傾向率為-4.5 mm/10a，呈顯著性減少趨勢(p<0.05)，其中降水量最大值出現在1983年(962.94 mm)，最小值出現在1997年(368.96 mm)，極值比1.61。從降水量累積距平曲線(圖1(b))可以看出：花園口站在1979-2015年間降水量變化呈現出先波動下降、后波動上升的趨勢，即降水量變化分為兩個階段：第1階段1979-2002年，該時期降水量平均值為596.24 mm；第2階段2003-2015年，該時期降水量平均值為639.39 mm，第1階段到第2階段的變化呈現了降水枯-豐轉變的過程，表明降水時間序列已發生變化。

1979-2015年，花園口站多年平均氣溫為15.14℃，年平均氣溫傾向率為0.6℃/10a。由線性擬合(圖2(a))可看出，花園口1979-2015年38 a的年平均氣溫呈現顯著上升趨勢(p<0.05)。其中最高值和最低值分別出現在2015和1984年，分別為16.67和13.68℃，極值比為0.22。與降水量的變化相比，氣溫的變化較為平緩。從圖2(b)的累積距平曲線可以看出，花園口站年平均氣溫表現出先下降后上升的發展趨勢，自1997年起，花園口站氣溫開始由低溫期向高溫期過渡，表明氣溫時間序列也已發生變化。

圖1 花園口站1979-2015年年降水量時間序列變化趨勢

圖2 花園口站1979-2015年年均氣溫時間序列變化趨勢

降水量的滑動樣本熵(SampEn)呈現逐漸增加的趨勢(圖3(a))，降水量變化分為3個階段： 第1階段(1979-1986年)、第2階段(1987-1993年)、第3階段(1994-2015年)，即發生了兩次突變，對應的突變點年份分別為1986和1993年，突變點前后的SampEn序列差異性顯著。氣溫的SampEn值在0.6～0.9之間(圖3(b)，其變化也呈現出3個階段：第1階段(1979-1991年)、第2階段(1992-1994年)、第3階段(1995-2015年)，對應的突變點年份分別為1991和1994年，氣溫時間序列的SampEn值均大于降水量序列，表明氣溫系統的復雜程度高于降水系統。

圖3 花園口站1979-2015年年降水量、氣溫時間序列SampEn分析

整體而言，從發生突變的年份來看，年降水量和年平均氣溫發生突變的時間有較高的對應關系。由以上水文序列發生突變可知，水文序列的統計規律并非一致，因此水文序列不滿足一致性。非一致性在徑流、降水和氣溫中已被廣泛認定。

3.3 潛在蒸散發量估算

利用Hargreaves方法在黃河下游4個水文站點(花園口、夾河灘、艾山、利津)中進行潛在蒸散發量估算，估算結果見圖4。從圖4可知，具有相似的規律特征，年平均潛在蒸散發量從850～1 200 mm不等，潛在蒸散發量均呈現隨時間減少的顯著趨勢(p<0.05)，4個站點的氣候傾向率分別為-20.8、-34.8、-76.3和-38.2 mm/10a，其中利津站在1990年前后出現了小幅回升，1979-1993年潛在蒸散發量呈現下降的趨勢，氣候傾向率為-34.60 mm/10a,2010-2014年潛在蒸散發量呈現上升的趨勢，氣候傾向率為51.79 mm/10a，這些變化符合近50年來中國潛在蒸散量的變化趨勢[29]。

3.4 Budyko假設徑流模擬精度分析

由于Budyko假設沒有考慮人類活動和氣候變化對徑流量的影響，因此利用降水量與潛在蒸散發量，同時基于Budyko假設改進的Zhang模型求得蒸發量，通過水量平衡法反算徑流量。本文以1982年為例，模擬計算了逐月徑流量并與實測值進行對比，如圖5所示。分析圖5可知：花園口站在8月模擬結果誤差最大(NSE<0)，在2月(NSE=0.81)和5月(NSE=0.97)模擬結果較為理想，夾河灘站和艾山站分別在7、8月模擬結果誤差最大(NSE=0)，利津站在4月(NSE=0.97)、5月(NSE=0.99)、6月(NSE=0.86)模擬結果較為理想，模擬精度較高。表3給出了黃河下游4個典型區1982年月均徑流模擬NSE檢驗結果，4個典型區的NSE均小于0。綜合以上研究結果，在月尺度徑流模擬中各個月模擬值與實測值之間誤差較大，模擬精度過低，因此Zhang公式不適用于月尺度徑流模擬中。

利用Zhang公式對4個水文站的1979-2015年徑流過程進行模擬，結果如圖6所示，模擬的徑流變化過程與實際的徑流變化過程基本吻合，能較好地反映徑流的變化趨勢和數值，但在2000年后個別年份模擬值與實際值差別較大，如花園口站2005年(NSE=-0.62)、2011年(NSE=-2.27)，艾山站2010年(NSE=-4.77)。但整體來看，4個典型區徑流模擬的納西系數均大于0.9(如表4所示)，說明模型精度較高，模擬結果可以接受，Zhang公式在年徑流量模擬方面具有一定的適用性。

圖4 黃河下游4個水文站1979-2015年潛在蒸散發量變化

圖5 1982年黃河下游4個水文站月徑流量觀測值與模擬值變化分析

圖6 1979-2015年黃河下游4個水文站年徑流量觀測值與模擬值變化分析

表3 1982年黃河下游4個水文站月徑流量NSE對比

表4 1979-2015年黃河下游4個水文站年徑流量NSE對比

徑流量模擬的差異性受諸多因素影響，水文氣象序列的選擇、人類活動和氣候變化等可直接影響徑流量模擬精度。本文利用基于Budyko假設改進的Zhang模型對黃河下游進行徑流量模擬，經過驗證表明模型精度較高，適應于黃河下游流域非一致性水文氣象序列的徑流量模擬，可為黃河下游徑流量分析、預報提供可靠的數據支持。

4 結 論

(1)黃河下游降水量呈現出顯著減少趨勢(p<0.05)，氣溫呈現出顯著增加趨勢(p<0.05)，利用滑動樣本熵從動力學的角度判別出降水量、氣溫均在20世紀90年代前后出現突變，且突變年份具有較高的同步性。

(2)黃河下游徑流變化過程受多因素影響，其中γ1和γ2均大于0.5，表明降水量和氣溫與徑流量的關聯程度均較高，是影響徑流量變化的主要因素。

(3)基于Budyko假設改進的Zhang模型模擬了黃河下游非一致性水文氣象序列下的年徑流量，NSE達到0.9以上，模擬效果良好，表明Zhang模型在黃河下游非一致性水文氣象序列中模擬徑流量具有可適性。