黄土高原现状产沙情势评价

刘晓燕

(黄河水利委员会，河南 郑州 450003)

2000—2021年，黄河潼关水文站年均输沙量2.41亿t/a、汛期6—9月含沙量18.3 kg/m3，分别较1919—1959年降低85%、67%。在此背景下，黄河水沙变化成为近年广受关注的热点问题[1-4]。考虑到水库和淤地坝拦沙是导致近20多年黄河沙量减少的重要原因[5-7]，而坝库拦沙量一般被视为“非可持续减沙量”[6]，因此本文重点关注黄河流域潼关以上黄土高原地区的现状产沙情势。

1 数据与方法

1.1 研究区范围本文研究范围包括黄河循化-青铜峡区间的黄土丘陵沟壑区(简称黄丘区)、内蒙古十大孔兑、黄河河口镇至龙门区间(简称河龙间)、汾河流域、渭河咸阳以上、泾河张家山以上和北洛河状头以上地区(图1)，总面积约39万km2。为论述方便，以下将其简称为“黄土高原”。不过，考虑到该区有大片水土流失轻微的土石山区、平原、风沙区和次生林区，故本文重点关注图1中的黄色区域，面积21.5万km2，称其为“黄河主要产沙区”。据1950—1969年实测数据推算，黄河主要产沙区的产沙量占黄土高原总产沙量的94%。

图1 研究区范围

1.2 数据采集与处理为真实反映林草植被和梯田状况，利用分辨率为30 m的遥感影像，配合大量野外样方核对，解译了1998年前后、2010年、2013年、2016年和2018年的土地利用和植被盖度；利用分辨率为2.1 m的遥感影像，提取了2012年和2017年梯田面积；利用分辨率为250 m的遥感影像，提取了2000—2021年植被盖度。

为科学描述植被梯田对流域水土流失的影响，本文采用刘晓燕等创建的易侵蚀区、林草有效覆盖率和林草梯田有效覆盖率等概念[8]。其中“易侵蚀区”是指流域剔除平原、城镇、河川地和石山区等轻微产沙区的其它区域，其面积用Ae表示；“林草有效覆盖率”指林草叶茎正投影面积占易侵蚀区的比例，计算方法详见文献[8]。将林草有效覆盖率Ve与“梯田面积占易侵蚀区面积的比例”之和称为林草梯田有效覆盖率，以综合反映林草植被和梯田对流域易侵蚀区的保护程度，记为Vet。

淤地坝拦沙所形成的坝地，也是影响流域产沙的关键因素。不过，至2019年，黄土高原坝地面积仅1016 km2、其中的60%集中在陕北，故未将其作为独立因素，而是计入梯田。

采用的降雨和水沙数据来自《中华人民共和国水文年鉴》和国家气象中心共享网站，包括研究区145个水文站设站至今的逐月水沙数据，黄河主要产沙区763个雨量站设站至今的逐日降雨数据；提取了日降雨大于25 mm的年降雨总量，简称有效降雨，用P25表示。

流域产沙量Ws是流域把口断面实测输沙量、流域内淤地坝和水库的拦沙量、灌溉引沙量和河道淤积量之和，计算公式为：

Ws=Sdischange+Ws-reservoir+Ws-warpingdam+Ws-irrigation+Wdeposit

(1)

式中：Sdischange为流域把口断面的实测输沙量；Ws-reservoir、Ws-warpingdam分别为流域把口断面以上区域水库和淤地坝的拦沙量；Ws-irrigation为流域把口断面以上区域自河道或水库取水的灌溉工程引沙量；Wdeposit为河道淤积量。

为定量比较不同流域在相同降雨情况下的产沙能力，引入产沙指数Si(t/(mm·km2))，它是指流域易侵蚀区内单位有效降雨在单位面积上的产沙量Ws，即：

(2)

