一种基于贝叶斯理论的天然河道水位流量关系确定方法与流程

本发明涉及水利工程水文测验技术领域，更具体地，涉及一种基于贝叶斯理论的天然河道水位流量关系确定方法。

背景技术：

水位流量关系是指水文测站断面水位与相应流量之间的关系。由于流量测验技术复杂、耗资昂贵，难以连续进行，在水文资料整编时通常将连续的水位资料，通过水位流量关系转换为连续的流量资料，因此，水位流量关系具有重要的实用意义。

传统的确定水位流量关系曲线方法需要先根据测站特性确定拟合线型，然后根据测站断面多次实测水位和其对应流量数据确定拟合参数，由此确定水位流量关系的数学方程式。拟合线型常选用幂函数型、多项式型和对数函数型，拟合参数则是在确定一定拟合准则后，通过相应算法求解，常用的拟合准则包括残差平方和最小准则、绝对残差绝对值和最小准则以及相对残差绝对值和最小准则，常用的算法包括最小二乘法、遗传算法、蚁群算法、粒子群算法、人工神经网络等。然而上述方法多关注算法原理及不同方法拟合效果对比，拟合线型实际上是根据实测水位流量数据点对的分布形式依靠经验选取的，运用所需数据量较大，线型选择依据不充分，误差来源不清晰，由于未基于水文测站固有特性充分考虑水流运动规律，拟合得到的水位流量关系物理意义不强。众多水文测站具有不同的测站特性，缺乏客观完整的反映水流运动规律的水位流量关系确定方法给生产实践中的定线推流工作带来了较大的不便。

技术实现要素：

本发明为克服上述现有技术所述的至少一种缺陷，提供一种基于贝叶斯理论的天然河道水位流量关系确定方法，仅考虑水位流量为单一对应关系的情况，旨在使水位流量关系的确定能充分利用现有资料和已知信息，参数物理意义明确，误差来源清晰，对特性不同的水文测站均具有较强的适用性。

为解决上述技术问题，本发明采用的技术方案是：一种基于贝叶斯理论的天然河道水位流量关系确定方法，包括以下步骤：

步骤1：水位流量关系模型建立

天然河道水文测站控制断面的过流流量通常用断面平均流速乘以断面面积求得，将断面平均流速及断面面积分别用水位表示，即可得到流量关于水位表达式的一般形式：

q＝a(h-b)^c(1.1)

式中q为流量，h为水位，a为系数，c为指数，b为零流水位。

本发明的步骤3是根据测站特性等已有信息构造水位流量关系表达式中各参数的先验分布，而系数a的表达式综合了多个物理量，无法直接构造先验分布，需要推导其更加详细的表达式。由于断面平均流速和断面面积是决定测站过流流量的两个关键因素，前者取决于测站控制类型，后者取决于断面形状。本发明步骤1将分四部分进一步说明：先识别测站控制类型，分析测站水力学属性推导断面流量关于平均流速和断面面积的表达式，再分析测站几何属性推导断面面积关于水位的表达式，最后联立方程得到流量关于水位的表达式。

步骤1又具体包括以下步骤：

步骤s11：测站控制识别

测站控制是水文测站对水流运动起控制作用的断面或河段的水力因素的总和，可分为断面控制和河槽控制两类：

(1)断面控制：多发生在河流上游，利用测站下游石梁、急滩、卡口、弯道等断面控制，形成临界流，控制断面以上水面线较为稳定且几近水平，控制断面以下水面线下降较快。常在稳定水面线测量断面水位，用堰流或孔口出流的水力学公式确定断面流量，由此建立断面控制的水位流量关系。断面控制的水位流量关系主要受断面形式及其几何形状影响。

(2)河槽控制：多发生在河流下游，依靠具有一定长度的顺直河段来实现。当该河段的水流为均匀流时，在控制断面测量水位，用曼宁公式确定断面流量，由此建立河槽控制的水位流量关系。根据曼宁公式，河槽控制的水位流量关系主要受测量河段过水断面几何形状、平均流动阻力及河底坡降影响。

步骤s12：测站水力学属性分析

