20世纪中期以来不同时段黄河年输沙量对水土保持的响应

李 敏， 朱清科

(1.水利部 黄河水利委员会 黄河上中游管理局， 710021，西安； 2.北京林业大学林业生态工程教育部工程研究中心，100083，北京)

黄土高原是世界上水土流失最严重的地区，流经黄土高原的黄河是世界上年输沙量最大的河流之一。20世纪中期黄土高原地区的水土保持就纳入国家重点治理。经过近70年的大规模治理，水土保持工程对遏制黄土高原水土流失产生了明显的作用，进而对减少黄河年输沙量发挥了巨大的“治本”效果。

近期，国家主要领导人视察黄河，提出：治理黄河，重在保护，要在治理。要坚持山水林田湖草综合治理、系统治理、源头治理。深入研究黄河年输沙量对水土保持的响应特性与效应，对于评判黄河治理效益，确定黄河治理重点具有重要意义。

黄河输沙量的变化与降雨和水土保持紧密相关。以往在研究水土保持与黄河年输沙量的关系中，多采用水文法和水保法、按照年代划分时段进行计算[1-5]。

理论和实践表明，在自然和人类活动的“干扰”下，河流输沙量对水土保持的响应遵从统计学原理。笔者在以往分析50多年整体系列的基础上[11]，进一步研究了不同时段黄河年输沙量对水土保持的响应特性。

1 研究区域与时段

黄河泥沙基本上来自流经黄土高原地区的河口镇(头道拐)到龙门区间的各支流和渭河、北洛河、汾河，其中主要产沙区域为黄河的河口镇(头道拐)到龙门区间(以下简称“河龙区间”)。本研究区域为渭河(华县站)、北洛河(状头站)、汾河(河津站)控制区域和河龙区间(以下简称“4站区域”)。4站区域总面积30.05万km2(图1)。

图1 研究区域示意图Fig.1 Schematic diagram of the study area

黄河自1919年设站观测，距今不到100年。本研究时段为1954—2015年，共62年。

1954—2015年(不计头道拐以上来沙)黄河年输沙量波动变化，进入21世纪平均值仅2.6亿t(图2)[6-10]，其中2015年仅0.55亿t。

黄河年输沙量承载了黄土高原生态环境状况、水土流失与水土保持的信息。

从20世纪50年代初期开始，黄河中游黄土高原地区的水土保持就开始了大规模的实施。20世纪70年代以前，由于受技术和资金条件的限制，工程质量普遍较低，水毁严重，期间还经历过毁林开荒等人为破坏，水土保持设施的面积和数量增长缓慢。进入20世纪80年代，黄河中游水土保持工程质量得到提升，资金投入得到保障，法规监督得以落实，水土流失治理面积基本以线性趋势持续增加(图3)。

2 研究方法

2.1 原理与方法

根据水文学和水土保持学基本原理，河流年输沙量对降雨和水土保持的响应遵从统计学规律。因此，笔者将研究区域黄河年输沙量的产生过程视为一个灰箱或者黑箱，采用逐步回归分析方法，筛选具有物理成因概念和含义的解释因子(自变量)，获得符合水文学和水土保持学基本原理、仿真模拟效果好的计算结果，以此研究黄河年输沙量对水土保持响应的特性。

图2 1954—2015年黄河4站区域年输沙量变化图Fig.2 Change of annual sediment transport in the area of 4 stations in the Yellow River in 1954 to 2015

图3 1980—2015年黄河中游黄土高原地区水土流失治理面积变化(资料来源：黄河上中游管理局)Fig.3 Change of the controlled area of soil erosion in the Loess Plateau of the middle reaches of the Yellow River in 1980 to 2015(Source: Upper and Middle Yellow River Bureau)

2.2 降雨因子选择

已知黄河泥沙主要源自河口镇到龙门区间(以下简称“河龙区间”)、源自汛期降雨(暴雨)；因此，对于降雨因子，初步选取4站区域年降雨量，即河龙区间(包括河龙区间5个亚区)年降雨量、7—8月降雨量、最大3日降雨量、最大一日降雨量等。以上共初步选择22个降雨因子。

2.3 水土保持因子选择与预处理

计算中对水土保持坡面措施面积、大型坝(骨干坝)坝控面积进行“耦合”处理[11-12]，以使坡面措施与淤地坝相互协调，既不重复计算，又不产生矛盾。

为了使今后对不同时间和空间尺度的计算结果具有可比性，笔者将“水土流失治理度”作为参选自变量。水土流失治理度的定义为

