计算天然河道水面曲线的新方法
金菊良1，杨晓华2，金保明3，丁晶1
(1.四川大学；2. 河海大学；3. 福建省南平市水电局)
摘要：在天然河道水面曲线计算中，传统的试算法或迭代法存在计算误差的累积问题。为此，文中指出，天然河道水面曲线计算的优化问题，实质上就是全河段总计算误差最小化问题，据此同时确定各断面的水位值；从而提出了用加速遗传算法来处理这一仅给出隐含表达式的复杂非线性优化问题。实例计算的结果说明了这一新方法是可行的、通用的和简便的，有效控制了全河段的总计算误差。 关键词：天然河道；水面曲线；计算误差；遗传算法
收稿日期：1999-09-20
基金项目：国家自然科学基金(49871018)、中国博士后科学基金和四川大学高速水力学国家重点实验室开放基金(9904)资助项目。
作者简介：金菊良(1966-)，男，江苏吴江人，四川大学水电学院副教授，博士后。
天然河道水面曲线的计算，就是根据河道地形、纵横断面资料和河道糙率，推求河段在某一流量下各横断面处的水位值，据此即可连出一条对应于该流量的水面曲线。在河道上修建拦河坝、水闸或桥梁等水利工程，会使上游水位抬高，可能会造成部分城市、工矿及农村淹没，为了估计淹没范围需计算壅水水面沿河段的高度，即计算水面曲线；河道疏浚、截湾、分流等工程设计也需进行天然河道水面曲线的计算；在沿河防洪工程规划设计中需计算各设计流量对应的洪水水面曲线；根据河道的预报流量推求相应的水面曲线，为堤防的防洪抗洪措施提供重要的决策支持。由此可见，河道水面曲线的计算不仅是河道水力学计算中的重要问题，而且具有广泛的工程应用背景。目前，计算河道水面曲线的常用方法是试算法、迭代法和图解法[1～3]，这类方法沿河道顺序逐段推求各计算河段横断面的水位值，存在各计算河段的误差累积问题，也不便控制全河段总计算误差。为此，本文提出了同时确定河道计算河段各断面的水位值的思想，这实质上就是全河段总计算误差最小化问题，可用作者研制的加速遗传算法(Accelerating Genetic Algorithm，简称AGA) [4]来处理这一复杂的非线性优化问题。
1 计算天然河道水面曲线的方法
根据河道地形及纵横剖面，把研究的整个河道划分成(m-1)个计算河段，共有m 个横断面。天然河道横断面一般很不规则，糙率及底坡沿程都有变化，岸线曲直相间，而且河底高低不平，兼有河床冲淤的影响。因此，可用分段求和法来研究水位在不规则河道中的沿程变化。对第i-1计算河段，第i-1断面和第i 断面的伯诺里(Bernoulli)能量方程为[1～3]。
(1)
其中z i-1、z i 分别为第i-1断面(上游断面) 和第i 断面(下游断面) 的水位(m)，s i-1、αi-1、ξi-1、n i-1分别为第i-1计算河段的纵向长度(m)、动能修正系数、局部阻力系数和河道糙率，Q 为河道流量(m3/s)，g 为重力加速度(m/s2) ，A i-1、A i 、R i-1、R i 分别为第i-1断面和第i 断面的过水断面面积(m2) 和水力半径(m)，i=2,3,…,m. 因为断面水力要素A 、R 是水位z 的函数，所以计算过程中需要建立各断面水力要素与水位的关系，通常采用拉格朗日(Lagrange)二次插值公式，即[2]
y=y1(x-x2)(x-x3)/[(x1-x 2)(x1-x 3)]+y2(x-x1)(x-x3)/[(x2-x 1)(x2-x 3)]+y3(x-x1) (2
(x-x2)/[(x3-x 1)(x3-x 2)] ) 其中x 1、x 2、x 3为插值基点，y 1、y 2、y 3为插值基点对应的函数值，它们都是已知
的。利用上式可求插值点为x 的插值函数值y. 将水位z 当作x ，断面水力要素如A 、R 等作为y ，就可求得相应于断面不同水位的水力要素值。于是，当流量及其他条件不变时，式(1)左右两端分别转化为第i-1断面水位和第i 断面水位的函数，即
E i-1(zi-1)=zi-1+(αi-1+ξE i (zi )=zi +(αi-1+ξ
)Q 2/(2gA2i-1)-Q 2n 2i-1s i-1/(2A2i-1R 4/3i-1) 222224/3
) i-1)Q /(2gAi )+Qn i-1s i-1/(2Ai R i
i-1
(3)
(4)
目前推求河道水面曲线的常用方法是试算法和迭代法。试算法的计算过程是：1) 已知计算河段的水位(例如z i ) ，按式(4)求函数E i (zi ) 值；2) 假定计算河段另一端水位(例如z i-1) ，按式(3)求出其对应的函数E i-1(zi-1) 值，比较E i (zi ) 与E i-1(zi-1) ，直至这两函数值之差小于允许误差d 时，所假设的值即为所求的水位值；3) 将本河段求出的水位值(例如z i-1) 作为下一计算河段已知断面的水位(例如z i ) 转步骤1) ，依次逐段进行，i=m,m-1,…,2，从而得到对应于某河道流量值的全河段各计算断面的水位值。迭代法的计算过程是：由式(3)和式(4)建立如下迭代公式[3] z i-1=Ei (zi )-(α
i-1
+ξ
i-1
)Q 2/(2gA2i-1)+Q2n 2i-1/si-1/(2A2i-1R 4/3i-1） (5)
与试算法类似，先假定z i-1值，然后按式(5)得到z i-1计算值，若计算值与假定值之差大于允许误差d ，则把计算值作为新的假定值，再按式(5)得到z i-1计算值，若计算值与假定值之差不大于d ，则把计算值作为最终结果；将本河段求出的z i-1值作为下一计算河段已知断面的z i 值，依次逐段进行，i=m,m-1，…,2, 从而得到对应于某河道流量值的全河段各计算断面的水位值。
可见，在试算法和迭代法中，每一计算河段都存在计算误差(允许误差为d) ，而计算又是逐段依次进行的，前面计算河段的计算误差将传递给后面计算河段的计算结果，因此存在各计算河段的误差累积问题，从而无法有效控制全河段的总计算误差。为此，在水利工程设计中常常需要预先确定多个控制断面对应于某河道流量的水位值(例如利用全河段内或其附近的水文站或水位站的水位流量关系曲线) ，以控制计算误差，而这在实际工程设计中时常无法满足。
在这里，本文提出同时确定河道各计算河段断面的水位值，以控制全河段总计算误差。由式(1)至式(4)，可得全河段的总计算误差为
f=｜E 1(z1)-E 2(z2) ｜+｜E 2(z2)-E 3(z3) ｜+…+｜
E m-1(zm-1)-E m (zm ) ｜
(6)
其中｜｜为取绝对值。因此，全河段水面曲线的计算问题实质上就是求如下优化问题
minf=｜E 1(z1)-E 2(z2) ｜+｜E 2(z2)-E 3(z3) ｜+…+｜
E m-1(zm-1)-E m (zm ) ｜
(zi ∈[ai ,b i ],i=1,m-1，假定全河段的最后断面的水位z m 可知)
其中[ai ,b i ]为第i 断面水位z i 的变化范围。该问题是关于所求断面水位｛z i ｜i=m-1,m-2,…,1｝作为变量的一复杂的非线性优化问题，它只给出隐含表达式，目前常规的优化方法很难处理。作为一种通用的优化方法，AGA 可方便地求解上述复杂的优化问题。AGA 的具体算法参见文献[4]。
(7)
2 算例
某河道上修建一座拦河闸，形成闸前全段缓流壅水[1]。已知河道流量为500m 3/s，闸前最高壅水位z 4=477.15m，在需要绘制水面曲线的河段范围，自闸前最高壅水位置将上游壅水段划分成3个计算河段，这3个河段(从上游到下游) 的长度分别为1500m 、1020m 、570m ，动能修正系数均为1.15，局部阻力系数分别取-0.33、0、0，河道糙率均为0.03，这4个断面(断面号按从上游到下游顺序递增) 的水力要素资料如表1所示(根据文献[1]图2～5河道水面曲线计算图查得).
表1 各断面水力要素资料
过水
断面号 i
水位 z/m
断 水力 断面面面半径 号 积 R/m A/m
2
过水水位 z/m
断 水力 断面面面半径 号 积 R/m A/m
2
水位 断面
i i
水力 过水断 水力
断面号 水位 半径 面面积 半径
z/m 面积 i z/m 2
R/m A/m R/m 2
A/m
过水
1 477.00 320 0.714 2 476.00 280 0.857 3 476.00 280 1.714 4 476.00 240 1.571 1 478.00 800 1.429 2 477.00 560 1.714 3 477.00 440 2.143 4 477.00 380 2.143 1 478.50 1120 1.857 2 478.00 960 2.429 3 478.00 640 2.714 4 478.00 560 2.857
现用AGA 同时计算该例各断面的水位值。目标函数取式(7).AGA的群体规模是遗传算法演化过程中每一父代个体的总数目，也就是优化问题的试探解的个数，在这里就等于AGA 每次演化迭代中所试探的水面曲线的条数，优秀个体是指每一父代中相对于目标函数最佳的几个个体。根据作者的经验[4]，群体规模和优秀个体数目可分别取300和10，用AGA 加速循环11次(586微机运行不到1分钟) ，即得各断面的水位值，参见表2，在该表中也列出了试算法和迭代法的计算结果(各计算河段的允许误差均取为0.01m).
表2 用AGA 同时计算各断面的水位值 单位：m
加速次数
z 1 477.0000，478.5000 477.9000，477.9332 477.9186，477.9188 477.9187 477.88 477.9217
[1]
优秀个体的变化区
间 z 2 476.0000，478.0000 477.6073，477.6515 477.6301，477.6303 477.6302 477.60 477.6306
z 3 476.0000，478.0000 477.3821，477.4139 477.3953，477.3955 477.3954 477.35 477.3954
最小目标 函数值f(1)
1 5 11 AGA 结果 试算法结果 迭代法结果
0.36679 0.00391 0.00000 0.00000 0.08188 0.00439
表2说明：1) 从全河段总计算误差(即目标函数) 值来看，AGA 的计算结果(0.00000)明显好于试算法和迭代法的计算结果，而迭代法的结果(0.00439)又好于试算法的结果(0.08188).2)从各断面水位的计算结果看，虽然迭代法能控制每个计算河段的计算误差，但由于误差的累积，沿计算方向(从断面3至断面1) 各断面的水位计算误差存在逐渐增大的趋势；在试算法中各断面的水位计算误差相对都较大；而AGA 则能有效控制全河段总计算误差，解决了误差累积问题，提高了计算精度。
3 结语
在天然河道水面曲线计算中，目前常用的基于逐河段计算的图解法、试算法和迭代法都存在计算误差的累积问题，不能有效控制全河段的总计算误差。为此，文中给出了在全河段计算误差最小化下同时推求各断面的水位值这一优化问题，并用作者提出的加速遗传算法成功地解决了这一仅给出隐含表达式的复杂非线性优化问题。实例计算的结果说明了，这一新方法是可行的、通用的和简便的，有效控制了全河段总计算误差，与常用方法相比，明显提高了计算精度。
参 考 文 献：
[1] 许念曾。河道水力学[M]。北京：中国建材工业出版社,1994.38-55. [2] 周雪漪。计算水力学[M]。北京：清华大学出版社,1995.226-233.
[3] 杨景芳。微机计算水力学[M]。大连：大连理工大学出版社,1991.211-218.
[4] 金菊良，杨晓华，储开凤，等。基因算法的应用及改进[J]。河海大学学报,1998，26(2)：75-79.
计算天然河道水面曲线的新方法
金菊良1，杨晓华2，金保明3，丁晶1
(1.四川大学；2. 河海大学；3. 福建省南平市水电局)
摘要：在天然河道水面曲线计算中，传统的试算法或迭代法存在计算误差的累积问题。为此，文中指出，天然河道水面曲线计算的优化问题，实质上就是全河段总计算误差最小化问题，据此同时确定各断面的水位值；从而提出了用加速遗传算法来处理这一仅给出隐含表达式的复杂非线性优化问题。实例计算的结果说明了这一新方法是可行的、通用的和简便的，有效控制了全河段的总计算误差。 关键词：天然河道；水面曲线；计算误差；遗传算法
收稿日期：1999-09-20
基金项目：国家自然科学基金(49871018)、中国博士后科学基金和四川大学高速水力学国家重点实验室开放基金(9904)资助项目。
作者简介：金菊良(1966-)，男，江苏吴江人，四川大学水电学院副教授，博士后。
天然河道水面曲线的计算，就是根据河道地形、纵横断面资料和河道糙率，推求河段在某一流量下各横断面处的水位值，据此即可连出一条对应于该流量的水面曲线。在河道上修建拦河坝、水闸或桥梁等水利工程，会使上游水位抬高，可能会造成部分城市、工矿及农村淹没，为了估计淹没范围需计算壅水水面沿河段的高度，即计算水面曲线；河道疏浚、截湾、分流等工程设计也需进行天然河道水面曲线的计算；在沿河防洪工程规划设计中需计算各设计流量对应的洪水水面曲线；根据河道的预报流量推求相应的水面曲线，为堤防的防洪抗洪措施提供重要的决策支持。由此可见，河道水面曲线的计算不仅是河道水力学计算中的重要问题，而且具有广泛的工程应用背景。目前，计算河道水面曲线的常用方法是试算法、迭代法和图解法[1～3]，这类方法沿河道顺序逐段推求各计算河段横断面的水位值，存在各计算河段的误差累积问题，也不便控制全河段总计算误差。为此，本文提出了同时确定河道计算河段各断面的水位值的思想，这实质上就是全河段总计算误差最小化问题，可用作者研制的加速遗传算法(Accelerating Genetic Algorithm，简称AGA) [4]来处理这一复杂的非线性优化问题。
1 计算天然河道水面曲线的方法
根据河道地形及纵横剖面，把研究的整个河道划分成(m-1)个计算河段，共有m 个横断面。天然河道横断面一般很不规则，糙率及底坡沿程都有变化，岸线曲直相间，而且河底高低不平，兼有河床冲淤的影响。因此，可用分段求和法来研究水位在不规则河道中的沿程变化。对第i-1计算河段，第i-1断面和第i 断面的伯诺里(Bernoulli)能量方程为[1～3]。
(1)
其中z i-1、z i 分别为第i-1断面(上游断面) 和第i 断面(下游断面) 的水位(m)，s i-1、αi-1、ξi-1、n i-1分别为第i-1计算河段的纵向长度(m)、动能修正系数、局部阻力系数和河道糙率，Q 为河道流量(m3/s)，g 为重力加速度(m/s2) ，A i-1、A i 、R i-1、R i 分别为第i-1断面和第i 断面的过水断面面积(m2) 和水力半径(m)，i=2,3,…,m. 因为断面水力要素A 、R 是水位z 的函数，所以计算过程中需要建立各断面水力要素与水位的关系，通常采用拉格朗日(Lagrange)二次插值公式，即[2]
y=y1(x-x2)(x-x3)/[(x1-x 2)(x1-x 3)]+y2(x-x1)(x-x3)/[(x2-x 1)(x2-x 3)]+y3(x-x1) (2
(x-x2)/[(x3-x 1)(x3-x 2)] ) 其中x 1、x 2、x 3为插值基点，y 1、y 2、y 3为插值基点对应的函数值，它们都是已知
的。利用上式可求插值点为x 的插值函数值y. 将水位z 当作x ，断面水力要素如A 、R 等作为y ，就可求得相应于断面不同水位的水力要素值。于是，当流量及其他条件不变时，式(1)左右两端分别转化为第i-1断面水位和第i 断面水位的函数，即
E i-1(zi-1)=zi-1+(αi-1+ξE i (zi )=zi +(αi-1+ξ
)Q 2/(2gA2i-1)-Q 2n 2i-1s i-1/(2A2i-1R 4/3i-1) 222224/3
) i-1)Q /(2gAi )+Qn i-1s i-1/(2Ai R i
i-1
(3)
(4)
目前推求河道水面曲线的常用方法是试算法和迭代法。试算法的计算过程是：1) 已知计算河段的水位(例如z i ) ，按式(4)求函数E i (zi ) 值；2) 假定计算河段另一端水位(例如z i-1) ，按式(3)求出其对应的函数E i-1(zi-1) 值，比较E i (zi ) 与E i-1(zi-1) ，直至这两函数值之差小于允许误差d 时，所假设的值即为所求的水位值；3) 将本河段求出的水位值(例如z i-1) 作为下一计算河段已知断面的水位(例如z i ) 转步骤1) ，依次逐段进行，i=m,m-1,…,2，从而得到对应于某河道流量值的全河段各计算断面的水位值。迭代法的计算过程是：由式(3)和式(4)建立如下迭代公式[3] z i-1=Ei (zi )-(α
i-1
+ξ
i-1
)Q 2/(2gA2i-1)+Q2n 2i-1/si-1/(2A2i-1R 4/3i-1） (5)
与试算法类似，先假定z i-1值，然后按式(5)得到z i-1计算值，若计算值与假定值之差大于允许误差d ，则把计算值作为新的假定值，再按式(5)得到z i-1计算值，若计算值与假定值之差不大于d ，则把计算值作为最终结果；将本河段求出的z i-1值作为下一计算河段已知断面的z i 值，依次逐段进行，i=m,m-1，…,2, 从而得到对应于某河道流量值的全河段各计算断面的水位值。
可见，在试算法和迭代法中，每一计算河段都存在计算误差(允许误差为d) ，而计算又是逐段依次进行的，前面计算河段的计算误差将传递给后面计算河段的计算结果，因此存在各计算河段的误差累积问题，从而无法有效控制全河段的总计算误差。为此，在水利工程设计中常常需要预先确定多个控制断面对应于某河道流量的水位值(例如利用全河段内或其附近的水文站或水位站的水位流量关系曲线) ，以控制计算误差，而这在实际工程设计中时常无法满足。
在这里，本文提出同时确定河道各计算河段断面的水位值，以控制全河段总计算误差。由式(1)至式(4)，可得全河段的总计算误差为
f=｜E 1(z1)-E 2(z2) ｜+｜E 2(z2)-E 3(z3) ｜+…+｜
E m-1(zm-1)-E m (zm ) ｜
(6)
其中｜｜为取绝对值。因此，全河段水面曲线的计算问题实质上就是求如下优化问题
minf=｜E 1(z1)-E 2(z2) ｜+｜E 2(z2)-E 3(z3) ｜+…+｜
E m-1(zm-1)-E m (zm ) ｜
(zi ∈[ai ,b i ],i=1,m-1，假定全河段的最后断面的水位z m 可知)
其中[ai ,b i ]为第i 断面水位z i 的变化范围。该问题是关于所求断面水位｛z i ｜i=m-1,m-2,…,1｝作为变量的一复杂的非线性优化问题，它只给出隐含表达式，目前常规的优化方法很难处理。作为一种通用的优化方法，AGA 可方便地求解上述复杂的优化问题。AGA 的具体算法参见文献[4]。
(7)
2 算例
某河道上修建一座拦河闸，形成闸前全段缓流壅水[1]。已知河道流量为500m 3/s，闸前最高壅水位z 4=477.15m，在需要绘制水面曲线的河段范围，自闸前最高壅水位置将上游壅水段划分成3个计算河段，这3个河段(从上游到下游) 的长度分别为1500m 、1020m 、570m ，动能修正系数均为1.15，局部阻力系数分别取-0.33、0、0，河道糙率均为0.03，这4个断面(断面号按从上游到下游顺序递增) 的水力要素资料如表1所示(根据文献[1]图2～5河道水面曲线计算图查得).
表1 各断面水力要素资料
过水
断面号 i
水位 z/m
断 水力 断面面面半径 号 积 R/m A/m
2
过水水位 z/m
断 水力 断面面面半径 号 积 R/m A/m
2
水位 断面
i i
水力 过水断 水力
断面号 水位 半径 面面积 半径
z/m 面积 i z/m 2
R/m A/m R/m 2
A/m
过水
1 477.00 320 0.714 2 476.00 280 0.857 3 476.00 280 1.714 4 476.00 240 1.571 1 478.00 800 1.429 2 477.00 560 1.714 3 477.00 440 2.143 4 477.00 380 2.143 1 478.50 1120 1.857 2 478.00 960 2.429 3 478.00 640 2.714 4 478.00 560 2.857
现用AGA 同时计算该例各断面的水位值。目标函数取式(7).AGA的群体规模是遗传算法演化过程中每一父代个体的总数目，也就是优化问题的试探解的个数，在这里就等于AGA 每次演化迭代中所试探的水面曲线的条数，优秀个体是指每一父代中相对于目标函数最佳的几个个体。根据作者的经验[4]，群体规模和优秀个体数目可分别取300和10，用AGA 加速循环11次(586微机运行不到1分钟) ，即得各断面的水位值，参见表2，在该表中也列出了试算法和迭代法的计算结果(各计算河段的允许误差均取为0.01m).
表2 用AGA 同时计算各断面的水位值 单位：m
加速次数
z 1 477.0000，478.5000 477.9000，477.9332 477.9186，477.9188 477.9187 477.88 477.9217
[1]
优秀个体的变化区
间 z 2 476.0000，478.0000 477.6073，477.6515 477.6301，477.6303 477.6302 477.60 477.6306
z 3 476.0000，478.0000 477.3821，477.4139 477.3953，477.3955 477.3954 477.35 477.3954
最小目标 函数值f(1)
1 5 11 AGA 结果 试算法结果 迭代法结果
0.36679 0.00391 0.00000 0.00000 0.08188 0.00439
表2说明：1) 从全河段总计算误差(即目标函数) 值来看，AGA 的计算结果(0.00000)明显好于试算法和迭代法的计算结果，而迭代法的结果(0.00439)又好于试算法的结果(0.08188).2)从各断面水位的计算结果看，虽然迭代法能控制每个计算河段的计算误差，但由于误差的累积，沿计算方向(从断面3至断面1) 各断面的水位计算误差存在逐渐增大的趋势；在试算法中各断面的水位计算误差相对都较大；而AGA 则能有效控制全河段总计算误差，解决了误差累积问题，提高了计算精度。
3 结语
在天然河道水面曲线计算中，目前常用的基于逐河段计算的图解法、试算法和迭代法都存在计算误差的累积问题，不能有效控制全河段的总计算误差。为此，文中给出了在全河段计算误差最小化下同时推求各断面的水位值这一优化问题，并用作者提出的加速遗传算法成功地解决了这一仅给出隐含表达式的复杂非线性优化问题。实例计算的结果说明了，这一新方法是可行的、通用的和简便的，有效控制了全河段总计算误差，与常用方法相比，明显提高了计算精度。
参 考 文 献：
[1] 许念曾。河道水力学[M]。北京：中国建材工业出版社,1994.38-55. [2] 周雪漪。计算水力学[M]。北京：清华大学出版社,1995.226-233.
[3] 杨景芳。微机计算水力学[M]。大连：大连理工大学出版社,1991.211-218.
[4] 金菊良，杨晓华，储开凤，等。基因算法的应用及改进[J]。河海大学学报,1998，26(2)：75-79.