(1)断面控制

①堰流

堰流中的薄壁堰具有稳定的水头与流量关系，一般多用于实验室及小河渠的流量测量，薄壁堰流量计算公式为：

式中q为过堰流量；aw为堰顶收缩断面面积；为堰顶收缩断面平均流速，其中为流速系数，k为与堰进口形式和过水断面改变有关的系数，有k＝hco/h0，hco为堰上收缩断面水深；g为重力加速度；h0为堰前全水头。

②孔口出流

本发明仅考虑水体经孔口流入大气即自由出流情况，式中q为孔口出流流量；ao为孔口断面面积；为孔口出流后收缩断面平均流速，其中φ为流速系数；g为重力加速度；h0为孔口前全水头。

(2)河槽控制

测量河段一般为明渠均匀流，用谢才公式确定过水断面平均流速，谢才系数c用曼宁公式计算，由此可得流量：

式中q为测量河段平均流量；a为过水断面面积；n为曼宁糙率系数；r为水力半径；i为河底坡降；水力半径由式r＝a/χ计算，其中χ为湿周。

步骤s13：测站几何属性分析

(1)断面控制

堰前全水头等于堰前压力水头与行近流速水头之和，行近流速水头量级较小，常可忽略不计，故堰上全水头h0可近似等于堰上压力水头h-b。堰顶收缩断面面积与断面形状有关，这里仅建立矩形、抛物线形、三角形3种标准形式常见堰的断面面积关于压力水头(即关于水位)的表达式：

①矩形堰

ar1＝bw1hco＝bw1kh0＝bw1k(h-b)(1.5)

式中ar1为矩形堰顶收缩断面过水面积，bw1为矩形堰净宽，h为堰前水位，b为堰顶高程。

②抛物线形堰

式中ap1为抛物线形堰顶收缩断面过水面积，c1为抛物线形状系数，满足y＝c1x²，带入坐标有其中hp1为漫滩水位，bp1为漫滩水位下的堰宽；h为堰前水位，b为堰顶高程。

③三角形堰

at1＝hco²tan(v1/2)＝k²(h-b)²tan(v1/2)(1.7)

式中at1为三角形堰顶收缩断面过水面积，v1为三角形堰夹角，h为堰前水位，b为堰顶高程。

(2)河槽控制

用曼宁公式计算河道流量时，式中过水断面面积a和水力半径r均与断面形状有关，这里仅建立矩形、抛物线形、三角形三种标准形式常见河道断面面积及水力半径关于平均水深的表达式：

①矩形河道断面

式中，ar2为矩形河道过水断面面积；bw2为矩形河道断面宽度；h为河道水位，b为零流水位；χr为矩形河道湿周；rr矩形河道水力半径。当河道宽浅时，水力半径近似等于平均水深，即rr≈h-b。

②抛物线形河道断面

式中，ap2为抛物线形河道过水断面面积；c2为抛物线形状系数，这里其中hp2为漫滩水位，bp2为漫滩水位下河道宽度；h为河道水位，bh为水位是h时的河宽，b为零流水位；χp为抛物线形河道湿周；rp为抛物线形河道水力半径，消去bh后，当河道宽浅时，即有

③三角形河道断面

式中，at2为三角形河道过水断面面积；v2为三角形河道夹角；h为河道水位，b为零流水位；rt为三角形河道水力半径，χt为三角线形河道湿周。

步骤s14：水位流量关系模型建立

基于测站水力学属性分析和几何属性分析，建立天然河道水文测站断面控制和河槽控制下流量q关于水位h的表达式分别如下：

(1)断面控制

根据测站水力学属性分析，由于行近流速水头量级较小，将其忽略不计的情况下，堰上全水头h0可近似等于堰上压力水头h-b。故式1.2可表达为：

①矩形断面堰流控制

将式1.5代入式1.11中，消去断面面积ar1可得：

式中cr为矩形断面堰流流量系数，

②抛物线形断面堰流控制

将式1.6代入式1.11中，消去断面面积ap1可得：

式中cp为抛物线形断面堰流流量系数，

③三角形断面堰流控制

