基于滑动分析法的黄河下游造床流量变化研究

(1.黄河水利委员会黄河水利科学研究院水利部黄河泥沙重点实验室，河南郑州，450003；2.中国科学院地理科学与资源研究所，北京，100101；3.武汉大学水资源与水电工程科学国家重点实验室，湖北武汉，430072)

泥沙淤积是黄河下游洪水危害的本质原因。黄河下游河床宽浅、散乱，主流摆动不定，河势变化剧烈。泥沙不断淤积及复杂多变的河势使得防洪负担日益加剧。黄河下游河道冲淤演变特性很大程度上取决于来水来沙条件，但是天然河道的水沙过程随时间发生变化，因此，在进行河床演变分析和河道整治时，往往需要确定与其多年流量过程综合造床作用相当的特征流量，即造床流量。目前确定造床流量的理论并不成熟，常用的方法有平滩流量法[1]、输沙率法[2]、输沙能力法等。其中，输沙能力法在工程实践中应用最多，最经典的是马卡维也夫提出的计算方法[3]，通常也称为马氏造床流量方法。该方法通常会得到2个造床流量，其中，第一造床流量相当于多年平均最大洪水流量，第二造床流量稍大于多年平均流量。张红武等[4]采用马氏造床流量方法分析了黄河“92·8”高含沙洪水异常现象成因，并针对多沙河流造床特点，引入水流挟沙能力修正了该计算方法[5]。随后，韩其为等[6-9]定义塑造河床纵剖面的流量为第一造床流量，塑造横断面的流量为第二造床流量，并将其应用于黄河下游调水调沙对造床流量影响研究中。江恩惠等[10]针对黄河下游来水来沙人为调控特征与变化趋势，提出了洪、中、枯兼容的“三级流路”整治模式，即整治导线按中水(相当于马氏第一造床流量)流路设计。基于黄河下游各水文站1950—2003年实测资料，陈建国等[11]系统分析了50 多年来黄河下游平滩流量和造床流量变化过程，提出增大河道的造床流量是治理河道萎缩的重要措施之一。陈绪坚等[12]通过计算黄河下游第一和第二造床流量，分析了黄河下游河槽的萎缩过程。胡春宏等[13]研究了1950—2003年黄河下游平滩流量、造床流量与洪水特征值的响应关系，阐明了黄河水沙过程调控与下游河道中水河槽塑造内在联系。孙东坡等[14]提出了针对黄河下游水沙变异特点的造床流量确定方法，指出维持黄河下游河道功能的河床塑造必须考虑适当的中常洪水，而不是长时间小流量。还有一些学者利用水文站长期观测数据探讨了造床流量与河床启冲流速[15]、平滩河宽[16]和洪峰流量重现期[17]之间的关系，以及通过物理试验研究了砾石河床[18]、砂质河床[19]演变规律与不同流量级冲刷之间的内在联系。造床流量对于分析天然河道水沙过程及河床演变特性十分重要。目前，大部分造床流量确定方法均以马氏造床流量方法为基础，流量级的人为界定往往带来误差。近10年来，随着黄河水沙调控工程体系的运用，黄河下游水沙关系发生了很大变化。因此，本文作者基于黄河下游花园口水文站1950—2015年实测数据，通过滑动分析法对马氏造床能力表达式进行修正，避免人为划分流量级导致的计算误差，研究近年来黄河下游造床流量变化规律，旨在探明人为调控作用下黄河下游水沙过程与河床的适应性。

1 黄河下游造床流量确定方法

1.1 马氏造床流量

经典马氏造床流量认为某个流量的造床作用既与该流量的输沙能力有关，也与该流量所持续时间有关。与流量Q输沙能力有关的造床作用与Q的m次方及比降J的乘积成正比，其中m为双对数坐标上实测断面日均输沙率Gs与相应日均流量Q关系曲线的斜率，与流量Q所持续时间有关的造床作用可用该流量出现的频率P来表示[1,3]，P的求解方式为将河段断面随时间变化的流量过程分为若干流量级，统计各流量级在各流量级内出现的频率。马氏造床能力的表达式为QmJP。诸多学者[20-27]在运用马氏造床流量计算方法时发现，人为划分流量级差会给造床流量的确定带来很大影响。

1.2 滑动分析方法

滑动分析方法通过顺序逐期增减新旧数据求取移动平均值或累加值，借以消除偶然变动和人为主观因素的影响，近年来该方法已被广泛地应用于黄河来沙来水量的计算分析[28-29]。本文基于滑动分析方法原理，对马氏造床流量进行修正，以消除人为划分流量级(参数P)的不良影响，具体计算步骤如下。

1)定义一系列流量级为

