3.4 黄河下游平滩流量对水沙条件变化的响应与模拟

本节选择花园口、高村、孙口、艾山、利津5个测站（见图3.2）作为研究对象，研究黄河下游平滩流量对来水来沙条件的响应关系，建立平滩流量的计算方法，对黄河下游主要测站平滩流量随水沙条件的变化过程进行模拟。花园口站、高村站、艾山站和利津站水沙资料为1956—2002年系列，平滩流量资料为1960—2002年；孙口站水沙资料为1964—2002年系列，平滩流量资料为1968—2002年。上述测站的资料统计见表3.4。

3.4.1 平滩流量与来水来沙条件的关系

黄河下游属于冲积性河流，河床形态与上游来水来沙关系密切。河床通过不断的泥沙冲淤来适应不断变化的来水来沙条件，结果，不同的水沙条件会塑造出不同的河槽形态，对应不同的平滩流量。考虑到河道的调整变化主要发生在流量较大的汛期，图3.25点绘了花园口站1950—2002年的平滩流量与汛期平均来流量的关系。由图3.25（a）可以看到，平滩流量随汛期平均流量的增大而增大，两者之间成正比关系，只不过相关程度不高，其回归关系为：

图3.25（b）为平滩流量与4年滑动平均汛期来流量的关系，其回归关系为：

图3.25 花园口站平滩流量与汛期平均流量的关系

比较式（3.44）和式（3.43）可以看到，平滩流量与4年滑动平均汛期流量之间的关系，较其与当年汛期平均流量之间的相关程度有较大提高，决定系数R²由0.41提高到0.62。

前述图1.11（a）显示了平滩流量和当年汛期平均流量的变化过程，可以看到，二者的总体上升和下降趋势一致，但平滩流量无法完全跟随当年汛期平均流量的波动变化。当采用4年滑动平均汛期流量时，两者的变化趋势变得比较一致［见图1.11（b）］，说明采用前期若干年的汛期滑动平均流量，基本上反映了来水条件的累积作用，同时还说明了这样一个事实，即河床的调整需要一定的时间，通常要滞后于水流的变化，目前的河槽形态是河道多年累积冲淤演变的结果，不仅与当年的水沙条件有关，还与前期的河床边界条件有关。为了说明平滩流量对前期河床形态的依赖性。图3.26点绘了当年平滩流量与往年（即上一年）平滩流量的关系，两者之间的相关程度很高（R²=0.77）。考虑到往年平滩流量本身也是过去水沙条件塑造的结果，以滑动平均流量来反映前期水流条件的影响也就在情理之中，图1.11（b）的结果说明了这种处理方法的合理性。

图3.26 花园口站当年汛期平滩流量与往年（即上一年）汛期平滩流量的关系

图3.27给出了花园口平滩流量与滑动平均汛期流量的相关程度随包括年数的变化，可以看到，一开始决定系数随着包括年数的增加而增加，并在包括年数为4左右达到最大；之后包括更多年数的相关程度反而降低。河床调整的这种滞后响应，在其他河床演变现象中也普遍存在（吴保生和邓玥，2005；Zheng等，2014a）。

关于平滩流量对来沙条件的响应，可以用来沙系数ξ（=S/Q）来反映。以往研究表明，黄河下游河道的冲淤变化与该系数的关系密切，来沙系数小时冲刷，平滩流量加大；来沙系数大时淤积，平滩流量减小。图3.28点绘了平滩流量与当年汛期平均来沙系数和4年滑动平均来沙系数的关系，依据图中数据得到的相关关系分别为：

需要说明的是，为了应用的方便，式（3.45）中的来沙系数是由汛期平均含沙量与汛期平均流量计算的，即，也可以先由日均含沙量和日均流量计算得到每天的来沙系数，由此得到汛期来沙系数的平均值，结果差别不大。此外，图3.27中给出的决定系数随滑动平均来沙系数包括年数变化的结果表明，当包括年数为5年时，决定系数最大。说明平滩流量的调整对来水量和来沙系数的响应方式有所不同，对来沙系数变化的响应与调整所需要的时间更长，但为了与图3.25和式 （3.44）一致，图3.28和式 （3.46）给出了采用4年滑动平均来沙系数的结果。

图3.27 花园口站平滩流量与滑动平均汛期流量的相关程度随包括年数的变化

