变化环境下大清河山区典型流域洪水演变规律研究

王思奇，朱宏鹏，鲁 帆*，江 明，周毓彦

(1.中国水利水电科学研究院，北京 100038；2.河北省保定市水政监察支队，河北 保定 071000)

随着人口增加和社会经济的快速发展，人类活动对流域下垫面及产汇流条件变化的影响不断加剧[1-2]。已有学者针对北方流域暴雨洪水演变特征问题开展了系列研究，梁艳琴[3]、陈旭等[4]的研究表明海河流域不同时段洪水总量、洪峰流量和单位线峰值普遍呈显著的减少趋势；常远远等[5]对黄河湫水河流域2015年稀遇暴雨洪水与20世纪60—80年代的13场洪水进行对比，得出相同降雨条件下的洪峰和次洪洪量较历史洪水均减小；雷成茂等[6]将2016年西柳沟流域的洪水过程与洪峰流量居历史前3位的大洪水相关特征平均值进行对比，发现其洪峰流量、径流系数均偏小；金双彦等[7]分析皇甫川流域次洪特征值变化特点，结果表明：近些年洪水发生次数和年最大洪峰流量高值出现次数均有减少趋势。人类活动对洪水演变的影响主要来源于土地开发利用、水土保持、植树造林、水利工程、城镇化建设等，这些措施不仅使土地利用等环境要素发生变化，而且能使流域产洪次数、地表径流模数和径流系数减小[8-11]。对于变化环境下的频率变化，已有学者采用时变参数模型分析了水文要素的非一致性，结果表明变化环境下应当充分考虑径流的非一致性特点，时变参数模型能够反映水文序列的非一致性特征[12-13]。针对海河流域大清河山区洪水特征变化的研究[8-10,14]表明：历年最大洪峰、次洪水量呈减少趋势，相同降雨量产生的次洪量减少，下垫面对流域径流系数的影响较大。为更深入地揭示该流域暴雨洪水的演变规律，本文利用秩次相关检验法研究了大清河山区典型站点洪水过程中洪峰流量值的整体趋势，对比流域特大洪水的特征变化，在此基础上，采用多时段对比法分析次洪洪峰流量、次洪水量、径流系数、产流阈值的变化情况，利用时变参数模型GAMLSS对洪水频率的演变特征进行分析，研究结果可为流域防洪规划与管理提供支撑。

1 研究区域及数据

大清河水系位于海河流域中部，东西长约275 km，南北宽约200 km，流域面积约4.3万km2，主要由南、北两支组成，流入白洋淀的支流为南支，流入东淀的支流为北支，北支主要为拒马河。其洪水主要由汛期暴雨造成，暴雨中心常出现在阜平、司仓及紫荆关一带，自有水文记载以来，1963年8月暴雨为最大，实测6天洪量为64.3亿m3，南支占80%，北支占20%，全水系30 d洪水总量79.98亿m3，其中南支占78%，北支占22%。紫荆关站是北支拒马河上游的主要控制站，流域面积为1 760 km2。阜平站是南支沙河上游的主要控制站，控制面积2 210 km2。该区域属大陆性季风气候，四季分明，春季干旱多风，夏季炎热多雨，土地利用类型主要以耕地、林地为主，是暴雨洪水多发地区之一，洪水暴涨暴落。

水文气象数据主要来源于《水文年鉴海河流域水文资料》第4册大清河水系以及中国气象数据网，选定1960—2018年的洪水水文要素摘录表及降水量表进行统计，进而得到次洪过程及次洪降水量。其中紫荆关水文站控制流域对应的6个雨量站为艾河村、团圆村、东团堡、石门、王安镇、平顶山，阜平水文站控制流域对应的10个雨量站为银厂、庄旺、站上、神堂堡、不老台、西弯、砂窳、龙泉关、桥南沟、阜平，大清河山区水系及上述水文气象站点的位置见图1。

2 研究方法

