CN115526074A - 基于河流表面流速的截面流速分布反演方法 - Google Patents

Abstract

本发明公开了一种基于河流表面流速的截面流速分布反演方法，包括：对河流截面进行网格化分割，将各网格点分别划分为内区和外区；计算外区涡流粘度并根据计算得到的外区涡流粘度结合雷达测量的表面流速分布得到内区和外区分界点处流速；根据内区和外区分界点处流速计算各内区网格点流速；根据雷达测量的表面流速分布以及内区和外区分界点处流速，计算各外区网格点流速；将反演的各内区网格点流速、内区和外区分界点处流速、外区网格点流速以及雷达测量的表面流速相组合得到截面流速分布。本发明通过分析截面内区流速分布和外区流速分布的差异和联系，分别采用对数亏损律和RANS方程反演内区流速分布和外区流速分布，从而获取河流截面的流速分布。

Description

基于河流表面流速的截面流速分布反演方法

技术领域

本发明属于微波雷达河流遥感的技术领域，具体涉及一种基于河流表面流速的截面流速分布反演方法。

背景技术

河流截面流速分布是一个关键的水文信息，它对流量估计、河道演化以及水沙输送均有重要意义。传统流速分布的测量工具如测流表、声学多普勒流速剖面仪仅能提供局部有效的流速数据，在工作时需要浸入水中，导致仪器的安装和维护成本较高。且这些仪器操作繁琐，单次截面流速分布测量的工作量偏大，实时性较差。随着电子技术的快速发展，以微波雷达为代表的非接触式测流技术逐渐应用于河流水文监测领域。微波雷达可以在不接触水体的情况下获得实时、准确、高分辨率的表面流速分布，然而微波雷达只能探测河流的表面流速，无法直接测量水面以下的流速，如何基于表面流速反演截面流速分布是微波河流遥感中一个亟需解决的问题。

发明内容

本发明的目的在于针对现有技术的不足之处，提供一种基于河流表面流速的截面流速分布反演方法，该方法通过分析截面内区流速分布和外区流速分布的差异和联系获取河流截面流速分布，解决了微波雷达无法测量水面以下流速的问题。

为解决上述技术问题，本发明采用如下技术方案：

一种基于河流表面流速的截面流速分布反演方法，其特征在于，包括如下步骤：

步骤1：对河流截面进行网格化分割，将各网格点分别划分为内区和外区；

步骤2：计算外区涡流粘度并根据计算得到的外区涡流粘度结合雷达测量的表面流速分布得到内区和外区分界点处流速；

步骤3：根据步骤2得到的内区和外区分界点处流速计算各内区网格点流速；

步骤4：根据雷达测量的表面流速分布以及内区和外区分界点处流速，计算各外区网格点流速；

步骤5：将反演的各内区网格点流速、内区和外区分界点处流速、外区网格点流速以及雷达测量的表面流速相组合得到截面流速分布。

进一步地，对河流截面进行网格化分割具体包括：

根据微波雷达测量的表面流速分布对应设置N根垂线，其中，N为微波雷达测量的表面流速点数，结合设定的垂向距离分辨率对各垂线进行垂向等分，使河流截面被划分为多个离散的网格点，按各网格点距河床的垂直距离划分为内区和外区，其中，靠近河床的网格点被划分为内区，远离河床的网格点被划分为外区。

进一步地，内区和外区的划分满足如下关系：

式中，z为网格点距河床的垂直距离，z₁为内区和外区分界点处的网格点距河床的垂直距离，h为网格点所在垂线的水深。

进一步地，步骤2具体包括：

建立流向方向上的RANS方程表示如下：

式中，v为液体的运动粘滞系数，u为流向方向的流速，y和z分别表示侧向方向和垂向方向，ε为涡流粘度，g为重力加速度，S_f表示河床坡降；

将上述RANS方程简化为：

式中，ε_outer为外区涡流粘度，u_outer表示外区流速，z为网格点距河床的垂直距离；

将RANS方程沿垂向方向从外区的下边界z₁积分到外区的上边界z₂得到：

式中，

为剪切流速，u_z2和u_z1分别为上边界流速和下边界流速；

根据Prandtl混合长度理论，内区涡流粘度ε_inner为：

ε_inner＝κzu_*；

式中，κ为karman常数，z为内区网格点距河床的垂直距离；

ε_outer和ε_inner在内区和外区的分界点处即下边界处相同：

即得到外区涡流粘度为：ε_outer＝κz₁u_*；