将式1.7代入式1.11中，消去断面面积at1可得：

式中ct为三角形断面堰流流量系数，

④圆形断面孔口出流控制

同理，忽略孔口出流行近流速水头，将式1.3中的孔口前全水头h0用孔口前压力水头h-b代替可得：

式中co为孔口出流流量系数，co＝φ。

(2)河槽控制

①矩形河道河槽控制

将式1.8中ar2的表达式代入式1.4中，消去过水断面面积a，考虑宽浅河道情况，再将水力半径rr近似表示为平均水深h-b，即可得到：

②抛物线形河道河槽控制

将式1.9中ap2的表达式代入式1.4中，消去过水断面面积a，考虑满足宽浅河道判断条件，再将水力半径公式带入式1.4中，即可得到：

③三角形河道河槽控制

将式1.10中at2和rt的表达式带入式1.4中，消去过水断面面积a及水力半径r，即可得到：

上述步骤推导得出了不同测站控制类型与不同断面几何形状下共7种流量q关于水位h的表达式，发现天然河道水位流量关系模型均符合q＝a(h-b)^c的幂函数形式，只是不同情况下系数a的表达式不同，指数c的值也不同。

步骤2：似然函数计算

水位流量数据为成对出现的离散型随机变量，本发明考虑其为单一对应关系，不考虑复杂绳套关系，且点对相互独立，则似然函数表示参数一定时水位流量点对数据被观测到的概率。由此可建立似然函数的一般表达式：

式中d为水位流量实测数据，共有n对，其中为实测水位，为实测流量；θ为参数集，θ＝(θrc,γ)，包括水力学参数θrc和误差参数γ，其中水力学参数为幂函数水位流量关系模型中的参数，有θrc＝(a,b,c)。

由于水位流量关系模型建立的目的是借此用实测水位推求监测成本较高的流量，即得到的是估计流量，需进一步评估流量估计值的准确度，故将式2.1的似然函数表示成实测流量关于实测水位与参数的联合概率形式，而建立实测流量与实测水位之间的函数关系，需要考虑误差来源。本发明步骤2先根据误差来源建立误差模型，再推导得到实测流量的似然函数表达式。

步骤s21.误差来源

本发明考虑测量误差和模型拟合误差，并将其用概率分布进行描述：

(1)测量误差

水位测量是直接从水尺读数，当测站水流稳定时，可认为水位测量值即为真实值，不存在测量误差；而受限于测流方式和测流精度，流量测量值往往不等于真实值，存在测量误差，可将上述关系概括为如下表达式：

式中为水位流量的测量值；(hi,qi)为水位流量的真实值；假设流量测量误差相互独立且服从正太分布，即εq,i～n(0,μq,i)，其中μq,i为标准差，流量测量误差可依据《河流流量测验规范》获取。

(2)拟合误差

拟合水位流量关系的幂函数模型是基于水流运动规律由水力学公式推导而来的。由于计算断面平均流速时，断面控制忽略了行近流速水头，河槽控制用的曼宁公式为经验公式，且计算断面面积时对过水断面水力半径也做了一定简化，因此用幂函数模型转换实测水位为相应流量估计值的过程中会由于上述假设条件的存在及经验公式的应用而引起误差，称为拟合误差，表示如下：

qi＝f(hi|θrc)+εf,i(2.3)

式中为通过水位流量关系模型转换得到的流量估计值；εf,i为拟合误差，假设拟合误差相互独立且服从正太分布，即εf,i～n(0,σf,i)，σf,i为标准差；基于实践经验，高流量的拟合误差通常较低流量的拟合误差大，这里用线性关系将σf,i表示为随流量估计值变化的函数，即γ为拟合误差参数，其中γ1为常数项，γ2为系数。

将式2.2和式2.3联立可得实测流量关于实测水位的表达式：

步骤s22.似然函数

将式2.4代入式2.1可得到实测流量关于实测水位与参数的联合分布：

为方便说明联合分布的表达方式，当pnorm[z|m,s]表示正太分布的概率密度函数时，z为样本序列，m为均值，s为标准差，式2.5中各项则与之相对应。