首先选择年最大值取样法选择年最大洪峰流量作为研究对象、利用斜线分割法计算其次洪的洪水量、利用算术平均法计算选择的洪水过程相应的面雨量，然后通过Kendall秩次相关检验法对年最大洪峰流量序列的趋势进行判别，通过多时段对比法得到次洪降水量、次洪水量、径流系数及产流阈值的变化情况，采用线性趋势法分析次洪过程中径流系数与降水量的相关关系，采用GAMLSS模型对变化环境下的洪水频率演变规律进行分析，计算流程见图2。

利用斜线分割法计算年最大洪峰对应场次洪水的洪量时，将洪水起涨点和地面径流终止点连成一条直线，直线以下为地下径流，直线以上的地表径流作为次洪洪水总量[15]，计算见式(1)。选择Kendall秩次相关检验法对研究的水文序列构建统计量U，见式(2)。对比计算的统计量值与一定显著性水平下的临界值，可判断水文序列的变化趋势[16]。径流系数指选定时段内的径流深度R与同时段内降水深度P的比值[17]，它代表降水量中有多少水转化成径流[18-19]，径流系数越高，表明地区降雨产流能力越强。产流是指降雨量扣除损失形成净雨的过程，其中降雨损失包括植物截留、下渗、填洼与蒸发，且以下渗为主，典型的产流机制包括蓄满产流和超渗产流，某些地区水文模型的构建兼具2种产流机制[20]。降雨产流阈值是指接受降雨的下垫面能够产生地表径流的最小降雨量[21-22]。因本文采用最大值取样法选取年最大洪峰对应的场次洪水，因此确定所选择洪水过程线中的流量起涨点tQ，根据降水量摘录表中tQ时刻之前各雨量站的时段降雨，求均值得到的区域降水量P即为产流阈值。

(1)

(2)

式中W——次洪洪水总量；S——地表径流与地下径流之和；QQ、QZ——起涨点tQ、终止点tZ所对应的流量值，m3/s；M——径流序列(Qt)中(Qi,Qj,j>i)中所有对偶值(Qi,Qj,j>i)中Qi>Qj出现的个数；N——径流序列的总长度。

在一定显著性水平α下的双尾检验临界值记为Uα/2，若|U|>Uα/2，则说明径流序列的变化趋势显著，否则变化不显著，且U>0时序列呈上升趋势，反之呈下降趋势。

(3)

(4)

本文选用Gumbel(GU)、Weibull(WEI)、Gamma(GA)、Logistic(LO)、Log Normal(LOGNO)、Normal(NO) 6种两参数连续分布作为备用分布函数，并考虑位置参数(μ)、尺度参数(σ)为常数和其随协变量(时间t)呈现线性变化的情况，合计共24种组合方式。计算不同拟合方案对应的AIC值，选定AIC值最小的分布函数及参数变化特征作为站点洪峰流量序列的最优拟合模型，对最优拟合模型进行模型残差评价验证拟合效果，根据确定的最优模型及其参数，计算站点指定时间及指定百分数情况下对应的设计洪峰流量值，并与一致性条件下的设计值对比。

3 洪水演变规律

3.1 洪水过程变化

对紫荆关、阜平水文站1960—2018年实测年最大洪峰流量序列的变化趋势用Kendall秩次相关检验法分析。根据式(2)计算得U分别为-5.22、-3.87，2个值均大于显著性水平α=0.1下的双尾检验临界值1.64，且均大于显著性水平α=0.05下的双尾检验临界值1.96，说明在显著性水平α=0.1及α=0.05下径流序列的变化趋势显著，且呈下降趋势。

