三门峡水库时空冲淤与滞后响应

郑 珊，吴保生，侯素珍，吕宜卫

（1.武汉大学 水资源与水电工程科学国家重点实验室，湖北 武汉 430072；2.清华大学 水沙科学与水利水电工程国家重点实验室，北京 100084；3.黄河水利委员会 黄河水利科学研究院，河南 郑州 450003）

1 研究背景

修建大坝是人类对河流最强烈的干扰之一，其影响既涉及下游河段[1-2]，也涉及大坝上游河道[3-4]。大坝修建后上游河段受到蓄水及水库运用的影响发生溯源冲淤[5]，同时上游来水来沙条件不断变化，上游水沙与下游水库运用协同作用，使库区河道的时空冲淤规律十分复杂。研究大坝上游河道对来水来沙与水库运用的响应过程与时空演变规律，有助于深化对库区河道演变规律的认识，为水库运用及泥沙管理提供科学参考。

冲积河流对外部扰动（如水沙变化、大坝修建等）具有滞后响应的特征[6-9]，这种时间滞后性可通过弛豫时间（relaxation time）[10-11]或前期水沙等条件的影响反映[6]；同时，河流系统内部上、下游河段之间具有联动性，反应在同一冲淤阶段或演变规律在上、下游河段的传播影响。大坝上游河道不仅在时间上对水库运用及水沙变化等具有滞后响应的特征，而且在空间上还存在沿程与溯源冲淤的传播及联动。例如，郑珊等[12]的研究认为渭河下游及小北干流河道演变滞后于来水来沙变化及三门峡水库运用；Wu等[6，13]研究发现潼关高程的变化滞后于来水来沙条件约4～6年；王兆印等[14]发现潼关高程大幅抬升和下降引起溯源淤积和冲刷行波，其向上游传播速度约10 km/a，传播过程中冲淤幅度不断衰减；Zheng等[15]对渭河下游时空冲淤规律研究表明，该河道受到来水来沙和水库运用的协同作用，在不同时段具有不同的时空冲淤分布特征。综上所述，关于大坝上游河床演变的时间滞后规律开展较早并已取得一定的认识，但关于库区河道冲淤的时空冲淤联动规律尚认识不足，仍处于定性研究阶段，且其数学计算方法仍具有改进空间[6，16]。

三门峡水库自1960年建成运用后，因泥沙问题先后采用了多种运用方式，同时上游来水来沙受自然和人为因素影响发生变异[17]，库区河道的冲淤演变规律十分复杂；黄河水利委员会对库区水沙和河道冲淤变化进行了长期系统的观测，因此，该水库是研究坝上游河道对扰动的滞后响应及时空冲淤演变规律的理想对象。目前已积累了较多关于三门峡水库的时空冲淤演变规律的研究成果，例如，韩其为[18]指出水库冲刷与淤积并不是完全可逆的，在研究三门峡运用方式时应考虑这一不可逆性；胡春宏等[19]及林秀芝等[20]分析了不同时期引起库区淤积的主要影响因素，并对水库的运用方式进行了研究。然而，在较长时空尺度上对三门峡库区冲淤演变的时间滞后与空间联动规律的研究仍有待深入。

本文以1960—2016年长系列水沙及库区冲淤资料为基础，深入分析库区河道的时空冲淤演变规律，研究河道对上、下游扰动的时间滞后与冲淤在时空上的传播特性，通过改进现有河床演变的滞后响应模型建立长时段库区冲淤的计算方法，为三门峡水库及多沙河流上水库运用与泥沙管理提供科学参考。

2 三门峡水库运用与水沙条件

三门峡水库运用后库区（大坝至潼关段）迅速淤积，同时引起渭河下游和黄河小北干流（潼关至龙门河段，图1）回水淤积严重。为减轻上游淤积、控制潼关高程抬升、减小上游河道的防洪压力，水库先后采用了不同的运用方式，包括1960—1961年蓄水拦沙，1962—1973年滞洪排沙，1974年至今蓄清排浑运用，2003年后在蓄清排浑的基础上进一步开展了“318”控制运用，即保持非汛期库水位不超过318 m。

图1 三门峡水库平面位置和测量断面布置示意图[21]

