基于特征河长和扩散波动力方程对Muskingum法的理论探讨

丁富平 姚志宏 徐敬华

摘要：为了探索特征河长、扩散波动力方程和Mmuskingum法槽蓄方程三者间的关系，根据水力学和水文学原理，以扩散波动力方程为基础，通过建立河段槽蓄方程，证明了Muskingum法槽蓄方程是扩散波动力方程的一阶近似以及流量比重系数X无量纲性的重要意义；在对特征河长概念外延的基础上，提出了利用AutoCAD绘图软件确定特征河长的新方法，并以两条河流的相关数据为依据进行了实际操作，结果表明，该方法具有操作简单，直观，适应性强，精度好等优点；同时推导出特征河长与X值之间的关系，并分析了不同X值对应的水位流量关系特点及河道断面特征，证明了特征河长、Mmuskingum法槽蓄方程与扩散波动力方程三者联系的纽带是其皆能反映附加比降对河道流量的作用；最后基于不同的假设，给出了计算X的3种公式，得到了不同河流的X值与其河段长L无关，但随计算空间步长、河底比降的增大而增大，随河道糙率和稳定流流量的增大而减小等重要结论。

Abstract: It is proved that the storage equation of Muskingum method was one-order approximation of the diffusion wave equation by establishing the storage equation of river based on the principle of hydraulics and Hydrology . A new approach according to the extension of concept of characteristic river length was proposed to determine the characteristics of long river by the AtuoCAD , and applied to two rivers cases, the results indicated that the approach is simple, visable，high-accuracy, adaptable and high efficiency. In addation, the discharge cross section and stage discharge curve were analyed within the limits of X ,which was calculated by using the relationships between characteristic river length and practical river length, the correlations between the characteristic river length，storage equation of the Muskingum method and the diffusion wave equation were given.Three X formulas were derived by different hypotheses,the results showed that X increased obviously with increacing of slope and space interval，and decreased with increacing of coefficient of roughness and constant flow rate.

关键词：特征河长；扩散波动力方程；Muskingum槽蓄方程；X；理论探讨

Key words： characteristic river length；diffusion wave equation；storage equation of Muskingum method；parameter X；theoretical studies

中图分类号：TV131.3文献标识码：A 文章编号：

早在17世纪，人们就开始了对水流演算数学方法的研究，并不断地发展着，1871年，Bare′de Saint-Venant建立了非恒定流一维分析的理论基础，由于其求解的复杂性，后来出现了许多简化的水流演算方法[1]，Muskingum法水流演算就是其中之一，它具有算法简单，便于程序化，易于稳定，适宜大尺度等优点[2]，进而得到了很多学者的广泛使用，并对其理论基础[3-6]和参数估计方法[7-8]等问题做了大量的研究，但是目前对特征河长、扩散波动力方程、Muskingum法三者联系方面的研究较少，因此，开展这方面的基础研究有助于加深了解Muskingum法水流演算的机理,提高河道水流演算的精度。鉴于此，本文在继承前人研究成果的基础上，根据水力学和水文学原理，借助特征河长的概念，以扩散波动力方程为基础，找到了特征河长、扩散波动力方程和Muskingum槽蓄方程之间的内在联系，拓展了特征河长和X的意义，提出了确定特征河长的绘图法，得到了计算X的3种公式及相关结论，并用相关数据加以验证，结果表明，由理论推导出的结论与计算数据一致，可见本文的推导是正确的。

1特征河长概念的外延

扩散波动力方程可表示为[4]：

（1）

式中：为稳定流流量；x为河段长度；为水深；为河底比降；

（2）

式中：为附加比降；

由式（1）和式（2）得：

（3）

一般地，即，符合二项展开式，对上式展开得：

（4）

对上式取一阶近似为：

（5）

所以附加比降引起的断面流量改变量：

（6）

水位变化引起的断面流量改变量可表示为：

（7）

对、取绝对值，由式（6）~（7）得：

（8）

在特征河长l内有：

（9）

（10）

将式（9）和式（10）代入式（8）得：

（11）

可见，上式等号两边分子部分的附加比降可消去，剩余的均为流量因素；分母部分为比降因素，所以在特征河长内，当时，有：

（12）

由式（8）~（10）得：

（13）

对于任意河段长L，由式（13）得：

（14）

当时，由式（13）和式（14）知：

（15）

其物理意义是，附加比降引起的流量改变量小于水位变化引起的流量改变量。即当河段长L大于特征河长时，由式（9）可得，L引起的水位改变量大于,由式（7）可得，引起的流量改变量大于，又因,所以。

同理，当时，

（16）

其物理意义是，附加比降引起的流量改变量大于引起的流量改变量。

2绘图法确定特征河长

