王思奇，朱宏鵬，魯 帆*，江 明，周毓彥

(1.中國水利水電科學研究院，北京 100038；2.河北省保定市水政監(jiān)察支隊，河北 保定 071000)

隨著人口增加和社會經(jīng)濟的快速發(fā)展，人類活動對流域下墊面及產(chǎn)匯流條件變化的影響不斷加劇[1-2]。已有學者針對北方流域暴雨洪水演變特征問題開展了系列研究，梁艷琴[3]、陳旭等[4]的研究表明海河流域不同時段洪水總量、洪峰流量和單位線峰值普遍呈顯著的減少趨勢；常遠遠等[5]對黃河湫水河流域2015年稀遇暴雨洪水與20世紀60—80年代的13場洪水進行對比，得出相同降雨條件下的洪峰和次洪洪量較歷史洪水均減??；雷成茂等[6]將2016年西柳溝流域的洪水過程與洪峰流量居歷史前3位的大洪水相關(guān)特征平均值進行對比，發(fā)現(xiàn)其洪峰流量、徑流系數(shù)均偏小；金雙彥等[7]分析皇甫川流域次洪特征值變化特點，結(jié)果表明：近些年洪水發(fā)生次數(shù)和年最大洪峰流量高值出現(xiàn)次數(shù)均有減少趨勢。人類活動對洪水演變的影響主要來源于土地開發(fā)利用、水土保持、植樹造林、水利工程、城鎮(zhèn)化建設(shè)等，這些措施不僅使土地利用等環(huán)境要素發(fā)生變化，而且能使流域產(chǎn)洪次數(shù)、地表徑流模數(shù)和徑流系數(shù)減小[8-11]。對于變化環(huán)境下的頻率變化，已有學者采用時變參數(shù)模型分析了水文要素的非一致性，結(jié)果表明變化環(huán)境下應當充分考慮徑流的非一致性特點，時變參數(shù)模型能夠反映水文序列的非一致性特征[12-13]。針對海河流域大清河山區(qū)洪水特征變化的研究[8-10,14]表明：歷年最大洪峰、次洪水量呈減少趨勢，相同降雨量產(chǎn)生的次洪量減少，下墊面對流域徑流系數(shù)的影響較大。為更深入地揭示該流域暴雨洪水的演變規(guī)律，本文利用秩次相關(guān)檢驗法研究了大清河山區(qū)典型站點洪水過程中洪峰流量值的整體趨勢，對比流域特大洪水的特征變化，在此基礎(chǔ)上，采用多時段對比法分析次洪洪峰流量、次洪水量、徑流系數(shù)、產(chǎn)流閾值的變化情況，利用時變參數(shù)模型GAMLSS對洪水頻率的演變特征進行分析，研究結(jié)果可為流域防洪規(guī)劃與管理提供支撐。

1 研究區(qū)域及數(shù)據(jù)

大清河水系位于海河流域中部，東西長約275 km，南北寬約200 km，流域面積約4.3萬km2，主要由南、北兩支組成，流入白洋淀的支流為南支，流入東淀的支流為北支，北支主要為拒馬河。其洪水主要由汛期暴雨造成，暴雨中心常出現(xiàn)在阜平、司倉及紫荊關(guān)一帶，自有水文記載以來，1963年8月暴雨為最大，實測6天洪量為64.3億m3，南支占80%，北支占20%，全水系30 d洪水總量79.98億m3，其中南支占78%，北支占22%。紫荊關(guān)站是北支拒馬河上游的主要控制站，流域面積為1 760 km2。阜平站是南支沙河上游的主要控制站，控制面積2 210 km2。該區(qū)域?qū)俅箨懶约撅L氣候，四季分明，春季干旱多風，夏季炎熱多雨，土地利用類型主要以耕地、林地為主，是暴雨洪水多發(fā)地區(qū)之一，洪水暴漲暴落。