将外区涡流粘度代入到上述将RANS方程中得到内区和外区分界点处流速：

式中，κ为karman常数，h为网格点所在垂线的水深，

为剪切流速，g表示重力加速度，s_f表示河床坡降，u_z1表示内区和外区分界点处的流速，u_z2为外区上边界处的流速，即为该网格点的表面流速，从雷达测量的表面流速分布中获得。

进一步地，步骤3中根据步骤2得到的内区和外区分界点处流速使用对数亏损律计算内区网格点流速，内区网格点流速计算公式如下：

式中，u_inner表示距河床垂直距离为z_inner的内区网格点流速，z₁为内区和外区分界点处的网格点距河床的垂直距离，u_z1表示内区和外区分界点处流速。

进一步地，步骤4中具体包括以下子步骤：

步骤4.1、建立在外区流向方向上的RANS方程，并采用有限差分法将其离散成各根垂线上的各个外区网格点流速对应的线性方程；

步骤4.2、将每个外区内网格点流速转换成步骤4-1中的线性方程，所有线性方程形成一个线性方程组，对于线性方程组用矩阵A、U、C、B分别对应表示线性方程组中的系数、外区网格点流速、各线性方程中包含的外区边界流速项以及右侧项，线性方程组表示为：

AU+C＝B；

步骤4-3：根据各垂线的表面流速、内区和外区分界点处流速计算矩阵C中各元素的值，进而对线性方程组进行求解，获得所有外区网格点流速，得到的外区网格点流速的计算公式如下：

U＝(B-C)A^-1。

与现有技术相比，本发明的有益效果为：

1.通过分析截面内区流速分布和外区流速分布的差异和联系，分别采用对数亏损律和RANS方程反演内区流速分布和外区流速分布，从而获取河流截面的流速分布，从而解决了微波雷达无法直接测量河流表面以下的流速，为后续的流量监测奠定了良好的基础，对河流水文监测、环境保护具有重要的应用价值；

2.本发明利用微波雷达测量的表面流速，在不接触水体的情况下完成了截面流速分布的反演，显著减少了流速分布测量时所需的工作量，对推进水文测验现代化具有重要意义；

3.本发明使用微波雷达作为探测工具，不受雨雾影响，能实现流速的全天候监测，适合在高流速、大流量情况下工作；此外，本发明不仅适用于岸基微波雷达，还可以拓展到其他的非接触式测流方法中，如大粒子图像测速法，空时图像测速法；

4.本发明计算简单快捷，能够满足雷达实时工作的要求，有助于掌握实时、准确的截面流速分布信息。

附图说明

图1为本发明实施例截面流速分布的反演方法的流程图；

图2为本发明实施例有限差分法原理图；

图3为本发明实施例反演的截面流速分布；

图4为本发明实施例反演的截面流速分布和声学仪器测量结果间的绝对误差分布；

图5为本发明实施例反演的截面流速分布和声学仪器测量结果的散点图。

具体实施方式

下面将结合本发明实施例对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动的前提下所获得的所有其他实施例，都属于本发明保护的范围。

需要说明的是，在不冲突的情况下，本发明中的实施例及实施例中的特征可以相互组合。

下面结合具体实施例对本发明作进一步说明，但不作为本发明的限定。

如图1所示，本发明实施例提供一种基于河流表面流速的截面流速分布反演方法，其特征在于，包括如下步骤：

步骤1：对河流截面进行网格化分割，将各网格点划分为内区和外区；

根据微波雷达测量的表面流速分布设置对应的N根垂线，其中N为微波雷达测量的表面流速点数，在本实施例中，结合设定的垂向距离分辨率0.02m对各垂线进行垂向等分，使河流截面被划分为多个离散的网格点。按各网格点距河床的垂直距离划分为内区和外区，内区和外区的划分满足如下关系：

式中，z为网格点距河床的垂直距离，h为网格点所在垂线的水深；

将所有网格点按照上述方法划分为内区和外区。当然在别的实施例中，也可以选择0.15h等其他垂直距离作为内区和外区的分界。

步骤2：计算外区涡流粘度并根据计算得到的外区涡流粘度结合雷达测量的表面流速分布得到内区和外区分界点处流速；

根据Boussinesq假设，对于平稳，均匀，充分成长的河流，建立流向方向上的RANS方程表示如下：

式中，v为液体的运动粘滞系数，u为流向方向的流速，y和z分别表示侧向方向和垂向方向，ε为涡流粘度，g为重力加速度，S_f表示河床坡降；