根据和式（12），本文提出了利用AutoCAD绘图功能确定特征河长的方法，此方法简单、方便，基本不需要数学计算，具体步骤如下：

（1）如果上、下断面的水位、流量资料是数据表格形式，则将流量数据输入到Excel的A列，将水位数据输入到B列，在C1单元格中输入“＝A1&”,”&B1”，即把A1和B1中的数据组成了一个坐标点，选定C1单元格当其右下角有“＋”时，拖动鼠标直到与最后的水位、流量数据平齐，这样单独的水位、流量数据在就在C列中形成了一组坐标点；如果上、下断面的水位、流量资料是经验公式格式，则根据计算需要选定相同间隔的一系列流量值，将其输入到Excel的A列中，在B列中编写该经验公式得到对应的水位数据，如前所述，将水位、流量数据在C列中转换成坐标点。

（2）利用AutoCAD中的“SPLINE”绘图功能绘制曲线，此曲线为三次样条曲线，其优点[9]为：在加密时，它和它的导数能在整个插值区间上充分靠近被插函数，能够保证计算精度和曲线的光滑性[10]。其绘图方法为：①从Excel的C列中复制所有的坐标点。②在AutoCAD的命令行中输入 “SPLINE” 命令后回车，将鼠标定位在命令提示“指定第一个点或 [对象(O)]:”的后面，然后按Ctrl+V，即将复制的坐标点输入到AutoCAD的命令行中，再按回车，然后选择拟合残差为零，即在命令行中输入“F”，回车后在命令行中再输入“0”，然后回车，即得到水位—流量的三次样条曲线。③给绘制出的曲线添加X轴和Y轴，X轴的起点必须从零开始。

（3）自计算流量处作X轴的垂线，与水位—流量曲线交于一点，自此点分别作Y轴的垂线和水位—流量曲线的切线，分别与Y轴交于上下两点，这两点间的距离为，相应的 ；自Y轴上较大的点做边坡系数为的斜线，自较小的点作水平线，两条直线的交点到Y轴的距离即为该流量对应的特征河长的长度（见式（12））。其中在绘制水位—流量曲线的切线时，由于此曲线是样条曲线，同AutoCAD中的其它弧线不同，所以本文采用先将该样条曲线偏移，得到一条与其形状相同的曲线，然后连接两曲线的端点得到一条线段，自原样条曲线端点做该线段的垂线，此垂线即为端点处样条曲线的切线（见图1）；绘制其它流量处的切线，要先将样条曲线在该点处打断形成一个端点，然后再按照上述方法确定其切线。

本文选取文献[3,11]中的水位和流量数据进行实际操作，以永定河雁翅—三家店上游断面的数据为例，在绘制过程中，选取103 m3/s为流量单位，使水位、流量的数值变化在用一个数量级上，以得到比例协调的样条曲线；在绘制的斜线时将斜率扩大103倍（即103），得到的单位为1m：1km的比例，即特征河长的长度单位为km，通过上述比例转换可得到坡度和长度协调的图形（见图1）。

图1绘图法确定特征河长示意图

Fig1 Determine the characteristic river length by the method of drawing

根据上下断面的水位流量数据，用绘图法可得到某流量下对应的两个特征河长值，取平均后作为此流量下的特征河长大小，并与文献[3,11]中的计算结果比较，结果见表1。

表1特征河长计算表

从表（1）可知，相对误差的最大值出现在沅陵—王家河河段的最小流量点，其值为17.5%；相对误差的最小值出现在雁翅-三家店河段的最大流量点，其值为0；除相对误差最大值外，其它各点的相对误差都在10%以内。此外，绘图法操作简单，直观，适应性强，人为干扰少，计算少，需时少，至于出现最大相对误差的原因分析如下：①选取的资料较少，根据样条函数的优点，我们可以用增加计算点个数的方法提高样条曲线的拟合程度，以精确地确定水位—流量曲线的斜率，提高计算精度。②流量最小点是水位—流量曲线的端点，其趋势由较大的水位流量坐标点确定，而无法反应较小坐标点对其变化趋势的影响，所以在计算某点的特征河长时，至少要选取一个比其小的坐标点，以更好的确定水位—流量曲线的切线方向；如：雁翅-三家店河段流量3×103m3/s和7×103m3/s得到的取平均后的相对误差就较小，其值仅为6.37%。

3对的推导

稳定流水深—流量关系可表示为：

（17）

式中，、分别为经验参数；其它符号的意义同前。

根据文献[4]，将式（17）代入式（1），求得并对其进行二项展开得：

（18）

对上式取一阶近似得：

（19）

河段长L的槽蓄量W可表示为：

（20）

式中，为河段的平均河宽； 为河段平均流量；为河段内的平均水深；其它符号的意义同前。

