河道糙率反问题在实时洪水预报中的应用

第二十一届全国水动力学研讨会暨第八届全国水动力学学术会议 暨两岸船舶与海洋水动力学研讨会文集 - 902 - 河道糙率反问题在实时洪水预报中的应用 吴晓玲，王船海 (河海大学，水文水资源与水利工程科学国家重点实验室，南京 210098) 摘要： 针对洪水预报中糙率参数的实时修正问题，采用卡尔曼滤波技术，实现动量方程 实时校正糙率参数，以提高水动力学模型的洪水预报精度。对复杂地形所引起的局部动量损失 进行处理，使反演的糙率参数更显客观。 选择长江干流寸滩至万县河段进行试验，应用 2004年该河段洪水资料作验证。结果表明， 反演的糙率参数对提高实时洪水预报精度有效。 关键词： 实时洪水预报；卡尔曼滤波；糙率参数；反问题 1 前言 对于水动力学计算中糙率参数的选择问题，国内外的学者倾注了大量的研究精力。杨世孝 等[1]依据反问题方法研究河网的非恒定流计算，利用各时段各河段的糙率平均值动态调整河段 的平均糙率；金忠青等[2]根据河道规模及水力特性，把河道划分为若干等级，采用最优化中的 复合形法对不同等级的河道糙率分别求解；董文军等[3]利用最小二乘逼近思想建立曼宁糙率的 辨识模型；李光炽[4]以曼宁及动量方程为基础，依据实测水位资料借助卡尔曼滤波技术反 推河道糙率。葛守西等[5, 6]针对应用卡尔曼滤波技术解决河段糙率问题，提出以水位、流量的 涨差作为观测量，建立以糙率为状态量的线性滤波器方程，通过实测水位、流量的变化对糙率 进行反推计算。 Fread等人[7]通过流量分级，分段拟和平均流量与糙率的关系函数，运用牛顿־拉斐森算法 得到当水位误差函数取最小时糙率与平均流量的关系；Khatibi等人[8, 9]应用最小二乘法，依据 实测水位解决潮汐河道糙率的最优估计问题，通过比较三种以水位作为“验证项”的最小二乘构 造，得出相对高潮水位而言，低潮位对糙率的反推估计更加有效的结论。 受上述文献的启发，本研究采用四点线性隐格式[10]求解圣维南方程组，分离动量方程中 的水位及糙率项，依据卡尔曼滤波量测方程形式，建立水位与糙率的对应关系，运用滤波器的 最优估计原理，反推河段糙率；同时针对河底高程变化剧烈的河段考虑局部动量损失的影响， 把反推得到的糙率应用到实时洪水预报中，通过实时动态调整河道糙率参数提高洪水预报精 度。 国家支撑项目(2006BAC05B02)和国家自然科学基金资助项目(50679025)资助. 第二十一届全国水动力学研讨会暨第八届全国水动力学学术会议 暨两岸船舶与海洋工程水动力学研讨会文集 - 903 - 2 模型的建立 2.1 河道糙率反问题-离散连续方程和动量方程 描述河道水流运动的圣维南方程组为 2 2 4/3 up ( ) 0 ( , , ) 0 ( , , ) 0 Z QB q t x u QQ Q ZgA gn t x A x R f Q Z t f Q Z t α ∂ ∂ + = ∂ ∂∂ ∂ ∂ + + + = ∂ ∂ ∂ = = down 上边界条件 下边界条件 ： ： (1) 式中：q为旁侧入流；Q、A、B、Z为河道断面流量、过水面积、河宽和水位；n为河道糙率； R为水力半径；α为动量校正系数。设河道被（L+1）个断面分成 L个子河段，对任一子河段 （j，j+1）应用四点线性隐格式离散式(1)，得到如下差分方程，    +=−++ ++=++− + ++ + + + + + ++ + ++ + )( )( 1 11 1 1 1 1 1 11 1 11 1 k j k jj k jj k jj k jj k jj k j k jjj k jj k jj k j k j QQHZFZFQGQE ZZCDZCZCQQ (2) 式中 Cj，Dj，Ej，Fj，Gj，Hj与断面形状，断面间距，旁侧入流及糙率等有关，当他们都 已知时，即由式(2)结合初始条件求水位、流量称之为正问题；若各河段的糙率参数预先未知， 则需要通过验证计算来确定糙率，使由式(2)求出的水位、流量的某些分量与相应实测值吻合 得最好，此为反问题[1]。反问题是一个参数辨识的过程，根据计算域中测点观测值，通过对所 建立的目标函数求最小来确定控制微分方程中的参数，因此在参数辩识中往往含有数学和 最优控制等技术[3]。 影响水力计算的因素包括河道的糙率、旁侧入流、上下游边界条件，对于非恒定流计算还 包括初始条件，一般情况下，初始条件随时间推移对流态的影响减弱。对于不同场洪水，旁侧 入流和边界条件均各异，而河道糙率与水位之间存在一定的联系，本研究针对河道糙率的反问 题进行研究，在实时校正中动态修正糙率，提高预报精度。 2.2卡尔曼滤波基本原理 一个动态系统的状态变化，根据经典物理学可知，在没有外部干扰时，系统未来状态可以 根据现时状态从已知的运动方程中确定。但总存在某种外部干扰或者对运动方程的描述不十分 精确的情况，卡尔曼滤波理论的基本思想就是采用信号与噪声的状态空间模型，利用前一时刻 的估计和观测值来更新对状态变量的估计，求现时刻的估计值[11, 12]。卡尔曼滤波系统的基本 方程由状态方程和量测方程组成，通常可表示为： 1 1| 1|k k k k k k kX X wΦ Γ+ + += + (3) 1 1 1 1 1k k k k kY H X U v+ + + + += + + （ 0≥k ） (4) 第二十一届全国水动力学研讨会暨第八届全国水动力学学术会议 暨两岸船舶与海洋工程水动力学研讨会文集 - 904 - 滤波公式按(5)~ (8)依次执行，便可实现卡尔曼滤波计算。 1| 1| ˆ ˆ k k k k kX XΦ+ += (5) 1 1| 1 1 1 1 1| ˆ ˆ ˆ( )k k k k k k k k kX X K Y U H X+ + + + + + += + − − (6) 1 1 1 1 1 1 1 1| |( )k k k k k k k k kK P H H P H R − + + + + + + += +T T (7) 1| 1 1 1 1|( )k k k k k kP I K H P+ + + + += − (8) 可见卡尔曼滤波的递推公式既可以得到滤波估计值，又可以得到误差的方差阵，即产生它 自身的误差，此为卡尔曼滤波器的一个重要优点。 2.3 卡尔曼滤波求解糙率反问题 由于在动量方程中参与计算的糙率是以平方的形式出现，作为状态量不太合适，因此，引 用参数θ，使糙率 new oldn n θ= × ，其中 nold为上一个时刻的糙率（计算过程中认为已知）； θ作为反问题中的研究对象，反映糙率的变化。 构建以θ作为状态量，水位作为观测量的卡尔曼滤波方程，对于内部有 m 个水位观测站 的河道，依据水位观测站划分为(m+1)个河段，不同的河段对应不同的糙率，状态方程为： 1 1| 1|k k k k k k kwθ Φ θ Γ+ + += + (9) 以上边界为流量条件( upQ )，下边界为水位条件( downZ )为例，建立θ与水位的关系： 1 1 1 1k k k kAZ B Cθ+ + + += + (10) 1 1 0 1 0 1 1 A − −    =     O O ； 1 0 1 1 0 0 k k m B Φ Φ Φ + +    =     O ； 1 0 1 1 m down k k B B C B Z + + −    =     M ； 1 1i i S i j j S Φ ϕ+ − = = ∑ ； 1 1i i S i j j S B β+ − = = ∑ ； 24 3 1 2 ( ) k j j j old jR u u x nϕ +  = ∆   ； 1 1 1 1 1 1 1 1 1 1 2 2 1( ) ( ) 2 ( ) ( ) j k k k k k k k k j j j j j j j j jk k j j Q Q Q Q Q Q gA x u u t gA β + + + ++ + + + + + ∆= − + − + −∆ 其中：0≤i≤m； Si≤j＜Si+1；Si为第 i个河段的首断面编号，0≤Si≤L。 依据一般情况，nold 的初值为初始水位对应的糙率，θ =1.0， { } 0kE ν = 应用式(5)至式(8) 第二十一届全国水动力学研讨会暨第八届全国水动力学学术会议 暨两岸船舶与海洋工程水动力学研讨会文集 - 905 - 以观测水位为依据校正糙率的θ参数。 2.4 局部动量变化在糙率反演中的应用 对于河流上游地形复杂的河道中，水的流态对糙率反推的影响也十分重要。本研究在上述 考虑的基础上对原有的动量方程作了稍许修改，主要针对断面之间地形突变情况。动量方程是 依据动量守恒定律建立，通过控制面流进控制体内的动量与作用于控制体外力的冲量之和，即 为控制体内动量冲量的增量[10]。对于流态急剧变化的局部动量损失一般不考虑，但是实际工 程中往往会遇到复杂地形所造成的复杂流态。相邻断面流速差距很大的河段应增加局部水头动 量损失，使其更完整地表达实际情况，则式(10)变为： 1 1 1 1 1k k k k kAZ B C Dθ+ + + + += + + ，相应的 式(10)写成： 1 11 1 1 1 1 1( ) ( )k k k k k kZ A B A C D vθ− −+ + + + + += + + + (11) 简化为： 1 1 1 1 1k k k k kZ H U vθ+ + + + += + + (12) 其中： 11 1k kH A B−+ += ； 11 1 1( )k k kU A C D−+ + += + ； [ ] 11 0 1 kk mD ϒ ϒ ϒ ++ ′= L ； 1 1i i S i s s S ϒ γ+ − = = ∑ ， 21( )2 k k s s s u u g γ ζ +−= ，s 为流速陡变的河断编号，Si≤s＜Si+1；Si为第 i 个河段的 首断面编号，0≤Si≤L。ζ 为局部阻力系数，取 0.5。 在实时洪水预报中，首先依据卡尔曼滤波技术对水位以及流量进行校正，此时 Qk+1认为 已知。在此基础上，构建式(9)和式(12)作为由水位实测值校正糙率的反问题解决方法，形成河 道水流与糙率的联合校正。在洪水预报中考虑局部动量损失的方式进行实测水位反推河段糙 率，实时动态调整不同水位条件下的糙率参数 new oldn n θ= × 。 3应用实例 3.1水文资料情况 选用长江干流寸滩至万县段为试验河段验证本文思路，该河段位于长江中上游三峡区间流 域内，总长约 321km，其间有寸滩、 清溪场、万县三个水文站；长寿，忠 县两个水位站（图 1）。 试验资料为 2004年 8月 18日 0 时至 10月 17日 0时共约 60d洪水过 程。以寸滩入流流量为上边界，武隆 支流作旁侧入流处理，万县水位为下 边界。河段除下边界以外，共有寸滩、 长寿、清溪场、忠县四个水位实测站 点，把整个河段分成四段。原始糙率参数是通过 2003 年三峡蓄水 135m后的洪水资料，应用 图 1 河道测站分布 第二十一届全国水动力学研讨会暨第八届全国水动力学学术会议 暨两岸船舶与海洋工程水动力学研讨会文集 - 906 - 一维水动力学模型[13]率定而得，把原始糙率的一维水动力学模型计算结果作为本次研究对比 的依据，称“原水动力学模型”。 该河段地形复杂，相邻断面之间河底高程起伏较大如图 2（黄海基面）所示，因此在动量 方程反推糙率时对流态变化较大的断面之间考虑局部动量变化。 3.2 水模型计算及结果 分析 本研究模型通过前 一时刻各站的实测水位 资料反推θ ，对模型糙 率参数进行实时调整， 用调整后的 nnew参与水 动力学预报（“糙率反推 预报模型”）。模型采取 “滚动预报”的方式进行，以寸滩，长寿，清溪场，忠县实测水位作为卡尔曼滤波计算校正的依 据。 与“原水动力学模型”相比，“糙率反推预报模型” 6 h预报结果在洪峰水位及低水位时都能 较接近于实测过程；且各站洪峰水位误差在“原水动力学模型”基础上降低 0.05~0.15m，同时， RMSE（误差均方根）也降低 0.18~0.32m，确定性系数普遍在 0.98以上（表 1），可见应用“糙 率反推预报模型”实时修正水动力学预报模型中的糙率参数，对于改善实时预报精度有所帮助。 表 1 2004年不同模式实时外推 6h预报效果分析 洪峰水位误差/m 确定性系数 RMSE/m 站名 方法 1 方法 2 方法 1 方法 2 方法 1 方法 2 寸滩 -0.46 -0.41 0.98 0.99 0.54 0.35 长寿 -0.30 -0.23 0.96 0.99 0.74 0.47 清溪场 0.19 0.14 0.95 0.98 0.79 0.47 忠县 0.47 0.32 0.98 1.00 0.38 0.20 注：方法 1表示“原水动力学模型”；方法 2表示“糙率反推预报模型”。 4 结论 针对水动力学实时预报系统中糙率参数估计问题，本研究提出糙率参数反推与卡尔曼滤波 技术相结合的方法，并对于地形变化突然的山区形河道考虑局部动量损失。将水位反推得到的 糙率应用到实时预报中，动态调整预报期的河道糙率参数。长江上游素来以险著称，河道地形 复杂多变，以三峡上游寸滩至万县段为例验证本文思路，寸滩、长寿、清溪场、忠县四个水位 实测站 6 h预报水位过程与实测值吻合，初步显示本方法可行。 由于本研究的试验河段取自山区性明显的长江上游干流河道，因此对于本研究方法是否适 应于平原河网地区有待继续深入研究。 40 55 70 85 100 115 130 145 160 175 0 50 100 150 200 250 300距离/km 高 程 /m 河底高程 2004-10-15 21:00水面线 水位观测站 寸滩 清溪场长寿 万县忠县 图 2河底高程及水面线示意图 第二十一届全国水动力学研讨会暨第八届全国水动力学学术会议 暨两岸船舶与海洋工程水动力学研讨会文集 - 907 - 参考文献 1杨世孝，肖子良．反求糙率的一种数值方法[J ]．数值计算与计算机应用，1994，12：247–260． 2 金忠青，韩龙喜，张健．复杂河网的水力计算及参数反问题[J]．水动力学研究与进展，1998：280–285． 3董文军，杨则燊．一维圣维南方程的反问题研究与计算方法[J]．水利学报，2002，9：61–65． 4 李光炽，周晶晏，张贵寿．用卡尔曼滤波求解河道糙率参数反问题[J]．河海大学学报( 自然科学版)，2003， 31（5）：490–493． 5 葛守西，程海云，李玉荣．水动力学模型卡尔曼滤波实时校正技术．水利学报．2005，6（36）：687–693． 6 程海云,葛守西,李玉荣．水动力学模型实时校正方法研究[J]．人民长江，2005，36(2)：6–8． 7 Fread D L, Smith G F. Calibration techniques for 1-D unsteady flow models [J]．Journal of Hydraulic Div, ASCE,1977, 104(7): 1 027–1 043. 8 Rahman H Khatibi, John J R Williams, Peter R Wormleaton. Identification Problem of Open-Channel Friction Parameters [J]．J Hydr Engrg, 1997, 123(12): 1 078–1 088. 9 Rahman H Khatibi, John J R Williams, Peter R Wormleaton. Friction Parameters for Flows in Nearly Flat Tidal Channels [J]．J Hydr Engrg, 2000, 126(10): 741–749. 10 王船海，李光炽．实用河网水流计算[M]．南京：河海大学，2003：9． 11 宋文尧，张牙．卡尔曼滤波[M]．北京：科学出版社，1991． 12 朱海，莫军．水下导航信息融合技术[M]．北京：国防工业出版社，2002：80–81． 13 王船海，李光炽．流域洪水模拟[J]．水利学报，1996 (3）：44–50． Applying inverse analysis of river channel roughness parameter in real time flood forecast WU Xiao-ling, WANG Chuan-hai (State Key Laboratory of Hydrology-Water Resources and Hydraulic Engineering, Hohai University, Nanjing 210098, E-mail: freebir7237@hhu.edu.cn ) Abstract: To detect the effect of roughness parameters in real time correction for flood forecasting, and improve the accuracy of hydraulic forecasting model, the Kalman filter was applied to realize the real time correction of the roughness parameters by the momentum equation. In order to reflect the roughness parameters impersonality, the local loss of momentum item was added to momentum equation where complex flow happened between adjacent cross sections of river channel by complex terrain. The calculation of flood forecasting for river section from Cuntan to Wanxian of Yangtze River verifies that the inverse analysis of channel parameters can be resultful in promoting the accuracy of forecasting by adopting the flood information of 2004. Key words: real-time flood forecast; Kalman filter; roughness parameter; inverse analysis