《1 前言》

1 前言

南水北调中线工程是缓解我国北方水资源严重短缺、优化水资源配置、改善生态环境的重大战略性基础设施，其渠道断面尺寸大，形状多样，有梯形、矩形等，是宝贵的原型试验平台。 目前，已经对南水北调中线京石段过水建筑物典型断面糙率进行了原型测试，按照常规水力学方法率定的糙率在 0.013 37~0.015 7 之间，大部分偏高，应用价值受到怀疑。目前，中线工程正在建设，需要准确可靠的实测糙率进行设计复核，同时，随着工程的建成，其运行水力控制也需要这一关键水力学参数。 因此，开展中线工程渠道糙率的系统辨识研究，不仅具有重要的理论价值，而且具有十分重要的实用价值。

系统辨识是 20 世纪 60 年代从现代系统理论发展出来的科学分支，国外对河道糙率的辨识研究开展较早。 Becker 等（1972）[1]首先构造了一个关于明渠非恒定流参数辨识的影响系数方法，结合改进的单纯形法，最小化误差平方或最小化最大误差绝对值来反演糙率；Chiu 等（1978）[2]以卡尔曼滤波配合观测的水深推求曼宁糙率； Wormleaton 等（1984）[3]将水深和流量的相对误差函数作为优化函数，利用影响系数方法和非线性最小二乘算法求解最适合的糙率以缩减水深、流量的相对误差；Crissman 等（1993）[4]以圣维南方程为基 础建立了河床糙率随时间变化的预测模式； Wasantha Lal（1995）[5]采用奇异值分解方法反演糙率；Khatibi 等（1997）[6]探讨了准则函数在不同噪声水平中的特性和样本量对计算结果可靠性的影响； Atanov 等（1999）[7]基于伴随方程方法，采用最小二乘原理求解目标泛函，反演计算梯型明渠糙率； Ramesh 等（2000）[8]采用连续二次规划方法反演糙率； Yan Ding 等（2004）[9]基于最优控制理论研究了二维浅水方程中糙率的辨识问题；Wong 等（2003）[10]建立了水深及糙率值之间的关系，并在进行运动波模式演算过程中，由水深变化即时修正糙率；Hsu Minghis等（2006）[11]建立了实际观测水深与计算值的目标函数，利用 Gauss -Newton 方法求解非线性的最小二乘问题。

在国内，金忠青等（1998）[12]根据实测资料，选择水位过程或流量过程作为目标函数的变量，构造各河段误差平方和这样一个目标函数，采用复合形法求解目标函数直接优化；董文军等（2002）[13]根据参数辨识理论对一维水流方程中的糙率进行理论推导，使用最小二乘法建立最优模型的目标函数；李光炽等（2003）[14]提出利用卡尔曼滤波来求解河道的糙率；程伟平等（2005）[15]引入控制论理论，应用带参数的卡尔曼滤波法进行河道糙率反演分析。

上述糙率的辨识研究主要针对河道，是在假设水力学参数，如流量、水深、水位等测量值准确无误差的基础上进行的。 由于原型水位、流量等测量数据包括很多不确定性因素，即随机误差，如测量仪器的误差；测量环境的影响，如安装仪器测点的选择、安装质量、气温和水温等因素，所以在一般情况下，按照这种方法获取的糙率受测量误差和计算误差的影响，不能推广应用到其他渠道系统，特别是人工渠道。

本文的目的是根据水力学原理，考虑渠道断面形状、底坡、渠长变化的影响，建立渠道沿程糙率与粗糙高度 ks 和水力半径 R 的函数关系，然后通过数学变换提出系统辨识的线性模型，最后以南水北调中线工程原型观测资料为基础，应用系统辨识的方法消除水力测量随机误差的干扰，得到通用的渠道沿程糙率经验公式。

《2 数学模型》

2 数学模型

《2.1 渠道沿程糙率的计算》

2.1 渠道沿程糙率的计算

对于恒定非均匀流渠道，渠道沿程糙率计算方程为

式（1） 中， Q 为流量, m3/s； y 为水深； A 为断面面积，m2； R 为水力半径，m； L 为渠道长度，m； n 为渠道糙率： s0 为底坡； g 为重力加速度，m/s2；下标“1”为渠道进口；下标“2” 为渠道出口； ζ 为渠道局部阻力系数；  为渠道过水断面平均值，即

在工程设计中，国内工程常常利用式（3）计算渠道沿程糙率，即

