具体实施方式

下面结合附图和具体实施例对本发明作进一步的详细描述。

一种天然河道推移质输沙率的估算方法，它包括以下步骤：

步骤1)：获取边界条件：利用声学多普勒流速剖面仪测量河道断面地形及流场分布，采样得到泥沙的粒径级配及其沿河宽方向的变化；

步骤2)：划分子断面：将整个河道断面沿河宽方向分为若干子断面；

步骤3)：计算各子断面的河床阻力系数：根据声学多普勒流速剖面仪的测量结果，明确各子断面的河床阻力系数；

步骤4)：计算各子断面的平均流速：

步骤5)：计算各子断面的沙波波高、运动速度：根据各子断面的水深、流速、泥沙粒径级配，确定推移质泥沙的中值粒径，并计算各子断面沙波的波高、运动速度；

步骤6)：计算整个断面的推移质输沙率：根据沙波的波高、运动速度确定各子断面的推移质输沙率，进而得到整个断面的推移质输沙率，即天然河道推移质输沙率。

本发明中，步骤4)计算各子断面的平均流速的方法有两种：

第一种是以流量Q作为计算输入条件：即针对河道断面的流量Q，通过水位与流量的关系得到河道水位Z，结合河道地形确定各子断面的水深、宽度，结合河床阻力系数得到各子断面的平均流速。

第二种是以河道水位Z为计算输入条件：即针对河道水位Z，通过水位与流量关系得到河道断面的流量Q，结合河道地形确定各子断面的水深、宽度，并结合河床阻力系数得到各子断面的平均流速。

在本发明中，只要完成一次从步骤1)到步骤3)的操作，后续针对任意流量Q或任意水位Z的输入，直接从步骤4)开始计算，即可确定推移质输沙率。

优选地，步骤1)中，通过底泥采样器采集泥沙，若样品以粒径小于50微米的粉沙或粘粒为主，则使用激光粒度仪测量泥沙粒径级配；若样品以粒径大于50微米的沙粒或砾石为主，则使用筛分法测量泥沙级配。

优选地，步骤2)中，根据河道地形、水流运动沿横向变化的剧烈程度，将整个河道断面沿河宽方向分为N个子断面，而且，任意子断面均应满足宽深比大于5。

优选地，步骤3)中，根据声学多普勒流速剖面仪关于流场、地形测量的结果，确定各子断面的过流量Q_0i、水深H_0i、宽度W_0i，结合曼宁公式，计算各子断面的河床阻力系数n_i如下；其中，角标i表示断面序号，i＝1、2…N；

其中，J₀为水力坡降，为未知量且不随子断面i发生变化，H_0i、W_0i、Q_0i分别为各子断面的水深、宽度、过流量；

从而，得到各子断面对应n_i值的比例，即，n₁：n₂：…n_i：…n_N。

优选地，步骤4)中，针对河道断面的流量Q，通过水位与流量关系得到相应的河道水位Z；将河道水位与河道地形高程结合，得到各子断面的水深H_i、宽度W_i；根据子断面的河床阻力系数n_i，通过如下公式(2)的曼宁公式，计算各子断面的流量分配比，即Q₁：Q₂：…Q_i…Q_N；

结合河道断面的流量Q，得到子断面的过流量Q_i和平均流速U_i＝Q_i/(H_i·W_i)。

优选地，步骤5)中，在已知U_i、H_i的基础上，结合泥沙起动流速公式，确定推移质泥沙的上限粒径D_i ^上；

其中，γ_s为泥沙容重，γ为水的容重，g为重力加速度。

进一步地，结合Shields曲线，推求摩阻流速u_*i，即，

其中，为Shields数Θ，为沙粒雷诺数Re_*，ν为水的粘滞系数，f()即Shields曲线；通过改变u_*i进行试算，使得以Re_*为横坐标、以Θ为纵坐标的点落在Shields曲线上；此时的u_*i即公式(4)的计算值。

结合泥沙悬浮指标，确定推移质泥沙的下限粒径D_i ^下，即

其中，卡门常数κ＝0.4，ω_i为粒径D_i ^下泥沙的沉降速度，采用张瑞瑾公式进行求解，

粒径小于D_i ^上且大于D_i ^下的非均匀沙即推移质；根据泥沙级配曲线，可以得到推移质的粒径组成，进而计算推移质的中值粒径D_mi；

结合下式，计算该子断面沙波的波高Δ_i，

沙波的运动速度c_i可通过如下公式得到，

优选地，步骤6)中，根据沙波的波高Δ_i、运动速度c_i，计算各子断面的推移质单宽输沙率p_i，即

p_i＝0.5·(1-e)·ρ_s·△_i·c_i (9)

其中，e为泥沙孔隙率e＝0.45,ρ_s为泥沙的密度ρ_s＝2650kg/m³；

结合子断面宽度W_i，得到相应的推移质输沙率；将所有子断面的输沙率求和，得到整个河道断面的推移质输沙率P，即

优选地，每隔一段时间进行河道断面、泥沙级配测量，用以修正边界条件，保证持续获得准确的推移质输沙率估算结果。

实施例1

如图1所示，天然河道推移质输沙率的估算方法，该方法包括以下步骤：