由水力学和水文学知识可得：,，将其代入式（20）得：

（21）

式中，I、O分别为河段L的上、下断面流量，其它符号意义同前。

取，并将代入上式得

（22）

Muskingum法槽蓄方程为：

（23）

式中，为蓄量常数；为流量比重因素。

比较式（22）和式（23）得：

（24）

（25）

对于同一个过水断面的流量过程有：

（26）

由式（17）可知，当时，，即水深—流量关系呈线性关系，在河段槽蓄量的计算中又假设了水深沿程呈线性关系，由上述两种线性假设可知，平均流量等于中断面的稳定流流量，即：

（27）

式中，、分别为河段上下断面的水深；其它符号的意义同前。

将式（26）、式（27）代入式（25）得：

（28）

将式（14）代入式（28）得：

（29）

由式（13）和式（28）可得：

（30）

式中，符号的意义同前。

对上述推导过程分析，可以得到如下认识：首先，在对X推导过程中，以扩散波动力方程为基础，通过对稳定流水位流量关系和水深沿程线性变化的假设，建立了河段槽蓄量计算公式（式（20）），通过理论推导得到了与Muskingum法槽蓄方程相同的形式，由于水深是由扩散波动力方程取一阶近似得到的(式（19）)，因此Muskingum法槽蓄方程实际上是扩散波动力方程的一阶近似；也可从特征河长角度加以证明:式（29）中的是由扩散波动力方程经过二项式展开取一阶近似得到的(式（5）)，将式（29）回代到Muskingum法槽蓄方程，也可以说明Muskingum法槽蓄方程是扩散波动力方程的一阶近似，而且是在线性假设条件下扩散波动力方程的一种简化形式，这就找到了Muskingum法水流演算与扩散波动力方程的关系，以及Muskingum法水流演算具有良好计算精度和普适性的理论基础。

其次，Muskingum法槽蓄方程首先是作为经验公式被提出的[5]，它之所以能够被实践证明是正确的，一个重要的原因是比重系数的成功引入，其引入的过程虽然带有很大的偶然性，但是的无量纲性是客观存在的，即推导得出的式（29），对此式进行量纲分析，可以很容易得到为无量纲数，如果不是无量纲数，就会使式（23）中和的含义更加模糊，就很难表示上、下断面流量对河段槽蓄量的影响，因此，具有无量纲性是其能被成功引入到Muskingum法槽蓄方程的立足点，另外，式（29）也是特征河长、Muskingum法槽蓄方程和扩散波动力方程三者联系的纽带。

最后，根据特征河长概念在任意河长L上的外延，可以很简单地计算X的变化范围以及代表的含义。当时，将式（13）代入式（28）得，，表示附加比降引起的流量改变量与水位变化引起的流量改变量相等；当时，将式（15）代入式（26）得，，表示附加比降引起的流量改变量小于水位变化引起的流量改变量，说明水位变化对流量的影响较大，由此可以推断河道断面应为宽浅型，反映在河道断面水位—流量曲线上，说明其斜率较小，式（17）的导数较大，即较大且（因为式（17）的导数应是水位的增函数）；当时，将式（16）代入式（28）得，，表示附加比降引起的流量改变量大于水位变化引起的流量改变量，说明水位变化对流量改变量的影响较小，由此推断河道断面应为窄深或狭小的“V”型断面，如长江万县到宜昌段，[3]，反映了河道断面水位—流量曲线的斜率较大，式（17）的导数较小，即较小且；此外，由于式（11）等号两边的附加比降可削去，因此等号左边附加比降引起的流量改变量的值是确定的，当计算河长一定时，则附加比降引起的流量改变量仅取决于水位—流量曲线的斜率（见式（14）），由此可知，附加比降对流量的影响可用水位流量关系来表示，恰恰是一个包含水位流量关系的参数（见式（28）），所以X能够反映附加比降对流量的作用，由此可见，特征河长、X和扩散波动力方程中的附加比降项皆能表示附加比降对流量的影响，这就是特征河长、Mmuskingum法槽蓄方程与扩散波动力方程三者间的交叉点和本质联系。

4X计算公式的不同形式

4.1基于式（17）的X计算公式

由式（26）和式（28），X的计算公式为：

（31）

当水深流量关系符合式（17）时，则上式可写为：

（32）

由上式可知，随河底比降的增大而增大，随稳定流流量的增大而减小；由于，所以X随水深的增大而增大。

从上式看，X也应随计算河段长度L的增大而增大，但从表2的第③和④列可以看出，X与L并没有明显规律，现对式（32）做进一步推导：

河底比降可表示为：

（33）

式中，为上断面河底高程；为下断面河底高程；其它符号意义同前。

将上式代入式（32）得：

（34）