图3.28说明平滩流量随来沙系数的增大而减小，两者之间成反比关系，并且采用4年滑动平均汛期来沙系数能够在一定程度上反映前期水沙搭配的影响，决定系数R²由0.36提高到0.71，相关程度有较大提高。假设花园口站汛期来沙系数由0.02增加到0.04（增加一倍），根据式（3.46）可以得到，花园口站平滩流量将由5296m³/s减小到4243m³/s，减少1054m³/s，即减少20%。

前述图1.12（a）点绘了平滩流量和汛期来沙系数的历年变化过程。可以看到，当来沙系数连续上升时，平滩流量减少；当来沙系数连续下降时，平滩流量增大，但平滩流量没有当年汛期来沙系数的波动大，反映了平滩流量较来沙条件变化的滞后性。图1.12（b）显示当采用4年滑动平均来沙系数时，两者的对应变化趋势大体相当，显示了水沙条件的累积作用。

图3.28 花园口站平滩流量与汛期来沙系数的关系

3.4.2 平滩流量的滑动平均计算方法

我们知道，河道的主槽断面形态和平滩流量是水沙条件综合作用的结果。对平滩流量Q_bf与滑动平均汛期流量和滑动平均汛期来沙系数的单变量相关分析已经表明，来水量连续的偏丰或偏枯，以及来沙量连续的偏少或偏多，均会造成河道主槽面积趋势性的增大或减小，采用基本上能够反映前期水沙条件的影响效果。此外，在滑动平均年数为4～5年时，决定系数达到最大值，说明在来水来沙条件发生变化时，河道主槽断面通过冲淤调整，达到新的平衡状态的时间约为4～5年，所以，以4年滑动平均汛期流量和4年滑动平均汛期来沙系数为基础建立Q_bf的计算公式，得到如下表达式：

根据1960—2002年实测平滩流量资料，回归得到相应各站的系数和指数见表3.5。式（3.47）不仅物理意义明确，反映了来水来沙条件决定河槽形态的河床演变学概念，而且结构简单，应用方便。对于给定的来水来沙条件，式（3.47）加上表3.5的系数，可以用来估算黄河下游平滩流量的变化。图3.29以花园口站和利津站为例，点绘了平滩流量与汛期水沙综合参数的关系，可以看到，平滩流量与的相关性较好。

图3.29 平滩流量与汛期水沙综合参数的关系

由式（3.47）和表3.5可以得到，平滩流量与汛期平均流量的0.18～0.41次方成正比，与汛期平均含沙量的0.11～0.32次方成反比。定性上讲，平滩流量随流量的增加而增加，随含沙量的增加而减小，这与过去的认识是一致的。各个测站系数和指数的不同，主要是由于各站来沙组成、河床边界条件及河型等的不同所致，其影响因素和变化规律有待深入研究。此外，式（3.47）还说明，汛期平均流量（或汛期水量）和汛期平均来沙系数（来沙量）的大小决定了平滩流量的大小，洪水过程的影响应该处于次要地位，其影响程度有待深入研究。

图3.30给出了花园口站和利津站计算与实测平滩流量的历年变化过程。可以看到，如果来水流量持续偏丰，来沙系数连续减少，平滩流量就会持续增大；相反，如果来流量持续偏枯，来沙系数连续增加，平滩流量就会持续降低。式（3.47）较好地反映了平滩流量长时期随水沙条件的持续上升或持续下降的过程。

图3.30 计算与实测平滩流量的历年变化过程

3.4.3 平滩流量滞后响应模型

（1）平滩流量滞后响应模型的建立

以平滩流量Q_b作为河床演变的特征变量，根据滞后响应模型的单步解析模式［式（2.8）］和多步递推模式［式（2.15）］，可以分别得到如下平滩流量的单步解析计算公式和多步递推计算公式：

式中：Q_b为t时刻的平滩流量；Q_b0为t=0时的Q_b值，即水沙条件变化前的初始平滩流量；Q_e为相对某一给定水沙条件下河床调整达到相对平衡时的平滩流量；β为平滩流量的变化速率，是一个与水流能量大小和河床边界可动性有关的参数；Δt为计算时间步长，取为一个运用年。

式（3.48）和式（3.49）都能够考虑前期来水来沙条件的累积作用和滞后影响，具有普遍的适用性，关键是如何针对具体情况给出平衡状态Q_e的表达式。从水沙特性决定河流特性的基本原理出发，对应平衡状态的平滩流量可以看作是来水来沙条件的函数，考虑到河道的调整变化主要发生在流量较大的汛期，Q_e可以表示为：