步骤1：获取边界条件。针对长江湖口站下游的某断面，利用声学多普勒流速剖面仪(ADCP)测量断面的地形及流场分布。将相关设备安装在配套的单体船上，控制船速2m/s左右，沿河宽方向往返测量四次。测量过程中，保证ADCP罗盘的参数在纵摇±10°和横摇±5°的范围内，保证探头入水深度0.1m。测量结果如图2所示。沿河宽方向布置左、中、右三个泥沙采样点，采集泥沙并测量粒径级配。泥沙级配见图3。

步骤2：划分子断面。根据河道地形、水流运动等沿横向(即河宽方向)变化的剧烈程度，将河道断面沿横向分为左、中、右三个子断面，见图2。此处分别称为1、2、3号断面。

步骤3：计算各子断面的河床阻力系数。根据图2关于ADCP测量的断面流场，1、2、3号断面分别对应：流量Q₀₁＝1055、Q₀₂＝18325、Q₀₃＝2917m³/s，宽度W₀₁＝200、W₀₂＝550、W₀₃＝200m，水深H₀₁＝17.59、H₀₂＝27.99、H₀₃＝12.56m。

结合如下的曼宁公式，可知n₁：n₂：n₃＝1.28:0.28:0.37。

步骤4：计算各子断面的平均流速。针对河道断面的流量Q，本实施例选择Q＝18679m³/s，根据“水位-流量”关系，得到此时的河道水位Z＝13.472m。结合图2的断面地形高程，明确各子断面水深H_i、宽度W_i。根据子断面的河床阻力系数n_i，通过如下的曼宁公式，计算各子断面的流量分配比，即Q₁：Q₂：Q₃＝1：20.21：2.93。

结合整个断面的过流量Q＝18679m³/s，得到子断面的过流量Q₁＝774、Q₂＝15638、Q₃＝2267m³/s。

子断面平均流速U_i＝Q_i/(H_i·W_i)，分别对应U₁＝0.254、U₂＝1.022、U₃＝0.924m/s。

上述结果见表1。

表1 Q＝18679m³/s条件下，各子断面的流动参数

步骤5：计算各子断面的推移质中值粒径以及沙波的波高、运动速度。在U_i、H_i已知的基础上，根据沙莫夫公式，计算推移质泥沙的上限粒径D_i ^上，

其中，γ_s为泥沙容重，γ为水的容重，g为重力加速度。

进一步地，结合Shields曲线(见图4)，推求相应的摩阻流速u_*i，即，

其中，为Shields数Θ，为沙粒雷诺数Re_*，ν为水的粘滞系数，f()即Shields曲线。通过改变u_*i进行试算，使得以Re_*为横坐标、以Θ为纵坐标的点落在Shields曲线上。此时的u_*i即公式(4)的计算值。

结合泥沙悬浮指标，确定推移质泥沙的下限粒径D_i ^下，即

其中，卡门常数κ＝0.4。ω_i为粒径D_i ^下泥沙的沉降速度。张瑞瑾公式建立了ω_i与D_i ^下之间的关系

相关计算结果见表2。

表2 Q＝18679m³/s条件下，各子断面的泥沙运动参数

由表2可知，1号断面对应D₁ ^上<D₁ ^下，说明不存在推移质运动，即所有发生运动的泥沙均以悬移质的方式向下游输移。2号和3号断面的推移质粒径范围(理论值)分别是0.455～2.102mm和0.448～2.025mm。

根据图3关于河床表面泥沙的粒径级配，2号和3号断面对应的最大粒径分别为0.606、0.742mm。因此，2号和3号断面的推移质粒径范围为(实际值)0.455～0.606mm和0.448～0.742mm。进一步结合图3关于上述粒径范围的泥沙级配，可知2号和3号断面的推移质中值粒径分别为D_m2＝0.501、0.514mm。

结合下式，分别计算各子断面沙波的波高Δ_i、运动速度c_i，

步骤6：计算整个断面的推移质输沙率与粒径级配。根据Δ_i、c_i，结合下式计算各子断面的推移质单宽输沙率p_i。相关结果见表3。

p_i＝0.5·(1-e)·ρ_s·△_i·c_i (9)

其中，e为泥沙孔隙率e＝0.45,ρ_s为泥沙的密度ρ_s＝2650kg/m³。

表3 Q＝18679m³/s条件下，各子断面的沙波运动参数

相应的，各子断面的输沙率P_i＝p_i·W_i,见表3。按下式计算整体断面的推移质输沙率P，

本实施例对应P＝51.553+10.658＝62.211kg/s。

实施例2(技术效果对比实施例)

本发明考虑了水流运动、泥沙粒径沿河宽方向的变化。倘若不考虑上述因素并将整个河道断面作为整体计算，得到的推移质输沙率为42.697kg/s。而实施例1的计算值为62.211kg/s。针对实施例1的河段，在Z＝13.524m条件下，利用推移质取样仪进行现场测量，得到的推移质输沙率为68.373kg/s。在现场测量过程中，沿河宽方向间隔80m布置一个采样点，将样品烘干、称重，除以取样时间，得到取样点位置的单宽输沙率。将各点的单宽输沙率取均值，再乘以河宽，得到整个断面的推移质输沙率。

对比上述结果，在河道水位与实施例1水位相近的前提下，本发明方法考虑了水沙运动沿河宽方向的变化，计算值与实测值的相对误差为9.01％。而不考虑水沙运动及沿河宽方向变化的方法，相应的误差为37.55％。