步骤3：先验分布构造

先验分布构造是指用概率密度函数对参数集θ进行表示，包括对水力学参数θrc和拟合误差参数γi进行先验分布构造，构造的关键在于利用实测水位流量数据以外的已有信息确定各参数分布形式。假定各参数彼此独立，则参数集先验分布的联合分布形式可表示为：

式中p(θ)为参数集的联合分布，其中p(a)、p(b)、p(c)为3个水力学参数的先验分布，p(γi)为误差参数的先验分布，m为误差参数的数量。先验分布构造主要是确定p(a)、p(b)、p(c)及p(γi)的分布形式及其超参数(以正太分布为例，正太分布的超参数为均值和标准差)的取值。本发明步骤3将从水力学参数先验分布构造及拟合误差参数先验分布构造两部分进行说明：

s31.水力学参数先验分布

由于天然河道测站控制条件和几何形状并不总是如步骤1中标准形式那样规则，因此需要先对天然河道水文站的测站控制类别进行识别，确定是断面控制还是河槽控制，然后用步骤1中建立的标准几何形状(矩形、抛物线形、三角形)去类比天然河道几何形状，参照标准几何形状测站控制下系数a的表达式、零流水位b的物理意义及指数c的值来构造各自的先验分布。本发明用正太分布描述3个水力学参数的先验分布，尚需进一步确定先验分布的超参数，即均值和标准差。

(1)系数a先验分布p(a)

根据步骤1中水位流量关系模型的最终表达式，系数a与多个物理量有关，需要根据相应物理量实际意义构造物理量参数的分布形式，再合成a的先验分布。a的均值可用各物理量参数的均值按a的表达式直接求解，a的标准差则需要在确定各物理量参数的标准差后，根据gum(jcgm2008)不确定度合成公式(式3.2)的一阶泰勒级数展开法计算。限于篇幅，本发明仅对矩形堰流断面控制和矩形河道河槽控制的不确定度合成公式进行推导，分别如式3.3和式3.4所示。

式中xi(i∈[1,n])为物理量参数，μ²(a)为a的方差，μ²(xi)为物理量参数方差。

矩形堰流断面控制：

矩形河道河槽控制：

式中ks为曼宁糙率系数n的倒数，即ks＝1/n，s为河底坡降，这里假设水面比降近似等于河底坡降，即i≈s。

(2)零流水位b先验分布p(b)

零流水位b的物理意义，断面控制时表示堰顶高程，河槽控制时表示平均河底高程，可直接根据测站断面特性确定b先验分布的均值和标准差。

(3)指数c先验分布p(c)

根据步骤1水位流量关系模型表示的推导过程，特定测站控制类型及断面几何形状下指数c的值恒定，当选定标准几何形状断面近似表示天然河道断面后，可将该标准情况下指数c值作为天然河道测站水位流量关系模型指数c值先验分布的均值，再根据精度要求设定标准差即可完成先验分布的构造。

s32.拟合误差参数先验分布

由于流量的真实值未知，因此较难确定拟合误差εf,i，对其标准差项中的超参数γ1和γ2采用宽均匀分布u(0,10000)进行表示。

步骤4.后验分布求解

水位流量关系模型可用详细的物理量参数表达(步骤1)，根据测站已有信息建立物理量参数的先验分布(步骤3)，将实测水位流量数据信息用似然函数表示(步骤2)，基于贝叶斯公式，即可得到参数集θ后验分布的表达式：

式中为后验分布的概率密度函数；为似然函数，同式2.5；p(θ)为先验分布，同式3.1；为边际分布，与参数无关；由于参数维度为5维，边际分布难以用常规数值积分方法求解。而马尔科夫链蒙特卡罗(mcmc)方法用于贝叶斯后验分布求解过程中，在接受率计算时可将边际分布抵消掉，求解效率较高，本发明选用该方模拟求解参数集θ的后验分布，边际分布视为一个正则化因子，由此可将式(4.1)简化成如下形式：

