基于滞后响应模型的三门峡水库冲淤计算方法

沈 逸，郑 珊，吴保生

(1.清华大学 水沙科学与水利水电工程国家重点实验室，北京 100084；2.武汉大学 水资源与水电工程科学国家重点实验室，湖北 武汉 430072)

1 研究背景

冲积河流与流域不断进行着物质和能量交换。由于河流因素的复杂性，上游输入的泥沙一般难以与当前河段的水流挟沙力相等，导致冲积河流的河床时刻处于冲淤变形中。然而，冲积河流的自动调整作用又使得河床演变朝着变形停止的方向发展，显示出平衡倾向性[1]。

河床受到外界扰动后的自动调整，其初始速度通常较快，但随着河床不断趋近于新的平衡状态，调整速率逐渐减缓，并最终趋近于零，如图1所示。宏观尺度上的河床变形总是滞后于来水来沙条件等外部因素的变化，称为河床演变的滞后响应[2]。对于图1所示的滞后响应模式，吴保生等[3-4]基于变率公式建立了河床滞后响应模型的理论框架，包括通用积分、单步解析、多步递推三种计算模式，适用于模拟不同条件下的河床演变过程。其中，滞后响应模型的单步解析公式可表示为：

y=(1-e-βt)ye+e-βty0

(1)

式中：y为特征变量当前值；ye为特征变量平衡值；y0为特征变量初始值；β为调整速率参数；t为调整时间。

图1 滞后响应模式及特征变量准平衡时间示意

特征变量平衡值ye和调整速率参数β是滞后响应模型的两个关键参数，平衡值ye表示河道受扰动后，特征变量调整的目标值；调整速率参数β则表示特征变量调整的快慢，β越大表示特征变量调整速率越快。平衡值ye往往根据具体研究对象的特点给出计算方法。例如，张为等[5]基于水流挟沙力公式推导得到挟沙强度的计算式，根据挟沙强度与汛期平均含沙量的比值反映河道的冲淤过程，得到三峡水库泥沙淤积量平衡值的计算公式。吕宜卫等[6]在推导下荆江河段单位河长累计冲刷量平衡值计算公式时，在考虑挟沙强度和汛期来沙量影响的基础上，引入了无量纲水面比降来减弱洞庭湖顶托效应对模型计算的影响。郑珊[7]将扰动后河道淤积体的平衡纵剖面形态概化为三角形或梯形，为累计淤积量平衡值的计算提供了一种概化方法。调整速率参数β多取为常数，仅少数研究在模拟黄河下游平滩流量时，把调整速率参数β表示为汛期平均流量的幂函数[8]。

此外，关于河床调整的滞后响应时间，以往研究主要基于资料相关分析，缺乏严格的数学或理论推导。例如，张原锋等[9]发现潼关高程与6年滑动平均年水量具有较好的相关关系。廖治棋等[10]发现监利站平滩面积受到前期3年内来水来沙条件的影响较大。唐小娅等[11]通过计算发现三峡水库汛期泥沙累计淤积量与5年线性叠加坝前水位关系最密切。王彦君等[12]的研究结果表明，黄河下游主槽断面形态受当年在内的前8年水沙条件的累积影响。

河床演变的滞后响应特征在水库淤积过程中体现明显。三门峡水库自1960年9月蓄水运用以来，泥沙淤积一直是制约水库防洪和发电等综合效益发挥的关键问题。过去围绕水库淤积和潼关高程变化与控制开展了大量研究[13-15]。吴保生等[16]发现三门峡库区淤积不仅与当年的流量和运行条件有关，还受到前3～4年的流量和运行条件的影响。近期吴保生等[17]采用水沙参数的多年平均值计算参数β，难以体现当前年河床冲淤状态对调整速率参数的影响。

本文将首先基于水库冲淤平衡纵剖面几何关系及其物理图景，推导得到潼关高程和库区累计淤积量平衡值的计算模式。然后引入调整速率参数β随河道冲淤及水沙条件发生变化的计算方法，采用式(1)表示的滞后响应模型单步解析模式，建立库区累计淤积量及潼关高程的计算方法。进而计算并分析河道达到准平衡状态所需要的时间，对比河道实际完成调整量与目标调整量的差距，揭示河床演变的滞后响应规律。

2 潼关高程与累计淤积量滞后响应模型公式推导

