林草和梯田措施对小流域降雨—径流—输沙过程的影响分析

白 娟，张 亦 弛，杨 胜 天 ，刘 晓 燕，甘 甫 平，闫 柏 琨，郭 艺

(1.中国自然资源航空物探遥感中心，北京 100083；2.北京师范大学地理科学学部，北京100875；3.北京师范大学水科学研究院，北京100875；4.水利部黄河水利委员会，河南 郑州 450004)

0 引言

水土流失是当前世界面临的严重环境问题之一[1-4]，如何有效控制水土流失的发生和发展是流域管理的重点问题。坡改梯和林草种植是坡面易侵蚀区常用的水土保持措施，不仅可削减自身产水产沙，还可拦截上方坡面的来水来沙，对坡面乃至流域尺度的水沙过程具有显著影响[5-9]。但目前对流域尺度林草、梯田措施减水减沙作用的研究仍有待完善：1)野外实验观测方法有助于揭示坡面尺度林草、梯田措施减水减沙的过程和机理[4,10-13]，但局限于小尺度研究[14]，受下垫面空间异质性影响较大[15]；2)水保法仅考虑林草、梯田措施自身的减水减沙作用，忽略了其对上方坡面来水来沙的削减，难以代表流域全局范围内的水沙过程；3)水文法可得到下垫面变化导致减水减沙的整体效果，但无法给出单项措施的减水减沙贡献，且计算结果受建模期雨量数据的代表性影响[16]；4)水文模型法侧重于模拟坡面措施自身的减水减沙作用，忽视了对林草、梯田措施在汇流过程中拦水拦沙作用的参数化表达，且计算单一措施的减水减沙效率时，较难区分其他措施的影响，计算不同措施耦合的减水减沙效率时，会产生重复计算。

黄河中游位于黄土高原，水土流失严重，直接影响黄河的生态安全和区域社会经济可持续发展[17-22]。为控制该区水土流失，国家实施了大规模的退耕还林还草和梯田建设工程，定量评价林草、梯田的减水减沙作用符合流域水沙调控的发展要求，对黄河流域水土流失防治和生态保护具有重要意义。Luo等应用修正通用土壤流失方程(Modified Universal Soil Loss Equation，MUSLE)，耦合分布式场次暴雨径流LCM模型，构建了分布式LCM-MUSLE坡面水沙联动模型，并对MUSLE方程中的作物管理因子(C)、地形因子(LS)和保持措施因子(P)进行了修订，实现其在黄土高原应用的本土化，提高了产沙模拟精度[23]，但该模型忽略了汇流过程中坡面措施减水减沙作用；Bai 等对该模型的汇流方法进行改进，将等流时面内林草、梯田控制区域视为最小减水减沙单元，构建“控制区域—等流时面—子流域”一体化坡面汇流系统，考虑林草和梯田措施在汇流过程中对径流和输沙的削减作用，诠释了水保措施导致的下垫面变化对流域水沙的影响[24]。在此基础上，本研究以黄河中游多沙粗沙区为研究对象，分别应用原分布式LCM-MUSLE坡面水沙联动模型和汇流方法改进后模型进行场次“降雨—径流—输沙”过程模拟，定量评估林草、梯田及其耦合措施在流域坡面汇流过程中减水减沙的贡献率，并讨论其主要影响因素，为黄河中游林草植被和梯田建设对流域水沙影响的定量评估提供科学依据。

1 研究区与数据

1.1 研究区概况

黄河中游多沙粗沙区暴雨多、洪水大、含沙量高、泥沙颗粒粗，是黄土高原水土保持治理的重点区域[25-27]。为控制水土流失，该区开展了造林、种草、封禁、修梯田等坡面措施以及修淤地坝等沟道措施[16]，致使黄河中游下垫面特征发生显著变化[25]。自2000年黄河来沙量明显减少，2000-2019年潼关来沙量为2.45亿t/a，2000-2009年和2010-2019年的下垫面产沙能力相比1919-1959年均值分别下降了57%和80%[28]。综合考虑研究区代表性、已有数据积累以及林草、梯田面积和空间分布的巨大变化等因素，本文选取偏关河偏关水文站以上流域和清涧河子长水文站以上流域为研究区(图1)。

图1 研究区位置Fig.1 Location of the study area

偏关河和清涧河均为黄河一级支流，均位于黄土高原丘陵沟壑区第1副区。偏关河流域面积为2 089.36 km2，偏关水文站控制面积1 922.81 km2，海拔984～2 162 m，属温带半干旱气候，多年平均降水量为429 mm，5-9月降水量占全年降水量的80.9%，产流类型为超渗产流，多年平均径流量为3 948万m3，多年平均输沙量为1 258万t，洪水输沙模数为6 523 t/(km2·a)，水土流失严重[29]。清涧河流域面积约4 080 km2，其子长水文站控制面积约889.65 km2，海拔1 047～1 585 m，属暖温带半干旱气候，流域多年平均降水量为482.4 mm，汛期6-9月降水量占年降水量的70%以上，产流类型为超渗产流，多年平均径流量为13 700万m3，多年平均输沙量为3 300万t，汛期输沙量占90%以上[30]。