由此可以看出，糙率 n 不仅与反映渠道表面平整度的 ks 值有关，还与水力半径 R 的大小有关，n 将随 R 的加大而加大。ks 的取值在工程初期一般取 0.000 61。

美国垦务局根据已建渠道的实测资料和室内试验分析，推荐采用下述方法确定混凝土渠道的糙率（王光谦等，2008）[16]

从式（3）和（4） 可知，渠道糙率是等效粗糙高度ks 和水力半径 R 的函数，可以写成统一形式

当令系数 c1 =18 和 c1 lgc2=19.55 -18lgks 时，则式（5）与国内现有经验公式（3） 相同。当令系数c1 =1 /0.056 5 =17.7 和 c1 lgc2 =17.7lg9 711 ＝70.6 时，则式（5） 与美国垦务局经验公式（4）中的第二式相同。 美国垦务局与我国经验公式系数的不同，主要反映了模型试验和原型观测的差别。南水北调中线工程干渠断面尺寸在我国人工渠道中是最大的，根据原型实测观测数据重新率定系数 c1 和 c2 的值，具有十分重要的实际意义。

《2.2 渠道沿程糙率系统辨识线性模型》

2.2 渠道沿程糙率系统辨识线性模型

为了便于应用系统辨识的最小二乘法确定系数c1 和 c2 ，式（5）可以改写为

式（6） 中

由式（1）和（6） 可得

整理得

式（8）中， 

在系统参数辨识过程中，线性方程式（8）的系数 和 b 是已知量， 和 是未知量，即需要辨识的参数。

对于 m 条渠道的系统，对渠道 i 式（8）可改写为

其矩阵形式为

式（11）中，

《2.3 系统辨识的最小二乘法》

2.3 系统辨识的最小二乘法

实测参数， 如流量、水深，存在各种噪声 （误差），系数矩阵 A 和向量 B 是依赖于实测参数，所以按照式（11）只能得出辨识参数 W 的估计值 ，即

式（13） 中，，其中上标 Ｔ 为转置符号。 ε 称为拟合误差，或者残差。

为了得到 W 的最优估计 ，最小二乘辨识算法的准则是 的平方和最小，即

为了使 J 极小化，可令 =0 ，即

得

式（16）中，矩阵 AＴA 是对称矩阵，当它为非奇异（正则矩阵）时，有

上述结果是在认为观测数据具有相同可信度的基础上推导出来的，即对每个 残差 给予相同的权，当对每个残差项加 以不同的权，并令 P 为加权矩阵时，则广义最小二乘法的准则函数为

式（18） 中， P 为对称正定矩阵，一般为对角线矩阵。同理，令，可得

式（19） 中，为加权最优估计，矩阵AT PA 也是对称矩阵，当它为非奇异时，有唯一解。

本文采用了全主元消去法直接求解式（16） 和式（19） 。 可以证明，最小二乘估计是无偏估计、一致估计和有效估计（ 徐枋同等，1999）[17]，可以有效消除随机误差的干扰。

《3 渠道沿程糙率的系统辨识》

3 渠道沿程糙率的系统辨识

下面将以南水北调中线京石段应急供水工程实测资料为依据，考虑渠道断面形状、底坡、渠长变化的影响，应用系统辨识的最小二乘法确定渠道沿程糙率与等效粗糙高度 ks 和水力半径 R 的函数关系。

《3.1 实测数据》

3.1 实测数据

表 1 列出了从唐河倒虹吸出口至漕河渡槽进口渠道特征参数一览表，渠道数 为 8，其 中渠道 1 ~5 是同形梯形渠道，它们的断面形状、尺寸、底坡完全相同；渠道 6 和 8 也是同形梯形渠道，它们的断面形状、尺寸、底坡完全相同；渠道 7 是隧洞，过水断面是矩形。 需要说明的是，渠道 5 和渠道 6 不是连续渠道，并且所有渠道中的流动都是非均匀流，水面线是雍水曲线。

《表1》

表1 渠道特征参数

Table 1 Characteristic parameters of channels

表 2 列出各渠段平均水力半径 R 和计算曼宁糙率 n 一览表，其渠道编号与表 1 是一一对应的，根据该表可以画出如图 1 所示 n - R 曲线。

《表2》

表2 渠道参数

Table 2 Parameters of channels