三门峡库水位随着水库运用方式的调整发生变化（图2）。可以看到，1960年水库蓄水使库水位迅速上升（9—10月上升约30 m），1961年水库高水位运用，1962—1973年滞洪排沙运用后库水位明显降低；1974—2002年蓄清排浑运用期间汛期与非汛期水位分别平均约304和316 m；自2003年“318”控制运用以来，汛期平均水位约306 m，非汛期平均水位维持在317 m且波动较小，该段时间内的汛期和非汛期平均库水位缓慢抬升，且汛期抬升更大，2016年汛期和非汛期库水位分别比2003年分别抬高约3.2 m和2.2 m（图2（b））。从最高日均库水位来看，蓄清排浑运用后，最高日均水位缓慢下降，从1974年的324.8 m下降至2003年的317.9 m。

基于日均库水位资料，计算三门峡日均库水位达到非汛期平均水位的日期，该日期一般出现在非汛期（汛期结束后水库开始蓄水），该日期出现越早，说明非汛期高水位蓄水开始越早，越不利于汛末或汛后库区冲刷和潼关高程降低。图2（c）显示2003年前该日期出现的时间波动较大，但具有一定的提前趋势，2003年后该日期明显提前，往往在汛期还没有结束库水位就已蓄至非汛期平均水平。

图2 三门峡库水位变化

1960—2016年三门峡入库水沙量具有明显的减小趋势（图3），表1 对比了蓄清排浑运用（1974年）前后、龙羊峡水库建成（1986年）前后及“318”控制运用前后的水沙量，1986年后入库水量明显减少，汛期水量占全年水量之比减小，出现非汛期总水量大于汛期水量的情况；来沙量持续锐减，2003—2016年来沙量仅为1960—1973年的16%。

图3 1960—2016年潼关站实测水沙量变化

表1 年均入库水沙条件

图4 三门峡水库月排沙比变化

定义三门峡水库排沙比为出库三门峡站沙量除以入库潼关站沙量，排沙比大于100%表示库区冲刷，反之库区淤积。不同月份的排沙比（图4）显示，1962—1973年滞洪排沙期间非汛期排沙效果较好，1974年蓄清排浑运用后，主要排沙期在6—10月，2005年后6—7月排沙比明显增加，排沙时间比以往提前且更为集中；2002年后8—10月份排沙比减小，部分年份甚至小于100%（即库区淤积）。

3 三门峡水库冲淤规律

3.1 河道纵向冲淤变化根据库区黄淤1 至黄淤41 断面（断面位置见图1）的实测资料，得到库区河床各年深泓纵剖面，对其进行线性拟合，假设拟合线的决定系数R2＞0.5 时，其斜率可近似为河床纵比降，得到河床深泓纵剖面及纵比降的变化（图5—图6），据此将河床冲淤变化分为如下4 个阶段：

（1）1960—1969年淤积期。该时段内来水来沙量较大，多年平均来水量450亿m3、来沙量14.2亿t，库水位较高，多年平均311 m，库区年均水面比降（定义为潼关与三门峡平均水位差除以河长）较小（约0.14‰），同时水库泄流能力不足，导致河床不断淤积抬高。1960—1964年河道以溯源淤积为主，靠近大坝淤积厚度最大，黄淤1断面深泓累计淤积21m，往上游淤积厚度逐渐减少（图5（a）），纵比降不断减小（图6），1964—1969年沿程淤积占主导，纵比降有所增大。

（2）1970—1973年冲刷期。该时段内入库水沙量仍较大，多年平均来水量302亿m3、来沙量12.4亿t，库水位在4个时段内最低，多年平均约305 m，该时段内大坝完成二期改建，泄流能力增大，库区年均水面比降最大（达0.2‰），库区产生明显的溯源冲刷（图5（b）），黄淤1 断面累计冲深约8 m，河道深泓纵比降明显增大（图6）。

（3）1974—2002年淤积期。该时段内水库采用蓄清排浑运用，多年平均来水311 亿m3、来沙8.5亿t，多年平均库水位312.1 m，库区多年平均水面比降约0.14‰，库区发生一定的淤积。1974—1978年间以溯源淤积为主，1978—2002年库区中段淤积加大，黄淤33断面以上淤积较少，黄淤1断面深泓累计淤积抬高约12 m。

图5 三门峡水库河床深泓纵剖面

图6 库区河床纵比降

（4）2003—2016年冲刷期。2003年后在蓄清排浑的基础上采用“318”控制运用，但冲刷从2002年即开始，如前所述，2003—2016年汛期和非汛期水位均稍高于1974—2002年，但2003年后来沙量锐减，约为上一时段来沙量的27%（表1），库区深泓普遍冲刷，黄淤1断面深泓累计冲深约11m，河道比降有所减小。