1.2 数据来源

研究数据包括：1)流域内雨量站观测的逐小时降雨数据以及水文卡口站观测的径流量和含沙量数据，对其进行时间重采样和空间插值[31]；2)流域泥沙粒径数据(D50)，由水利部门公布的各水文站粒径数据经反距离加权插值后获取；3)ASTER GDEM 数据，用于提取地形指数、等流时线、河网和子流域等空间信息；4)土地利用类型和植被覆盖度数据，基于Landsat影像获取，前者参照《土地利用现状分类(GB/T21010-2007)》标准，后者基于归一化植被指数(NDVI)，采用像元二分模型提取得到，研究中分别用1978年、2010年的土地利用和植被覆盖度数据代表1980s和2010s的土地利用和植被覆盖度状况；5)梯田数据，采用水利部黄河水利委员会(YRCC)基于2012年ZY-3卫星影像解译的结果；6)土壤数据，包括土壤类型、土壤机械组成、土壤有机碳含量和土壤饱和含水量，土壤类型来源于中国1∶100万土壤图，土壤机械组成和有机碳含量等属性通过查询世界土壤数据库HWSD (Harmonized World Soil Database version 1.1)获取，并运用SPAW(Soil-Plant-Atmosphere- Water Field & Pond Hydrology)模型计算各土壤类型的饱和含水量。

2 研究方法

2.1 改进的LCM-MUSLE坡面水沙联动模型

图2 改进的LCM-MUSLE模型结构Fig.2 Framework of modified LCM-MUSLE model

(1)

(2)

Vui=ΔAVcontrol,i/ΔAi

(3)

(4)

ΔTi,j=Ti,j-Ti,j-1

(5)

(6)

Tui=ΔATcontrol,i/ΔAi

(7)

Tci=ΔAterrace,i×di

(8)

采用等流时线法计算各子流域坡面汇流后，利用马斯京根方法计算河道汇流[24]，汇沙过程还需考虑子流域泥沙颗粒大小和流域水力学特性。通过等流时线法计算各子流域的产沙量后，根据汇沙模型[36]计算河道汇沙过程前子流域i的输沙量SYi：

(9)

式中：Yi为等流时线法坡面汇沙后计算出的子流域产沙量(t)；B为汇沙系数；Ti为子流域到流域出口的汇流时间(h)；D50i为子流域泥沙粒径中值(mm)。

2.2 模型精度评价

基于小时尺度径流量、输沙量观测数据，运用纳西效率系数(Nash)对模型模拟结果进行精度评价，并分析其在研究区的适用性，计算公式为：

(10)

2.3 林草、梯田措施减水减沙作用评估

林草、梯田措施的减水减沙作用是指在相同降雨条件下实施措施后流域较天然时期减少的产水产沙量[37]。将流域植被基本稳定、水保措施相对较少、人类活动影响较小的1980s作为天然时期，将流域植被恢复明显、水保措施建设密集、人类活动影响显著的2010s作为现状年。根据研究区在不同时期有无林草、梯田措施，提出O1(原模型)、R1(林草措施)、R2(梯田措施)、R3(林草+梯田措施)4种模拟情景。由于1980s的梯田数据不完整，且数据质量较差，故对于1980s的场次降雨仅模拟O1和R1情景。

通过构建削减率(RR)和单位面积削减量(AR)两个指标，定量分析林草、梯田在汇流过程中的减水减沙作用。以径流为例，将场次小时尺度径流模拟值加总，获取场次径流总量，然后统计场次降雨在4种情景下的模拟径流总量，分别计算R1、R2和R3相对O1的径流削减率(式(11))，以及R1和R2的单位面积径流削减量(式(12))；同理可计算R1、R2和R3相对O1的输沙削减率，以及R1和R2的单位面积输沙削减量。

RRi=(y-yi)/y×100%

(11)

ARi=(y-yi)/Ai

(12)

式中：RRi、ARi分别为流域场次降雨在i(i=1,2,3)措施下的径流削减率(%)和单位面积径流削减量；y为O1情景的模拟径流总量(m3)；yi为R1、R2 或 R3情景的模拟径流总量(m3)；A为流域林草或梯田的面积(km2)。

3 结果分析

3.1 场次“降雨—径流—输沙”模拟与精度验证