考虑资料序列长度以及全国水资源评价工作的时间节点，将资料系列划分为3个阶段(1960—1979、1980—1999、2000—2018年)，对3个阶段的洪水进行统计，见表1。紫荆关以上流域洪峰的时期均值先下降后上升，洪量均值一直下降。同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个流域均出现2次，紫荆关的特大洪水发生在1963、2012年，阜平发生在1963、2016年。相较于1963年的特大洪水，紫荆关2012年、阜平2016年的特大洪水均属于单峰洪水，具有短历时强降雨的特点，且峰值、面雨量均减少，见图3。紫荆关站在1963、2012年的洪峰值分别为4 490、2 156.5 m3/s，阜平站1963、2016年的洪峰值分别为3 380、2 020 m3/s。紫荆关站1963年特大洪水峰现时间历时79 h，面雨量达到380 mm，次洪水量折合径流深为213 mm。紫荆关2012年特大洪水峰现时间历时9 h，面雨量为166 mm，次洪水量折合径流深为39.6 mm，与1963年特大洪水相比，峰值减少约52%，面雨量均值下降214.5 mm。阜平站1963年特大洪水峰现时间历时109 h，面雨量达到541 mm，次洪水量折合径流深为205.6 mm。阜平站2016年特大洪水峰现时间历时4.6 h，面雨量为76 mm，洪量折合径流深为16.5 mm。与1963年特大洪水相比，峰值减少40%，面雨量下降465 mm。

3.2 径流系数变化

根据年最大洪峰选取的场次洪水整理的次洪降水量以及次洪洪量，求得2个流域1960—2018年逐年的径流系数值，并计算3个阶段径流系数的均值，见图4、5。2个代表区域的径流系数阶段均值均下降，点据下移趋势明显。紫荆关以上流域1960—1979年径流系数均值为0.126，相比于1960—1979年，1980—1999年径流系数下降显著，减少量近0.07，下降幅度达到55.9%；同1980—1999年相比，2000—2018年径流系数减少0.019，下降幅度为34.4%。阜平以上流域1960—1979年径流系数均值为0.138，相比于1960—1979年，1980—1999年径流系数减少近0.03，减少比例为22.8%；同1980—1999年相比，2000—2018年径流系数减少约0.05，减少比例达到45.6%，下降显著。

郑江坤等[26]对于蒲江县朝阳水库的4个坡耕地径流小区表明，各小区径流系数均随降雨等级增加呈增加趋势，因此对于径流系数的变化，除对比阶段均值的变化情况外，对区域次洪径流系数与降水量的相关关系进一步分析。以降水量作为横坐标，点绘2个区域整个序列以及3个阶段的径流系数值，分析径流系数的变化趋势，得到径流系数与降水量、时期的相关关系，见图6、7。根据径流系数随降水量的变化图可以看出，线性趋势线的斜率均大于0，即次洪径流系数与降水量呈正相关，降水初期发生蓄渗，不会立即形成径流，随着降水量增加，满足截留、填洼等条件后，进入产流过程，但是当降雨等级较小时，转化为地表径流部分较小，而大暴雨和特大暴雨在一定程度上削弱了土壤前期含水量的影响，因此次洪径流系数随降水量的增加而增加。紫荆关以上流域的年序列变化线性趋势线的斜率为0.001 2，大于0，说明随着降水量的增加，径流系数呈上升趋势；对不同阶段进行分析，3个阶段的线性趋势线的斜率也均大于0，分别为0.001 4、0.000 9和0.000 2，说明径流系数随降水量的增加而增加。阜平以上流域的年序列变化线性趋势线的斜率为0.000 9，与紫荆关以上流域类似，随着降水量的增加，径流系数呈上升趋势；对不同阶段进行分析，3个阶段的线性趋势线也均大于0，分别为0.002 4、0.002 3和0.001 3，进一步验证相同降水量的条件下，径流系数值下降，与紫荆关以上流域类似，在降水量相同的情况下，径流深减少幅度大。

3.3 产流阈值变化

根据年最大洪峰值对应的次洪场次，对研究区域产流阈值的年序列变化情况进行统计，见图8、9。根据图8，紫荆关以上流域产流阈值整体呈上升趋势；3个时期的均值分别为8.4、12.7、13.7 mm，同1960—1979年相比，1980—1999年的产流阈值增加4.3 mm，增加幅度达到50.6%；2000—2018年与1980—1999年相比，增加约1 mm，增幅仅为8.1%。阜平以上流域的产流阈值变化情况与紫荆关以上流域类似，阶段均值均在增加，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个区域的产流阈值变化情况进行比较，产流阈值均呈上升趋势，但是不同区域阶段均值变化不一致，说明人类活动能使区域下垫面条件有显著改变，从而影响降水下渗，导致产流阈值的增加，但是不同区域人类活动的影响有所不同。