水文氣象數(shù)據(jù)主要來源于《水文年鑒海河流域水文資料》第4冊大清河水系以及中國氣象數(shù)據(jù)網(wǎng)，選定1960—2018年的洪水水文要素摘錄表及降水量表進行統(tǒng)計，進而得到次洪過程及次洪降水量。其中紫荊關(guān)水文站控制流域?qū)?個雨量站為艾河村、團圓村、東團堡、石門、王安鎮(zhèn)、平頂山，阜平水文站控制流域?qū)?0個雨量站為銀廠、莊旺、站上、神堂堡、不老臺、西彎、砂窳、龍泉關(guān)、橋南溝、阜平，大清河山區(qū)水系及上述水文氣象站點的位置見圖1。

2 研究方法

首先選擇年最大值取樣法選擇年最大洪峰流量作為研究對象、利用斜線分割法計算其次洪的洪水量、利用算術(shù)平均法計算選擇的洪水過程相應的面雨量，然后通過Kendall秩次相關(guān)檢驗法對年最大洪峰流量序列的趨勢進行判別，通過多時段對比法得到次洪降水量、次洪水量、徑流系數(shù)及產(chǎn)流閾值的變化情況，采用線性趨勢法分析次洪過程中徑流系數(shù)與降水量的相關(guān)關(guān)系，采用GAMLSS模型對變化環(huán)境下的洪水頻率演變規(guī)律進行分析，計算流程見圖2。

利用斜線分割法計算年最大洪峰對應場次洪水的洪量時，將洪水起漲點和地面徑流終止點連成一條直線，直線以下為地下徑流，直線以上的地表徑流作為次洪洪水總量[15]，計算見式(1)。選擇Kendall秩次相關(guān)檢驗法對研究的水文序列構(gòu)建統(tǒng)計量U，見式(2)。對比計算的統(tǒng)計量值與一定顯著性水平下的臨界值，可判斷水文序列的變化趨勢[16]。徑流系數(shù)指選定時段內(nèi)的徑流深度R與同時段內(nèi)降水深度P的比值[17]，它代表降水量中有多少水轉(zhuǎn)化成徑流[18-19]，徑流系數(shù)越高，表明地區(qū)降雨產(chǎn)流能力越強。產(chǎn)流是指降雨量扣除損失形成凈雨的過程，其中降雨損失包括植物截留、下滲、填洼與蒸發(fā)，且以下滲為主，典型的產(chǎn)流機制包括蓄滿產(chǎn)流和超滲產(chǎn)流，某些地區(qū)水文模型的構(gòu)建兼具2種產(chǎn)流機制[20]。降雨產(chǎn)流閾值是指接受降雨的下墊面能夠產(chǎn)生地表徑流的最小降雨量[21-22]。因本文采用最大值取樣法選取年最大洪峰對應的場次洪水，因此確定所選擇洪水過程線中的流量起漲點tQ，根據(jù)降水量摘錄表中tQ時刻之前各雨量站的時段降雨，求均值得到的區(qū)域降水量P即為產(chǎn)流閾值。

(1)

(2)

式中W——次洪洪水總量；S——地表徑流與地下徑流之和；QQ、QZ——起漲點tQ、終止點tZ所對應的流量值，m3/s；M——徑流序列(Qt)中(Qi,Qj,j>i)中所有對偶值(Qi,Qj,j>i)中Qi>Qj出現(xiàn)的個數(shù)；N——徑流序列的總長度。

在一定顯著性水平α下的雙尾檢驗臨界值記為Uα/2，若|U|>Uα/2，則說明徑流序列的變化趨勢顯著，否則變化不顯著，且U>0時序列呈上升趨勢，反之呈下降趨勢。

(3)

(4)

本文選用Gumbel(GU)、Weibull(WEI)、Gamma(GA)、Logistic(LO)、Log Normal(LOGNO)、Normal(NO) 6種兩參數(shù)連續(xù)分布作為備用分布函數(shù)，并考慮位置參數(shù)(μ)、尺度參數(shù)(σ)為常數(shù)和其隨協(xié)變量(時間t)呈現(xiàn)線性變化的情況，合計共24種組合方式。計算不同擬合方案對應的AIC值，選定AIC值最小的分布函數(shù)及參數(shù)變化特征作為站點洪峰流量序列的最優(yōu)擬合模型，對最優(yōu)擬合模型進行模型殘差評價驗證擬合效果，根據(jù)確定的最優(yōu)模型及其參數(shù)，計算站點指定時間及指定百分數(shù)情況下對應的設(shè)計洪峰流量值，并與一致性條件下的設(shè)計值對比。