式中：D为水土流失治理度，%；A1为水土保持综合治理面积，万km2；A2为水土流失面积，万km2 [13-15]。

2.4 研究时段划分

对于黄河年输沙量，一般认为：在20世纪60年代前后为“天然状态”的输沙量，即流域下垫面破坏严重，水土保持等治理工程措施薄弱阶段；在70—80年代为水土保持措施减沙效果显现时期，该时期减沙效果波动变化；80年代后期以来水土保持工程项目持续实施，减沙效果越亦显著。据此，本文将1954—2015年划分为年输沙量“高位运行，剧烈波动” “波动递减” “低位运行”3个时段分别进行计算分析。

3 结果与分析

3.1 1954—1966年计算结果分析

观察1954—1966年黄河4站年输沙量变化(图2中1954—1966年时段)可以看出：该时段黄河年输沙量均值达到15.7亿t；波动剧烈，年际间波动幅度最大超过20亿t。简单回归分析表明，尽管该时段年输沙量波动幅度很大，但是时段变化趋势微弱，基本没有统计学意义上的增减变化趋势。

图5 1954—1966年实测年输沙量与计算年输沙量比较Fig.5 Comparison of measured annual sediment transport with calculated annual sediment transport in 1954 to 1966

统计分析表明，由于1954—1966年研究区域黄河年输沙量较大，年际间的波动剧烈；同时该时段水土保持工程面积较小，治理薄弱，水土流失治理度偏低，黄河4站年输沙量对水土流失治理度没有统计学意义上的响应(图4)。

图4 1954—1966年黄河4站年输沙量与水土流失治理度散点图Fig.4 Scatter plot of annual sediment load with controlled percentage of soil erosion from 1954 to 1966 in the area of four stations in the Yellow River

经过对22个降雨因子和水土流失治理度因子的逐步回归计算结果表明，仅“河吴区间7—8月降雨量”1个因子可以“解释”该时段年输沙量近78%(R=0.892 330,调整的R2=0.777 731)的变化。线性回归显著性水平α=0.001。

比较回归模型计算值与实测值(图5)看出，计算值与实测值拟合较好，反映了黄河年输沙量年际变化来自于对河吴区间7—8月降雨量年际变化的响应。

其中，1966年回归标准化残差为2.184 6，该年的计算值明显小于实测值。初步分析，造成这一现象的原因是该年黄河流域发生较大的局部暴雨，冲毁了部分水土保持工程措施，造成以往多年拦淤减少的泥沙在当年被冲入黄河，使“河吴区间7—8月降雨量”不能够“解释”这些输沙量，因而计算年输沙量明显小于实测年输沙量[16-18]。

对1954—1966年时段的回归计算结果说明，这一时段黄河年输沙量的变化主要来自汛期降雨的变化，由于该时段水土流失治理措施规模较小，输沙量对其没有响应。

3.2 1967—1987年计算结果分析

3.2.1 计算结果的总体分析 观察1967—1987年黄河4站年输沙量变化(图2 1967—1987年时段)可以看出：该时段黄河年输沙量波动下降，简单回归分析表明，该时段年输沙量呈指数减少趋势。

经过对1967—1987年22个降雨因子和水土流失治理度因子的逐步回归筛选计算，建立了由“河吴区间7—8月降雨量”“水土流失治理度”“河龙区间最大一日降雨量”“ln区域年降雨量”等4个自变量组成的非线性回归方程。回归方程R=0.954 689，显著性水平α=0.001。

根据回归方程，对计算年输沙量与实测年输沙量比较(图6)看出，计算值与实测值拟合很好，不但描述了输沙量年际间的波动变化，而且反映了长期递减的趋势。

图6 1967—1987年实测年输沙量与计算年输沙量比较Fig.6 Comparison of measured annual sediment transport with calculated annual sediment load in 1967 to 1987

回归分析表明，1967—1987时段，黄河年输沙量的主要影响因子仍然是汛期降雨，年际间输沙量的波动变化响应来自汛期降雨，“河吴区间7—8月降雨量”因子“解释”了该时段黄河年输沙量接近63%的变化，其次“水土流失治理度”1个因子可以“解释”黄河年输沙量接近50%的变化，黄河年输沙量波动下降的现象主要是对该因子的响应结果。

3.2.2 回归方程的可视性分析 由于是非线性关系，不同水土流失治理度对年输沙量的影响程度不同，影响的过程(趋势)也不同。根据回归方程做出在相同降雨条件下不同水土流失治理度不同视角的三维曲面图形，直观对比分析不同水土流失治理度对年输沙量的影响程度和影响趋势。为了做出三维图形，对贡献率最小的“ln区域年降雨量”固定为时段的平均值，选择时段起始年(1967年)和时段末尾年(1987年)，分别做出该2个年份年输沙量对水土流失治理度的响应曲面进行比较分析(图7)。