1.3 产沙情势内涵及评价方法众所周知，通常所说的黄河天然时期沙量16亿t/a，是黄河陕县水文站1919—1959年实测输沙量的平均值。分析表明，1919—1959年黄土高原汛期降雨量平均为327 mm，与1919—2020年均值“329 mm”接近，即能大体反映长系列降雨情况。因此，黄河沙量“16亿t/a”一直是官方公布的天然时期黄河沙量特征值，相应时期的黄土高原年均产沙量(16.5亿t/a)反映了天然时期的产沙情势。

黄土高原现状产沙情势是指在不考虑坝库拦沙情况下，现状下垫面在长系列降雨情况下的多年平均产沙量。本文采用两种途径认识黄土高原的现状产沙情势，一是利用实际产沙量推算，二是理论推算。

途径1：基于现状年实际产沙量，推算黄土高原在长系列降雨情景下的年均产沙量。黄土高原2000年以来或2010年以来的年均产沙量，是认识现状产沙情势的重要途径，但存在明显弊端：一是该时段的降雨条件不仅与1919—1959年不同，也与人们预测的未来降雨条件不同；二是黄土高原各地的丰枯程度相差很大。为规避以上问题，采用以下评价方法：首先，获取各支流现状年的实测输沙量、水库和淤地坝拦沙量、灌溉工程引沙量，采用式(1)计算其年均产沙量Ws；第二，利用现状年实测降雨量和年均产沙量，采用式(2)计算各支流的现状产沙指数；第三，基于对黄土高原降雨变化趋势的认识，设计多种情景的长系列降雨；最后，基于现状产沙指数和设计降雨系列的有效降雨P25，利用式(2)计算相应降雨情景的年均产沙量，即Ws=Si×P25×Ae。

该方法的关键环节是，获取的坝库工程拦沙量、灌溉引沙量和河道冲淤量是否符合实际。

途径2：利用遥感水文统计模型，推算黄土高原在长系列降雨情景下的年均产沙量。研究发现[8]，在流域尺度上，林草有效覆盖率与产沙指数之间呈现明显的指数关系。进一步考虑不同规模梯田加入的影响，刘晓燕等创建了遥感水文统计模型[9]，用于评价流域在一定植被、地形和土壤条件下的产沙能力，基本公式如下：

Si=a×e-b×Vet

(3)

式中a和b是与雨强、地形和土壤有关的参数。为克服雨型和降雨落区(涉及地形和土壤的差异)的影响，构建该模型采用的数据多为2～3个雨强相似年份的均值，因此采用式(3)计算得到的产沙量单位应为t/a。此外，模型构建时采用的林草有效覆盖率，是基于分辨率为30～56 m遥感影像提取的信息，故应用该模型时，提取信息的遥感影像分辨率也应与之一致，以尽可能减少因分辨率不同而带来的植被信息差异。

采用该模型，首先将遥感影像提取的流域现状林草梯田有效覆盖率代入式(3)，计算相应的产沙指数Si；然后，将降雨量P25和流域的易侵蚀区面积Ae带入式(2)，即可得到流域的年均产沙量。

2 现状年界定及下垫面特征值

从认识黄河水沙变化角度，过去哪个时段的下垫面可以作为现状下垫面？这是一个看似简单的问题，但由于研究者掌握的数据资源、研究项目启动时间、成果服务对象等有所差异，因此，目前被称为现状下垫面的时段差异很大，包括2000—2012年、2007—2014年和2000—2016年等。

分析黄土高原各支流汛期含沙量和汛期径流系数变化过程可见(图2和图3)，2008年以来，绝大多数支流汛期含沙量和汛期径流系数逐渐趋稳。洪水含沙量和洪量是流域产沙量的函数，汛期含沙量和汛期径流系数趋于稳定，意味着流域产沙能力趋于稳定。鉴于此，并考虑时间系列的长度，本文将黄土高原2010年以来的下垫面称为现状下垫面。

图2 汛期含沙量变化

图3 汛期径流系数变化