式中：Q_f为汛期平均流量，m³/s；ξ_f为汛期来沙系数，kg·s/m⁶；K、b、c分别为待定系数和指数。

将式（3.50）代入式（3.48）和式（3.49），分别得到：

以上两式便是冲积河流平滩流量的滞后影响模型。当已知初始平滩流量Q_b0时，可以根据本时段的来水来沙情况，用式（3.51）来推求时段末的的平滩流量；当初始平滩流量Q_b0未知时，根据前期若干时段的水沙条件，用式（3.52）来计算任何水沙系列条件下平滩流量的响应变化。需要注意的是，计算中实际包括了n+1年的水沙条件。

（2）平滩流量变化过程的模拟

选择黄河下游花园口、高村、孙口、艾山、利津5个测站，收集整理了1950—2003年的水沙量及1960年以来的平滩流量资料，用以检验式（3.51）和式（3.52）平滩流量滞后响应模型的合理性和适用性。作为代表，表3.6给出了花园口站的流量、来沙系数及平滩流量。

对于已知Q_b0的情形，根据表3.6给出的花园口站资料，取Δt=1年，得到式（3.51）中的系数和指数K、b、c和β的值分别为488.9、-0.25、0.18和0.315，即花园口站平滩流量的单步解析模式计算公式为：

或

上两式中，Q_b、Q_b0、Q_f的单位均为m³/s，ξ_f的单位为kg·s/m⁶。采用同样的方法，可以得到对于其余各站平滩流量的计算公式，表3.7给出了所有各站依据实测资料得到的系数和指数值。

由式（3.53）或（3.54）可以看到，初始平滩流量Q_bo的权重为73%（e^-0.315Δt=0.73），而本时段来水来沙的权重仅为27%（1-e^-0.315Δt=0.27），这一结果定量说明了前期河床边界条件对本时段平滩流量的作用是巨大的。

表3.7中的系数同样适用于滞后响应模型的递推关系式（3.52），可得花园口站平滩流量的多步递推模式计算公式为：

或

表3.7给出了各站计算与实测平滩流量之间的决定系数R²值，包括采用单步解析模式和多步递推模式的计算结果。作为例子，表3.6给出了采用式（3.54）和式（3.56）计算的花园口站历年平滩流量与实测值之间的相对误差，式（3.54）的平均相对误差为7.2%（相对误差绝对值的算术平均值），最大为-19.10%，式（3.56）的平均相对误差为7.9%，最大为19.45%。图3.31和图3.32还对花园口站和利津站的计算与实测平滩流量进行了比较，可以看到，计算与实测平滩流量吻合较好，虽然个别年份存在一定误差，但计算与实测平滩流量的总体趋势是基本一致的。

图3.31 计算与实测平滩流量的比较[式（3.5 4）]

图3.32 计算与实测平滩流量的比较[式（3.5 6），n=4]

由表3.7、图3.31和图3.32的结果不难发现，单步解析模式［式（3.54）］的计算精度要高于多步递推模式［式（3.56）］的计算精度，这是因为单步解析模式在每一个计算时段的计算中采用了初始平滩流量Q_b0的实测值，包含了更多的已知条件。在实际应用中可以根据研究的目的和已知条件来选择模型，当已知初始平滩流量Q_b0且只预报未来的一个时段时，应采用单步解析模式；而多步递推模式可以用来研究长时间水沙系列对平滩流量的作用，进行不同设计方案的比较。

3.4.4 考虑与不考虑初始条件的模型比较

为了探讨由Q_e0代替Q_b0可能带来的误差，采用滞后响应模型考虑初始条件的多步递推模式［式（2.14）］，并采用式（3.50）计算平滩流量的平衡值，使用表3.7中的系数和指数可得：

或

式（3.58）代表的多步递推模式Ⅲa为考虑初始条件的计算公式，而式（3.56）代表的多步递推模式Ⅲb为不考虑初始条件的计算公式。为了对这两种计算公式进行比较，图3.33点绘了采用两个公式计算的花园口站和利津站平滩流量的决定系数随包括年数n的变化。可以看到，开始时式（3.58）的R²值远大于式（3.56）的值，但在n=4后两者的R²值趋于相同，原因是当n逐渐增大时，式（3.58）中已知Q_b0的影响变弱，以致于最后与不包括Q_b0的式（3.56）的计算结果一样。图3.33还表明，在来水来沙条件发生变化时，黄河下游河道冲淤调整达到新的平衡状态的时间约为5年，这与以往认识是一致的（吴保生等，2007）。