经筛选，1981-2010年两流域均有8个场次“降雨—径流—输沙”观测数据用于模型模拟和精度验证，每场次降雨以其起始时间(年/日/时)进行编号。偏关站以上流域以1981/203/17、1983/215/22和1983/235/16共3个场次作为参数率定期，1988/199/13、1989/203/19、2006/195/05、2006/224/08和2010/263/20共5个场次作为验证期；子长站以上流域以1986/187/14、1987/238/00和1988/218/14共3个场次作为参数率定期，1988/197/04、1988/237/18、2006/188/20、2006/237/02和2006/263/20共5个场次作为验证期。

由两流域径流量(图3)和输沙量(图4)模拟值与观测值对比可知：总体上，对于1980s的各场次降雨，R1的峰值流量和径流总量比O1均有明显削减，对于2010s的各场次降雨，径流量峰值排序为O1>R2>R1>R3，表明林草措施对峰值流量和径流总量的削减作用大于梯田措施，综合考虑两种措施对峰值流量和径流总量的削减作用最强。对于场次输沙量而言，2010s O1与R1两种情景下输沙量的差异小于二者径流量的差异，这是由于在MUSLE产沙模型中，Rs、C和P因子的设定已体现对产沙过程的削减作用。对于1980s的各场次降雨，R1的输沙量峰值均低于O1的峰值；对于2010s的各场次降雨，输沙量峰值排序为O1>R2>R1>R3，表明林草措施对输沙量峰值削减幅度比梯田措施大，两种措施耦合对输沙量的削减作用最强。

图3 偏关站以上流域和子长站以上流域径流量模拟值与观测值对比Fig.3 Comparison of observed and simulated runoff in upstream watershed of Pianguan station and upstream watershed of Zichang station

图4 偏关站以上流域和子长站以上流域输沙量模拟值与观测值对比Fig.4 Comparison of observed and simulated sediment discharge in upstream watershed of Pianguan station and upstream watershed of Zichang station

由两流域场次“降雨—径流—输沙”模拟精度(表1)可知，偏关站以上流域在1980s和2010s的场次径流量平均Nash系数分别从O情景的0.46、-15.29 提高到R3情景的0.62、0.39，场次输沙量平均Nash系数分别从O情景的0.08、-2.47提高到R情景的0.22、0.37；子长站以上流域在1980s和2010s的场次径流量平均Nash系数分别从O情景的-0.86、-5.28提高到R3情景的0.11、0.43，场次输沙量平均Nash系数分别从O情景的-1.40、-0.34提高到R3情景的-0.24、0.34，表明原模型在现状年模拟精度不高，改进后模型适用于现状年的径流(输沙)量模拟。

表1 两流域场次径流量/输沙量平均Nash系数Table 1 Average Nash coefficients of runoff and sediment discharge in upstream watershed of Pianguan station and upstream watershed of Zichang station

3.2 林草、梯田措施减水减沙作用分析

统计1980-2010年各场次降雨在不同情景下的径流(输沙)量，计算林草措施、梯田措施及二者耦合的径流(输沙)削减率，以及林草、梯田措施的单位面积径流(输沙)削减量(图5、图6)。可以看出：

图5 林草、梯田措施径流削减作用Fig.5 Runoff reduction effects of vegetation and terraces

图6 林草、梯田措施输沙削减作用Fig.6 Sediment reduction effects of vegetation and terraces

对于偏关站以上流域而言，1980s场次降雨林草措施的平均径流削减率为17.86%，平均单位面积径流削减量为1 378.84 m3/km2；2010s场次降雨林草措施的平均径流削减率RR1为48.08%，受植被覆盖度增加影响，比1980s增加了30.22%，平均单位面积径流削减量为4 063.59 m3/km2；梯田措施的平均径流削减率RR2为26.64%(小于林草措施)，平均单位面积径流削减量为15 827.36 m3/km2(远大于林草措施)；林草梯田耦合措施的平均径流削减率RR3为69.26%。1980s场次降雨林草措施的平均输沙削减率为12.02%，平均单位面积输沙削减量为341.19 t/km2；2010s场次降雨林草措施的平均输沙削减率RR1为32.59%，受植被覆盖度增加影响，比1980s增加了20.57%，平均单位面积输沙削减量为663.61 t/km2；梯田措施的平均输沙削减率RR2为24.50%(小于林草措施)，平均单位面积输沙削减量为3 384.71 t/km2(远大于林草措施)；林草梯田耦合措施的平均输沙削减率RR3为54.18%。