3 洪水演變規(guī)律

3.1 洪水過程變化

對紫荊關(guān)、阜平水文站1960—2018年實測年最大洪峰流量序列的變化趨勢用Kendall秩次相關(guān)檢驗法分析。根據(jù)式(2)計算得U分別為-5.22、-3.87，2個值均大于顯著性水平α=0.1下的雙尾檢驗臨界值1.64，且均大于顯著性水平α=0.05下的雙尾檢驗臨界值1.96，說明在顯著性水平α=0.1及α=0.05下徑流序列的變化趨勢顯著，且呈下降趨勢。

考慮資料序列長度以及全國水資源評價工作的時間節(jié)點，將資料系列劃分為3個階段(1960—1979、1980—1999、2000—2018年)，對3個階段的洪水進行統(tǒng)計，見表1。紫荊關(guān)以上流域洪峰的時期均值先下降后上升，洪量均值一直下降。同1960—1979年比較，1980—1999年的洪峰、洪量分別下降71.1%、68.7%；同1980—1999年相比，2000—2018年的洪峰上升9.5%，洪量反而下降23.8%。阜平以上流域洪峰和洪量均值基本呈下降趨勢。同1960—1979年比較，1980—1999年的洪峰、洪量分別下降了45.2%、48.7%；同1980—1999年相比，2000—2018年的洪峰下降35.1%，洪量下降58.7%。

表1 不同時期洪峰流量、洪水量變化均值

對于洪峰值在2 000 m3/s以上的特大洪水，2個流域均出現(xiàn)2次，紫荊關(guān)的特大洪水發(fā)生在1963、2012年，阜平發(fā)生在1963、2016年。相較于1963年的特大洪水，紫荊關(guān)2012年、阜平2016年的特大洪水均屬于單峰洪水，具有短歷時強降雨的特點，且峰值、面雨量均減少，見圖3。紫荊關(guān)站在1963、2012年的洪峰值分別為4 490、2 156.5 m3/s，阜平站1963、2016年的洪峰值分別為3 380、2 020 m3/s。紫荊關(guān)站1963年特大洪水峰現(xiàn)時間歷時79 h，面雨量達到380 mm，次洪水量折合徑流深為213 mm。紫荊關(guān)2012年特大洪水峰現(xiàn)時間歷時9 h，面雨量為166 mm，次洪水量折合徑流深為39.6 mm，與1963年特大洪水相比，峰值減少約52%，面雨量均值下降214.5 mm。阜平站1963年特大洪水峰現(xiàn)時間歷時109 h，面雨量達到541 mm，次洪水量折合徑流深為205.6 mm。阜平站2016年特大洪水峰現(xiàn)時間歷時4.6 h，面雨量為76 mm，洪量折合徑流深為16.5 mm。與1963年特大洪水相比，峰值減少40%，面雨量下降465 mm。

3.2 徑流系數(shù)變化

根據(jù)年最大洪峰選取的場次洪水整理的次洪降水量以及次洪洪量，求得2個流域1960—2018年逐年的徑流系數(shù)值，并計算3個階段徑流系數(shù)的均值，見圖4、5。2個代表區(qū)域的徑流系數(shù)階段均值均下降，點據(jù)下移趨勢明顯。紫荊關(guān)以上流域1960—1979年徑流系數(shù)均值為0.126，相比于1960—1979年，1980—1999年徑流系數(shù)下降顯著，減少量近0.07，下降幅度達到55.9%；同1980—1999年相比，2000—2018年徑流系數(shù)減少0.019，下降幅度為34.4%。阜平以上流域1960—1979年徑流系數(shù)均值為0.138，相比于1960—1979年，1980—1999年徑流系數(shù)減少近0.03，減少比例為22.8%；同1980—1999年相比，2000—2018年徑流系數(shù)減少約0.05，減少比例達到45.6%，下降顯著。

