郑金德，粟 忠，范玉洁，许景璇，李勇军，苏 钊，吴树诚，方荣杰

(1.桂林理工大学，广西环境污染控制理论与技术重点实验室，广西 桂林 541004；2.桂林理工大学，岩溶地区水污染控制与用水安全保障协同创新中心，广西 桂林 541004；3.桂林水文中心，广西 桂林 541001)

1 引言

桂林市以其独特的人文和山水自然风光而闻名于世，已被评为中外著名的国家级风景和旅游城市。但是桂林市每年汛期的暴雨洪水几乎都会给这座秀美的城市带来不同程度的洪涝灾害，而且在枯水期漓江的水量又严重不足，给桂林市的城乡饮用水、漓江旅游用水等造成一定的困扰。随着旅游产业以及经济的发展，桂林市对防汛的要求也越来越高。桂林市急需准确的雨情、水情信息，为桂林市的防汛工作提供相应的参考依据。

目前，关于漓江流域河段流量演变规律的研究相对较少。马斯京根法主要应用于河段流量演算中，由McCarthy于1938年第一次提出。该方法从提出至今，一直被国内外广大从事水文工作的人员广泛使用。马斯京根模型的主要参数为K、x，对演算参数求解的方法大多采用人工试算法。在运用马斯京根模型进行流量演算的过程中，发现演算出来的结果精度不是很高，为此许多科研工作者对马斯京根法进行了改进。在对K、x求解的方面，JyS.Wu等提出了改进图解法， O，Donnell提出了一种直接三参数马斯京根法，运用矩阵计算方法直接推求K、x的值，且矩阵法较图解法推求出来的结果更合理[1,2]。最小二乘法[3,4]是基于数学层面的计算方法，其概念清晰，在水文统计过程中被普遍应用。例如，董晓华等[5]采用最小二乘法，对马斯京根法河道流量演算精度的影响做了探究，并提出了采用最小二乘法有助于提高河道流量计算的精度；卫小蓉[6]、薛梦楠[7]等也基于最小二乘法对不同河道的流量进行模拟，二者均取得了较好的演算结果。

本文针对漓江流域上游三街水文站至桂林水文站段干支流，选取三街水文站、灵川水文站和桂林水文站的逐日实测流量数据，在有约束条件下推导出流量演进系数的最优估计值，客观地建立适用于漓江流域上游的马斯京根流量演算模型，并探究该河段的流量演进规律，进而为漓江流域防洪工作地进行及水资源的合理调度提供相应的理论依据。

2 材料与方法

2.1 流域概况

漓江流域上游桂林水文站以上的集水区总面积为2762 km2，平均海拔约150 m，河段主要干流长105 km，流经桂林市区的河段长4.93 km，整个流域的地势自西北向东南倾斜[8]。且桂林水文站以上流域内雨水充足，水资源丰富；流域多呈羽状，支流短，河道类型为峡谷式河槽，河床的坡降大，漓江流域上游存在多个暴雨区。桂林市气候属亚热带季风气候，多年平均降雨量为1889.4 mm，汛期大多集中在3～8 月份[9]，汛期降雨量约占年降水量的75%，在汛期内漓江流域降水较多，频繁发生水涝灾害。桂林市区三街水文站至桂林水文站主要河流基本特征见表1，河流及水文站网分布见图1。

表1 漓江桂林水文站以上主要河流基本特征

2.2 数据来源

本次模型探究所使用的资料由桂林水文中心所提供，主要研究河段为漓江上游的三街水文站至桂林水文站，该河段内甘棠江为漓江的一级支流，甘棠江上布设有灵川水文站，其它支流并未布设有相应的水文测站。三街水文站处于漓江上游河段的干流上，该水文站以上的流域内建有斧子口水库、川江水库和小溶江水库，各水库所处位置见图1，且流域以上支流众多，地形地势复杂。进行流量演进的前提是演算河道上下游的流量资料要相互对应，三街水文站为新建的水文站，所具有的流量监测数据年份较少，并且由前文可知桂林市的汛期主要集中在每年的3～8月份。因此，本次研究所使用的数据为三街水文站、灵川水文站和桂林水文站2019、2020年3～8月份的逐日实测流量数据。以2019年3～8月份的逐日实测流量数据进行模型参数的率定，以2020年3～8月份的逐日实测流量数据对模型进行了验证。

图1 桂林市区河流及水文站网分布

2.3 马斯京根流量演算模型

马斯京根模型主要是一种在河道流量(洪水)演算中被广泛应用的模型，它利用水量平衡方程式代替连续方程，以槽蓄方程代替复杂的水动力方程，极大地简化了计算过程[10]。本文依据漓江流域上游河段干支流洪水汇集的特性，采用先合后演[11]的流量演算方法，将干、支流断面上三街水文站与灵川水文站的实测流量相互叠加在一起，以此为基础作为漓江流域上游河段的总入流量，按照无支流河段的演算方法进行连续马斯京根演算。

2.3.1 马斯京根方程

在第1、2时段内对水量平衡方程和槽蓄方程[12,13]进行差分求解，可推导出流量演算方程式：

Q2=C0I2+C1I1+C2Q1

(1)

式(1)中：

(2)

(3)

(4)

C0+C1+C2=1.0

(5)

因此，在使用马斯京根方程时，只要确定了参数K、x与△t的值，就可根据式(2)、(3)、(4)分别推导得出C0、C1、C2的值，然后再根据流量演算方程式(1)演算出下断面的最终流量过程Q(t)。