式中符号∝表示正比于，根据参数后验分布模拟结果，取参数均值可建立最大出现概率水位流量关系模型，取参数均值加减2倍标准差可建立具有95％置信水平的水位流量关系模型的上下包络线。为提高模拟计算效率，本发明设计了自适应mcmc算法求解参数集θ的后验分布，算法设计如下：

step1设定初值：由先验分布生成初始参数集由此计算初始协方差矩阵v＝(v1,...,vp)，设定参数计数变量初值k＝0；

step2迭代求解：

fori＝1:ncycles(自适应循环求解结构)

forj＝1:nddapt(非自适应循环求解结构)

k＝k+1；

ford＝1:p

①构造提议分布进行状态转移生成

②计算状态转移概率

③状态转移接受概率为αd＝min(1；τ)，当接受状态转移时，

当拒绝状态转移时，则

ford＝1:p

计算第d个参数的接受率αd；

如果αd≤αmin，vd＝φ^-×vd；(接受概率过低时减小协方差)

如果αd≥αmax，vd＝φ⁺×vd；(接受概率过高时增大协方差)

如果αmin≤αd≤αmax，vd＝vd；(接受概率适中时不做改变)

算法设计nddapt＝100，ncycles＝100，表示先进行100次普通计算，得到参数接受概率αd后经过判断再进行100次自适应计算；自适应计算时设计αmin＝0.1，αmax＝0.5，φ^-＝0.9，φ⁺＝1.1，认为αd落在[0.1,0.5]区间内时，模拟效果良好，当αd小于区间下限0.1时，需要减小协方差以提高接受概率，减小系数为0.9，当αd大于区间上限0.5时，需要增大协方差以减小接受概率，增大系数为1.1。

与现有技术相比，有益效果是：

1.确定稳定水位流量关系的影响因素为测站控制属性，建立的幂函数模型充分考虑了水文测站水流运动规律，参数的物理意义明确；

2.传统方法都将参数视为固定值，本发明所用的贝叶斯方法将参数视为服从一定分布的随机变量，在得到参数计算结果的同时还可直接得到其出现概率，便于进行不确定度分析；

3.传统方法认为水位流量实测值即为真实值，误差仅为拟合误差，本发明引入了测量误差，构造包含测量误差和拟合误差的误差模型计算似然函数，且将拟合误差设计成随流量增大而增大的线性形式，更加符合实际；

4.传统方法水位流量关系的确定全部取决于实测数据，本发明应用贝叶斯方法考虑参数物理意义，根据物理量参数已有信息即可初步建立水位流量关系先验模型，将实测数据耦合到似然函数中对先验模型进行修正，可最大化利用已有信息，减少实测数据获取，降低测量成本；

5.传统方法基于拟合误差最小准则用优化算法求解参数，目标函数为复杂的非线性形式，存在计算量大，收敛困难等问题，本发明基于贝叶斯理论用mcmc方法直接求解参数的后验分布，且设计了自适应算法提高了计算效率。

附图说明

图1是本发明方法流程图。

图2是本发明测站控制示意图，其中(a)为断面控制，(b)为河槽控制。

图3是本发明断面控制示意图，其中(a)为矩形断面堰流控制，(b)为抛物线形断面堰流控制，(c)为三角形断面堰流控制，(d)为孔口出流断面控制。

图4是本发明河槽控制示意图，其中(a)为矩形河道河槽控制，(b)为抛物线形河道河槽控制，(c)三角形河道河槽控制。

图5是本发明测流断面横剖图，虚线框为标准形式矩形河道河槽控制几何形状。

图6是本发明水位流量关系先验模型，图中qmaxprior是参数先验分布均值对应的曲线，qhigh_prior和qlow_prior是参数后验分布95％置信区间对应的曲线。

图7是本发明水位流量关系后验模型，图中qmaxpost是参数后验分布均值对应的曲线，qhigh_post和qlow_post是参数后验分布95％置信区间对应的曲线。

具体实施方式

附图仅用于示例性说明，不能理解为对本发明的限制；为了更好说明本实施例，附图某些部件会有省略、放大或缩小，并不代表实际产品的尺寸；对于本领域技术人员来说，附图中某些公知结构及其说明可能省略是可以理解的。附图中描述位置关系仅用于示例性说明，不能理解为对本发明的限制。

