联合径流侵蚀功率与连通性指数识别流域侵蚀对植被恢复的响应

2023-11-09郭宗俊吴磊张慧勇刘帅

农业环境科学学报 2023年10期

郭宗俊，吴磊，张慧勇，刘帅

（1.西北农林科技大学旱区农业水土工程教育部重点实验室，陕西杨凌 712100；2.西北农林科技大学黄土高原土壤侵蚀与旱地农业国家重点实验室，陕西杨凌 712100；3.西北农林科技大学水利与建筑工程学院，陕西杨凌 712100）

流域侵蚀是最受关注的环境问题之一[1-2]，大多数研究运用模型方法探究侵蚀及泥沙输移过程的规律[3-4]。其中，RUSLE（Revised Universal Soil Loss Equation）模型被广泛应用于模拟流域多年侵蚀情况，且已在黄土高原较多流域中得到检验；而在泥沙输移的研究中，Borselli 等[5]提出了泥沙连通性指数（Index of Connectivity，IC），该指数可用于识别泥沙连通路径并评估坡面泥沙输移的可能性；径流侵蚀功率则主要描述降雨径流对流域侵蚀过程的驱动，并在不同流域中得到运用[6-7]，其代表了径流能量对泥沙输移的影响，能够明确流域不同位置的侵蚀潜能差异。IC 和土壤侵蚀模型是研究流域泥沙输移的有效工具[8]，两者的结合可以克服仅依靠侵蚀模型而不能充分表达输沙潜力时空特征的缺点[9]，同时借助径流侵蚀功率能够表达径流所造成流域侵蚀的能力，将三者相结合来确定侵蚀泥沙来源、范围和流域内关键的侵蚀产沙区域以及分析泥沙连通性和侵蚀之间的联系是探索流域侵蚀产沙模式的关键[10]，因此，本研究将结合RUSLE、IC 和基于SWAT（Soil and Water Assessment Tool）模型的径流侵蚀功率开展。

植被恢复改变了流域的下垫面情况，进而改变了流域土壤侵蚀和泥沙输移的空间格局[11]。黄河流域经过长时间的治理和恢复，包括退耕还林还草、粮食换绿和三北防护林等措施与计划的实施[12-13]，其生态环境明显好转，流域侵蚀和泥沙输移的环境发生了较大变化。以往的一些研究评估了植被恢复所引起的土壤水文性质响应[14]，也明确了植被状况对泥沙的形成和迁移具有显著影响[15-17]。然而，就流域尺度而言，植被恢复如何影响流域侵蚀和泥沙连通性以及它们之间的关系仍不够明确，识别流域中侵蚀和泥沙输移过程对植被变化的响应是有必要的，可有效解释流域尺度中水文要素与植被的关联[18-20]。

延河流域水土流失严重，以往在延河流域开展的水沙研究十分广泛[21-23]，但对该流域侵蚀产沙过程的整体识别和探究还不充分；另外，由于近年来开展的生态措施，延河流域植被条件持续向好，探究其植被变化所产生的影响可为流域进一步管理提供建议。因此，本研究结合RUSLE、IC 和基于SWAT 模型的径流侵蚀功率方程，建立了综合模型方法框架，进行长序列的侵蚀过程模拟以识别流域的侵蚀过程及其规律，并评估植被变化对土壤侵蚀和泥沙输移的影响，本研究的方法和结果可为进一步的生态措施实施和生态文明建设提供理论依据。

1 材料与方法

1.1 研究区概况及数据来源

延河位于黄河流域中部地区（图1），是黄河的一级支流，也是黄河主要泥沙来源之一，流域总面积为7 725 km2，其中黄土丘陵沟壑区占流域的94.6%以上，主要土壤类型为黄绵土（约85%），还包括部分冲积土、黑垆土、黏土。流域处于半干旱地区，属于大陆性季风气候，夏季气温较高，降雨较多，冬季寒冷干燥，降水稀少[24]。流域多年平均降水量为435 mm，多集中于6—9 月[25]。流域内沟道密布，地形破碎，梁峁起伏，沟壑密度在2.1～4.6 km·km-2之间[26]。流域内植被具有明显的地带过渡性，从南向北依次为森林区、森林草原区和草原区。由于恶劣的自然条件和不合理的耕作方式，延河流域内土壤侵蚀严重。目前，流域内开展了大规模水土保持工作，包括退耕还林还草、淤地坝等措施，侵蚀问题得到了一定的控制。

图1 延河流域DEM及子流域划分Figure 1 The DEM and subbasin division of Yanhe watershed