鄭江坤等[26]對于蒲江縣朝陽水庫的4個坡耕地徑流小區(qū)表明，各小區(qū)徑流系數(shù)均隨降雨等級增加呈增加趨勢，因此對于徑流系數(shù)的變化，除對比階段均值的變化情況外，對區(qū)域次洪徑流系數(shù)與降水量的相關(guān)關(guān)系進一步分析。以降水量作為橫坐標，點繪2個區(qū)域整個序列以及3個階段的徑流系數(shù)值，分析徑流系數(shù)的變化趨勢，得到徑流系數(shù)與降水量、時期的相關(guān)關(guān)系，見圖6、7。根據(jù)徑流系數(shù)隨降水量的變化圖可以看出，線性趨勢線的斜率均大于0，即次洪徑流系數(shù)與降水量呈正相關(guān)，降水初期發(fā)生蓄滲，不會立即形成徑流，隨著降水量增加，滿足截留、填洼等條件后，進入產(chǎn)流過程，但是當降雨等級較小時，轉(zhuǎn)化為地表徑流部分較小，而大暴雨和特大暴雨在一定程度上削弱了土壤前期含水量的影響，因此次洪徑流系數(shù)隨降水量的增加而增加。紫荊關(guān)以上流域的年序列變化線性趨勢線的斜率為0.001 2，大于0，說明隨著降水量的增加，徑流系數(shù)呈上升趨勢；對不同階段進行分析，3個階段的線性趨勢線的斜率也均大于0，分別為0.001 4、0.000 9和0.000 2，說明徑流系數(shù)隨降水量的增加而增加。阜平以上流域的年序列變化線性趨勢線的斜率為0.000 9，與紫荊關(guān)以上流域類似，隨著降水量的增加，徑流系數(shù)呈上升趨勢；對不同階段進行分析，3個階段的線性趨勢線也均大于0，分別為0.002 4、0.002 3和0.001 3，進一步驗證相同降水量的條件下，徑流系數(shù)值下降，與紫荊關(guān)以上流域類似，在降水量相同的情況下，徑流深減少幅度大。

3.3 產(chǎn)流閾值變化

根據(jù)年最大洪峰值對應的次洪場次，對研究區(qū)域產(chǎn)流閾值的年序列變化情況進行統(tǒng)計，見圖8、9。根據(jù)圖8，紫荊關(guān)以上流域產(chǎn)流閾值整體呈上升趨勢；3個時期的均值分別為8.4、12.7、13.7 mm，同1960—1979年相比，1980—1999年的產(chǎn)流閾值增加4.3 mm，增加幅度達到50.6%；2000—2018年與1980—1999年相比，增加約1 mm，增幅僅為8.1%。阜平以上流域的產(chǎn)流閾值變化情況與紫荊關(guān)以上流域類似，階段均值均在增加，3個時期的均值分別為10.3、12.5、16.0 mm，1980—2000年同1960—1979年相比，增加2.2 mm，增加幅度為21.7%；2000—2018年與1980—1999年相比，增加量為3.5 mm，增加比例達到27.9%。對2個區(qū)域的產(chǎn)流閾值變化情況進行比較，產(chǎn)流閾值均呈上升趨勢，但是不同區(qū)域階段均值變化不一致，說明人類活動能使區(qū)域下墊面條件有顯著改變，從而影響降水下滲，導致產(chǎn)流閾值的增加，但是不同區(qū)域人類活動的影響有所不同。

3.4 洪水頻率演變

根據(jù)紫荊關(guān)站、阜平站1960—2018年的年最大洪峰流量序列，利用備選的6種分布函數(shù)類型，并確定位置參數(shù)(μ)、尺度參數(shù)(σ)為常數(shù)或參數(shù)隨時間t呈現(xiàn)線性變化，分別簡記為s、l，進而對洪峰流量序列進行擬合，選擇AIC值最小者為最佳模型，具體AIC值見表2。

表2 水文站點不同分布函數(shù)和變化趨勢模型的AIC值