在河流截面的中心部分，u在z方向的变化远大于y方向上的变化，

和

可以忽略，在外区，紊流应力占主导，粘性应力项

和

可以忽略，RANS可以简化为：

式中，ε_outer为外区涡流粘度，u_outer表示外区流速；

将RANS方程沿垂向方向从z(z≥0.2h)积分到水深h可得：

在河流表面，

RANS方程可以表示为：

将RANS方程沿垂向方向从外区的下边界z₁＝0.2h积分到外区的上边界z₂＝h可得：

式中，

为剪切流速，u_s为表面流速，u_0.2h为内区和外区分界点处流速；

根据Prandtl混合长度理论，内区涡流粘度ε_inner为：

ε_inner＝κzu_*，0＜z≤0.2h；

式中，κ为karman常数，z为内区网格点距河床的垂直距离；

由于ε_outer和ε_inner在内区和外区的分界点z＝0.2h处相同，故ε_outer＝ε_inner|_z＝0.2h；

因此，在任何一根垂线处，将内区涡流粘度代入上述公式即可得单根垂线的外区涡流粘度以及内区和外区分界点处流速：

ε_{outer_i}＝0.2κh_iu_*i

式中，ε_{outer_i}为第i根垂线的外区涡流粘度，κ为karman常数，h_i为第i根垂线的水深，

为第i根垂线的剪切流速，g表示重力加速度，s_f表示河床坡降，u_0.2hi表示第i根垂线内区和外区分界点处的流速，u_si为第i根垂线的表面流速。

步骤3：根据步骤2得到的内区和外区分界点处流速计算内区网格点流速；

在该步骤中，根据步骤2得到的内区和外区分界点处流速使用对数亏损律计算各垂线的内区网格点流速；其中，单根垂线的内区网格点流速计算公式如下：

式中，i为垂线号，u_{inner_i}表示距河床垂直距离为z_{inner_i}处的内区网格点流速；

对N根垂线重复该步骤，直至完成所有垂线的内区网格点流速计算。

步骤4：根据雷达测量的表面流速分布以及内区和外区分界点处流速，计算各外区网格点流速；

步骤4.1，建立在外区流向方向上的RANS方程，可以表示为：

上式中与步骤2相同的符号的含义也相同；根据图2所示的有限差分法进行离散，

其中，i为垂线号，j为第i根垂线内外区的网格点号，Δy为相邻两根垂线间的距离，Δz为设定的垂向距离分辨率。第i根垂线内第j个外区网格点流速对应的RANS方程可以转换为如下式的线性方程：

步骤4-2：将每个外区内网格点流速转换成步骤4-1中的线性方程，所有线性方程形成一个线性方程组；假设截面内包含M个外区网格点，对外区内每个网格点设置一个从步骤4-1中得到的线性方程，M个线性方程组成的线性方程组可以表示为：

式中，a_ij为第i个线性方程中u_j的系数，u_j表示第j个外区网格点流速，c_i是第i个线性方程中包含的外区边界流速项，

为线性方程的右侧项；

用矩阵A、U、C、B分别对应表示线性方程组中的系数、外区网格点流速、各线性方程中包含的外区边界流速项以及右侧项，得到：

线性方程组可以表示为：

AU+C＝B；

步骤4-3：根据各垂线的表面流速、内区和外区分界点处流速计算矩阵C中各元素的值，进而对线性方程组进行求解，获得所有外区网格点流速，得到的外区网格点流速的计算公式如下：

U＝(B-C)A^-1。

步骤5：将反演的内区流速分布、内区和外区分界点处流速、外区流速分布、以及雷达测量的表面流速相组合得到截面流速分布。图3为应用本实施例方法反演的截面流速分布结果。从图3中可以看出，较大的流速主要集中在河流中心区以及接近河流表面的位置，较小的流速集中在侧壁区以及靠近河床的位置。

为说明本发明的有益效果，本申请人将将本实施例反演的截面流速分布与声学仪器测量结果进行比较，见图4和图5，从图4中可以看出，采用本实施例的方法反演的整个河流截面，80％的流速反演误差都低于0.1m/s。图5为应用本实施例反演的截面流速分布和声学仪器测量结果的散点图，两者间的相关系数为0.9，上述结果说明本发明反演的截面流速分布具有较高的准确性。