延河子流域划分见图1。研究区数据包括数字高程数据（DEM）、100 m 分辨率土地利用数据、1∶10万土壤数据、250 m 分辨率NDVI（Normalized Difference Vegetation Index）数据、气象数据包括多站点1985—2020 年日数据（降雨、气温、风速、日照），以及甘谷驿水文站点的实测日流量数据，其详细来源可见表1。土地利用类型包括耕地、林地、草地、水域、城镇和裸地，气象数据用于构建SWAT 模型的气象数据库以及RUSLE 模型的降雨侵蚀力空间分布，NDVI 则作为重要的影响因素被考虑在各模型的构建中。

表1 数据及来源Table 1 Date and source

综合模型方法框架包括RUSLE、IC和基于SWAT的径流侵蚀功率方程，其构建流程和分析思路见图2，旨在识别流域侵蚀过程中的两个重要关系：（1）泥沙连通性和径流侵蚀功率与侵蚀的变化关系；（2）泥沙连通性和径流侵蚀功率与植被的互联关系。其中，借助嵌套流域理论对泥沙连通性和径流侵蚀功率进行了表达。嵌套流域理论是在流域空间分级基础上，体现水文过程尺度效应的重新划分。本研究借助SWAT 模型划分了水文响应单元和子流域，并根据子流域分布进行了嵌套流域的划分，以研究径流侵蚀功率和泥沙连通性的尺度效应。

图2 模型方法框架流程图Figure 2 Flow chart of the model methodology framework

1.2 基于RUSLE 模型的侵蚀计算方法

采用通用土壤流失方程RUSLE 计算流域年侵蚀模数，其表达式如下：

式中：A是土壤侵蚀模数，t·hm-2·a-1；R是降雨侵蚀力因子，MJ·mm·hm-2·h-1·a-1；K是土壤可蚀性因子，t·hm2·h·hm-2·MJ-1·mm-1；L是坡长因子，无量纲；S是坡度因子，无量纲；C是植被覆盖与管理因子，无量纲；P是保持措施因子，无量纲。

（1）降雨侵蚀力（R）采用侵蚀性日降雨的计算方法[27]，该方法在黄土高原运用广泛。包括延河流域内和周边共12 个气象站点的日降雨数据被用于获取逐站点R值，并进行整个流域的空间插值得到逐年的R值栅格图层。

（2）土壤可蚀性因子（K）用于评价土壤受侵蚀的难易程度，该因子的大小与土壤质地、内部结构、有机质含量等多因素有关。K值根据各类型土壤的不同组分进行计算[28]。

（3）坡度和坡长（LS）是用于衡量地形对土壤侵蚀影响的重要指标，利用流域的DEM 数据可以获取坡度因子和坡长因子[29]。

（4）植被覆盖因子（C）的计算方法很多，包括土地利用赋值法、标准小区法和遥感数据法等。本研究采用了基于NDVI 的遥感数据法，首先由NDVI 值计算出流域的植被覆盖度（f），再根据植被覆盖度计算C值[30]，其公式如下：

式中：Ni为第i个单元的NDVI 值；Nmax为流域NDVI 的最大值，一般指纯植被像元的NDVI 值；Nmin为流域NDVI的最小值，一般指全裸土像元的NDVI值。

（5）水土保持措施因子（P）根据土地利用类型赋值。水体和城镇的P值设定为1 和0.01，森林和草地的P值分别设定为0.7 和1；在耕地中，通常认为坡度越大，保持措施效果越好，根据已有研究对不同坡度的耕地进行了P因子赋值[31]，范围在0.1～0.8之间。

1.3 基于IC的泥沙连通性计算方法

IC 被定义为单位量的泥沙通过坡沟系统从流域上部坡面到达河道的潜能，能够反映泥沙路径的连续性和连通程度[32]，目前该指数已被开发为可靠的指标，用于表达泥沙从流域上游坡面到达下游河道的可能性，IC（IC）的计算公式如下：

式中：Dup是上坡组分；Ddn是下坡组分是上坡区域的平均权重因子；-S是上坡区域的平均坡度；A是上坡区域面积；di是第i个单元沿着汇流路径到出口的距离；Wi是第i个单元的阻抗因子；Si是第i个单元的坡度。其中，W因子被考虑为流域逐年的C因子，从而表达植被变化对IC 的影响；汇流路径的出口设置为逐个嵌套的流域出口，以表达不同嵌套流域的IC。

泥沙输移比（Sediment Delivery Ratio，SDR）作为流域泥沙研究的重要指标，可由IC 计算，其表达式可以呈现整个流域的泥沙输移比情况[33]：

式中：SDRi是第i个单元的泥沙输移比；SDRmax是流域理论上最大的泥沙输移比，通常取为1；IC0和k是计算泥沙输移比的修正系数。

