基于奇异矩阵分解法的河道糙率反演计算方法

王玲玲,钟 娜,成高峰

(1.河海大学水利水电学院,江苏南京 210098;2.华东勘测设计研究院,浙江杭州 310014)

河道糙率不仅与表面粗糙程度有关,还受水力要素及水流特性的影响,是一个综合水力摩阻系数,对河道水流及其冲淤变化的计算成果有很大影响.文献[1-3]提出了几种河道糙率的计算公式和方法,但这些计算公式和方法[1-3]大多属于半经验性的,所得结果往往与实际偏差较大.工程实践中常根据河道的实测水文资料通过试错法推求糙率.若河网规模较大,用试错法调试每一河段的糙率,工作量将十分巨大.文献[4-5]提出的系统参数反演的脉冲谱、离散-优化等方法,虽具有理论基础,但不能应用于大规模欠定的河网糙率反分析.本文借助文献[6]的基本思路,采用建立在对输出误差不断进行迭代校正基础上的反分析方法和奇异矩阵分解(singular value decomposition,SVD)法[7-8]建立了河道糙率反演计算方法,并通过编程实现了东江108 km河道的糙率反演,得到了较为精确的糙率分布.

1 反演计算方法的建立

1.1 理论基础

假设待率定参数与有效空间内的变量存在连续函数关系,给待率定参数一个微小的变化量,建立影响系数矩阵,若计算域中各物理量的输出值与实测值之间误差不断减小,则可以求得待率定参数的空间分布.

设ykj表示k(k=1,2,…,l)时刻、j(j=1,2,…,m)位置通过数学模型计算的物理量(如测站水位),Ykj是同一变量的实际观测值.目标则是找出参数xi(糙率),i=1,2,…,n,使误差εkj=ykj-Ykj在ykj的非恒定变化过程中逐渐减小.其中:n为率定糙率参数总数;m为水文测站总数;l为非恒定计算的时刻点总数.

一维河道沿程水位ykj与待率定糙率xi之间的关系满足圣维南方程组,其具体数值求解方法详见文献[9].

用矩阵表示,则式(1)可写成

式中:A k,m×n——影响系数矩阵;k——计算的时刻点.

Ak,m×n的确定过程如下:给每个河道断面糙率参数xi一个初始值并将其代入河道计算模型,计算出初始糙率条件下的各个测站的非恒定水位过程ykj(xi);按照1到n的顺序给每一个待率定参数xi一个很小的增量(如Δxi=0.001)就得出新的参数组合,将其代入河道计算模型,计算出测站非恒定过程的水位值ykj(xi+Δx);将每次计算得出的ykj(xi+Δx)和ykj(xi)代入式(1)进行计算,就可以得到参数xi对j测站物理量的影响系数akji.当有n个参数需要率定时,必须进行n+1次河道模型的计算.

设有m个测站,根据各个测站每一时刻计算值与实测值的误差,可得出各个测站在整个计算时段内的平均误差 εsj:

写成矩阵形式,得

其中

求解方程(5),即可得到满足误差 εs=O的率定参数xi的修正值矩阵Δx.对糙率初始分布进行修正,可获得下一轮计算的新的糙率分布:

式中:r——迭代次数;ρr——控制增量幅度的欠松弛参数.修正过程中,对 xi采用约束条件

予以控制.其中 xub,xlb为指定的上、下边界.

反复进行上述迭代与修正过程,直到各测站水位计算值与实测值误差的2范数‖εs‖2在给定的误差允许范围内,且糙率参数修正值收敛为止.

1.2 SVD方法的应用

上述数理问题求解结果的精度,很大程度上依赖于所形成大型代数方程组的求解方法.对于测站数m远小于待率定断面数n的情况,所形成的代数方程组将严重病态,其求解成为反演计算的关键.SVD方法可直接用于求解m=n,m＞n,m＜n,特别是系数矩阵为奇异矩阵的问题[11].方程(5)的求解可通过将系数矩阵 A s分解成 A s=U∑VT的方法来实现.将式(5)写成

图1 东江河道走势及部分水位测点分布Fig.1 Trend of Dongjiang River and distribution of measuring points

其中 Δz=VTΔx,d=U-1εs.由此可以得到