从图7总体上看出，1967年和1987年的水土流失治理度对年输沙量影响显著不同。在曲面的最高点，即2个降雨因子均处于最大值时，1987年输沙量响应曲面比1967年输沙量响应曲面降低超过20亿t；在曲面的最低点，即2个降雨因子均处于最小值时，2个年输沙量响应曲面相差幅度约5亿t。对比分析说明，随着水土流失治理度的增加，一方面，总体上年输沙量显著减少，另一方面，年输沙量响应曲面在高降雨条件下的减少幅度大于低降雨条件[11]。

这一现象说明水土流失治理不仅大幅度减少了黄河年输沙量，而且使年输沙量对降雨变化响应的“敏感度”降低。

A为1967年水土流失治理度水平下年输沙量响应曲面,B为1987年水土流失治理度水平下年输沙量响应曲面。A is the annual sediment transport response surface for the controlled percentage of soil erosion in 1967, and B is the annual sediment transport response surface for the controlled percentage of soil erosion in 1987.图7 年输沙量对不同水土流失治理度的响应曲面Fig.7 Response surface of annual sediment transport with different controlled percentage of soil erosion

图8 1988—2015年实测年输沙量与计算年输沙量比较Fig.8 Comparison of measured annual sediment transport and calculated annual sediment transport in 1988 to 2015

3.3 1988—2015年计算结果分析

3.3.1 回归计算结果分析 在1988—2015年时段，黄河4站年输沙量仍然以指数变化趋势急剧波动下降，显著减少(图2 1988—2015年时段)。进入21世纪，年输沙量在6亿t以下继续波动减少，其中2015年仅0.55亿t。

1988—2015年时段，黄河年输沙量对水土保持的响应十分显著。回归分析表明，仅“水土流失治理度”1个因子就能够“解释”黄河年输沙量接近80%的变化(图8)，回归方程的相关系数R达到0.881 389，显著性水平α=0.001。

如前分析，黄河年输沙量的年际间波动来自于降雨量的年际变化，而长期递减趋势则主要来自于对水土流失治理度的响应，所以“水土流失治理度”的变化特性较好地拟合了黄河4站年输沙量的变化趋势。

图9 1988—2015年回归结果预测分析Fig.9 Prediction and analysis of regression results in 1988 to 2015

3.3.2 1988—2015年时段减沙分析 进一步分析1988—2015年回归结果，计算在80%和95%置信度时的回归计算值上限和下限，分析黄河4站年输沙量变化趋势(图9)。2008年时水土流失治理度超过50%，计算年输沙量不到2亿t，其80%置信度上限5.9亿t；2015年水土流失治理度超过60%，计算年输沙量仅0.7亿t，其80%置信度上限2.7亿t，95%置信度上限5.6亿t。

这一结果表明，根据1988—2015年数据建立的回归模型分析，在当前水土流失治理程度条件下，以95%的置信度预测，黄河年输沙量将小于6亿t。

与研究区域人为破坏严重、水土流失治理薄弱状态下的年输沙量比较，当前黄河中游黄土高原地区的水土保持治理措施年减沙量超过10亿t。

4 结论与讨论

1)依据水土保持和降雨与黄河年输沙量间的统计学关系，划分3个时段，采用回归分析方法对观测数据进行计算。结果表明：1954—1966时段流域下垫面破坏严重、水土流失治理度偏低，黄河四站年输沙量对水土保持没有响应；1967—1987时段水土流失治理度逐渐增加，黄河四站年输沙量开始对其有所响应；1988—2015时段水土流失治理度持续增加，黄河四站年输沙量对其响应强烈，仅水土流失治理度1个因子就较好地“解释”了20世纪80年代以来黄河年输沙量锐减的趋势。

2)在计算时段内(以1967—1987年时段为例)，相同降雨条件下，水土流失治理度的增加一方面显著减少黄河年输沙量，另一方面，随水土流失治理度的提高，年输沙量对降雨的响应趋于“迟缓”。由此，在1988—2015年时段，仅由水土流失治理度1个因子就能够很好“解释” 黄河年输沙量的变化。

3)根据现有数据分析，在水土流失治理度超过60%(2015年治理度水平)时，以95%置信度预测黄河年输沙量小于6亿t。当前由于实施水土保持治理工程，黄河中游黄土高原地区的年度水土流失量减少了10亿t以上。