南水北调中线工程渠道糙率的系统辨识

杨开林，汪易森

(1.中国水利水电科学研究院流域水循环模拟与调控国家重点实验室，北京 100038;2.国务院南水北调工程建设委员会专家委员会，北京 100038)

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.1 渠道沿程糙率的计算

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

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

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

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

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

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

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

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

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

式(6)中

由式(1)和(6)可得

整理得

式(8)中，

在系统参数辨识过程中，线性方程式(8)的系数a1、a2和b是已知量，w1和w2是未知量，即需要辨识的参数。

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

其矩阵形式为

式(11)中，

2.3 系统辨识的最小二乘法

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

式(13)中，ε = [ε1，ε2，…，εM]T，其中上标T为转置符号。ε称为拟合误差，或者残差。

式(16)中，矩阵ATA是对称矩阵，当它为非奇异(正则矩阵)时，有

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

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

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

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

3.1 实测数据

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

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

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

表1 渠道特征参数Table 1 Characteristic parameters of channels

表2 渠道参数Table 2 Parameters of channels

图1 实测糙率n与水力半径RFig.1 Calibrated channel roughness n，versus hydraulic radius R

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

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

把上述系数矩阵和向量代人式(16)，可得辨识参数W的最优估计为W^LS=[ ]67.6 28.0T，因为w1=c1lgc2，w2=c1，代人式(5)可得

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

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

或者

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

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

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

表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 渠道沿程糙率与水力半径关系曲线Fig.2 Curves of channel roughnessversus hydraulic radius

5 结语

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

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

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

[1]Becker L，Yeh W W G.Identification of parameters in unsteady open channel flows[J].Water Resources Research，1972，8(4):956－965.

[2]Chiu C L，Emmanuel O I.Kalman filtering in open channel flow estimation[J].Journal of the Hydraulics Division，1978，104(8):1137－1151.

[3]Wormleaton P R， Karmegam M.Parameter optimization in flood routing[J].Journal of Hydraulic Engineering，1983，110(12):1799－1814.

[4]Crissman R D，Chiu C L.Uncertainties in flow modeling and forecasting for Niagara River[J].Journal of Hydraulic Engineering，1993，119(11):1231－1250.

[5]Wasantha Lal A M.Calibration of riverbed roughness[J].Journal of Hydraulic Engineering，1995，121(9):664 －671.

[6]Khatibi R H，Williams J J R，Wormleaton P R.Identification problem of open － channel friction parameters[J].Journal of Hydraulic Engineering，1997，123(12):1078 －1088.

[7]AtanovG A， Evseeva E G，Meselhe E A.Estimation of roughness profile in trapezoidal open channel[J].Journal of Hydraulic Engineering，1999，125(3):309 －312.

[8]Ramesh R，Datta B，Bhallamudi S M，et al.Optimal estimation of roughness in open－channel flows[J].Journal of Hydraulic Engineering，2000，126(4):299 －303.

[9]Yan Ding， Jia Yafei，Sam S Y Wang，et al..Identification of manning’s roughness in shallow water flows[J].Journal of Hydraulic Engineering，2004，130(6):501 －510.

[10]WongT S W， Zhou M C.Kinematic wave parameters and time of travel in circular channel revisited[J].Advances in Water Resources，2003，26:417 －425.

[11]Hsu Minghis， Fu Jincheng，Liu Wencheng.Dynamic routing model with real－time roughness updating for flood forecasting[J].Journal of Hydraulic Engineering，2006，132(6):605 －619.

[12]金忠青，韩龙喜.复杂河网的水力计算及参数反问题[J].水动力学研究与进展，1998，13(3):280－285.

[13]董文军， 杨则燊.一维圣维南方程的反问题研究与计算方法[J].水利学报，2002(9):61－65.

[14]李光炽，周晶晏，张贵寿.用卡尔曼滤波求解河道糙率参数反问题[J].河海大学学报，2003，31(5):490－493.

[15]程伟平，刘国华.基于广义逆理论的河网糙率反演研究[J].浙江大学学报:工学版，2005，39(10):1603－1608.

[16]王光谦，黄跃飞，魏加华，等.南水北调中线工程总干渠糙率综合论证[J].南水北调与水利科技，2006，4(1):8－14.

[17]徐枋同，李永华.系统辨识理论与实践:在水电控制工程中的应用[M].北京:中国电力出版社，1999.