陈学秋，瞿思敏,2，李代华，石 朋,2，王轶凡，单 帅，勾建峰

(1.河海大学水文水资源学院，南京 210098；2.河海大学水文水资源与水利工程国家重点实验室，南京 210098；3.云南省水文水资源局文山分局，云南 文山 663099)

富春江上游兰溪入库站水位受富春江水库水位顶托、开机开闸回水波、动库容和上游橡胶坝等人类活动影响，水位变化剧烈，传统的水位流量关系线无法应用。为应用水位比降法准确估计兰溪站入库流量，为洪水调度提供科学依据，首先需获取兰溪站及其下游女埠站水位资料。在调查过程中发现由于实际条件局限，无法实时获得女埠站的水位数据，因此本文采用水动力学方法和经验相关法分别对女埠站水位进行模拟计算。

以圣维南方程组为理论基础的水动力学方法是河道洪水演算的常用方法之一[1]。圣维南方程组是描述河道一维非恒定流水力要素随时空变化的一阶双曲线型拟线性偏微分方程组[2]，随着计算机技术的快速发展，高精度的数值求解方法被广泛应用于河道洪水模拟，并发展到多维河道水流[3]运动及河网计算[4]，以此为基础建立的水动力学模型具有一定物理基础，对水文资料要求较低，计算精度高，适用于平原地区河道洪水演算。此外，水动力学模型受到初始条件[5]、边界条件以及河道糙率[6]等影响，对水文地形资料的选取有一定要求。经验相关法是水位计算研究中出现最早、应用最广泛的一种数理统计方法，其基于大量实测水文资料的统计分析建立水文要素间的相关关系，根据已知站点水位过程推求未知站点水位过程。

1 研究方法与模型原理

1.1 水动力学模型

圣维南方程组是描述河道一维非恒定水流运动的基本方程，它刻画了各断面水力要素随时间、空间位置的变化规律。圣维南方程组依据质量守恒和动量守恒2大基本定律[7]，由连续方程和动量方程组成：

(2)

式中：x为距离，m；t为时间，s；B为过水断面底宽，m；A为过水断面面积，m2；Q为流量，m3/s；Z为水位，m；qL为旁侧入流，m2/s，入流为正，出流为负；vx为旁侧入流沿水流方向的速度，m/s，若旁侧入流方向垂直于水流方向，则vx=0；α为动量修正系数；g为重力加速度，m/s2；C为谢才系数，m1/2/s；R为水力半径，m。

圣维南方程组是一组一阶双曲线型拟线性偏微分方程，本文拟采用简化四点线性隐格式对圣维南方程组进行求解[8]，差分格式如下：

(3)

式中：f为变量，代表水位、流量等；θ为加权系数，0≤θ≤1。

1.2 经验相关法

经验相关法是水位计算研究中出现时间最早、应用最广泛的方法。该方法对大量实测资料进行统计分析，找出具有相关关系的水文要素作为自变量和因变量，利用数理统计方法对其进行拟合，得到各变量之间的经验公式或相关关系图[9]。经验相关法依据变量多少分为单因素和多因素相关关系。依据自变量和因变量的关系分为线性和非线性相关关系。因此要根据研究河段影响水位变化的因素不同，因地制宜地选择合适的水文要素拟定经验公式，进行水位计算。

多因素线性相关基本模型一般表达为：

Y=b0+b1X1+b2X2+…+bnXn+ε

(4)

式中：X1,X2,…,Xn为自变量；Y为是因变量；b0,b1,…,bn为系数；ε为随机误差，均值为0。

2 模型应用

2.1 流域概况

富春江水库位于浙江省钱塘江中游富春江七里垅峡谷的出口处。钱塘江位于东经118°～121°，北纬28°～31°，流域面积49 930 km2，是浙江省第1大河。富春江流域处于亚热带地区，气候温润多雨，多年平均气温为17.1 ℃。流域多年平均降水量为1 659.2 mm，最大年降水量为2 271.9 mm。流域降水量时空分配不均，年内分配以6月份最多，12月份最少，空间分布西多东少。富春江流域山地占61.5%，丘陵占23.7%，平原占14.8%，地形自西向东倾斜，河流的坡降以上游较大，衢江以上约为0.5%，衢县至梅城为0.25%～0.5%，梅城至坝址约为0.15%，坝址以下坡降平缓，已受钱塘江强潮影响。