图3.33 计算与实测平滩流量之间决定系数随包括年数n的变化

作为例子，表3.8列出了用式（3.58）对花园口站1961年平滩流量的计算步骤。这里把第4个时段，即1961年，看作是当前时段，Q_b4是前期5个时段（1957—1961年）水沙条件累积作用的结果。对于其他年份或采用不同n值的情形，其计算过程与表3.8所示相同，这里不再赘述。

需要说明以下几点：

1）式（3.50）计算平滩流量平衡值时只采用汛期平均参数，忽略了非汛期来水来沙的影响，原则上讲有所不足。我们曾采用年平均径流量和年平均来沙系数进行计算，发现结果与只采用汛期的结果基本相同。我们知道，洪水是塑造河槽形态的决定因素，因此，忽略非汛期来水来沙的影响也是可以理解的。

2）对于不同的测站，式（3.51）和式（3.52）中的系数和指数不同，影响了公式的通用性，原因可能是模型中还有一些水文泥沙和河床边界条件参数没有得到考虑，如洪水过程、来沙组成、床沙黏粒含量、河道整治工程、河床形态、河型等。这些参数及非汛期水沙条件是造成公式中的指数和参数不同的主要原因，需要进一步深入研究。

3）式（3.50）中用Q_e0替代Q_b0的处理方法有不足之处。Q_b0是一个水沙系列变化前的初始平滩流量，而Q_e0是该水沙系列对应平衡状态的平滩流量。从物理概念上讲，前者是一个状态的初始值，后者是该状态对应的潜在平衡值，将这两个值取相等的处理方法，虽然作为近似处理方法是可以接受的，但更合理的方法还有待进一步探讨。

3.4.5 平滩流量滞后响应模型的改进方法

从式（3.52）可以看出，平滩流量是当年以及前期n年内的平滩流量平衡值的加权函数，其权重取值根据对应的年份和体现河道调整速率的参数β共同确定；式（3.50）表明，平滩流量平衡值Q_e取决于当年的上游来水来沙条件，并以汛期平均流量和汛期平均来沙系数来体现。虽然模型应用到黄河下游取得了较好的效果，但仍然有部分测站部分年份的计算值与实测值的符合程度不够理想，原因来自两个方面：一是原模型的结构形式在一些情况下会导致Q_e0的权重分布不合理；二是模型着重考虑了汛期来水来沙量对平滩流量调整的影响，但未考虑来水过程的作用以及不同来流情况下河床调整速率的差别。因此，本节围绕以上两方面对模型进行修正。

（1）公式结构的改进

由第2.2节滞后响应模型的不同模式可知，式（3.52）是基于模式Ⅲb提出的，该模式为消除模型对Q_b0的依赖，以Q_e0代替Q_b0，此时可能出现第0年Q_e0的权重大于第i年Q_ei的权重的情况（i>1）。式（2.17）则解决了这一问题，因此，可采用式（2.17）计算平滩流量：

式（3.59）是在假定各个时段的参数β不变的前提下得出的。但实际上，参数β反映了河道通过自我调整趋近于平衡状态的速率，β值越大河道调整越快，反之调整越慢，因此，不同年份的河道形态以及上游来水来沙条件有所不同，河道调整的速率也应该存在一定的区别。为此对式（3.59）的结构形式做进一步的改进，记第i年（i=0，1，…，n）对应的参数β为β_i，经过计算时段Δt后，式（3.59）记为：

下一时段末的平滩流量可以表达为：

将式（3.60）代入式（3.61）可得：

依此类推，经过n个计算时段，平滩流量Q_bn为Q_b0、Q_e1、…、Q_en的函数，同时以Q_e0代替Q_b0，并参考式（3.59）调整Q_e0的影响权重，得到Q_bn的表达式如下：

显然，当β₀=β₁=…=β_n=β时，式（3.63）即简化为式（3.59）的形式。式（3.63）即为对模型结构进行修正后的计算公式，该式与式（2.28）考虑不同时段β取不同值的模式类似，不同之处是采用初始值的平衡值代替了实际值，并对初始权重进行了调整。

（2）平滩流量平衡值Q_e及参数β计算方法的改进

平滩流量平衡值Q_e是指作用于冲积河流的来水来沙条件维持足够长时间不变，所塑造出的达到冲淤平衡状态的河道形态对应的平滩流量。确定的来水来沙条件对应着唯一的Q_e。但实际河流的历年水沙条件总是不断变化的，河道形态的调整难以在短时间内达到平衡状态，因此实际中很难直接观测到Q_e。