由上式可知，X与计算河段长度L无关，而随上、下断面高程差的增大而增大，见表2第⑦列；

表2是以X值从大到小对河流进行排序得到的。从表中可以看出，X值与河段长度L之间并没有明显的规律，与式（34）中得到的X与L无关的结论相同；由于河段比降，水位流量关系等其它水力因素的影响，表2中X值与河流上、下断面高程差之间的规律不完全明显，因此，按照X值的范围将全部河流分为四组，即：0.5~0.4，0.4~0.3，0.3~0.1和0.1~0，然后对每组河流上下断面的高程差求平均值，其值分别为：75.35m，28.68m，12.78m和5.80m，由此可以看出，X值随着河流上下断面高程差的减小而减小，与式（34）得到的结论相符（需要解释的是：0.1~0组中钱塘江的河底比降与其它两条河流的相近，可消除河底比降对X值的影响，且与伊洛河都有支流汇入的影响，所以，尽管钱塘江的X值较大但仍将其放入0.1~0组，这种处理能尽量保证该组河流水利特征的一致性，从计算结果来看，这种处理也是合理的）。

表2 中国部分河流的Muskingum法参数X及河段基本特征[3,4]

4.2基于宽浅矩形河槽的X计算公式

式（34）中虽然没有包括河底比降，但是不能排除河底比降对X的影响，因为经验系数和都受到过水断面处河底比降的影响，因此对X的计算公式再做进一步的探讨。

对于宽浅矩形河槽稳定流流量符合下式：

（35）

式中，为糙率，其它符号意义同前。

将式（33）和式（35）代入式（31）得：

（36）

由上式可知，宽浅矩形河槽的X随河道糙率的增大而减小，随河底比降的增大而增大；由表2第⑤列可知，河底比降与X值之间的规律性较好，与上式的结论相符。

基于空间步长的X计算公式

在洪水波流量演算时，需将连续性方程和动力方程化为差分方程来求解，在河长L内确定空间步长后，需要计算的X值，当，由式（32）可得：

（37）

由上式和式（33）可知， 不能如式（32）中一样被削去，所以X值随着空间步长的增大而增大。

5结论

从特征河长的概念入手，分析了附加比降和水位变化对流量的影响，并将特征河长的概念外延到任意河长L（式（15）和式（16））；以扩散波动力方程为基础，通过建立河段槽蓄方程，证明了muskingum法槽蓄方程是扩散波动力方程的一阶近似，由式（29）证明了X的无量纲性是Muskingum法槽蓄方程成立的根本立足点；揭示了特征河长的概念与X值的关系（式（28）~（30）），确定了X的变化范围，分析了不同X值对应的河道的断面特点和水位流量关系特征，能够弥补 “楔蓄”概念对X出现负值等现象无法解释的缺点[11]；最后指出水位流量关系、X值和附加比降对流量改变量的影响是等价的，即找到了特征河长、Mmuskingum法槽蓄方程与扩散波动力方程三者的本质联系。

通过对特征河长的重新认识，结合式（11）、式（12）提出了利用AutoCAD软件确定特征河长的新方法，此方法有操作简单，直观，适应性强，精度好，人为干扰少，计算少，需时少等优点。

基于不同的假设，给出了X的3种新的计算公式，得到了X值与计算河长L无关，随附加比降、计算空间步长的增大而增大，随河道糙率和计算流量的增大而减小等结论，探讨了X值与水深的关系。

参考文献:

[1] 〔美〕David R. Maidment，张建云，李纪生译.水文学手册[M].北京：科学出版社，2002，387-388.

[2]赵人俊，王佩兰.二维洪水演算方法[J].河海大学学报（自然科学版），1992， （03）：1-6.

[3]芮孝芳.水文学原理[M].北京：中国水利水电出版社，2004，191-195.

[4]芮孝芳.Muskingum法及其分段演算的若干理论探讨[J].水科学进展，2002，13

[5]芮孝芳.水文学的研究方法与理论创新[J].水利水电科技进展，2003，25（2）：46-50.

[6]袁晓辉,张双全,张勇传等.非线性马斯京根模型参数率定的新方法[J].水利学报,2001,(5):77-81.

[7]鲁帆,蒋云钟,王浩等.多智能体遗传算法用于马斯京根模型参数估计[J].水利学报，2007，38（3）：289-294.

[8]邓建中，刘之行.计算方法（第二版）[M].西安：西安交通大学出版社，2001，82-83.

[9] 黄振平,顾巍巍,葛慧等.基于三次样条函数的马斯京根参数估计方法[J].河海大学学报(自然科学版),2008，36（1）：11-13.

[10] 包卫民.水文预报（第三版）[M].北京：中国水利水电出版社，2006，98-99，322.

注：文章内所有公式及图表请以PDF形式查看。