2.1 潼关高程平衡值三门峡水库自1973年11月以来采用蓄清排浑的运用方式，非汛期抬高水位运用，汛期平水期控制水位305 m发电，遇洪水时降低水位泄洪排沙。胡春宏等[18]指出三门峡水库蓄清排浑的运用方式是一个根据来水来沙条件和实际需求不断调整的过程。自蓄清排浑运用开始到1985年期间，汛期4个月全部敞泄排沙，后又发展为当潼关站出现3000 m3/s流量的洪水期时水库敞泄排沙。1986—2000年间，由于汛期大于3000 m3/s流量的天数减少，敞泄排沙的入库流量标准下调至2500 m3/s。2000年之后，敞泄排沙的入库流量标准又进一步降至2000～1500 m3/s。

图2 三门峡水库深泓纵剖面示意

图2显示1974年汛后河道深泓纵剖面可用直线拟合，拟合直线的斜率可表示为深泓纵剖面比降Js，截距即为坝前断面深泓高程Z*s。据此，可以得到图3所示水库冲淤平衡状态下的水库纵剖面示意图，潼关高程平衡值Ze(m)及潼关站深泓高程Ztgs(m)分别表示为：

Ze=Ztgs+htg

(2)

Ztgs=Z*es+JesL

(3)

式中：htg为潼关断面在1000 m3/s流量时的水深，m；Z*es为水库冲淤平衡时河道深泓纵剖面起始高程，即坝前断面深泓高程，m；Jes为水库冲淤平衡时河道深泓纵剖面比降；L为潼关断面至大坝的河流长度，取113.5 km。

(4)

(5)

式中：Qtgi为日均入库流量；Si为日均入库含沙量；Zdi为日均坝前水位；n为加权计算的天数，平年n=365，闰年n=366；ε和σ分别为流量和含沙量的指数。

根据以往研究[19-20]，河道深泓纵剖面比降Jes与水沙条件及水库运用水位相关，可表示为：

(6)

将式(4)和式(6)代入式(3)，得到Ztgs的表达式：

(7)

进而将式(7)代入式(2)，即可得到潼关高程平衡值Ze的完整表达式：

(8)

式(8)可进一步简化为：

(9)

(10)

图3 三门峡水库冲淤平衡纵剖面概化示意

2.2 累计淤积量平衡值为了简化计算，将河道断面形态概化为宽度为B的矩形，进而可以得到由初始河床纵剖面和与当前河床纵剖面组成的概化水库平衡淤积体。考虑到淤积量计算中的概化河底高程与深泓高程存在一定差别，图3中把累计淤积量平衡值计算中涉及的概化河底高程和河床纵剖面的符号统一用下标m表示。

累计淤积量平衡值计算中涉及的初始河床边界参数，有概化河底纵剖面的初始比降Jm0、潼关断面的初始概化河底高程Ztgm0和坝前断面的初始概化河底高程Zdm0。

图3中的冲淤平衡时的概化河底纵剖面比降、概化河底纵剖面起始高程及潼关断面概化河底高程分别记为Jem、Z*em、Ztgm。类似式(4)和(7)，Z*em和Ztgm可分别表示为：

(11)

(12)

(13)

按照图3所示水库冲淤平衡纵剖面，累计淤积量平衡值Ve可表示为：

(14)

分别将Z*em和Ztgm的计算式(11)和(12)代入式(14)有：

(15)

进一步化简，可得如下累计淤积量平衡值计算公式：

(16)

2.3 滞后响应模型单步解析模式迭代计算方法当考虑β随时段变化时，应用滞后响应模型的多步递推模式计算较为复杂[17]，因此，本文采用单步解析模式(式(1))对潼关高程和累计淤积量进行计算。对于外部扰动阶梯状变化的情况，一个时段(水文年)的河床调整结果，无论是否已经达到平衡状态，都将作为下一个时段的初始条件对河床演变产生影响，并由此使前期的水沙条件对后期的河床演变产生影响[2]。

因此，在本文模型中，每个时段的特征变量平衡值分别根据对应时段内的年均入库流量、年均来沙系数、年加权水位和年均坝前水位确定，将上一时段特征变量的计算结果作为下一时段的初始条件，采用式(1)逐时段计算每个时段末特征变量的状态值，依次递推便可以得到经过多个时段后的特征变量状态值，如图4。

图4 单步解析模式迭代计算流程示意

单步递推公式的一般形式可表示为：

yi=(1-e-βΔt)ye，i+e-βΔtyi-1

(17)

式中：Δt为时段长度；i为时段序号。