以上仅为本发明较佳的实施例，并非因此限制本发明的实施方式及保护范围，对于本领域技术人员而言，应当能够意识到凡运用本发明说明书内容所作出的等同替换和显而易见的变化所得到的方案，均应当包含在本发明的保护范围内。

Claims (6)

1.一种基于河流表面流速的截面流速分布反演方法，其特征在于，包括如下步骤：

步骤1：对河流截面进行网格化分割，将各网格点分别划分为内区和外区；

步骤2：计算外区涡流粘度并根据计算得到的外区涡流粘度结合雷达测量的表面流速分布得到内区和外区分界点处流速；

步骤3：根据步骤2得到的内区和外区分界点处流速计算各内区网格点流速；

步骤4：根据雷达测量的表面流速分布以及内区和外区分界点处流速，计算各外区网格点流速；

步骤5：将反演的各内区网格点流速、内区和外区分界点处流速、外区网格点流速以及雷达测量的表面流速相组合得到截面流速分布。

2.根据权利要求1所述的基于河流表面流速的截面流速分布反演方法，其特征在于，步骤1中，对河流截面进行网格化分割具体包括：

根据微波雷达测量的表面流速分布对应设置N根垂线，其中，N为微波雷达测量的表面流速点数，结合设定的垂向距离分辨率对各垂线进行垂向等分，使河流截面被划分为多个离散的网格点，按各网格点距河床的垂直距离划分为内区和外区，其中，靠近河床的网格点被划分为内区，远离河床的网格点被划分为外区。

3.根据权利要求2所述的基于河流表面流速的截面流速分布反演方法，其特征在于，内区和外区的划分满足如下关系：

式中，z为网格点距河床的垂直距离，z₁为内区和外区分界点处的网格点距河床的垂直距离，h为网格点所在垂线的水深。

4.根据权利要求1所述的基于河流表面流速的截面流速分布反演方法，其特征在于，步骤2具体包括：

建立流向方向上的RANS方程表示如下：

式中，v为液体的运动粘滞系数，u为流向方向的流速，y和z分别表示侧向方向和垂向方向，ε为涡流粘度，g为重力加速度，S_f表示河床坡降；

将上述RANS方程简化为：

式中，ε_outer为外区涡流粘度，u_outer表示外区流速，z为网格点距河床的垂直距离；

将RANS方程沿垂向方向从外区的下边界z₁积分到外区的上边界z₂得到：

式中，

为剪切流速，u_z2和u_z1分别为上边界流速和下边界流速；

根据Prandtl混合长度理论，内区涡流粘度ε_inner为：

ε_inner＝κzu_*；

式中，κ为karman常数，z为内区网格点距河床的垂直距离；

ε_outer和ε_inner在内区和外区的分界点处即下边界处相同：

即得到外区涡流粘度为：_outer＝κz₁u_*；

将外区涡流粘度代入到上述将RANS方程中得到内区和外区分界点处流速：

式中，为karman常数，h为网格点所在垂线的水深，

为剪切流速，g表示重力加速度，s_f表示河床坡降，u_z1表示内区和外区分界点处的流速，u_z2为外区上边界处的流速，即为该网格点的表面流速，从雷达测量的表面流速分布中获得。

5.根据权利要求1所述的基于河流表面流速的截面流速分布反演方法，其特征在于，步骤3中根据步骤2得到的内区和外区分界点处流速使用对数亏损律计算内区网格点流速，内区网格点流速计算公式如下：

式中，u_inner表示距河床垂直距离为z_inner的内区网格点流速，z₁为内区和外区分界点处的网格点距河床的垂直距离，u_z1表示内区和外区分界点处流速。

6.根据权利要求1所述的基于河流表面流速的截面流速分布反演方法，其特征在于，步骤4中具体包括以下子步骤：

步骤4.1、建立在外区流向方向上的RANS方程，并采用有限差分法将其离散成各根垂线上的各个外区网格点流速对应的线性方程；

步骤4.2、将每个外区内网格点流速转换成步骤4-1中的线性方程，所有线性方程形成一个线性方程组，对于线性方程组用矩阵A、U、C、B分别对应表示线性方程组中的系数、外区网格点流速、各线性方程中包含的外区边界流速项以及右侧项，线性方程组表示为：

AU+C＝B；

步骤4-3：根据各垂线的表面流速、内区和外区分界点处流速计算矩阵C中各元素的值，进而对线性方程组进行求解，获得所有外区网格点流速，得到的外区网格点流速的计算公式如下：

U＝(B-C)A^-1。

Images