由表可知，采用LOGNO對2個站點洪峰流量序列的模擬效果均較好，LOGNO分布中，μ、σ常數(shù)、線性模型均為最優(yōu)。站點最優(yōu)模型的概率分布類型、位置參數(shù)、尺度參數(shù)及對應的GD、AIC、SBC 3個指標的值及站點最優(yōu)模型殘差序列的統(tǒng)計特征(均值、方差、偏態(tài)系數(shù)、峰態(tài)系數(shù)和Filliben系數(shù))均列于表3。與最優(yōu)模型對應的參數(shù)形式為：紫荊關(guān)站位置參數(shù)和尺度參數(shù)均為線性函數(shù)，阜平站的位置參數(shù)為線性函數(shù)，尺度參數(shù)為常數(shù)。圖10展示了各站最優(yōu)模型的殘差蠕蟲圖，中間的紅線是由圖中散點系列擬合的多項式曲線，圖中所有的散點基本位于上下2條曲線之間的置信區(qū)間內(nèi)，可認為各站點最優(yōu)模型的殘差序列均服從標準正態(tài)分布，從而判斷前面所構(gòu)建的優(yōu)選模型的分布類型和參數(shù)選擇是合理的。

表3 站點最優(yōu)模型分布類型指標及殘差特征

根據(jù)各優(yōu)選模型的分布和參數(shù)，可以計算站點年最大洪峰流量序列在2000—2018年指定百分位數(shù)下對應的設(shè)計洪峰流量平均值，結(jié)果列于表4中。對比發(fā)現(xiàn)考慮水文非一致性情形下的設(shè)計洪峰流量值均小于穩(wěn)態(tài)LOGNO分布的設(shè)計值，按非一致性序列分析方法計算的設(shè)計洪峰相對于一致性序列分析方法而言，都有明顯的下降，也說明了洪水序列的非一致性對設(shè)計洪水頻率分析結(jié)果的影響顯著。

表4 指定分位數(shù)對應的設(shè)計洪峰流量值 單位：m3/s

5 結(jié)論

本文分析了大清河山區(qū)典型流域的暴雨洪水基本特征，重點針對年最大洪峰對應的洪水過程和降水量資料，研究洪水水文要素、徑流系數(shù)、產(chǎn)流閾值及洪水頻率的變化情況。主要研究結(jié)論如下。

a)2個站點的洪峰流量下降趨勢顯著，同1960—1979年比較，紫荊關(guān)以上流域、阜平以上流域1980—1999年平均洪峰值分別下降71.1%、45.2%，階段平均洪量值分別下降68.7%、48.7%。同1980—1999年相比，紫荊關(guān)以上流域2000—2018年的階段平均洪峰流量上升9.5%，階段平均洪量下降23.8%；阜平以上流域階段平均洪峰、洪量分別下降35.1%、58.7%。

b)2個代表區(qū)域的徑流系數(shù)階段均值均下降，相比于1960—1979年，紫荊關(guān)以上流域1980—2000年徑流系數(shù)下降幅度為55.9%，2000—2018年與1980—1999年相比，下降34.4%；阜平以上流域徑流系數(shù)下降幅度分別為22.8%、58.0%；隨著降水量的增加，徑流系數(shù)呈上升趨勢，徑流系數(shù)與降水量呈正相關(guān)。白洋淀流域產(chǎn)流閾值整體呈現(xiàn)上升趨勢，2個流域2000—2018年的均值分別為14、16 mm，比1960—1979年增加5～6 mm。

c)對比GAMLSS模型中6種不同分布類型的模擬結(jié)果，帶時變參數(shù)的LOGNO分布更適用于紫荊關(guān)站和阜平站，且最優(yōu)模型對應的殘差序列均服從標準正態(tài)分布，說明優(yōu)選模型的分布類型和參數(shù)選擇是合理的。考慮水文非一致性情形下的設(shè)計洪峰流量值均小于穩(wěn)態(tài)LOGNO分布的設(shè)計值，按非一致性序列分析方法計算的設(shè)計洪峰相對于一致性序列分析方法而言，都有明顯的下降。