近年来富春江流域下垫面因子发生了重大改变，流域内部人类活动影响增加，特别是流域内小水电站较前几年大量增加，使富春江流域洪水预报的难度大大增大。富春江上游有兰溪入库站，控制流域面积的85%，是富春江入库估计有效的参照站。本文研究对象为兰溪及其下游三河、梅城、江南站之间河段。

2.2 水动力学模型应用

2.2.1 模型构建

兰溪断面为富春江水库入库控制断面，取兰江兰溪至三河段为研究河段，河段长21 km，其间女埠站距上游兰溪站6 km，河段示意见图1。

图1 研究河段站点分布Fig.1 Distribution of hydrometric stations

河道首、末断面采用85国家高程基准面作为绝对基面，兰溪站冻结基面与绝对基面高差为20.135 m，三河站冻结基面与绝对基面高差为20.096 m，2站相对高差为0.039 m。兰溪站至三河站之间通过线性内插得到20个断面，各断面之间相差1 km，计算断面示意见图2。内插断面及三河站无实测断面资料，考虑到计算河段较短，且两岸多为钢筋混凝土防洪堤坝，断面形状变化不大，因此所有断面均使用兰溪站平均实测大断面数据进行计算。

图2 研究河段计算Fig.2 Calculation structure of Lanxi-Sanhe river section

基于水位数据的一致性原则，本文选取2013年5场洪水过程，分别是20130323、20130423、20130601、20130626及20131007次洪水，其中前4场洪水用于参数率定，第5场洪水用于验证。

初始条件包括各断面计算开始时刻水位或流量[5]。本文已知上下边界初始时刻水位，对其进行线性内插可得各断面初始水位。各断面水位均取85国家高程基面。

根据已有水文资料，本文采用兰溪站水位过程作为上边界条件，三河站水位过程作为下边界条件，中间不考虑旁侧入流，计算女埠站水位过程。

2.2.2 参数率定及模型计算结果

对圣维南方程组采用4点隐式格式进行差分，需要率定的参数包括河道糙率和权重系数。权重系数θ在[0.5,1]内无条件稳定，率定结果为0.75。考虑到河道糙率的空间变异性，模型对所有计算河段糙率逐一率定，率定结果见表1。

表1 各河段糙率率定结果Tab.1 Calibration results of reach roughness

参数率定及验证过程水位计算结果见表2，计算过程见图3、图4。5场洪水洪峰误差均小于实测变幅的20%，峰现时间均不超过±1 h，确定性系数均大于0.8，说明此次模拟结果满足许可误差，精度较高。率定期洪水过程线出现计算洪峰水位大于实测洪峰水位，计算峰现时间晚于实测峰现时间，退水段洪水过程线拟合不佳的现象。总结为以下2点：①涨水段与退水段模拟不一致，总体看来涨水段模拟结果比较符合实测值，退水段模拟结果偏高；②计算过程线相比实测过程出现坦化，不能反映真实水位的动荡变化，验证期洪水过程线表现尤为明显。

图3 率定期洪水模拟水位过程Fig.3 Simulated water level process of floods in model calibration

图4 验证期洪水模拟水位过程(20131007次洪水)Fig.4 Simulated water level process of floods in model validation(flood 20131007)

类型洪号实测洪峰水位/m计算洪峰水位/m峰值绝对误差/m峰值相对误差/%峰现时间误差/h确定性系数率定2013032323.8623.880.020.0800.942013042323.9924.010.020.08-10.922013060125.1125.160.050.1900.972013062625.3925.400.010.0400.98验证2013100723.3823.400.020.0900.88

注：峰现时间误差一栏中，“-”代表洪峰滞后，“+”代表洪峰提前，以下表格同。