3.2 河道横向及断面形态变化河道断面形态变化与上述4 个冲淤演变阶段相对应，图7 给出了典型断面的变化过程。1970—1973年溯源冲刷后大部分河段形成“高滩深槽”的断面形态。1970—1973年河道纵剖面冲刷下降（图5（b）），图8（a）显示这一冲刷过程中河宽普遍减小，二者也综合反映了河道断面形态向“高滩深槽”发展。需要注意的是，三门峡库区河道宽浅不一，通过对比上述4个时段始末断面的形态变化，得到河床的冲淤河宽（即河床冲淤变化的横向宽度，1973年后近似主槽宽度，图8（a）），黄淤28—34河段（距大坝约60～85 km）明显较宽，如黄淤28断面1960—1969年河床淤积宽度约6 km，而部分窄河段如黄淤1断面的冲淤宽度不到1 km；1970—1973年溯源冲刷后河道冲淤宽度普遍减小，之后变化不大。

图7 典型断面形态变化

河道深泓或主槽的横向摆动速率（图8（b））不仅与河道宽度有关，而且与河道冲淤状态有关，宽浅段（如黄淤28—34）断面深泓的摆动速率大于窄河段，处于冲刷状态的河道深泓摆动幅度小于淤积时的摆动，尤其当1970—1973年河道发生溯源冲刷时深泓摆动速率最小，1960—1969年河道淤积时深泓摆动速率最大。此外，断面面积变化的空间分布（图8（c））明显地反映了库区“淤积一大片，冲刷一条线”的演变特征，1960—1969年宽浅段淤积较多，淤积量甚至超过近坝段，1970—1973年溯源冲刷时，冲刷量具有往上游逐渐减小的趋势，大约在黄淤27断面以上河道冲刷量较小。1974年蓄清排浑运用以来的冲淤速率远小于前两个时段，因此也有部分研究[12，20]认为库区河道达到了动态冲淤平衡状态。

3.3 库区时空冲淤变化蓄清排浑运用以来，三门峡水库年内呈现明显的“汛期冲刷非汛期淤积”或“洪冲枯淤”的演变特征（图9（a））；年际间长时段冲淤过程（图9（b））可分为1960—1969年快速淤积、1970—1973年快速冲刷、1974—2002缓慢淤积与2003—2016缓慢冲刷4个阶段。

对比三门峡水库库水位、潼关以下库区累计淤积量、潼关高程以及小北干流和渭河下游累计淤积量随时间的变化过程（图10），可以看出：1960—1961年库水位大幅抬升后，潼关以下库区快速淤积（1960—1964年），潼关高程抬升（1960—1969年），小北干流和渭河下游大幅淤积（1960—1973年），越往上游淤积持续的时间越长，其中库水位于1961年达到最大值，潼关以下累计淤积量于1964年达到最大值，而小北干流和渭河下游在1973年左右累计淤积量达到最大值，比库区潼关以下淤积量达到最大值滞后约9年，比库水位达到最大值晚了约12年，说明空间溯源淤积的影响传播时段较长，在库水位开始下降后，上游河道还在受到溯源淤积的影响。此外，1960—1964年潼关以下库区累计淤积37 亿m3，单位河长累计淤积0.33 亿m3，1960—1973年小北干流和渭河下游分别累计淤积约18 亿m3和10 亿m3，单位河长累计淤积0.14 亿m3和0.05 亿m3，说明溯源淤积向上游传播过程中影响时长增加的同时，淤积幅度在不断减小。

1961—1972年三门峡库水位不断下降，累计下降约28 m，引起上游河道发生溯源冲刷，潼关以下库区在1964—1973年间累计冲刷约9亿m3，潼关高程在1969—1975年持续下降，累计下降约2.6 m，小北干流和渭河下游在1973—1976年发生轻微的冲刷，冲刷量约1～1.5亿m3。由此可见，溯源冲刷向上游传播过程中的冲刷幅度减小，且影响时段较短，例如潼关以下库区共冲刷约9年，而小北干流和渭河下游冲刷持续3年左右，且冲刷幅度小很多。

图8 库区各断面沿程冲淤河宽、深泓摆动速率与面积的变化

1973—1979年库水位抬升，之后稍有下降，库区淤积较缓慢，1973—2002年累计淤积3.5亿m3，小于1976—2002年间小北干流和渭河下游的累计淤积量（分别为8.6 亿m3和3.7 亿m3），说明该时段库区淤积除了受到水库运用影响外，更多受到入库水沙条件的影响，小北干流和渭河下游在1985年后淤积明显加快，这与龙羊峡水库运用后来水量比来沙量减幅更大、含沙量明显增加有关（见表1）。

