无定河流域土壤保持服务供需关系及服务流模拟

2022-05-26王玉丹周自翔王一达

水土保持学报 2022年3期

王玉丹，李晶，周自翔，王一达

(1.陕西师范大学地理科学与旅游学院，西安 710119；2.西安科技大学测绘科学与技术学院，西安 710054)

生态系统服务指生态系统形成和所维持的人类赖以生存和发展的环境条件与效应。目前，对于生态系统服务的研究多从自然角度出发进行静态研究，以生物物理区域或行政区划作为生态系统边界，忽视了系统的开放性和景观单元的空间异质性导致的生态系统服务供需错位现象。生态系统服务流能够在供给区与需求区之间建立空间联系，是生态系统服务管理的前提。

土壤保持服务为自然生态防止水土流失和保持泥沙的能力，能够防治土壤流失、保育土壤生产肥力，为调节服务中的重要内容，也是土壤健康、植被固碳、水质净化等功能的重要保障，研究土壤保持服务供需匹配及其服务流能够为防治水土流失提供新视角。生态系统服务潜在供给为生态系统长期持续提供的服务且不以人的意志为转移，土壤保持服务潜在供给即为生态系统服务的控制侵蚀和拦截泥沙的能力；生态系统服务潜在需求则包括人类社会已获取和人类希望获取但受条件限制而未能获取的生态系统服务，而土壤保持服务的相关研究则对潜在需求没有进行明确的定义。目前，对于土壤保持服务供需匹配关系开展了较多研究，但对于土壤保持服务流的模拟研究则较少。肖玉等将生态系统服务流分为原位、全向和定向服务流3种类型；李双成等将其分为供给移动服务流和使用移动服务流。对于土壤保持服务，上游地区生态系统防治土壤流失能为下游提供侵蚀调控服务，此为沿水系进行自然传输的定向服务流。Zheng等尝试通过Dinf算法来模拟土壤保持服务空间流动路径，并对供给区与需求区进行甄别；王勇智、梁丽丽利用SPANs(service path attribution networks)模型对土壤保持服务空间流动路径进行模拟，但在对土壤保持服务流流量模拟时是从输沙量角度来研究，以输沙量的减少程度来研究对下游地区提供土壤保持服务的增加程度。

无定河流域位于黄土高原与毛乌素沙地过渡带，暴雨集中在夏季，水土流失严重，生态环境具有波动性和脆弱性。鉴于此，本文将土壤保持服务作为研究对象，基于SWAT(soil and water assessment tool)分布式水文模型，对无定河流域土壤保持服务供需时空变化进行分析，并从供给移动服务流角度出发，对土壤保持服务中基于“土壤侵蚀—运移—沉积”过程的定向服务流—减沙服务流进行模拟与探究，揭示其空间流转机制。

1 材料与方法

1.1 研究区概况

无定河流域位于黄土高原中部，跨陕西省、内蒙古自治区两省，发源于陕西省白于山北麓，流域面积30 261 km，海拔为566～1 824 m，总体呈西北高、东南低的趋势。流域地处半湿润、半干旱地区交界带，为温带大陆性季风气候，雨热同期，夏秋季的集中降雨及黄土土质疏松易受降水侵蚀的特性成为无定河流域土壤侵蚀严重的重要影响因素。流域西北部为毛乌素沙地南缘，属于风沙区，分布有流动、固定沙丘，植被稀疏，风蚀剧烈；东南部为黄土丘陵沟壑区，人类活动频繁，以农耕地为主，垦殖指数高，为泥沙的主要来源区。受恶劣的自然条件和强烈的人类活动影响，水土流失严重，土壤侵蚀速率是黄土高原地区最高的，年输沙量占黄河下游输沙量的14.3%。严重的土壤侵蚀、脆弱的生态环境成为制约无定河流域的经济、社会可持续发展的关键问题之一。

1.2 数据来源