1.4 基于SWAT模型的径流侵蚀功率计算方法

本研究构建了7 个基于不同年代土地利用的SWAT 模型（1985、1990、1995、2000、2005、2010、2015年，整体研究时段为1985—2020），逐个模型分别模拟并率定5 年，如1985 年土地利用模型模拟了1985—1989 年子流域的月径流情况，然后对每个子流域进行径流侵蚀功率的计算。计算方法参考龚珺夫等[34]开发的径流侵蚀功率计算方法，将子流域的年径流过程概化为以月为时间步长的暴雨，其公式如下：

式中：Eyear为年径流侵蚀功率，m4·s-1·km-2；Q′m，year为最大月流量模数，m3·s-1·km-2；Qm，year为年内最大月流量，m3·s-1；A′为子流域的控制面积，km2，即上游至下游的嵌套流域面积；Hyear为年平均径流深，m；Qyear为年平均流量，m3·s-1。

2 结果与分析

2.1 流域径流侵蚀功率时空分布

7个基于不同年份土地利用的SWAT模型的率定方式基于SWAT-CUP 中的SUFI-2 算法，率定的参数则通过敏感性分析共选取12 个相对敏感的参数，其敏感性排序如下：径流曲线数（CN2）、土壤有效含水量（SOL_AWC）、地表径流滞后系数（SURLAG）、植物蒸发补偿系数（EPCO）、土壤蒸发补偿系数（ESCO）、基流α 因子（ALPHA_BF）、主河道水力传导系数（CH_K2）、主河道的曼宁系数（CH_N2）、地下水延迟系数（GW_DELAY）、初始土壤蓄水量（FFCB）、浅层含水层初始水深（SHALLST）、土壤饱和导水率（SOL_BD）。模型都通过了率定要求（表2），并可以相互验证，模拟月值和实测值的对比见图3。

表2 模拟结果评价Table 2 Evaluation of simulation results

图3 月径流模拟值与实测值Figure 3 Monthly runoff simulation value and measured value

在建模过程中，子流域划分采取了相同的模式，只针对土地利用进行了不同的建模，最终流域被划分为50 个子流域，模型经过率定符合要求后，通过模型获取流域内各个子流域的月径流量，最后依据公式（6）～公式（8）计算出其年径流侵蚀功率。其中，依据面积大小共选取了12级嵌套流域进行分析（表3），径流侵蚀功率是基于这些嵌套流域进行计算，其计算过程是将每一级的嵌套流域作为一个单独的流域出口进行分析，将其上游的径流纳入计算，最后在ArcGIS平台绘制出逐年径流侵蚀功率空间分布图。

表3 嵌套流域划分Table 3 Nested watershed division

年径流侵蚀功率具有显著的时空变异性（图4）。年径流侵蚀功率在时间尺度上表现为整体的下降趋势，1985 年最大子流域径流侵蚀功率为13.28×10-4m4·s-1·km-2，2000 年最大径流侵蚀功率下降到10.06×10-4m4·s-1·km-2，但2005 年侵蚀功率空间分布图显示子流域侵蚀功率有一定的增大，为19.36×10-4m4·s-1·km-2，而到2020 年最大径流侵蚀功率仅有4.40×10-4m4·s-1·km-2。

图4 子流域侵蚀功率的时空变化Figure 4 Spatiotemporal variation of erosion power in subbasins

在空间尺度上，年径流侵蚀功率呈现出支流大而干流小的特点（图4），上游子流域的侵蚀功率最大，并且干流子流域的侵蚀功率很小，而支流子流域的侵蚀功率很大，即使它们位于中游或下游地区。流域内上游至下游一些子流域的侵蚀功率年际变化见图5，位于上游的7号和14号子流域侵蚀功率较大，其次是中游地区的31 号和40 号子流域，流域下游的43 号和50 号子流域侵蚀功率非常小，表现出上游到下游逐渐减小的规律。同时，侵蚀功率随着子流域控制面积的增大而减小（图6），两者满足幂函数关系，说明在一定范围内的汇流面积对径流侵蚀功率有较大影响，但超过特定阈值后，这种影响快速降低，如图显示延河流域的面积阈值应在1 000 km2以内。

图5 子流域侵蚀功率年际变化趋势Figure 5 Annual variation trend of erosion power in subbasins

图6 子流域多年平均侵蚀功率与控制面积的关系Figure 6 Relationship between average erosion power and control area in subbasins

2.2 嵌套流域IC空间分布