年均库水位由2003年的311.7 m 抬升至2016年的314.2 m，而在此期间潼关以下库区发生缓慢冲刷（约1.4 亿m3），潼关高程下降0.84 m，小北干流和渭河下游分别冲刷2.8 亿m3和2.4 亿m3，明显大于潼关以下冲刷量，说明2002年后库区冲刷主要受上游来水来沙条件影响，该时段内水量与1986—2002年相差不大，但来沙量减少69%，水流含沙量锐减（表1），引起库区河道发生冲刷。

图9 三门峡水库冲淤变化

4 三门峡水库冲淤过程模拟

4.1 滞后响应模型及其改进吴保生[7]基于变率原理，提出河床演变的滞后响应模型，用于模拟冲积河流受到扰动后的非平衡调整过程，其多步模式如下：

式中：yn为河床演变特征量y在第n个Δt时段末的值；ye，i为y在第i个时段对应的平衡值；i、n分别为时段序号和总时段数；y0为y在i=0时刻的初始值，当初始值未知时，可用初始平衡值ye0估算，从而取消模型对初始条件的依赖[7]，即：

式（2）显示河床演变特征量yn可看作是不同年份平衡值ye，i的加权平均，初始年（i=0）权重为e-nβΔt，第i年（i=1，2，…，n）的权重为（1-e-βΔt）e-（n-i）βΔt，距离当前年越近相应权重越大，当前年（i=n）对应的影响权重（1-e-βΔt）最大。式（2）中各年影响权重之和等于1，即：

李凌云等[16]考虑到初始年份（i=0）的影响权重e-nβΔt在β＞0.693 时大于其上一年的影响权重（1-e-βΔt）e-（n-1）βΔt（i=1），与扰动的影响随着时间衰减的普通认知相矛盾，进而对式（2）进行了改进，提出滞后响应模型的改进模式：

图10 三门峡水库冲淤的时间滞后与空间传播

式（4）初始年份的影响权重为（1-eβ）e-nβ，小于其他年份，但各年权重之和小于1：

这可能引起对河床演变特征量平衡值的偏大估计，基于此，本文将式（4）中各年影响权重进行归一化处理，得到权重归一化的滞后响应模型计算式：

式（6）中各年影响权重λi符合扰动距当前年越远其影响权重越小的一般规律，且各年权重之和为1。

根据滞后响应模型的应用[6，8]，一般取Δt=1，β在0到1之间，通过多元回归得到最佳的前期影响年份n，基于此，本文以Δt=1，β=0.3为例，图11对比式（2）、式（4）和式（6）中当n取不同值时河床演变特征量平衡值的权重，当n较小时（如n=4）式（2）初始年的权重明显偏大，且大于其后1年甚至几年的权重，当n较大时（如n=8）这种偏差虽小但仍存在。式（4）中各年权重随时间非线性衰减，但权重之和小于1。式（6）在式（4）基础上进行归一化处理后，整体增加了各年的影响权重，符合各年权重随时间非线性衰减的规律且权重之和为1。

图11 式（2）、式（4）、式（6）平衡时权重（取β=0.3，Δt=1）

4.2 库区累计冲淤量计算水流功率是河道输沙能力和造床作用的重要参数，假设三门峡库区累计淤积量可表示为水流功率的函数，水流功率可表示为：

式中：Ω为水流功率，N/s；γ为水的容重，N/m3；Q为潼关站流量，m3/s；Ztg、Zsjt分别为潼关和史家滩（坝前）水位，m；ΔL为潼关至史家滩的河长，m。

相对于运用水位Zsjt来说，潼关断面水位Ztg变化相对较小，可视Ztg为常数，例如可取Ztg=328 m作为特征水位或代表水位。在式（7）基础上进一步考虑年入库水流含沙量及汛期水量占全年比值的影响[11]，可得库区累计淤积量平衡值Ve（亿m3）的计算式：

式中：W、S分别为年入库水量（亿m3）和年均含沙量（kg/m3）；根据参考文献[13]，采用流量加权平均坝前水位Zˉ（m）代替库水位，其中Qtg和Qsjt分别为潼关与史家滩站的日均流量；α为汛期水量占全年水量百分比；K、a、b、c和d为待定系数和指数。