对于子长站以上流域而言，1980s场次降雨林草措施的平均径流削减率为21.69%，平均单位面积径流削减量为3 867.98 m3/km2；2010s场次降雨林草措施的平均径流削减率RR1为61.01%，受植被覆盖度增加影响，比1980s增加了39.32%，平均单位面积径流削减量为7 594.16 m3/km2；梯田措施的平均径流削减率RR2为8.55%(小于林草措施)，平均单位面积径流削减量为21 638.68 m3/km2(远大于林草措施)；林草梯田耦合措施的平均径流削减率RR3为65.54%。1980s场次降雨林草措施的平均输沙削减率为19.46%，平均单位面积输沙削减量为1 423.69 t/km2；2010s场次降雨林草措施的平均输沙削减率RR1为43.71%，受植被覆盖度增加影响，比1980s增加了24.25%，平均单位面积输沙削减量为541.07 t/km2；梯田措施的平均输沙削减率RR2为13.94%(小于林草措施)，单位面积输沙削减量为1 525.85 t/km2(远大于林草措施)；林草梯田耦合措施的平均输沙削减率RR3为55.84%。

总体上，两流域径流(输沙)削减率排序均为林草+梯田措施>林草措施>梯田措施，且林草和梯田措施的径流(输沙)削减率之和略大于二者耦合结果，可能是当梯田起到减水作用后，同一等流时面内植被的径流削减率会降低。此外，本文假设偏关站以上流域均为一类梯田，子长站以上流域梯田均无田埂，可能会影响梯田措施减水作用的估算结果。

3.3 林草、梯田减水减沙影响因素分析

流域林草植被面积和植被覆盖度变化直接影响流域场次径流量、输沙量，林草植被的空间分布也会对汇流过程中林草的减水减沙作用产生影响。为分析流域场次降雨径流(输沙)量变化与林草植被数量及其空间分布的关系，分别计算林草面积比、植被覆盖度和林草减水面积比3个指标与各场次降雨的林草径流削减率、输沙削减率的Pearson相关系数。其中，林草面积比为流域内林草地面积占易侵蚀面积比例，林草减水面积比为流域内林草减水面积(即林草面积及其集水区面积之和)占易侵蚀面积比例。结果显示：在α=0.05的显著性水平下，林草径流削减率与林草面积比、植被覆盖度和林草减水面积比的相关系数分别为0.81、0.80和0.34，林草输沙削减率与上述3个指标的相关系数分别为0.76、0.71和0.48，说明流域林草径流(输沙)削减率与植被数量明显相关，与植被空间分布的相关性较弱。

流域梯田面积和梯田田埂完整度会直接影响流域场次降雨径流(输沙)量，梯田的空间分布(用梯田集水区与流域面积比间接表征)也会对汇流过程中梯田的减水减沙作用产生影响。为分析流域场次降雨径流(输沙)量变化与梯田集水量及梯田空间分布的关系，分别计算梯田比、田埂完整度、梯田集水区与流域面积比3个指标与各场次降雨的梯田径流(输沙)削减率的Pearson相关系数。其中，梯田比为某地区水平梯田面积占区内轻度以上水土流失面积比例，该指标比梯田面积与流域面积比更能反映梯田对流域主要产沙区的控制程度[37,38]；研究区梯田均为水平梯田，故田埂完整度可反映梯田集水量。结果显示：在α=0.05的显著性水平下，梯田径流削减率与梯田比、田埂完整度、梯田集水区与流域面积比的相关系数分别为0.98、0.56和0.94，梯田输沙削减率与上述3个指标的相关系数分别为0.98、0.44和0.97，说明流域梯田径流(输沙)削减率与梯田面积和梯田空间分布明显相关，与梯田田埂完整度的相关性较弱。

4 结论

本文分别采用改进的分布式LCM-MUSLE坡面水沙联动模型和原模型在黄河中游偏关站以上流域和子长站以上流域开展场次降雨水沙过程模拟，计算林草、梯田及二者耦合措施的减水减沙作用并讨论其主要影响因素。结论如下：1)改进后的分布式LCM-MUSLE坡面水沙联动模型更适用于现状年的径流(输沙)量模拟，可反映下垫面显著变化(坡改梯和植被恢复)对径流(输沙)的影响；2)单位面积梯田的减水减沙量远大于林草，林草措施总的径流(输沙)削减率显著大于梯田措施；3)流域尺度梯田措施下径流(输沙)削减率影响因素排序为梯田比>梯田集水区与流域面积比>田埂完整度，林草措施下影响因素排序为林草面积比>植被覆盖度>林草减水面积比。

本文一定程度上剖析了近年来黄河中游来水来沙量锐减的原因，建立了径流量、输沙量对于梯田、林草变化的响应机制，但仍存在如下不足：文中田埂完整度的设定来源于已发表文献，与田埂真实完整程度可能存在偏差，造成径流、泥沙削减率略高于真实值；未考虑水库、淤地坝等工程措施以及开矿、修路等土地利用开发对河道及河床周围水文过程的影响。未来将开展下垫面变化对河道及河床周围水文过程的影响研究，以期更全面地刻画多沙粗沙区“降雨—径流—输沙”过程。