IC 的空间分布呈现出“坡面小沟谷大”的特点（图7），由于泥沙的传递路径在流域内的各处不尽相同，远离河流的坡面泥沙需要经过长距离输移，泥沙输移较困难，使得其IC 较低，减少了产沙的发生；在靠近河流的沟谷地区，泥沙很容易被输送，其区域内的IC 较高。整个流域的IC 变化范围非常大（-13.11～1.95），表现为中游较小、上游和下游较大的空间格局。IC 的计算考虑了侵蚀泥沙到达流域出口的输移距离，而即使针对不同的流域出口，位于上游的嵌套流域总是具有较大的IC，最大的IC 出现在嵌套流域2 中，说明上游地区的侵蚀泥沙较容易到达河道，随后通过水文过程输出流域形成产沙。因此，综合IC的分布和大小程度而言，嵌套流域6以上的区域应该是延河流域侵蚀产沙防治的重点。最小的IC 出现在嵌套流域7中（图7），位于中游地区；而下游地区的IC 仍然较大，因此需要根据流域不同区域的特征进行不同泥沙管理措施。

2.3 嵌套流域侵蚀规律

为更明确识别径流侵蚀功率和泥沙连通性对流域侵蚀的影响，绘制了嵌套流域的多关系曲线见图8。首先，在嵌套流域中，面积的尺度效应对径流侵蚀功率的影响非常大，侵蚀功率随嵌套流域控制面积的增大而减小，两者满足幂函数关系（R2=0.82），在小于某个阈值时，侵蚀功率较大，因此需要关注小尺度中的径流事件以防止严重侵蚀；同时IC 与嵌套流域控制面积也满足较好的线性相关性（R2=0.84），即流域的控制面积越大，IC 值通常越小，上游的泥沙越难以离开该流域，但这种变化并不剧烈。此外，IC 与侵蚀保持了更强烈的线性相关性（图8），即使这种关系以平均情况描述。我们可以得到一种观点是，高侵蚀现象通常伴随着高泥沙连通性。在实地研究中发现[32]，高连通性区域观察到泥石流，并贡献37%以上的侵蚀泥沙，而低连通性区域则成为泥沙的沉积区，识别这些区域对流域泥沙管理可能具有关键作用。另外，即使径流侵蚀功率和侵蚀过程间并没有表现出强烈的相关性，但径流仍然是流域泥沙运动的重要动能和媒介，未来研究可能需要关注更细致时间尺度或事件尺度的径流-侵蚀过程的研究，以揭示其过程的动态规律。

图8 径流侵蚀功率和泥沙连通性的尺度效应及与侵蚀模数的相关关系Figure 8 Scale effect of runoff erosion power and sediment connectivity and their correlation with erosion modulus

3 讨论

3.1 植被变化对侵蚀的影响

总体而言，1985—2020 年间流域的土地利用空间格局有较大变化（表4），包括耕地向林地、草地的大量转化。耕地占比在2000 年之前稳定在43%左右，而在2010 年前后快速减少；同时期内，草地占比稳定在45%，但在2010 后迅速增加到52%，面积增加高达53 480 hm2；林地则呈现出整体的上升趋势，这主要得益于退耕还草、植树造林等生态工程的开展。流域的NDVI 呈现总体上升趋势（图9a1～图9a6），但在2000 年前变化较小，维持在0.30～0.35 之间；而2000 年至今，NDVI 显著上升，达0.70 左右（图9b），表明流域植被在快速增加。

表4 各年份土地利用情况（%）Table 4 Land use in different years（%）

流域植被的持续恢复改变了侵蚀发生的环境，流域侵蚀情况得到控制且持续向好。模拟结果表明，在空间尺度上（图10），高侵蚀地区集中在上游，即使在多年平均条件下，仍然有9 个上游子流域侵蚀模数超过80 t·hm-2·a-1，属于极强烈侵蚀区，并且2 号子流域的侵蚀模数高达238 t·hm-2·a-1，属于剧烈侵蚀区；低侵蚀地区则集中在流域中南部。在时间尺度上，流域的侵蚀模数逐步减少，2000 年前流域的平均侵蚀模数较大，多个年份的侵蚀模数超过80 t·hm-2·a-1，最高可达144 t·hm-2·a-1（1985 年），而在2015 年左右大幅下降至10～30 t·hm-2·a-1（图11），包括一些早期为高侵蚀区的子流域其侵蚀模数下降率高达80%以上。此外，侵蚀和NDVI 在时间上呈强烈的负相关趋势（图11），两者的拟合关系较明确（R2=0.641 9），植被增加可为侵蚀的衰减提供合理的解释。另外，一些研究认为气候因素特别是降雨的减少是黄河流域近年来泥沙减少的重要原因[35-36]，但在构建延河流域RUSLE 模型时，我们发现降雨侵蚀力变化并不大，近10 年的降雨侵蚀力水平与研究时段的早期水平相当。在辽西地区的研究中发现对侵蚀影响最大的是地表径流，其次是降雨侵蚀力，最后是植被覆盖；但在降雨一定的前提下，植被则是最关键的因素[37]，因此植被恢复仍可作为流域泥沙管理的重要手段。