将潼关高程和累计淤积量作为特征变量代入式(17)，通过逐时段递推可以得到潼关高程和累计淤积量的逐年计算值：

Zi=(1-e-βΔt)Ze，i+e-βΔtZi-1

(18)

Vi=(1-e-βΔt)Ve，i+e-βΔtVi-1

(19)

式中：Zi-1和Zi分别为第i-1和i年末的潼关高程计算值；Ze，i为第i年的潼关高程平衡值；Z0为潼关高程初始值，采用1974年末潼关高程实测值；Vi-1和Vi分别为第i-1和i年末的累计淤积量计算值；Ve，i为第i年的累计淤积量平衡值；V0为累计淤积量初始值，采用1966年末累计淤积量实测值。在潼关高程和累计淤积量的计算中，Δt=1 a，调整速率参数β的单位为a-1。式(18)和式(19)中特征变量平衡值Ze，i和Ve，i分别采用式(9)和式(16)计算。

调整速率参数β需综合考虑来水来沙条件及河道冲淤状态的影响。河道调整速率与水流挟沙力及水体中含沙量大小相关，当水体中床沙质含沙量大于挟沙力时，泥沙将落淤，含沙量越大，河床淤积越显著；反之，当水体中床沙质含沙量小于挟沙力时，冲刷河床，冲刷与水动力强弱紧密相关[21]。此外，以往研究表明河流受扰动后的调整速率往往先快后慢[22]。综上所述，假设库区淤积或潼关高程抬升时(即平衡值大于初始值)，β为含沙量的指数衰减函数；当库区冲刷或潼关高程降低时(即平衡值小于初始值)，β为流量的指数衰减函数：

(20)

3 潼关高程和累计淤积量的模拟与验证

3.1 潼关高程和累计淤积量计算结果根据1975—2018年实测水沙资料和潼关高程数据，通过多元非线性回归方法，得到潼关高程单步递推公式(式(18))中平衡值Ze和调整速率参数βtg计算公式的参数值。然后，使用拟合得到的单步递推公式(式(18))计算1975—2018年的潼关高程。

采用多元回归拟合得到潼关高程平衡值Ze和调整速率参数βtg的计算公式分别为：

(21)

(22)

(23)

由式(22)可得潼关高程处于上升状态时调整速率参数与年均含沙量正相关；潼关高程处于下降状态时调整速率参数βtg与年均入库流量正相关，但随着流量的增大调整速率参数增加不明显。

潼关高程历年计算值与实测值的对比情况见图5(a)。潼关高程实测值与计算值在1975—2002年的决定系数R2为0.80，1980—2002年以及1980—2018年的决定系数R2均为0.88，在2003—2018年的决定系数R2为0.75。计算值与实测值在1975—2002年和2003—2018年相对误差的平均值分别为0.06%和0.02%，均处于较低水平。由图5(b)可知，潼关高程计算值与实测值接近，计算效果较以往有所提高。其中，1980—2018年和2003—2018年的决定系数R2分别从之前的最好计算结果0.85和0.74[17]提升至0.88和0.75。

根据1967—2018年实测水沙资料和累计淤积量数据，通过多元回归方法得到累计淤积量单步递推公式(式(19))中关键参数平衡值Ve和调整速率参数βV的计算公式。然后，使用拟合得到的单步递推公式(式(19))计算1967—2018年的累计淤积量。拟合得到累计淤积量平衡值Ve和调整速率参数βV的计算公式分别为：

(24)

(25)

(26)

由式(25)可知库区处于淤积状态时调整速率参数与年均含沙量正相关；库区处于冲刷状态时调整速率参数与年均入库流量正相关。

图5 潼关高程计算值与实测值对比

累计淤积量历年计算与实测值的对比情况见图6(a)。计算值与实测值在各个时段均保持着较高的相关性，决定系数R2在1967—1979年、1980—2002年和2003—2018年分别为0.96、0.92和0.84。计算值与实测值在1967—2002年和2003—2018年相对误差的平均值分别为0.83%和0.45%，计算效果较好。由图6(b)可知计算值与实测值在1967—2002年和2003—2018年均基本保持一致。

图6 累计淤积量计算值与实测值对比