实施例1：

如图1所示，一种基于贝叶斯理论的天然河道水位流量关系确定方法，包括以下步骤：

步骤1.水位流量关系模型建立。以某水文站为研究对象，收集该站水文测验报告等相关资料，得到某年份该站实测水位流量数据以及测流大断面数据，从工作人员处还得知该水文站测流河段顺直，测流断面附近无石梁、卡口等结构，该水文站的测站控制可视为河槽控制，水位流量关系由式1.4所示的曼宁公式确定。根据图5所示的测流断面横剖图，测站控制近似为矩形，分析各测次平均水深与平均河宽的比例，满足宽浅河道条件(b≥5(h-b))，水力半径近似等于平均水深，由此可建立式1.16所表达的水位流量关系模型。实测水位流量数据共19对，如图6和图7中guagings所示，这里不单独列出。

步骤2.似然函数计算。计算似然函数需要先建立误差模型，主要考虑测量误差和拟合误差两种误差来源。该测站在水位自记井内用浮子和传感器自动记录水位，测量结果精确，可将测量水位视为真实水位，测验断面设在基下8m的位置，用船载adcp测流，测流误差在实测流量5％～7％范围内服从正太分布，由此可建立式2.2所表达的测量误差模型。将拟合误差表达成估计流量的线性函数，即高流量的拟合误差大于低流量的拟合误差，由此建立式2.3所表达的拟合误差模型。联立式2.1、2.2、2.3即可得到式2.5所示的实测流量关于实测水位与参数集的似然函数。

步骤3.先验分布构造。3个水力学参数的先验分布用正太分布pnorm[z|m,s]表示，需要确定各参数先验分布的超参数(均值和标准差)，2个拟合误差参数的先验分布用宽均匀分布u(0,10000)表示，这里对水力学参数先验分布的构造做详细说明。图5给出了测站控制断面横剖图，观察其几何形状，用标准形式矩形河道河槽控制(如图4(a))近似表示该天然河道河槽控制。由式1.16所示的标准形式矩形河道河槽控制水位流量关系模型，系数包含3个物理量参数，其中n为曼宁糙率系数，i为水面比降，bw2为标准矩形宽度；b为零流水位；指数c为恒定值，等于1.67。构造系数a的先验分布要先构造出3个物理量参数的先验分布，a的先验分布超参数——均值根据其表达式直接由3个物理量均值计算得到，a的先验分布超参数——标准差则按式3.4所示的不确定度传播公式计算，为方便计算，将用ks表示，即水面比降i用河底坡降s表示，即i≈s，因此可得零流水位b的先验分布根据断面平均河底高程设定。将指数c先验分布的均值设成1.67，再设定合理的标准差即可。构造出的参数集先验分布如下：

系数a的物理量参数：ks～n(25,5)；bw～n(60,5)；s～n(0.003,0.001)；

水力学参数：a～n(82.16,22.46)；c～n(1.67,0.025)；b～n(-0.5,0.25)；

拟合误差参数：γ1～u(0,10000)；γ2～u(0,10000)；

基于上述参数集的先验分布即可初步建立水位流量关系先验模型，如图6所示，先验模型的95％置信区间包络线范围较大表示其不确定度较高，实测水位流量散点集中分布于95％置信区间下包络线的附近。

步骤4.后验分布求解。根据式4.1所示的贝叶斯公式，将实测数据耦合到似然函数中对先验模型进行修正，用设计的自适应mcmc算法对参数集求解，可得水位流量关系后验模型，如图7所示，发现经实测数据修正后的水位流量关系模型置信区间大大缩小，不确定度较低，实测水位流量散点均均匀分布在置信区间内，可利用qmaxpost所示的曲线进行流量的插补延长。

显然，本发明的上述实施例仅仅是为清楚地说明本发明所作的举例，而并非是对本发明的实施方式的限定。对于所属领域的普通技术人员来说，在上述说明的基础上还可以做出其它不同形式的变化或变动。这里无需也无法对所有的实施方式予以穷举。凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明权利要求的保护范围之内。