图10 子流域多年平均侵蚀模数Figure 10 Average annual erosion modulus of subbasins

图11 1985—2020流域平均侵蚀模数和NDVI年际变化趋势Figure 11 Annual variation trend of average erosion modulus and NDVI of watershed from 1985 to 2020

3.2 植被变化对泥沙输移比的影响

IC 和NDVI 具备更强的相关性，其线性相关系数高达0.98（图12）。本研究使用的IC 属于结构泥沙连通性，其计算式表明地形因素决定了IC 的上限，即考虑为泥沙输移主要受到地形因素的限制[38]。而IC 的低值则受到阻抗因素（如C因子）变化的影响，延河流域的平均IC 下降，这与流域内植被的增加密切相关，IC 的表征反映了植被对于泥沙路径的阻碍作用，如延河流域较平坦的中游地区，作为城市所在区域且较早开展植被恢复措施，植被覆盖度较高（图9b1～图9b3），其IC 较低；而在地势起伏、地形破碎的上游地区以及坡陡谷窄的下游地区，植被稀少导致IC 较大。IC 呈现“坡面小沟谷大”的分布特点，也与植被覆盖的时空分布相关联，实际上，流域早期的植被覆盖较低，侵蚀泥沙有较大概率直接从坡面和沟坡进入河道后进行输移，其IC 一般较高，然而随着植被增加和相关水土保持措施的实施，大部分侵蚀泥沙被截留，泥沙路径的连通程度受到阻碍，IC 逐渐减小（图12），但上述情况多发生于坡面，而沟谷自然成为IC 较高的区域。

图12 1985—2020流域平均IC及SDR与NDVI的年际变化关系Figure 12 The interannual variation relationship between average IC,SDR,and NDVI in the watershed from 1985 to 2020

通过式（5）获取了流域近40 年的泥沙输移比（图12），流域平均SDR有下降趋势，2015 年后的SDR普遍在0.5 以下，这说明流域的一系列植被恢复举措取得了成效，更多的泥沙被截留在流域内，即植被恢复对减少流域侵蚀产沙具有积极的影响[39]。IC 和SDR的减小也说明侵蚀泥沙的输移受到了更大的阻碍，泥沙路径的连通程度受到增加植被的阻碍或许是两者减小的重要原因。因此，植被恢复对于流域泥沙管理至关重要，研究植被恢复如何影响流域侵蚀和产沙能够帮助理解流域的泥沙演变规律，未来更应该着重于这种研究的定量识别。

3.3 适用性和局限性评价

总体上，构建的模型方法框架能够较好地应用于流域研究，首先在模型的构建上，RUSEL、SWAT 和IC模型都易于建立，且SWAT 模型通过了水文站实测数据验证。此外，上述模型已在延河流域得到一些运用[19，40-41]，RUSLE、IC和SDR结果则与已有研究保持一致[42-43]，三者的结合在该流域创新应用是可行的。

同时模型方法框架也存在一定局限：首先，联合IC 和径流侵蚀功率主要分析了坡面侵蚀泥沙的传递和输移过程，缺少对河道泥沙输移更详细的描述。其中，IC 虽在一定程度上反映了河道的泥沙路径，但这种泥沙路径通常被认为是连通的，因为河道中的泥沙输移主要受径流的动力驱动影响[44]，植被或地表的阻碍作用影响较小，因此构建的IC 在表达河道泥沙输移时缺乏对径流驱动的考虑。此外，研究所涉及的泥沙连通性只具备结构特征，其计算过程依赖于地形和景观指标，缺乏动态特征的表达，如径流和降雨等动态要素，为此本研究结合径流侵蚀功率表达了子流域径流对侵蚀过程的影响，补充了这一不足，但仍需重点关注泥沙连通性的动态特征。在后续流域泥沙研究中，还应关注更细致时间尺度或事件尺度的泥沙过程，因为研究结果发现径流侵蚀功率和侵蚀过程间没有较好的关联性，但在以往的事件尺度研究中发现径流影响了侵蚀的发展和衰减，并且径流的产与汇极大地影响了泥沙的输移和沉积[45]，大流域尺度或大时间尺度的研究可能难以表达这些过程。