3.4 洪水频率演变

根据紫荆关站、阜平站1960—2018年的年最大洪峰流量序列，利用备选的6种分布函数类型，并确定位置参数(μ)、尺度参数(σ)为常数或参数随时间t呈现线性变化，分别简记为s、l，进而对洪峰流量序列进行拟合，选择AIC值最小者为最佳模型，具体AIC值见表2。

表2 水文站点不同分布函数和变化趋势模型的AIC值

由表可知，采用LOGNO对2个站点洪峰流量序列的模拟效果均较好，LOGNO分布中，μ、σ常数、线性模型均为最优。站点最优模型的概率分布类型、位置参数、尺度参数及对应的GD、AIC、SBC 3个指标的值及站点最优模型残差序列的统计特征(均值、方差、偏态系数、峰态系数和Filliben系数)均列于表3。与最优模型对应的参数形式为：紫荆关站位置参数和尺度参数均为线性函数，阜平站的位置参数为线性函数，尺度参数为常数。图10展示了各站最优模型的残差蠕虫图，中间的红线是由图中散点系列拟合的多项式曲线，图中所有的散点基本位于上下2条曲线之间的置信区间内，可认为各站点最优模型的残差序列均服从标准正态分布，从而判断前面所构建的优选模型的分布类型和参数选择是合理的。

表3 站点最优模型分布类型指标及残差特征

根据各优选模型的分布和参数，可以计算站点年最大洪峰流量序列在2000—2018年指定百分位数下对应的设计洪峰流量平均值，结果列于表4中。对比发现考虑水文非一致性情形下的设计洪峰流量值均小于稳态LOGNO分布的设计值，按非一致性序列分析方法计算的设计洪峰相对于一致性序列分析方法而言，都有明显的下降，也说明了洪水序列的非一致性对设计洪水频率分析结果的影响显著。

表4 指定分位数对应的设计洪峰流量值 单位：m3/s

5 结论

本文分析了大清河山区典型流域的暴雨洪水基本特征，重点针对年最大洪峰对应的洪水过程和降水量资料，研究洪水水文要素、径流系数、产流阈值及洪水频率的变化情况。主要研究结论如下。

a)2个站点的洪峰流量下降趋势显著，同1960—1979年比较，紫荆关以上流域、阜平以上流域1980—1999年平均洪峰值分别下降71.1%、45.2%，阶段平均洪量值分别下降68.7%、48.7%。同1980—1999年相比，紫荆关以上流域2000—2018年的阶段平均洪峰流量上升9.5%，阶段平均洪量下降23.8%；阜平以上流域阶段平均洪峰、洪量分别下降35.1%、58.7%。

b)2个代表区域的径流系数阶段均值均下降，相比于1960—1979年，紫荆关以上流域1980—2000年径流系数下降幅度为55.9%，2000—2018年与1980—1999年相比，下降34.4%；阜平以上流域径流系数下降幅度分别为22.8%、58.0%；随着降水量的增加，径流系数呈上升趋势，径流系数与降水量呈正相关。白洋淀流域产流阈值整体呈现上升趋势，2个流域2000—2018年的均值分别为14、16 mm，比1960—1979年增加5～6 mm。

c)对比GAMLSS模型中6种不同分布类型的模拟结果，带时变参数的LOGNO分布更适用于紫荆关站和阜平站，且最优模型对应的残差序列均服从标准正态分布，说明优选模型的分布类型和参数选择是合理的。考虑水文非一致性情形下的设计洪峰流量值均小于稳态LOGNO分布的设计值，按非一致性序列分析方法计算的设计洪峰相对于一致性序列分析方法而言，都有明显的下降。