式中λi为∑对角线上的各元素.如果至少存在一个 λi的值为0或接近于0,则矩阵为奇异矩阵.此时,对应于该 λ值的1/λ均赋值为0.通过这种方式可以得到该奇异矩阵的唯一解.另外,控制式(8)中最小的奇异值 λmin的舍入误差对于减小结果误差也很有效.因为 λmin很小,则1/λ很大,结果的误差会被无限放大.对于非奇异矩阵,用λmax/λmin代替奇异矩阵中的控制条件 λmin,可防止该矩阵成为病态矩阵[12].

2 反演计算方法的应用

将上述反演计算方法应用于广东东江108km长河段的糙率计算.

2.1 水文资料及计算工况

东江是珠江水系三大河流之一,干流全长520 km,集水面积27040 km2.本次计算从三方塘水文站至岗头村水文站,全长108km,沿途设有18个水位测点,河道走势及部分测站位置如图1所示.

沿河划分64个计算断面,即63个河段,空间步长Δx=108.331m,时间步长Δt=20min;河床糙率初值为0.03,欠松弛参数为0.8,上、下边界取 x ub=0.1,x lb=0.001.根据2005年6月2日15:00—4日12:00约50h时间域内16个测点的实测水位资料反演63个河段的糙率值.本算例是个欠定问题,所形成的代数方程组将严重病态.

采用上游三方塘水位站和下游岗头村水位站的实测水位变化过程作为上、下游边界条件(图2).以50h恒定流的计算结果作为非恒定流计算的初始条件.

2.2 糙率反演结果及分析

2.2.1 糙率系数计算收敛过程

采用上述糙率反演计算方法,经过8次迭代计算,得到64个断面糙率系数的稳定结果.图3和图4所示为部分断面糙率系数计算收敛过程和相应各水位测站计算值与实测值误差的降低过程.

图2 上、下游边界水位过程线Fig.2 Hydrographs of upstream and downstream water levels

图3 沿程若干断面糙率系数的收敛过程Fig.3 Convergence process of roughness coefficients of various sections

图4 水位测点计算值与实测值误差的降低过程Fig.4 Reduction process of errors between calculated and measured values

由图3和图4可见,本文计算方法在得到稳定的糙率系数时,计算水位与实测水位误差也稳定减小.

2.2.2 糙率反演结果

图5所示为8次迭代计算后的糙率沿程分布.

图5表明,东江河床的糙率范围在0.02～0.045之间,在平原河流糙率的经验值范围内.结合地形资料分析可知:糙率在0.02～0.03之间的河段主要表现为单一河槽的河段,且河槽中植被较少;糙率大于0.03的河段则主要表现为滩槽共存的复式河段,滩地有较为稠密的杂草,阻水能力较强.因此,由糙率反演结果也可以大致判断河段的断面形态以及植被的稠密程度.

图5 8次迭代计算后的糙率沿程分布Fig.5 Distribution of roughness coefficients after 8 iterations

图6 糙率率定过程中缆西镇水位非恒定过程Fig.6 Unsteady p rocess of water levels in roughness calibration at Lanxi town

2.2.3 非恒定水位逼近过程

虽然沿程各水位测点的实测水位与计算水位的差值在整个非恒定过程中是稳定收敛的(图4),但不足以说明非恒定过程各个时刻计算值与实测值吻合程度的提高,因此需要对糙率率定过程各测站水位计算值与实测值的非恒定变化过程进行比较.图6为糙率率定过程缆西镇测站各迭代步水位计算值与实测值的比较结果.

由图6可以看出,随着糙率的迭代收敛,测站水位的计算值趋向于实测值的过程:由糙率初值(0.03)计算的水位非恒定过程,与实测值相差很大;采用糙率率定模型对初值进行修正后获得的糙率新值,计算的水位过程与实测值的吻合程度,明显好于修正前;随着率定过程的进行,水位结果逐渐趋向实测值;当糙率系数经过8次迭代时,该糙率计算的非恒定过程与实测资料基本吻合.可见,糙率的寻优过程,就是各个测站各个时刻的水位计算值逐渐趋向实测值的过程.同时说明,在糙率率定时,采用水位计算值与实测值的误差在整个非恒定过程中逐渐减小作为目标函数是可行且易于实现的.

3 反演计算方法的验证