为便于计算，Q_e可以表达为描述来水来沙条件的特征指标的函数形式。式（3.50）将Q_e表达为Q_f和ξ_f的函数，体现了汛期平均来水来沙量的影响。除此之外，来水来沙过程及来沙组成也影响河道形态调整的最终结果。选择洪峰流量（以汛期最大日均来流量来代表）从一定程度上反映来水过程的不同，在相同的汛期平均流量下，洪峰流量越大说明集中出现的大流量越大，反之洪峰流量越小说明来水过程越趋于均匀。利用汛期日均最大流量对Q_e的计算方法进行修正，将式（3.50）改写为如下形式：

式中：Q_m为汛期日均最大流量，m³/s，代表洪峰流量；K、a、b和c为相应的系数和指数；其他物理量意义与式（3.50）相同。

此外，河道调整的速率与河道边界及来水来沙条件有关。不同的河道边界及来水来沙条件对应的水流剪切力及河床组成物质抗拒运动的力不同，从而影响河道冲淤调整的速度。为简化计算，这里仅考虑来水量对参数β的影响。针对不同年份汛期平均流量Q_f的差异，按如下方法计算参数β：

式中：β₀和d为系数和指数。

式（3.64）和式（3.65）给出了模型中参数的表达形式，但相应的系数和指数需要根据实测资料确定。具体方法是：将式（3.64）和式（3.65）代入模型中，利用实测资料拟合相应的系数和指数的值，在兼顾参数取值合理的前提下，使得模型计算的平滩流量与实测值符合最好。此外，为方便应用，本节仅选择部分指标代表来水来沙条件，并不能完全描述来水来沙条件的所有特征。如来沙过程、组成等因素的影响，本节暂未作进一步的探讨。

（3）计算过程与结果

由前述可知，黄河下游平滩流量与包括当年在内的前期5年左右的水沙条件有关，因此，在本节平滩流量的计算中取n=4。采用修正后的模式（3.63），计算方法如下：

1）已知汛期最大日均流量Q_m、汛期平均流量Q_f和汛期平均含沙量C_f，计算汛期平均来沙系数ξ_f=C_f/Q_f，由式（3.64）计算平滩流量平衡值Q_ei。

2）利用式（3.65）计算每年的河床调整速率参数β_i。

3）将Q_ei和β_i代入式（3.63）计算平滩流量Q_bn。

计算所需参数及计算值与实测值决定系数见表3.9。

图3.34以花园口站为例给出了模型修正前后平滩流量计算值与实测值的对比情况，可以看出修正后模型的计算精度较修正前有较大幅度的提高。修正前［式（3.52）］花园口站平滩流量计算值与实测值的决定系数、最大相对误差和平均相对误差分别为：0.83、20%和8%，修正后［式（3.63）］计算值与实测值的决定系数、最大相对误差和平均相对误差分别为：0.88、20%和7%。说明通过对模型结构的修正，并同时考虑洪峰流量对Q_e的修正和汛期平均流量对β的修正，模型更加真实地反映了黄河下游平滩流量随水沙条件变化的调整规律。

图3.34 模型修正前后花园口站平滩流量计算值与实测值的对比

从表3.9可以看出，体现来水过程修正的指数c为正值，洪峰流量与平滩流量平衡值成正相关关系，说明在其他条件相同时，来水过程中集中出现的大流量越大越有利于河槽的塑造，相反趋于均匀的来水过程容易导致河槽萎缩；参数β的修正指数d为正值，说明汛期来流量越大，河道调整速度越快，这也符合河床演变学的一般认识。

此外，分析中还按上述修正过程对原模型逐一进行单独修正，结果发现，各种修正方法对平滩流量计算精度均有一定程度的提高，而综合修正后计算精度有进一步的小幅提高。以修正前后计算值与实测值决定系数R²提高的百分比为衡量标准，以花园口站为例，单独在计算Q_e时考虑洪峰流量的影响使模型计算精度在原来基础上提高5%；仅按式（3.59）对模型中Q_e0影响权重进行调整，模型计算精度可在原来基础上提高6%；而其他条件不变的情况下，仅按式（3.63）和式（3.65）考虑来流量对参数β的影响，模型计算精度在原来基础上提高5%。综合上述三项修正，模型计算精度提高7%。