伴随着水库控制方式的调整和来水来沙条件的变化，三门峡库区共经历了1960—1969年的快速淤积期、1970—1973年的快速冲刷期、1974—2002年的缓慢淤积期与2003年之后的缓慢冲刷期4个阶段[23]。4个时段中，2003—2018年淤积量的变化范围最小，为29.41亿～30.66亿m3，最大、最小值差值仅为1.25亿m3。从图6(b)中可以看出2003—2018年数据点(红点)的分布范围明显小于1967—2002年数据点(蓝点)。如何准确模拟2003年之后累计淤积量的小幅度变化是滞后响应模型研究中的难点，本文通过改进调整速率参数的计算公式提升了模型在2003年之后累计淤积量的计算精度，将决定系数R2从之前的最好计算结果0.77[17]升高至目前的0.84。

3.2 平衡值公式起始高程和比降变化范围及合理性检验在收集1967—2018年黄淤1—41断面实测汛后深泓高程数据的基础上，得到了黄淤1—41断面深泓高程的线性拟合公式，并对汛后深泓纵剖面比降、深泓起始高程以及汛后潼关站的深泓高程进行了统计，结果如表1所示。其中，1975—2018年为潼关高程的计算时段，1967—2018年为累计淤积量的计算时段。

1975—2018年汛后潼关高程的多年平均值为327.55 m，汛后潼关深泓高程的多年平均值为324.62 m，两者之差Δh为2.93 m。1975—2018年实测深泓线的起始高程变化范围为291.78～298.79 m，加上Δh，即为式(10)所示Zw的合理变化范围：

294.71 m≤Zw≤301.72 m

(27)

表1 三门峡水库实测汛后纵剖面参数统计结果

模型计算结果显示，1975—2018年期间，水库冲淤平衡时深泓纵剖面比降Jes的变化范围为0.250‰～0.262‰，在同时段汛后深泓比降Js的变化范围之内(表1)。Jes的时段平均值为0.255‰，与同一时段的汛后深泓比降Js多年平均值为0.254‰接近。模型计算得到Zw的变化范围为296.80～300.26 m，在式(27)所示的合理变化范围之内。

在累计淤积量的模拟中，由于概化河底高程与实测深泓高程存在差别，所以仅对概化河底纵剖面平衡比降的变化范围进行讨论。

(28)

1967—2018年实测深泓线比降的变化范围为0.144‰～0.326‰，多年平均值为0.253‰。考虑到在平衡条件下比降的变化范围较小，且为了保证概化河宽B率定结果的准确性。使Jem在实测深泓线比降多年平均值0.253‰上下约0.022‰的范围内变化：

0.231‰≤Jem≤0.275‰

(29)

根据实测断面统计，在1967—2016年，黄淤1—41断面冲淤河宽时段均值的变化范围为0.4～2.3 km，据此可对概化河宽B的变化范围进行如下约束：

0.4 km≤B≤2.3 km

(30)

模型率定得到概化河宽B=1.6 km，模型计算得到Jem的变化范围为0.237‰～0.273‰，在式(29)所规定的合理变化范围之内；时段平均值为0.243‰，同一时段汛后深泓比降Js的多年平均值为0.253‰，可以发现Jem和Js的时段均值接近。因此，本模型中比降Jem的变化范围是合理的。

4 关于滞后调整过程的讨论

4.1 达到准平衡所需时间的计算在滞后响应模型概化示意图1中，根据“扰动”和“准平衡态”两个临界点，将特征变量变化过程分为扰动前、调整阶段、准平衡阶段三个时段。据式(1)可将调整时间t表示为：

(31)

将河道受到扰动后特征变量从初始值调整到平衡值所需的变化量Δy称为目标变化量，Δy=(ye-y0)。因为平衡值ye只可以无限逼近而无法达到，所以还需要定义一个准平衡态。设定特征变量的变化量达到目标变化量Δy的95%时河流系统进入准平衡态，把特征变量y从初始值y0达到准平衡态y′e所需的时间T称为准平衡时间。即当y=y′e时，y′e-y0=0.95(ye-y0)或y′e=0.95ye+0.05y0，将y=y′e代入式(31)，可得：

(32)

由式(32)计算得到潼关高程和累计淤积量在不同年份的准平衡时间如图7所示。注意到潼关高程与累计淤积量准平衡时间均在2014年和2015年出现峰值，原因是这两年的入库含沙量严重偏低，分别只有3.18和2.70 kg/m3，分别按式(22)和(25)中βtg淤和βV淤的公式，计算所得2014年和2015年的潼关高程调整速率参数βtg分别只有0.22 a-1和0.20 a-1(图8(a))，累计淤积量调整速率参数βv分别只有0.09 a-1和0.06 a-1(图8(b))。根据式(32)可知，特征变量的准平衡时间与调整速率参数β成反比，所以潼关高程和累计淤积量准平衡时间均在2014年和2015年出现高峰。由此可见，严重偏小的含沙量是导致2014年和2015年调整时间偏长的原因。