分析图4可见，2012年以来黄土高原林草植被盖度的增幅明显变缓，2018年以来植被盖度没有增加。因此，采用分辨率为30 m的遥感影像，提取了2010年和2016年的林草地面积和易侵蚀区面积，取其平均值作为2010—2019年的代表值；取2010年、2013年、2016年和2018年4个典型年林草地植被盖度(遥感分辨率30 m)的平均值，作为2010—2019年林草植被盖度。将2018年的林草有效覆盖率，作为现状年林草有效覆盖率的最大值。进而计算了现状年不同时段的林草梯田有效覆盖率Vet，结果见表1。

图4 林草植被盖度变化(遥感影像分辨率250 m)

表1 现状年黄土高原林草梯田有效覆盖率(Vet)

3 降雨条件设计

为深入认识现状产沙情势，针对黄河主要产沙区，设计了4种降雨情景：

(1)情景A。分析表明，与1956—2019年平均值比较，2010—2019年黄河主要产沙区汛期降雨偏丰10.5%，有效降雨P25偏丰24.0%；有效降雨偏丰主要表现在中东部的河龙间和泾河流域，而六盘山以西大部分地区明显偏枯，见图5。情景A为2010—2019年实际降雨的放大情景，放大原则是：凡其有效降雨P25达到1966—2019年平均值的1.2倍者，采用2010—2019年实测降雨数据；凡不足1.2倍者，采用1966—2019年实测降雨量的1.2倍。该情景6—9月降雨量总体偏丰15.7%，其中有效降雨P25偏丰29%。

图5 2010—2019年有效降雨P25丰枯状况

(2)情景B。情景A虽消除了图5的降雨偏枯区域，但各地的丰枯程度不均，其中河龙间和泾河降雨更偏丰。为更合理地比较各支流的现状产沙能力，将1966—2019年实测降雨量P25直接放大1.28倍，作为情景B。情景B的汛期降雨偏丰12.5%，其中有效降雨P25偏丰23%。

(3)情景C。该情景是1956—2019年实测降雨系列，其平均降雨量与1919—2019年百年均值相近。

(4)情景D。该情景年均降雨量与1956—1967年相当，其汛期降雨偏丰17.3%，P25偏丰32%。

表2是各支流在设计降雨情景的有效降雨量(P25)，图6诠释了4种设计降雨情景的有效降雨(P25)丰枯程度，图中1919—1955年数据是利用李庆祥等构建的我国各地区1910—2009年汛期降雨网格化数据[10]、1956—2009年实测汛期降雨量和黄河陕县站1919—1955年实测输沙量等信息推算而成。不过，因1930年以前黄土高原只有太原、三门峡和平遥等3个雨量站，故其降雨数据可能与实际相差较大，仅供参考；1936年以后，数据的代表性逐渐提高。对比可见：

表2 各情景的有效降雨(P25)设计结果

图6 1919—2020年黄河主要产沙区有效降雨变化

(1)2010—2019年降雨量大体介于情景A和情景B之间，但前者在六盘山以西地区降雨偏枯，尤其是有利于产沙的暴雨偏少，见图5。值得注意的是，六盘山以西的祖厉河、清水河和马莲河上游，不仅是地形特殊、容易产沙的黄土丘陵区第5副区[9]，而且是目前林草梯田有效覆盖率最低的区域，见图7。

图7 2018年林草梯田有效覆盖率

(2)与1919—2019年百年均值相比，A、B、D三种情景的有效降雨分别偏丰29%、23%、32%，6—9月降雨偏丰15.7%、12.5%、17.3%。

4 结果分析

4.1 实际产沙情势通过对黄土高原所有水库和淤地坝地理位置、建成时间、控制面积、总库容和已淤积库容的深入调查和辨识，刘晓燕等计算了1950年代以来不同时期的坝库拦沙量和灌溉引沙量[7]。