从图 1 可见，受糙率噪声影响，糙率 n 随水力半径 R 变化波动较大，下面将应用系统辨识的最小二乘法减少噪声干扰，确定渠道沿程糙率与等效粗糙高度 ks 和水力半径 R 的函数关系。

《图1》

图1 实测糙率 n 与水力半径 R

Fig.1 Calibrated channel roughness n，versus hydraulic radius R 

《3.2 渠道沿程糙率参数的最优估计》

3.2 渠道沿程糙率参数的最优估计

当取每座阻水桥局部阻力系数 为 0.12 时，根据式（7）和式（9） 可得

把上述系数矩阵和向量代人式（16），可得辨识参数 W 的最优估计为，因为，代人式（5）可得

或者

当考虑每段渠道长度不同对每个残差 的影响时，可取加权矩阵 P 为对角线矩阵且对角线元素 为

式（21）中， 

采用广义最小二乘法参数估计，有

可得辨识参数 W 的最优估计为= ，将其代人式（5） 可得

或者

需要说明，虽然式（20）和式（22）中没有显示出现等效粗糙高度 ks ，但是公式中已经包含了 ks 的影响。

《4 几个重要渠道沿程糙率经验公式比较分析》

4 几个重要渠道沿程糙率经验公式比较分析

目前，国内外都有一些渠道沿程糙率计算经验公式，比较典型的是国内常用经验公式（3） 和美国垦务局经验公式（4），下面根据南水北调中线京石段应急供水工程实测资料把它们与本文两个经验公式（20） 和公式（22）进行比较分析。 应用上述 4 个经验公式，可以计算得到如表 ３ 所示结果，其中渠道编号与表 1 和表 2 一一对应，实测沿程糙率是南水北调中线京石段应急供水工程实测率定资料。

《表3》

表3 典型经验公式渠道糙率计算比较表（按照水力半径 R 的大小排列）

Table 3 Comparison of channel roughness under four formulas（according to the hydraulic radius R）

根据表 3 数据可以画出如图 2 所示曲线，其中曲线 1 和 2 分别对应本文经验公式（20)和(22）；曲线 3 对应国内常用经验公式（3） ， ks =0.000 61；曲线 4 对应美国垦务局经验公式（4） ；点 5 是南水北调中线京石段应急供水工程实测率定。

从图 2 可得如下结论。

1) 国内常用沿程糙率经验公式计算结果，除一段渠道与实测值吻合外，其余渠道均远远小于实测值。 国内经验公式 n 的平均值约为 0.013 6，实测平均值约为 0.014 9，两者平均相对偏差约为 9.5 %。造成误差大的原因主要是该公式是根据实验室资料得到的，并且 ks ＝0.000 61 与 实际工程可能差别较大。

2) 本文沿程糙率经验公式（20） 和公式（22）与美国垦务局经验公式计算结果非常接近。 从表 3 和图 2 可见，本文沿程糙率经验公式（22）和美国垦务局经验公式对应水力半径 R 的 n 值偏差小于0.000 1，相对偏差小于 0.7 %。 说明南水北调中线京石段应急供水工程渠道建设质量达到国外同等水平。

3) 本文沿程糙率经验公式（22） 计算南水北调中线京石段应急供水工程 n 的平均值约为 0.014 8，实测 n 的平均值约为 0.014 9，两者相对偏差小于0.7 %。

综上所述，本文经验公式（22） 考虑了渠道长度的影响，可以作为人工混凝土渠道工程的计算依据。

《图2》

图2 渠道沿程糙率与水力半径关系曲线

Fig.2 Curves of channel roughness versus hydraulic radius

《5 结语》

5 结语

本文的创新点是提出了渠道沿程糙率的系统辨识模型，该模型考虑了渠道几何参数，如断面形态、长度、底坡等的影响，将渠道沿程糙率与粗糙高度 ks 和水力半径 R 的关系用对数函数描述，依据水力学原理，然后通过数学变换提出了适合最小二乘法进行系统辨识的线性模型。

同时，以南水北调中线京石段应急供水工程实测资料为依据，假设每座阻水桥局部阻力系数均为0.12，考虑渠道断面形状、底坡、渠长变化的影响，应用最小二乘法得到的渠道沿程糙率计算公式为

研究也表明，现有国内常用沿程糙率经验公式计算结果除一段渠道与实测值吻合外，其余渠道均远远小于实测值。 本文沿程糙率经验公式与美国垦务局经验公式计算结果非常接近，可以作为其他人工混凝土渠道工程设计的依据。