图7 潼关高程和累计淤积量准平衡时间逐年分布

图8 潼关高程和累计淤积量调整速率参数逐年变化

由于2014年和2015年潼关高程和累计淤积量准平衡时间与前后年份差异较大，为避免个别年份的极端值掩盖准平衡时间的整体特性，在准平衡时间的统计中略去2014年和2015年的数据。同时，为了便于对比，将潼关高程与累计淤积量准平衡时间的研究时段取为一致，统计结果见表2。

表2 潼关高程与累计淤积量准平衡时间 (单位：a)

统计结果显示，潼关高程准平衡时间多年平均值约为10.2 a，累计淤积量准平衡时间多年平均值约为6.7 a，潼关高程的准平衡时间长于累计淤积量的准平衡时间表明潼关高程对水沙条件的响应要慢于累计淤积量。同时，注意到2003年之后，受水库调控方式改变及来沙量减小等因素的影响，潼关高程和累计淤积量的准平衡时间均有一定程度的延长。

如果水沙条件不变，自河道受扰动开始，经过时间T(准平衡时间)后，河道将进入准平衡态。但由于水沙条件总是处于不断变化之中，这一准平衡态一般难以达到。此外，库区河道达到淤积平衡与水库淤满并不是一个概念，淤积平衡是指一定水沙条件所对应的库区河道平衡形态，而水库淤满或达到淤积使用寿命则是指淤沙库容淤满而没有足够的库容来拦截泥沙。在三门峡水库蓄清排浑运用下，若水沙条件不变，达到淤积平衡同时保持一定的淤沙库容是可以实现的。

4.2 特征变量逐时段以平衡值为目标的调整过程由于外部条件的变化，特征变量每个时段的平衡值yei会随着外部条件的改变而处于动态调整之中。

根据滞后响应模型单步递推公式(式(17))，可将每个计算时段内特征变量平衡值与计算值的关系表示为：

(33)

式中特征变量计算值与平衡值的关系可进一步表达为：

yi-1

(34)

yi-1>yi>yei(冲刷时)

(35)

图9(a)(b)分别为潼关高程和累计淤积量平衡值与计算值的关系，可以看到在模型计算中，特征变量平衡值起到了“调整目标”的作用。

图9 特征变量平衡值与计算值关系

在每个计算时段内，特征变量会以时段初值yi-1(即上一时段末特征变量的计算值)为基础，以本时段的平衡值yei为目标，随着时间推进向yei靠近，经过一个时段的调整得到时段末计算值yi。当yei高于yi-1时，特征变量向增大方向调整；反之，特征变量向减小方向调整。按照这一模式，特征变量的计算值从初始年份开始，通过逐时段调整得到当前年的计算值。这也说明了当前时段的河床演变不仅受当前水沙条件的影响，而且通过边界条件，还受前期水沙条件的影响[2]。

5 结论

(1)基于水库冲淤平衡纵剖面几何关系和机理分析，推导了潼关高程和累计淤积量平衡值的计算公式；考虑调整速率随库区冲淤和水沙条件的变化，提出了滞后响应模型调整参数β的计算方法；进而采用滞后响应模型的单步解析模式，建立了潼关高程和累计淤积量的滞后响应计算方法，论证了模型的合理性，计算结果表明，潼关高程1975—2002年和2003—2018年计算值与实测值决定系数R2分别达到0.80和0.75，累计淤积量1967—2002年和2003—2018年计算值与实测值决定系数R2分别达到0.95和0.84。(2)定义特征变量的变化量达到目标变化量的95%时河流系统进入准平衡态，研究了潼关高程和累计淤积量达到准平衡所需时间的变化规律，结果显示，潼关高程达到准平衡的平均时间约为10.2年，累计淤积量准平衡时间约为6.7年。潼关高程的准平衡时间长于累计淤积量的准平衡时间，反映了潼关高程的变化滞后于累计淤积量的变化。(3)探讨了潼关高程和累计淤积量逐时段以平衡值为目标的调整过程与特点，在每个时段内，特征变量会以时段初值为基础，以本时段平衡值为调整目标，随着时间推进向时段平衡值靠近。一个时段的调整结果会作为下一个时段的初始条件使前期的水沙条件对后期的河床冲淤产生影响，显示了河床演变的自动调整过程与滞后响应特性。