显然，确定模型演算参数x、K的值，成为使用马斯京根流量演算模型的关键。模型参数x、K值的确定，一般根据马斯京根法W～Q`呈单一直线的基本假设对x、K的值进行确定，传统的确定x、K值的方法有试错法、最小二乘法、最小面积法等。利用传统的方法推求x、K值是使W～Q`呈单一直线的最优估计值，但是根据所推导出的x、K值计算的流量演算系数未必能使演算出流量与实测出流量的拟合误差最小[14]。而基于最小二乘法，直接对流量演算系数C0、C1、C2的最优估值进行推导，最后再根据相应公式进行x、K值的反算，再利用所推求出的系数进行河道流量演算，则可取得较好的演算结果。

2.3.2 模型参数优选依据

模型参数的优选以演算出流量与实测出流量的离差平方和最小为依据，构造目标函数[15]表达式(6)与约束方程表达式(7)：

(6)

(7)

式(6)中S为离差平方和，n为实测流量点数。

由式(2)、(3)、(4)可知C0、C1、C2必须满足式(5)的要求。因此可以根据拉格朗日乘数法构造一个非线性规划问题，由式(6)、式(7)构成等式约束方程，再构造拉格朗日函数，将等式所需要约束的极值问题转化为一个无线性约束的极值问题，其极值函数为式(8)：

L(λ,C0,C1,C2)=

(8)

式(8)中，λ为拉格朗日乘数。式(8)达到最小的C0、C1、C2值即为式(7)的解，根据极值存在的必要条件，则有式(9)、(10)：

(9)

(10)

经处理化简后得式(11)、(12)、(13)：

(11)

(12)

(13)

联立式(5)与式(11)、(12)、(13)即可求出流量演算系数C0、C1、C2的最优估值。然后根据式(14)、(15)即可反求出x、K的值。

(14)

(15)

2.4 模型的验证方法

本文主要以确定性系数[16]和相关系数[17]作为模型预测精度的评价标准，确定性系数DC按式(16)计算，相关系数R按式(17)计算：

(16)

(17)

3 结果与分析

3.1 马斯京根流量演算系数的推求

选取三街水文站、灵川水文站至桂林水文站2019年3～8月份的逐日实测流量资料，代入式(8)对马斯京根模型的流量演算系数C0、C1、C2进行推求，然后根据式(14)与式(15)推求出K、x的值；在对传统的马斯京根模型的流量演算系数进行推求，两者所推求出的流量演算系数值见表2。把推求得到的C0、C1、C2值代入式(1)进行模型率定期的流量演算，演算出来的流量值与实测流量值成果见图2。

表2 河段流量演算系数值

图2 演算流量与实测流量对比(率定期)

由图2可知，利用上文所介绍的方法与传统的方法对马斯京根法的流量演算系数C0、C1、C2进行最优取值，然后在进行河道流量演算，得出实测流量值与演算流量值具有大致相同的变化趋势，演算出的流量值与水文站实测的流量值拟合较好。

3.2 模型检验

把由2019年3～8月份的逐日实测流量数据推求出的河段流量演算系数值，代入到三街水文站至桂林水文2020年3～8月份的逐日流量数据进行模拟，模拟结果见图3。从图2、图3可以看出，无论是在模型的率定期或者验证期，所演算得出的流量值与实际监测得到的流量值之间是具有较大偏差的。在流量较小的阶段，演算值与实测值之间可以达到很好的拟合效果，而在流量较大的阶段，演算值与实测值之间尽管存在着大体一致的变化趋势，但是仍然存在着一定的偏差。产生误差的主要原因是，在进行流量演算的过程中并未加入桃花江与南溪河的流量进行运算，而由表1 可知桃花江的多年平均流量为12.2 m3/s，如若加上桃花江与南溪河的流量进行演算，则可取得更好的结果。

由图3和表3可以看出，把由2019 年的流量数据确定出来的马斯京根演算系数C0、C1、C2值，代入到三街水文站2020 年的流量监测数据进行演算，演算出的结果与桂林水文站实测的流量结果拟合程度较好，二者都达到了乙级标准，可为洪水预报工作提供相应的参考依据。

由表3可以明显看出，基于改进的马斯京根法所推求出来的流量值，不论是在率定期还是验证期，其确定性系数与相关系数都比由马斯京根法推求出来的高，但是对比二者的平均相对误差值，改进的马斯京根法却比马斯京根法的高。从图2、图3可以看出不论是在模型的率定期还是验证期，改进的马斯京根法与实测流量值之间具有更好的拟合效果，曲线的变化趋势更加一致。改进的马斯京根法是基于最小二乘法，直接推导出马斯京根模型的河道流量演算系数，其计算过程简单，取得的拟合效果更好。而由马斯京根法推求的流量演算系数，计算过程繁琐，演算系数的推求受计算者的影响较大，所取得的演算结果往往与实测的结果拟合效果较差。因此，改进的马斯京根法比马斯京根法更适合应用于漓江上游的流量演算。

图3 演算流量与实测流量对比(验证期)

表3 模型精度评定结果

4 结论

(1)基于改进的马斯京根流量演算法，推求出了适合漓江上游流域的流量演算系数C0、C1、C2(0.37；0.62；0.01)的值。改进的马斯京根流量演算法较传统的马斯京根法取得更好的拟合效果，在模型验证期内，改进的马斯京根法所推求出的结果确定性系数达到了0.812，相关系数达到了0.940，平均相对误差为-0.195，达到了乙级精度。可为桂林市的洪水预报作业提供相应的参考依据。

(2)改进的马斯京根法较传统的马斯京根法，减少了繁琐的参数试算过程，提高了参数的选取效率，同时也增强了演算结果的可靠性。

(3)本次模型探究所使用的数据是水文测站每日监测的流量数据，其精度可能会存在一定的问题，且使用的数据资料年份较少，会给参数的确定带来一定的影响，再加上有一些支流上并未布设水文测站，所以在今后的工作中应及时掌握流域的来水情况，及时地对模型参数进行相应的修正。