综合分析水动力学模型计算误差来源主要是模型结构的不足及断面资料缺乏：①水动力学模型用差分方程逼近微分方程时引入了误差，误差项使计算相位发生变化，导致洪峰滞后。②模型对河道底坡进行简化，可能与河底真实底坡相差较大，如果内插所得断面底坡高于真实值，则会导致水深起算点比实际偏高，计算过程坡降偏小，重力作用被削减，流速减小，计算水位偏高。③糙率是河道洪水坦化的内在因素[7]，涨退水过程中的水力要素往往不一致，本文模型概化中涨水和退水段采取同一套糙率值进行计算，不能反映其时间变异性，导致退水过程计算误差较大。

2.3 经验模型应用

2.3.1 断面资料

本文通过分析江南水位站、梅城水位站、三河水位站和兰溪水位站的水位资料，建立相关关系，推求女埠站水位资料。各水位站水文资料统计见表3，各水位站的地理位置见图1。

2.3.2 经验模型计算

首先通过分析各水位站水位资料，找出女埠站与其他水位站之间的函数关系。通过20130323次、20130423次、20130601次及20130626次洪水实测资料率定得到女埠站计算水位表达式(回归统计数据见表4，回归参数分析见表5)：

表3 各水位站水文资料统计Tab.3 Hydrological data of water stations

Z(女)=0.697 908Z(兰)+0.033 926Z(三)+

0.277 292Z(梅)+0.028 699Z(江)-1.007 58

(5)

表4 回归数据统计Tab.4 Regression data statistical table

表5 回归参数分析Tab.5 Regression parameters statistical table

注：表中Z(兰)、Z(三)、Z(梅)、Z(江)分别表示兰溪站、三河站、梅城站和江南站水位。

将20131007次洪水各水位站实测资料代入经验公式计算得到女埠站计算水位过程，与实测水位过程进行比较分析(见表6及图5)。本文总结的经验公式在进行水位过程计算时模拟程度较好，确定性系数较高，洪峰误差及峰现时间误差皆在许可误差范围内。由洪水过程线可知，涨水过程计算结果略低于实测过程，退水过程则相反，说明计算过程较实测过程较为滞后。

表6 经验相关法计算结果Tab.6 Calculation results of the empirical correlation method

图5 20131007次洪水计算过程Fig.5 Simulated water level process of flood 20131007

基于数据分析的经验方法对资料的数量要求较大，本文只选取4场洪水率定经验公式，1场洪水进行验证，虽然计算结果总体较好，但是所得经验公式是否具有普适性，能否满足不同年份、不同场次洪水过程对水位计算的要求，还需要进一步验证。

3 结 语

本文选择水动力学方法及经验相关法计算富春江入库控制断面水位，2种方法计算结果精度均满足要求。在洪峰水位计算上水动力学方法优于经验相关法，而在洪水过程的模拟上则是经验相关法的确定性系数更高。总体而言，本文建立的水动力学方法具有一定物理基础，对水文资料要求不高，能够较好地对富春江入库控制断面水位进行模拟计算，反映河道水流运动的内在规律。如果模型概化能够考虑糙率的时间变异性，则计算结果会更加准确，因此值得深入研究。

□

[1] 芮孝芳.水文学原理[M]. 1版.北京:中国水利水电出版社,2004:6-7.

[2] M 乔芬,A 雷麦斯,芮孝芳. 圣维南方程组的解析解和数值解[J]. 河海大学科技情报,1987,(2):155-163.

[3] 王船海,李光炽.流域洪水模拟[J].水利学报,1996,(3):44-50.

[4] 李光炽.流域洪水演进模型及其参数反问题研究[D]. 南京：河海大学,2001.

[5] 张小琴,包为民,王 涛.水动力学模型初始条件问题浅议[J].水电能源科学,2008,26(2):62-65.

[6] 张小琴,包为民.糙率对感潮河段水位预报的影响[J].水文,2009,29(3):19-23.

[7] 赵振兴,何建京.水力学[M].2版.北京:清华大学出版社,2010:108-109.

[8] 胡庆云,王船海.圣维南方程组4点线性隐格式的稳定性分析[J].河海大学学报(自然科学版),2011,39(4):397-401.

[9] 包为民.水文预报[M].4版.北京：中国水利水电出版社，2009:132-134.