河道冲淤和采砂也是影响流域产沙还原的重要因素[11]。调查表明，研究区采砂主要表现在河龙间和十大孔兑，其中河龙间黄河干流采砂主要集中在府谷至龙门区间，近20年平均约4000万t/a(粒径0.15～5 mm)；支流采砂主要在窟野河、皇甫川和十大孔兑等粗沙区河流。对比府谷至龙门区间“支流入黄沙量”与“干流上下断面的输沙量之差”发现(图8)：1996年以前，黄河干流该河段以“淤积为主”，但现在以“冲刷为主”，尤其是大量粗泥沙汇入的府谷至吴堡区间。1956—1996年，府谷-吴堡区间支流入黄沙量的15.4%(0.46亿t/a)会沉积在干流；2000—2020年，该区支流入黄沙量0.26亿t/a、黄河府谷以上流域来沙0.30亿t/a，而吴堡断面输沙量达0.94亿t/a，说明干流河道总体冲刷。故而推断，2000年以来，黄河干流采砂主要为前期淤积的特粗泥沙，对支流入黄泥沙向龙门以下输出影响极小。

图8 1956—2020年黄河府谷-龙门区间干流冲淤变化(正值表示淤积)

采砂对黄河减沙量的影响，极可能主要发生在窟野河、皇甫川和十大孔兑等粗泥沙区支流。实测数据表明，这些粗沙支流入黄断面悬移质粒径大于0.1 mm的粗泥沙比例，已由1980年代以前的20～40%降为5%左右，说明大量粗泥沙未输移至黄河干流。调研表明，近20年该区采砂量约2100万t/a。

基于以上分析，采用式(1)，还原得到黄土高原在不同时期的产沙量，见表3。由表3可见，潼关以上黄土高原2010—2019年实际产沙量为4.43亿t/a。

表3 2010—2019年黄土高原产沙量还原结果 单位：万t/a

基于表3的结果，利用同期实测降雨数据，可计算出2010—2019年各支流的产沙指数。进而，计算出该时期下垫面在不同设计降雨情景下黄土高原的可能产沙量，结果见表4。对应降雨情景A和降雨情景B，计算的黄土高原年均产沙量分别为4.60亿t/a和4.40亿t/a。

表4 2010—2019年下垫面在不同降雨情景下的产沙量计算结果 单位：万t/a

4.2 理论推算的产沙情势采用遥感水文统计模型，计算了2010—2019年下垫面在设计降雨情景的产沙量，结果见表4，表中兰州至青铜峡未控区间、苦水河、北洛河刘家河-状头区间、泾河景村-张家山区间、渭河拓石-咸阳区间和汾河兰村-河津区间，采用2010—2019年实际产沙量，合计0.43亿t。由表4可见，现状下垫面在四种设计降雨条件下的产沙量分别为4.77亿t/a、4.59亿t/a、3.90亿t/a、4.98亿t/a。

对比可见，基于2010—2019年下垫面背景，在情景A和情景B等两种降雨情况下，利用实际产沙量推算的黄土高原产沙量分别为4.60亿t/a、4.40亿t/a，而用遥感水文统计模型推算的理论产沙量分别为4.77亿t/a、4.59亿t/a，二者相差约4%，说明理论方法基本可靠。

4.3 现状产沙情势综合分析综上，黄土高原2010—2019年下垫面在不同降雨情景下的多年平均产沙量为3.9亿t/a～5.0亿t/a。毋庸置疑，“现状下垫面产沙量3.9亿t/a～5.0亿t/a”并不意味着极端暴雨年份没有大沙。1933年，黄河实测沙量39.1亿t，是过去百年的实测最大值、1919—1959年均值的2.44倍；1933年8月的特大暴雨落区有5个，其中3个位于泾河上中游和清水河部分地区，最大暴雨中心位于泾河西北部马莲河上游[12]。已有研究认为[12]，基于2010—2019年各支流的降雨量-输沙量关系，如果重现1933年汛期的全部降雨，并假定坝库不拦沙、不水毁，入黄沙量可能达到17亿t左右。利用遥感水文统计模型和1933年降雨数据，本文也粗略测算了1933年降雨重现在现状下垫面的可能产沙量，约12亿t～13亿t(不考虑1933年之前连续11年枯水期的重力侵蚀累积影响)。推算沙量如此之大的原因有三方面：一是该年的最大暴雨中心恰为极易产沙的黄土丘陵沟壑区第5副区，故1933年泾河张家山站输沙量高达11.7亿t、黄河上游青铜峡站达4亿t；二是最大暴雨落区目前几乎是研究区植被梯田覆盖最差的区域(图7)，2018年该区林草梯田有效覆盖率仅36%、相应产沙指数约49t/(mm·km2)；三是有效降雨大，黄河主要产沙区面平均有效降雨量P25约280 mm，达多年均值的2.2倍(其中泾河2.6倍)。