式中：i为流量级的序号；n为流量级的总数；和分别为第i级流量的上、下界；ΔQ1为每组流量级长度；ΔQ2为流量级滑动距离；具体定义如图1所示，其中和分别为第i+1 级流量的上、下界。

图1 流量级、ΔQ1和ΔQ2定义示意图Fig.1 Sketch map of definitions of discharge stages,ΔQ1 and ΔQ2

2)定义第1 流量级的初始值为所分析流量数据中最小值Qmin，则该系列流量级可表述为

3)绘制Qi-φi图，积分求取其围成的面积S=∫Qiφidi；对于某个水文站某段时间内水沙过程，S为常数；为了便于对比不同时间段流量的造床能力，对进行归一化处理，即Φi=φi/S；Qi-Φi变化规律与Qi-φi变化规律是一致的，两者横坐标相同，区别仅在于纵坐标的值不同(Φi纵坐标值是φi纵坐标值的1/S倍)；绘制Qi-Φi图，在图中查出Φi的最大值和次大峰值，分别为第一造床能力和第二造床能力。2个流量中，流量较大值为第一造床流量，较小值为第二造床流量。

2 黄河下游水沙关系研究

2.1 参数敏感分析

基于黄河下游花园口水文站1950—2015年实测数据(地理位置如图2所示)，采用滑动分析法修正马氏造床流量计算方法，分析黄河下游造床流量变化规律。首先在双对数坐标系中点绘花园口1950—2015年日均流量Q和相应日均输沙率Gs，如图3所示，对数据进行拟合，得到Q和Gs关系为Gs=5e-5Q2.1(e为自然底数；R2为函数拟合可决系数，R2=0.9，其越接近1，表明函数拟合程度越好)。

图2 黄河部分水文站和水库地理位置Fig.2 Geographic locations of some hydrological stations and reservoirs in the Yellow River

图3 花园口1950—2015年流量Q与输沙率Gs关系Fig.3 Relationship between discharge Q and sediment transport rate Gs based on data from the Huayuankou hydrological station from 1950 to 2015

马氏造床流量方法中流量级长度取值直接影响其组数，根据文献[20-27]，马氏方法选取流量级组数从十几组到近百组不等。根据实测资料，花园口水文站1986—2015年实测日均流量中最小值为7.58 m3·s-1，最大值为7 270 m3·s-1。当流量级取100 m3·s-1时，流量级组数为73 个；当其取800 m3·s-1时，流量级组数为10 个。因此，选取100和800 m3·s-1以分别代表流量级组数较多和较少时情况，进而对比分析2种流量级划分对造床流量的影响，计算结果如图4所示，其中比降J取花园口实测比降资料的平均值0.000 2。由图4(a)可知：当流量级为100 m3·s-1时，计算得到第一造床流量(粉色标记)和第二造床流量(红色标记)分别为1 260和860 m3·s-1。由图4(b)可知：当流量级为800 m3·s-1时，第一造床流量和第二造床流量分别为5 050和1 060 m3·s-1。可以看出，人为划分流量级对马氏造床流量计算结果影响较大；若流量级取得小，则会造成某一小流量级因出现频率过高而掩盖真实值；若流量级取得大，则计算精度不高。

图4 花园口1986—2015年不同流量级下马氏造床流量分析Fig.4 Analysis of dominant discharge under different man-made discharge stages using Makayev method based on data from the Huayuankou hydrological station from 1986 to 2015

分析滑动分析法中ΔQ1和ΔQ2这2个参数对计算结果的影响。以花园口水文站1986—2015年实测数据为研究对象，首先选定ΔQ2=1 m3·s-1(最小的正整数流量)，这样，流量级滑动时既可确保流量数据得到充分利用，又不致于计算量太大，然后分别计算当ΔQ1为200，300，400，500和600 m3·s-1时Qi-Φi关系曲线，如图5所示。从图5可以看出：随着ΔQ1增大，Φi整体分布略降低，但规律上几乎没有变化，第一和第二造床流量分布没有变化。

图5 花园口1986—2015年不同ΔQ1下Qi-Φi关系曲线Fig.5 Relationship curves between Qi and Φi with different ΔQ1 based on data from the Huayuankou hydrological station from 1986 to 2015

另一方面，选定ΔQ1=500 m3·s-1(由图5可知ΔQ1在200～600 m3·s-1范围内对造床流量的确定没有影响)，设定ΔQ2为1，50，100，300和600 m3·s-1，以分析ΔQ2对Qi-Φi关系曲线的影响，如图6所示。由图6可知：当ΔQ2增大时，由于数据点减少，Qi-Φi关系曲线变得越来越不光滑(特别是当ΔQ2大于ΔQ1时)，但整体分布规律几乎没有变化，第一和第二造床流量分布受其影响不大。因此，建议在使用该方法时ΔQ2小于ΔQ1，这样可确保Qi-Φi关系曲线中有效数据的数目，以更加真实地反映曲线的分布规律。