SWAT模型模拟的基础数据包括：无定河流域四、五级水系图(中国科学院资源与环境数据中心所，https://www.resdc.cn/)、DEM数字高程数据(JAXA全球ALOS—30 m地形数据，https://www.eorc.jaxa.jp/ALOS/en/aw3d30/data/index.htm)、土地利用(中国科学院资源与环境数据中心，30 m分辨率，https://www.resdc.cn/)、世界土壤数据(harmonized world soil database，HWSD，1∶100万)、气象数据(国家气象科学数据中心，http://data.cma.cn)和水文数据(黄河水利委员会，http://www.yrcc.gov.cn/)。

1.3 研究方法

1.3.1 SWAT模型构建 SWAT模型由美国农业部(USDA)开发，用于流域尺度水质与水量模拟，并对水质与水量受土地利用、土地管理等的影响情况进行定量分析和预测。本文根据无定河流域的地形、地貌、气候、水系等的特点设定子流域最小面积阈值为40 000 hm，把无定河流域划分为45个子流域，并利用土地利用类型、土壤数据、地形坡度等进行水文响应单元(HRU)分析，再输入逐日尺度气象数据库并由模型根据流域特征自动构建、模拟。

(1)SWAT模型算法。SWAT模型可以模拟流域内部多种地理过程。本文的主要研究内容为无定河流域土壤保持服务，通过SWAT模型得到无定河流域产沙量，在算法部分将SWAT模型的土壤侵蚀模拟算法作为重点。

SWAT模型使用通用土壤侵蚀方程MUSLE(modified universal soil loss equation)计算泥沙量，其计算公式为：

式中：SED为土壤侵蚀量(t)；为地表径流量(mm/h)；为洪峰流量(m/s)；为水文响应单元的面积(hm)；为土壤侵蚀因子(无量纲)；为植被覆盖和管理因子(无量纲)；为地形影响因子(无量纲)；为保持措施因子(无量纲)；CFRG为粗碎屑系数(无量纲)。

(2)基于SWAT模型的率定与验证。利用ArcSWAT对无定河流域作水文模拟后，使用SWAT—CUP(SWAT—calibration and uncertainty programs)对径流、泥沙进行率定与验证。根据前人研究及无定河流域的实际情况，最终选取27个调径流参数(表1中仅体现<0.5的18个)与7个调泥沙参数敏感性参数(表1)。

基于白家川水文站的月径流量与月泥沙量，利用SWAT—CUP并根据“先径流，后泥沙”的原则，对无定河流域2008—2013年径流量与2008—2018年泥沙量进行率定与验证，选用决定系数(correlation coefficient,)和系数(nash—sutcliffe efficiency)对结果进行适用性评价。当>0.6且>0.5时认为模拟结果可用。由图1可知，对于径流量，2008—2010年为率定期，=0.7，=0.64，2011—2013年为验证期，=0.78，=0.7；由图2可知，对于泥沙量，2008—2013年为率定期，=0.67，=0.51，2014—2018年为验证期，=0.74，=0.72。

表1 模型敏感性参数

图1 逐月径流量率定期与验证期的模拟值与观测值对比

1.3.2 土壤保持服务供需模拟

(1) 土壤保持服务供给模拟。本文中土壤保持服务的供给量为潜在土壤侵蚀量与实际土壤侵蚀量之差，计算公式为：

=SED-SED=11.8×(××)×××CFRG×(1-×)

式中：为土壤保持服务供给量(t)；SED为潜在土壤侵蚀量(t)；SED为实际土壤侵蚀量(t)。SWAT模型可以直接模拟得到实际土壤侵蚀量，但无法直接获取潜在土壤侵蚀量，因此将SWAT模型源代码中的因子和因子修改为1，模拟出假设不存在植被覆盖和保持措施下的侵蚀量，则为潜在土壤侵蚀量。

(2)土壤保持服务需求模拟。由于目前研究没有对土壤保持服务潜在需求进行明确定义，因此，将土壤保持服务潜在需求定义为不考虑技术、成本和可达性等限制因素、潜在土壤侵蚀量与通过土壤保持措施适宜分布区的布设获得的最小侵蚀模数之差，即将潜在土壤侵蚀量与最小可能土壤侵蚀模数之差作为土壤保持服务需求量。最小可能土壤侵蚀模数为某一地区容纳最大适宜水土保持措施量时的土壤侵蚀模数。

图2 逐月泥沙量率定期与验证期的模拟值与观测值对比

根据中国水利部发布的《土壤侵蚀分类标准》(SL190—2007)，黄土高原容许土壤流失量为1 000 t/(km·a)，但这一值于2007年针对西北黄土高原这一大面积区域提出，时效性较差，且对于无定河流域来说区域适用性较小。因此，拟采用Gao等布设不同土壤保持措施适宜分布区的方法，并考虑到无定河流域位于半湿润半干旱地区交界处，西北部为毛乌素沙地过渡带，生态环境脆弱，在水土保持措施布设方法基础上修改布设的优先级，将草地设为最优先的地位，建立了相应的土地利用类型，得到2000—2020年无定河流域最小侵蚀模数，求得无定河流域月平均最小侵蚀模数为45 t/(km·月)。

因此，土壤保持服务需求量具体计算公式为：

=-45

式中：为土壤保持服务需求量(t)；为潜在土壤侵蚀量(t)。

(3)土壤保持服务盈亏量。为进一步分析土壤保持服务供需关系，使用土壤保持服务盈亏量来做进一步探究，即土壤保持服务供给量与需求量之差。计算公式为：

式中：即土壤保持服务盈亏量(t)；为土壤保持服务供给量(t)；为土壤保持服务需求量(t)。当>0，则子流域土壤保持服务为盈余，当<0，则为亏损；当=0，则子流域的土壤保持服务为盈亏绝对平衡状态。

1.3.3 土壤保持服务空间流动模型减沙服务是基于“土壤侵蚀—运移—沉积”过程的定向服务流，因此以减沙服务量为无定河流域土壤保持服务流的研究对象。

(1)减沙服务流动模型。使用Tarboton提出的多流量模型Dinf算法来模拟土壤保持服务流动路径。在3像元×3像元窗口中，中心栅格单元与周围8个栅格单元形成三角形中的最大坡度为该栅格单元坡度，以该三角形的坡向为流向，且当前栅格单元流经的水流方向能够确定2个下游栅格单元。

(2)减沙服务流量及受益分析。减沙服务流量计算公式为：

号子流域减沙服务受益值为与其连通的上游各子流域减沙服务总流量与号子流域减沙服务量之和，计算公式为：

=∑+

式中：为第个子流域减沙服务受益值(t)；∑为与子流域相连通的所有上游子流域减沙服务流流量的累计值(t)；为第个子流域中减沙服务量(t)。

2 结果与分析

2.1 土壤保持服务供需时空变化分析

由图3可知，2000—2020年无定河流域土壤保持服务供给量与需求量变化趋势具有较强的一致性，且两者总体上都呈波动上升的趋势，而在2000—2020年降水量也大体呈上升趋势，黄土高原降水集中多为暴雨是影响水土流失的重要因素，降水量的增加导致土壤侵蚀强度增加，土壤侵蚀潜力也随之增加，表现出高供给高需求的情况。此外，整体上无定河流域土壤保持服务需求量大于供给量，且供给与需求均主要集中在每年的6—9月，与降雨量较多的月份相对吻合。其中，2001年8月、2013年7月、2017年7月供给量及需求量达到了2000—2020年所有月份中的较大值，供给量均超过2.5×10t，分别为2.92×10，2.85×10，3.14×10t，而这3个月的需求量则均超过了3.5×10t，分别为3.55×10，3.59×10，3.87×10t。相应地，月均降水量在这3个月也达到了较大值，分别为213，209，180 mm。无定河流域土壤保持服务对降水情况尤其是夏季强降雨情况反映迅速，具有较强的相互对应关系，进一步说明降水是影响无定河流域土壤保持服务的决定性因素。

图3 2000-2020年无定河流域月土壤保持服务供需量与降水量

空间上，为研究无定河流域土壤保持服务空间分布规律，选取2000年、2005年、2010年、2015年、2020年绘制空间分布图，由图4可知，整体上无定河流域的土壤保持服务供给量与需求量的空间分布格局相对具有一致性，具有明显的空间分异规律，均呈西北低、东南高的分布特点，且在2000—2020年均呈波动上升趋势。而在我国，降水量空间分布特点为东南向西北逐渐递减，与无定河流域内土壤保持服务供需量空间分布特点相符。在研究期间内，无定河中上游地区供给量与需求量相对较稳定，常年低于500 t/km。这部分区域以草地为主要土地利用类型，且降水少，大部分区域位于内蒙古自治区，地广人稀，侵蚀潜力小，潜在侵蚀量与实际侵蚀量均小，供给量与需求量较少。而中下游地区供给量与需求量上升趋势明显，这些地区以耕地为土地利用类型，降水量大，潜在侵蚀量与实际侵蚀量较多，供给量与需求量均高。到2020年，39，42，43号子流域供给量与需求量均超过了60 000 t/km，表明39，42，43号子流域在研究期间表现出了较大的不稳定性，2020年耕地面积占比分别高达97.35%，98.48%，89.36%，且子流域面积小，单位面积土壤保持服务供给量与需求量显著增高。

2.2 土壤保持服务盈亏时空变化分析

为了凸显无定河流域土壤保持服务的供需情况，本文使用盈亏量来表征土壤保持服务的供需关系。由图5可知，无定河流域土壤保持服务2000—2020年间盈余的月份较少，共有107个，占研究期的42%，且盈余量小，均低于9×10t；而占研究期58%的145个月份表现为亏损状态，亏损量大且呈波动性增加趋势，其中亏损量较大的月份主要集中在7—9月，这进一步体现出在无定河流域降水量尤其是强降雨对于土壤保持服务的盈亏量起着绝对性的作用。

由图4可知，空间上，无定河流域土壤保持服务盈余的区域主要为中上游地区，表现为盈余区域的同时盈余量较少，均低于500 t/km。这表明这些区域土壤保持服务的供需基本保持平衡，供略大于需。在2000—2020年间中上游地区实际土地利用类型以草地为主，尽管在2020年耕地、裸地、城市用地面积比2000年有所增加，但土地利用类型总体上变化不大，且坡度以5°以下为主，不易形成地表径流，侵蚀强度较低，因此多为盈余区域的同时盈余量少。无定河流域中下游地区则表现为亏损，2000—2020年亏损量大且整体为亏损加重趋势，较大面积亏损区域亏损量达到500 t/km以上。除了其位于半湿润地区降水量相对较多雨水冲刷量大外，这些地区坡度主要为15°以上，易受侵蚀，且城市用地扩张、植被占地面积降低导致土壤可侵蚀性增加，土壤实际侵蚀量随之增加。

2.3 土壤保持服务空间流动模拟

减沙服务流流量及受益值结果见图4。服务流的箭头指向为减沙服务在无定河流域各子流域内的流动路径，而符号的粗细则表示流的大小。无定河流域内减沙服务流的流动路径与河网水系具有较高的一致性，2000—2020年期间流量波动性增多，尤其2020年在下游地区增长明显。以研究期内流量变化较大的39号子流域为例(图6)，研究期内39号子流域减沙服务流流量呈波动上升趋势。表明自2000年黄土高原地区实施退耕还林还草措施以来，无定河流域保沙能力提升，产沙量减少，生态系统提供的减沙服务逐步增多。

图4 无定河流域土壤保持服务供给量、需求量、盈亏量、受益量时空分布

图5 2000-2020年无定河流域土壤保持服务月盈亏量

生态系统通过植被根系保育土壤，通过截留、吸收、下渗作用降低降水和径流对土地的冲刷力从而能够减少泥沙。当上游子流域泥沙量减少时能够减少径流中携带的泥沙量，从而为下游地区提供减沙服务，此传输过程具有累加效应。因此减沙服务受益量在传输过程中不受子流域的盈亏所影响，子流域的受益值包括其自身和与其联通的上游子流域减沙服务的总量。

鉴于此，选取3号子流域、21号子流域和45号子流域作对比。3号子流域为无定河流域的边缘子流域，其不与上游连通，受益值为自身减沙服务的值，研究期内年均受益值为6.75×10t；21号子流域与1，2，3，4，7，8号这6个子流域连通，其受益值为自身减沙服务的值与这6个子流域产生的减沙服务受益量之和，研究期内21号子流域年均受益值达到3.98×10t；而作为无定河流域出水口的45号子流域，在经过1～44号子流域减沙服务受益量的累加后，其最终的年受益值为8.44×10t。

注：以39号子流域为例。

3 讨论

土壤保持服务是生态系统调节服务之一，能够抑制侵蚀产沙和调控泥沙传输从而减少对当地和下游带来的负面生态效应。土壤保持服务流通过模拟生态系统阻滞泥沙传输路径，明确服务的供给区与受益区，量化供给与需求的空间失耦程度，建立起流域上中下游之间的空间联系，将上中下游视为一个整体，为流域水土流失综合治理和生态修复提供参考。土壤保持服务与人类福祉密切相关，未来应该关注土壤保持服务与人类福祉的关系，建立两者之间的互馈机制，进一步促进水土流失防治政策的科学制定。

本文在对土壤保持服务需求量化时，重构水土保持措施适宜分布区得到最小土壤侵蚀模数，避免了将实际侵蚀量直接作为土壤保持需求量所存在的问题。从供给移动服务流角度对减沙服务流进行模拟，一定程度解决了以输沙量减少程度作为服务供给量所导致的土壤保持服务与输沙产沙危害两者的概念冲突问题。路径模拟方法上选取多流量Dinf算法来模拟，根据李宗梅等的研究，Dinf算法提取水系的质量高于常用的单向流D8算法，对D8算法易产生的“平行”伪河道问题进行了改进。

使用SWAT水文模型模拟时存在不确定性的问题。SWAT模型能够基于物理过程较好地对水文过程进行长期连续模拟，特别是对月尺度水文过程模拟，但黄土高原地区河流泥沙通常由汛期的短历时高强度暴雨造成，对其水土流失具有决定性影响。无定河流域径流量与泥沙量对单日降雨量反映迅速，在SWAT模型中不能较好地体现。此外，已有研究发现，SWAT模型对夏季汛期径流量模拟较好，对于3—4月的春汛或降雨较少的枯水期模拟效果则较差，这与SWAT模型在国内其他干旱半干旱地区应用结果一致，径流量枯水期模拟结果或进一步影响泥沙模拟的效果。

对无定河流域土壤保持服务的研究是以SWAT模型中使用的MUSLE模型所得结果为基础的，其为水力侵蚀模型。无定河流域土壤侵蚀模式以水力侵蚀为主，约1/6的输沙量是由风力作用产生的，由风蚀造成的土壤侵蚀在本研究中没有考虑。