2010年以来，虽然其它省区的梯田增量很小，但甘肃省在黄土高原持续实施了大规模梯田建设。因此，2010—2019年的黄土高原产沙环境仍处于改善过程中，即该时期的下垫面一致性仍存在缺陷。为使研究更贴近当前产沙情势、并提高成果对认识未来产沙情势的参考价值，以下拟识别出下垫面一致性更强的“现状时段”，进一步诠释黄土高原的现状产沙情势。

文献[13]指出，由于水资源条件限制，黄土高原植被状况已趋稳定；基于云雾山自然保护区30多年的长期观测数据，程积民认为[14]，封禁15～20年，植被的生物量及盖度即达峰值。对比图4可见，以上结论与我们遥感跟踪监测的结果完全一致。未来，虽然林草植被的类型仍将继续演替，但受当地气候和土壤条件制约，可体现生物量的植被盖度很难继续增加。遥感调查表明，与1978年前后相比，黄河主要产沙区的2018年林草地面积仅增加3211 km2，其中，河龙间、北洛河上游和十大孔兑增加10 501 km2，其他区减少7290 km2，加之农村留守人口数量已稳定数年，因此未来退耕的潜力很小。由此可见，由于林草植被盖度和林草地面积的增加潜力极小，故林草有效覆盖率继续增长的潜力非常有限。

实地调查走访了解到，近十多年，虽然甘肃每年统计的梯田建设面积一直在增长，但在2018年以来的新增梯田面积中，大多为老旧梯田改造升级，这与前阶段情况差别很大。至“十四五”末，预计黄土高原研究区将新建梯田面积370万亩，由此可使黄河主要产沙区的林草梯田有效覆盖率增加0.8个百分点。之后，随着农村劳动力减少，对梯田的需求也会大幅降低。

由此可见，只要黄土高原气候不发生大波动，2018年以来的植被状况和梯田保有面积不仅最具一致性，而且基本反映其未来情景。因此，基于2018年林草梯田有效覆盖率(表2)，采用遥感水文统计模型，计算了4种设计降雨情景下的黄土高原产沙量，结果分别为4.31亿t/a、4.17亿t/a、3.55亿t/a、4.48亿t/a，其中最枯年份为1.6亿t、最丰年份为7.4亿t，该结果可基本反映黄土高原当前下垫面在长系列降雨情况下的产沙情势。

与天然时期黄土高原年均产沙量“16.5亿t/a”相比，黄土高原2010—2019年下垫面和2018年下垫面产沙能力分别降低约72%、75%，该减幅与采用不同技术路线得到的土壤侵蚀调查结论基本一致，见表5。

表5 历次土壤侵蚀遥感调查结果(本文界定的潼关以上黄土高原范围)

值得注意的是：(1)黄土高原产沙能力进一步降低的潜力很小。一方面，如前文所述，未来植被进一步改善和梯田大量增加的空间十分有限。另一方面，基于流域产沙与植被梯田覆盖状况的响应规律[8]，林草梯田有效覆盖率大于60%后流域产沙指数趋于稳定(图9)，而黄河主要产沙区的2018年林草梯田有效覆盖率已达60.3%。因此，“十四五”以后，若无遏制流域产沙的创新措施，该区产沙指数很难继续降低。

图9 流域尺度上植被梯田覆盖状况与产沙指数的关系