2.2 不同代表年段花园口造床流量计算

根据图5和图6计算结果可知：ΔQ1在200～600 m3·s-1范围内、ΔQ2在1～600 m3·s-1范围内对造床流量的确定没有影响。考虑黄河造床流量计算时经常采用的马氏流量级组数以及黄河日均流量分布，计算时选定ΔQ1为500 m3·s-1，选定ΔQ2为1 m3·s-1。近年黄河水沙变化较典型年份有1986年和1999年：1986年，黄河上游龙羊峡水库和刘家峡水库得到联合运用；1999年，黄河下游小浪底水库得到运用。本文将研究时段分为1950—1985年、1986—1999年、2000—2015年以及1950—1999年、1986—2015年，相应的Qi-Φi关系曲线如图7所示。在图7中分别找到Φi的最大峰值和次大峰值对应的特征流量，其中，较大流量为第一造床流量，较小流量为第二造床流量，如表1所示。

图6 花园口1986—2015年不同ΔQ2下Qi-Φi关系曲线Fig.6 Relationship curve between Qi and Φi with different ΔQ2 based on data from the Huayuankou hydrological station from 1986 to 2015

由图7和表1可知：1950—1985年花园口第一造床流量的造床能力大于第二造床流量的造床能力；1986—1999年水沙条件发生变化，花园口第一和第二造床流量均显著减小，且第一造床流量的造床能力小于第二造床流量的造床能力。可以看出，1986年以后黄河下游水沙关系发生了显著变化，第一造床流量明显降低，意味着洪水频次和洪峰流量显著减少；第二造床流量也显著降低，表明来水流量减少，这是造成下游河道淤积萎缩的主要原因。

另一方面，相比于1986—1999年，花园口2000—2015年第一造床流量降低，第二造床流量增加，且第一造床流量的造床能力小于第二造床流量的造床能力。通过分析可知，自1999年小浪底水库运用及调水调沙以来，下游河道发生持续冲刷，主槽过流输沙能力提高，第二造床流量增加。为了更加全面分析目前黄河下游河道造床流量分布现状，本文还计算了黄河下游高村水文站(距上游的花园口水文站近189 km，地理位置如图2所示)和泺口水文站(距上游的花园口水文站近469 km，地理位置图2所示)1986—2015年第一和第二造床流量。计算时，高村水文站比降为实测比降的平均值即0.000 15，拟合Gs与Q得到m为2.0；泺口水文站比降为实测比降的平均值即0.000 1，拟合Gs与Q得到m为1.8，计算结果如图8和表1所示。可以看出，两者造床能力现状均为第二造床流量大于第一造床流量。

2.3 滑动分析法与马氏方法对比

采用滑动分析法(ΔQ1和ΔQ2分别为图5和图6所示范围)和马氏方法(流量级为100和800 m3·s-1)对花园口1950—1985年、1986—2015年造床流量进行计算，所得造床流量如表2所示。由表2可知：马氏方法中流量级对造床流量影响很大，流量级取值不同，计算得到的造床流量确定值差异很大；而当ΔQ1在200～600 m3·s-1范围内、ΔQ2在1～600 m3·s-1范围内时，采用滑动分析法均可准确计算出第一和第二造床流量。

图7 花园口不同研究时段Qi-Φi关系曲线Fig.7 Relationship curves between Qi and Φi based on data from the Huayuankou hydrological station under different research periods

图8 高村和泺口水文站1986—2015年Qi-Φi关系曲线Fig.8 Relationship curves between Qi and Φi based on data from the Gaocun and Luokou hydrological stations from 1986 to 2015

表1 不同研究时段花园口第一、第二造床流量Table1 The first and second dominant discharge based on data from the Huayuankou hydrological station under different research periods

表2 基于滑动分析法和马氏方法计算得到的花园口造床流量Table2 Dominant discharge based on data from the Huayuankou hydrological station according to moving analysis method and the Makayev method

3 结论

1)滑动分析法的引入有效避免了马氏造床流量计算方法中人为划分流量级差带来的误差。当流量级长度参数ΔQ1在200～600 m3·s-1范围内、流量级滑动距离参数ΔQ2在1～600 m3·s-1范围内时，改进后的造床流量计算方法对第一和第二造床流量的确定没有影响，均可准确计算出造床流量。

2)1986—1999年，黄河下游水沙条件发生了显著的变化，第一造床流量和第二造床流量显著减少，下游河道淤积萎缩严重；自1999年小浪底水库运用及调水调沙以来，第二造床流量增加，下游河道主槽过流输沙能力提高。

3)目前黄河下游造床能力已由第一造床流量大于第二造床流量转变为第二造床流量大于第一造床流量。