需要说明的是，库区累计淤积量平衡值Ve是根据水流功率的定义进行一定简化、并考虑库区泥沙淤积的影响因素后得到，其与水流功率、挟沙力及输沙率具有一定的反比关系，水流功率越大，水流挟沙力和输沙率越大，库区淤积量越小，因此，式（8）中经验参数的正负值应取为：K＞0，a＜0，b＜0，c＞0和d＜0。

采用权重归一化的滞后响应模型式（6）计算三门峡库区累计淤积量V（亿m3），结合式（8），取Δt=1年，得：

参考文献[6]，三门峡水库冲淤主要受当前年和前期约4年的水沙和水库运用的影响，因此取n=4。采用1974—2002年的实测数据对式（9）中的参数进行率定，通过多元回归分析得：β=0.255，K=40.8，a=-0.097，b=-0.025，c=0.008和d=-0.135。采用2003—2016年的资料进行验证，计算与验证结果如图12所示。由图可知式（9）的计算结果较好，能够较好地模拟三门峡库区长时段的冲淤变化过程，1978—2002年、2003—2016年模型计算值与实测值之间的决定系数R2分别为0.88和0.55。

图12（c）对比了库区累计淤积量的实测值与其平衡值（式（8）计算）随时间的变化过程。由于入库水沙和坝前水位等条件不断变化，累计淤积量的平衡值也不断变化，并且其波动幅度一般大于实测累计淤积量。外界条件变化后，河道调整趋向使累计淤积量向其平衡值发展，二者虽难以完全相等，但在多数年份差别不大，反映了库区河道可能处于动态冲淤平衡状态[14，21]。

图12 库区累计淤积量的实测值与计算值的对比及其与平衡值的对比

可采用假设的水沙条件及库水位，利用式（9）分析不同水沙及水库运用条件下库区的冲淤变化，以对比不同水沙条件及水库运用对库区演变的影响。前文分析表明，2003年后库区缓慢的冲刷（2003—2016累计冲刷0.05亿m3）主要是由于来沙量减少；假设2003—2016年水沙与1986—2002年年均情况相同，即水量变化不大，沙量为实际多年平均沙量的3倍多，由式（9）计算得，在假设水沙条件下，库区将淤积约0.6亿m3，由此可见，2003年后来沙量减少是库区发生冲刷的主要原因。需要说明的是，自2006年以来黄河水利委员会连续实施了利用桃汛期洪水降低潼关高程的原型实验，对降低潼关高程和库区冲刷也起到积极作用[22]，限于文章篇幅，未对其作用进行具体的分析与计算，可在后续研究中进行深入探讨。

此外，式（9）虽然可较好地模拟累计淤积量的变化过程，但仍存在一定误差，这与模型的概化，如采用常数328m特征水位代替潼关站年均水位，未考虑日水沙过程、桃汛期降水冲刷、来沙粒径变化、库区河道形态等因素的影响，后续研究可进一步考虑这些因素改进式（9），提高其理论基础与计算精度。

5 结论

三门峡库区经历了1960—1969年快速淤积、1970—1973年快速冲刷、1974—2002年缓慢淤积和2003—2016缓慢冲刷4个阶段。根据库区及回水区河道的时空冲淤演变规律与模拟，得到主要结论如下：

（1）溯源冲刷向上游传播过程中，冲刷幅度逐渐减小，且影响时段较短。越位于上游的河道，发生溯源冲刷的时间越晚，冲刷幅度越小，持续时间越短。如1961—1972年库水位下降，库区潼关以下河道于1964年开始冲刷约9年，回水区小北干流和渭河下游于1973年开始累积冲刷3年，且冲刷幅度明显较小。

（2）溯源淤积的传播时段较长，越往上游，溯源淤积幅度越小，但淤积持续的时间越长，甚至在库水位开始下降后，上游河道可能还受溯源淤积的影响。如水库运用初期库水位上升引起上游河道发生溯源淤积，小北干流和渭河下游淤积量达到最大值的时间（1973年）比库区淤积达到最大（1964年）滞后约9年，比库水位达到最大值（1961年）滞后约12年。

（3）改进了已有河床演变的滞后响应模型，得到权重归一化的滞后响应模型计算式，修正了原模型由于特征量平衡值的权重之和小于1可能引起的对平衡值的偏大估计；基于水流能量的概念及相关假设，推导得到三门峡水库累计淤积量平衡值的计算式，同时考虑当前年及前期4年水沙及水库运用的影响，得到三门峡库区累计淤积量的滞后响应模型计算方法，较好模拟了水库自蓄清排浑运用以来的冲淤变化过程。