(2)黄土高原产沙量存在明显的反弹风险。如果扣除栽种在梯条田上的乔灌植物(产沙计算时，已将此类地块计入梯田)，目前大部分地区仍以草本植物为主、间有小灌木，抵御干旱和人畜干扰的能力较弱，尤其是降雨量300～450 mm的主要产沙区，经历了2021年春夏大旱后植被盖度明显下降，见图10。在现状人类活动强度下，研究区植被盖度主要取决于前期降雨条件，见图11。若经历1991～2002年连续干旱年，植被盖度可能会下降17%、退至2010年，随后黄土高原产沙量为5.5亿t/a～6.5亿t/a(梯田采用2021年数据)、至潼关5亿t/a～6亿t/a；若自然修复形成的草灌植被全部“牺牲”，黄土高原产沙量存在反弹至8亿t/a左右的风险。

图10 2000—2021年黄土高原不同气候带林草植被盖度变化

图11 现状人类活动强度下植被盖度与降雨量的关系

(3)潼关以上流域产沙量并不等于潼关断面输沙量。因水库和淤地坝等沟道工程拦截、灌溉引沙和河道冲淤调整等，能够入黄或输移至潼关断面的沙量，必然小于流域产沙量。正是由于一半左右的泥沙被层层拦截，2000年以来潼关年均输沙量只有2.41亿t/a。分析认为[7]，现状坝库仍具有一定的拦沙能力，“十四五”期间还将有4000余座淤地坝和泾河东庄水库等即将建成投运，故未来30年内潼关以上地区坝库拦沙作用仍可维持在1.6亿t/a～2.0亿t/a。由此可见，未来30年内，尽管黄土高原产沙量达3.6亿t/a～4.5亿t/a(若不补建新拦沙坝，假定植被维持2018年水平)，但潼关输沙量仍可能只有2亿t/a～3亿t/a。之后，随着东庄等现状坝库淤满，潼关沙量将逐渐反弹：考虑青藏高原产沙、灌溉引沙和河道冲淤，潼关输沙量最终可能与黄土高原的产沙量相当。

(4)未来潼关输沙量基本取决于黄土高原入黄沙量。青藏高原多年平均产沙量0.3亿t/a，天然时期也是潼关输沙量的组成部分，但因龙羊峡、李家峡、九甸峡和纳子峡等一系列大中型水库拦截，预计未来200年内很难有泥沙下泄。黄河青铜峡至三湖河口区间存在灌溉引沙和风沙入黄问题，但2010年以来的灌溉引沙量仅0.12亿t/a；而风沙粒径太粗，很难输移至潼关断面。

5 结论

(1)综合考虑影响流域产沙的植被和梯田因素变化、各支流实测洪水及其含沙量变化，认为黄土高原2010年以来的下垫面，可作为现状产沙情势评价的代表下垫面。尤其是2018年以来的黄土高原下垫面，其产沙环境更具一致性。

(2)利用1956—2019年实测降雨数据，设计了四种汛期降雨情景。通过实际产沙和理论产沙分析认为，如果汛期降雨较1956—2019年均值偏丰0～17.3%，黄土高原2010—2019年下垫面的多年平均产沙量为3.9亿t/a～5亿t/a。基于2018年以来下垫面水平的多年平均产沙量为3.6亿t/a～4.5亿t/a。

(3)综合分析未来植被梯田发展趋势、流域产沙与植被梯田覆盖程度的响应规律，认为黄土高原产沙能力进一步降低的潜力很小；而且，在降雨量450 mm以下的黄土高原，由于新生草被脆弱，未来存在产沙量反弹至5亿t/a～8亿t/a的风险。因此，服务于治黄重大工程布局的沙量设计必须有底线思维。

(4)由于已建和在建坝库拦截和灌溉引沙等因素影响，未来30年内，可输移至黄河潼关断面的沙量仍将远小于黄土高原产沙量。之后，随着已建和在建坝库淤满，若不补建新拦沙坝，潼关沙量将逐渐反弹。

由于未来潼关输沙量基本取决于黄土高原入黄沙量，故本文对黄土高原现状产沙情势的评价结果，对认识黄河来沙形势具有重要意义。