为了验证糙率反演计算方法的精度,采用上述糙率率定结果,对2005年6月10日17:00—12日18:00该河段各测点的水位和流量进行了计算,并与M IKE11软件计算结果以及实测结果进行了对比.图7所示为槟榔潭整个计算时段水位本文方法计算值与M IKE11软件计算值以及实测值的比较,图8所示为6月11日16:00该河段流量本文方法计算值与M IKE11软件计算值的比较以及本文方法计算值与观音阁1号、2号和3号3个测站实测值的比较.

图7 槟榔潭水位本文方法计算值与M IKE11软件计算值以及实测值的比较Fig.7 Comparison among values calculated by present method and software M IKE11 and measured ones of water levels at Binglangtan

图8 流量本文方法计算值与MIKE11软件计算值以及观音阁实测值的比较Fig.8 Comparison among values calculated by present method and software M IKE11 and measured ones of discharges at Guanyinge

图7表明,水位本文方法计算值与实测值的误差较小,绝对误差在10 cm以内,而与MIKE11软件计算值的误差较大.图8表明,2种方法的流量计算结果均与实测结果接近,而本文方法的流量计算结果更接近于实测结果.由此可以看出,本文糙率反演方法得到的结果满足精度要求,具有一定的可靠性.

4 结 语

为了获得精确的糙率系数,采用奇异矩阵分解法对河道糙率进行率定,形成了一套包括模型率定、水动力场预测的计算模型,并通过编制程序来实现算法.东江108km天然河道糙率率定结果表明,该糙率率定方法精度较高,可用于河道过流能力及复杂河网水动力模拟分析,可解决欠定的糙率率定问题并提高模型的实用性和洪水预报精度.

[1]何建京,王惠民.流动形态对曼宁糙率系数的影响研究[J].水利水电技术,2002,22(6):22-24.(HE Jian-jing,WANG Huim in.The effects of types of flow on manning's roughness coefficient[J].Hydrology,2002,22(6):22-24.(in Chinese))

[2]程进豪,安连华,王华,等.黄河山东段河床糙率分析[J].水利学报,1997,28(1):39-43.(CHENG Jin-hao,AN Lian-hua,WANG Hua,et al.An analysis of the bed roughness in Shandong reaches of the Yellow River[J].Journal of Hyd raulic Engineering,1997,28(1):39-43.(in Chinese))

[3]李榕.关于影响曼宁粗糙系数n值的水力因素探讨[J].水利学报,1989,20(12):62-66.(LI Rong.The study on hydraulic factors effect on manning's roughness coefficient n[J].Journal of Hydraulic Engineering,1989,20(12):62-66.(in Chinese))

[4]WANG Ling-ling.The algorithms and applications for determining permeability coefficients by pu lse spectrum technique[J].Journal of Hydrodynamics:Ser B,2002,14(2):112-115.

[5]戴会超,王玲玲.参数控制反问题的求解及其在坝工渗流中的应用[J].水力发电学报,2005,24(6):72-77.(DAI Hui-chao,WANG Ling-ling.The algorithms of parameter-control inverse problems and their applications in seepage of dam engineering[J].Journal of Hydroelectric Engineering,2005,24(6):72-77.(in Chinese))

[6]WASANTHA LAL A M.Calibration of riverbed roughness[J].Journal of Hydraulic Engineering,1995,121(9):664-671.

[7]杨克君,曹叔尤,刘兴年.复式河槽综合糙率计算方法比较与分析[J].水利学报,2005,36(7):780-786.(YANG Ke-jun,CAO Shu-you,LIU Xing-nian.Analysis on methods for predicting composite roughness of river channel with compound cross section[J].Journal of Hydraulic Engineering,2005,36(7):780-786.(in Chinese))

[8]潘光斌,陈光裕.基于SVD的测量系统建模方法研究[J].电子测量与仪器学报,2006,20(5):60-62.(PAN Guang-bin,CHEN Guang-yu.Design of indirectly measuring model based on SVD[J].Journal of Electronic Measurement and Instrument,2006,20(5):60-62.(in Chinese))

[9]李光炽,王船海.实用河网水流计算[M].南京:河海大学出版社,2002.

[10]BECKER L,YEH W W-G.Identification of parameters in unsteady open channel flows[J].Water Resour Res,1972,8(4):956-965.

[11]MENEKE W.Geophysical data analysis:discrete inverse theory[M].New York:Academic Press,Inc.,1984.

[12]钟娜.复式河道一维数值模拟[D].南京:河海大学,2007.