退耕还林背景下北方农牧交错带土壤保持功能时空变化

刘孟竹, 张红娟, 任贺宇, 裴宏伟

(1.河北建筑工程学院 市政与环境工程系, 河北 张家口 075000; 2.河北省水质工程与水资源综合利用重点实验室， 河北 张家口 075000)

生态系统提供给人类生存所需的各种惠益统称为生态系统服务[1]。土壤保持服务作为生态系统中土壤保育、水源涵养等生态服务功能的基础，对于区域内的土地安全以及生态可持续发展意义重大。土壤侵蚀与沉积是生态系统维系自身稳定发展的一种自然过程[2]，但其演变过度时会造成水土流失、土地退化、地力下降、河道淤塞等，严重时会导致一系列自然灾害发生，威胁人类生产生活[3]。有关土壤保持的研究一直以来都是国内外学者广泛关注的重要环境议题[4]。欧阳志云等[5]学者对中国生态系统服务进行的一项评估中发现，包括土壤保持服务在内的6项生态服务在2000—2010年中得到明显提升，而退耕还林等生态保护政策对此起到重要作用；Ren等[6]分析了长江上游流域水土流失的分布，研究指出森林生态系统对土壤保持的作用最大；同时，Lufafa等[7]通过预测模型评估了维多利亚湖流域的土壤损失分布，发现农田土壤保持能力最高，其次为牧场。综合上述研究来看，InVEST(Integrated Valuation of Ecosystem Services and Trade-offs)模型在对土壤保持的量化和绘图中受到了广泛使用，在土地利用变化对土壤侵蚀和保持的影响研究中发挥了重要作用。

20世纪60年代，由Wischmeier提出的可定量估算土壤保持量的土壤侵蚀方程(USLE)，被认为是目前应用最广泛的土壤侵蚀模型[8]。该模型计算土壤侵蚀的原理是通过潜在土壤侵蚀量减去实际土壤侵蚀量得到土壤保持量，然而该过程并未考虑土壤自身拦截上游沉积物的能力，使得计算结果有失偏颇[9]。这一问题在InVEST模型中得到了有效解决。InVEST模型全称为生态系统服务功能权衡综合评估模型，包括了水源涵养、生物多样性、碳储量、生境质量等子模型。其中，泥沙输移模块，在USLE的基础上考虑了地块拦截上游沉积物这一重要水文过程，使得结果更符合实际，目前已被应用于中国的山区、流域以及非洲桑塔尼亚、美国夏威夷等地区[10]。

北方农牧交错带处于中国的半干旱半湿润过渡区，生态环境极为脆弱[11]，该区同时也是我国中东部地区的“生态屏障”[12]，具有极其重要的生态区位。长期以来，北方农牧交错带 “重农轻牧”、“重生产轻生态”，产业结构严重失衡、超载放牧以及水资源无序超采等已经造成严重的草场退化、土地沙化、水土流失等生态环境问题[13]。对于本地区而言，水土资源的合理、持续利用，是实现农牧业可持续发展、经济稳定增长的重要保障。退耕还林工程实施以来，区域内的植被条件得到一定程度的改善。北方农牧交错带处在农耕区向草原牧区过渡的特殊地带，东西两侧地理跨度大，气候也不尽相同。目前，对于北方农牧交错带土壤保持的相关研究主要集中在该区西段[14]以及较为分散的市域[15]，针对整个片区的相关研究目前鲜有报道。基于此，本文收集该区2000—2018年土地利用、降雨、土壤、DEM (Digital Elevation Model)等多源数据，通过InVEST模型定量估算该区的土壤保持总量，并分析过去近20 a的变化特点，以期为该区生态治理以及可持续发展提供参考。

1 材料与方法

1.1 数据源及处理

本研究所用的土地利用数据来源于中国科学院资源环境科学数据中心，该数据源于Landsat影像解译制成，分辨率为30 m。土地利用数据时间节点选取了2000年、2018年2期，该数据精度在85%以上，能满足本文研究需要。土地利用分类根据《土地利用现状分类》标准并参考北方农牧交错带地区的实际情况将其分为6类：耕地、林地、草地、水域、建设用地、未利用地。DEM数据来源于地理空间数据云(http:∥www.gscloud.cn/)SRTMDEMUTM 90 M数字高程产品，裁剪后的研究区DEM栅格图在ArcMap软件中经过了填洼处理。土壤数据来自于国家冰川冻土沙漠科学数据中心的1∶100万中国土壤数据集(http:∥www.crensed.ac.cn/portal/)。降雨数据来自于中国气象科学数据共享服务网(http:∥data.cma.cn/)，主要为46个遍布北方农牧交错带地区的气象站研究期间的逐日降雨数据。所有空间数据均统一成WGS_1984坐标系，模型中所需的其余参数数据均参考模型数据库。

1.2 研究区概况

本文参考中国农业农村部于2017年发布的《农业部关于北方农牧交错带农业结构调整的指导意见》，该指导意见中详细划分的我国北方农牧交错带重要区域共涉及宁夏、甘肃、内蒙古、陕西、山西、河北、辽宁在内7省的146个县市(图1)。图1中具体划定的北方农牧交错带总面积约为4.7×105km2，经DEM数据裁剪后该区海拔在27～3 061 m。该区年均气温2～8℃，近30 a来年均增温0.32 ℃/10 a，年均降雨量400 mm左右，降水年际变化较大[16]。该区地貌单元主要以高原、丘陵结合为主，部分区域地貌为平原、沙地、山地[9]。一直以来，北方农牧交错带存在着过度开垦、超载过牧等现象，土地荒漠化、生物多样性丧失、水土流失现象严重[17]。目前研究区已实施的几个主要生态修复工程包括有“三北防护林”、“京津风沙源治理”、“退耕还林”等。

图1 研究区高程

1.3 研究方法

1.3.1 流域边界确定 InVEST模型需要输入研究区域内的子流域矢量图，本研究通过填洼后的DEM数据在ArcMap软件中采用水文分析工具得到。汇水累积量阈值是提取流域的关键参数，当值过低时，会导致提取的子流域过于密集细碎。如果该值偏高会造成提取的各子流域面积过大，研究细节则被忽略。流域汇水阈值选取500，1 000，2 000，…，15 000等一组参数进行调参，当该阈值设置为5 000时，获取的河网以及流域整个的边界与研究区地理图基本吻合，子流域大小合适，最后确定该阈值下研究区生成共计生成3 385个子流域。

1.3.2 土壤保持量计算 InVEST土壤保持模型是在USLE的理论基础上加以改进。在沉积物持留方面，该模型考虑了植物对土壤侵蚀的减缓以及对上坡沉积物具有拦截作用这一重要水文过程，使得模型结果更科学合理。本研究采用改进的USLE方程，具体计算方法参考以往研究[18]如下：

RKLSx=Rx×Kx×LSx

(1)

USLEx=Rx×Kx×LSx×Cx×Px

(2)

(3)

SEDRETx=RKLSx-USLEx+SEDRx

(4)

式中:SEDRETx为栅格x的土壤保持量(t);RKLSx为栅格x基于地貌和气候计算的的潜在土壤侵蚀量(t);USLEx为栅格x的土壤实际侵蚀量;SEDRx为栅格x的土壤持留量;Rx,Kx,LSx,Cx,Px为单元栅格x的降雨侵蚀力因子、土壤可蚀性因子、坡长坡度因子、植被覆盖与管理因子、土壤保持措施因子;SEx为x栅格泥沙持留率;SEz为上游栅格z泥沙持留率。USLEy为进入栅格x的上游栅格y产生的泥沙量。

1.3.3 降水侵蚀力因子R降水侵蚀力(R)是土壤侵蚀模型中最重要的参数，其反映了降雨条件下雨水对土壤的分离、运移、冲洗而导致水土流失的潜在能力。由于北方农牧交错带所处在干旱半干旱区，年降雨量偏低，经过降雨资料的筛选发现该区大部分气象站点日降雨量≥12 mm的天数一年中极少，故不采用日降雨量降雨侵蚀力算法模型。本文选择月尺度公式对研究区气象站点数据进行计算，公式参考以往研究[19]如下：

(5)

式中:p表示2000—2018年气象站点的年降雨量(mm);pi表示月降雨量(mm)。其中R需乘转换系数17.02得到国际单位(MJ·mm)/(hm2·h·a)。将每一个站点逐月降雨量与年平均降雨量代入上式计算,统计结果通过克里金法插值得到研究区降雨侵蚀力图。

1.3.4 土壤可蚀性因子K土壤可蚀性用来表征土壤对降雨对其侵蚀力的敏感性。Wischmeier等[20]定义土壤可蚀性因子K为标准小区(坡长22.12 m，坡度9%)上单位降雨侵蚀力引起的土壤侵蚀率。本文选择改进的USLE方程[21]，得到北方农牧交错带土壤可蚀性因子K图。公式如下：

KUSLE={0.2+0.3exp[-0.0256Ssand×(1-Ssilt/100)]}×(Ssilt/Sclay+Ssilt)0.3×

[1-0.25Corg/(Corg+exp(3.72-

2.95Corg)]×{1-0.75Sn/[Sn+

exp(-5.51+22.9Sn)]}

(6)

Sn=1-Ssand/100

(7)

式中:Ssand为沙粒含量;Ssilt为淤泥、细沙含量;Sclay为黏土含量;Corg为各图层有机碳含量;Sn为代数符号,无实义;KUSLE为改进USLE方程构建的土壤侵蚀力因子。

该模型K值与中国地区实测K值存在出入，故采用张科利等[22]研究构建的K值修正公式，公式如下：

Kreal=-0.01383+0.51575KUSLE

(8)

式中:Kreal为修正后的土壤侵蚀力因子;KUSLE为改进USLE方程构建的土壤侵蚀力因子。

1.3.5 坡长坡度因子LS 地形是土壤侵蚀中较为重要的参数，其中坡长、坡度因子能表征坡度与地表的关系，是影响坡面汇流以及产沙的重要地貌因素。将研究区的DEM数据导入到InVEST模型中，模型会根据输入的坡度阈值用不同公式分坡度来自动进行计算得到LS因子图层。本文中参考多数研究选择，将坡度阈值设置47%(约为25°)，计算公式由模型运行[23]。

1.3.6 植被覆盖与管理因子C覆盖和管理因子C是指植被覆盖或实施田间管理的区域土壤侵蚀总量与实施清耕的连续休闲地土壤侵蚀总量的比值，在0～1取值[24]。多数研究采用归一化植被指数NDVI进行计算提取得到C值。考虑到常用的NDVI产品分辨率较低，本文参考以往研究对土地利用类型进行赋值，以生物物理系数表导入模型自动生成覆盖和管理因子C图层，其中耕地、林地、草地、水域、建设用地、未利用地分别赋值为0.05，0.03，0.04，0，0，1[25]。

1.3.7 水土保持措施因子P水土保持措施因子P是指采取水土保持措施后的土壤侵蚀量与不采取任何措施的顺坡耕作时相应土壤侵蚀量的比值，介于0～1。0值时代表没有土壤侵蚀区，1值表示未采取水土保持措施的区域。其中耕地、林地、草地、水域、建设用地、未利用地分别赋值为1，1，1，0，0，1。

2 结果与分析

2.1 各省土壤保持量变化分析

2000年、2018年北方农牧交错带土壤保持及其变化由InVEST模型和ArcMap软件量化与绘制(图2)。根据统计，2000年、2018年研究区土壤保持总量分别为214.76亿t，214.93亿t，研究期间共计增长0.17亿t。结合表1分析，土壤保持量增长区域主要分布在河北、辽宁和山西，因为在研究期间上述3省的平均土壤保持量呈增长趋势，近20 a来分别增加了1.825，0.299，0.453 t/hm2；相反地，宁夏和内蒙古土壤保持量呈减少趋势，2000—2018年分别由274.259，258.234 t/hm2减少到274.092，258.227 t/hm2；陕西和甘肃在2000年、2018年土壤保持量相差不大。从分布来看，北方农牧交错带东南部(河北和山西)土壤保持量较高，并且呈增长状态；相反地，西北部(陕西、宁夏、甘肃、内蒙古)以及辽宁在土壤保持量能力上表现不突出。综上来看，北方农牧交错带7省土壤保持量变化的特点是“三增两减两不变”，分布上呈现“东南强，西北弱”的趋势。

表1 2000-2018年北方农牧交错带各省土壤保持量变化

图2 2000-2018年北方农牧交错带土壤保持变化

2.2 不同坡度等级土壤保持量分析

北方农牧交错带范围较广，但为了更好区分不同坡度区域的土壤保持水平，本文采用了坡度每隔15°一段的等距划分法，并赋予等级(表2)。实际上，表2中数据与2000年数据并无差异，本文选择了邻近年份2018年的统计数据进行分析。北方农牧交错带低坡度(0°～15°)土地面积居多，总面积3 908.96万hm2，占比约83.58%；其次是较低坡度(15°～30°)，737.94万hm2，占比15.78%。中等及中等以上坡度面积占比不到1%。尽管低坡度(0°～15°)区域土壤保持总量最多，约116.22亿t(54.09%)，但其平均土壤保持量最低，仅297.31 t/hm2。较低坡度(15°～30°)区域面积约为研究区总面积的15.78%，但占据了研究区42.02%的土壤保持总量，约92.28亿t。中等及中等坡度以上区域面积占比0.64%，土壤保持总量比重达到3.89%。综上，随着坡度的上升，土壤保持总量减少，但单位面积的土壤保持量增多，即坡度越陡，土壤保持能力越强。多数研究表明，大于15°以上的土地更易发生土壤侵蚀，随着坡度越陡且越发严重[26]，但地表植被会对土壤潜在侵蚀产生一定的阻控，土壤保持能力随之增强[27]。总的来看，北方农牧交错带整个区域处于低坡度区，坡度在0°～15°的土地面积占比83.58%，同时该区土壤保持量在坡度分布上呈现梯级上升式的“总量减、均值增”的趋势。

表2 2018年北方农牧交错带不同坡度土壤保持量

2.3 不同土地利用类型土壤保持量分析

北方农牧交错带不同土地利用类型的单位面积土壤保持量在2000年、2018年及其期间的变化统计见表3。2018年研究区林、草、耕地土壤保持总量分别为76.61，95.74，38.92亿t，分别占据土壤保持总量的36.14%，43.99%，17.50%。研究期间，林地相较于其他土地利用类型，平均土壤保持量在整个研究时期均保持最高，在2000年、2018年分别为923.28，863.62 t/hm2，下降了59.66 t/hm2；其次是草地，由2000年的208.87 t/hm2上升到2018年的520.39 t/hm2，近11.52 t/hm2的增长。耕地的平均土壤保持量变化不显著，在研究初、末期分别为237.96，238.44 t/hm2，未利用地的平均土壤保持量最低，2018年仅42.56 t/hm2，不足林地的1/20。综上来看，对研究区土壤保持功能提高贡献最大的是林地和草地

表3 2000-2018北方农牧交错带不同土地利用类型土壤保持量

2.4 “退耕还林”区域对土壤保持的影响

为定量评估退耕还林(草)工程对北方农牧交错带地区的2000—2018年的土壤保持功能影响，在InVEST模型中，分别以该区2000年、2018年的土地利用数据作为变量，其余参数均以2000年数据为基准代入到模型中分别进行了2组试验。同时，提取了研究区2000—2018年期间耕地转林地、耕地转草地、未利用地转林地、未利用转草地的的4个图层。通过ArcMap软件将2000—2018年土壤保持变化栅格图与耕地转林、草地图层叠加分析，得到结果为研究区退耕还林(草)区域土壤保持量变化栅格图(图3)。

图3 北方农牧交错带2000年-2018年“退耕还林”区土壤保持变化

北方农牧交错带2000—2018年各个省“退耕还林”区统计见表4。整个研究区由耕地、未利用地转变为林地和草地的总面积约为109.20万hm2，其中耕地转林地、耕地转草地、未利用地转林地、未利用地转为草地的面积分别为43.00(39.38%)，53.92(49.38%)，0.92(0.84%)，11.36(10.40%)万hm2。“退耕还林”实施的主体仍然是耕地，未利用地植被重建只占到11%。在各个省中，辽宁、陕西、内蒙古“退耕还林”面积最多，分别为24.05(22.02%)，22.81(20.88%)，22.59(20.69%)万hm2；河北和山西次之，分别为13.68(12.53%)，14.30(13.09%)万hm2；宁夏和甘肃最少，分别为4.94万(4.50%)hm2，6.83(6.25%)万hm2。陕西省和甘肃省是最早的第一批退耕还林试点区，而研究区内陕西省的林、草地恢复显然是显著于甘肃省的；内蒙古相较于其余省虽然在“退耕还林”面积中处于前列，但其区域面积最大，“退耕还林”面积相对并不显著。

表4 2000-2018年北方农牧交错带各省“退耕还林”区域面积 hm2

“退耕还林”区域的土壤保持量变化见表5。整个研究区耕地、未利用地转化成林、草地带来的土壤保持增量为1 014.50万t，占研究区土壤保持量总增量(约1 700万t)的59.68%。每1 hm2的耕地转为林、草地增长的土壤保持增量分别为1.38，0.80 t，但每1 hm2的未利用地转变为林、草地后土壤保持量分别增加67.72，75.32 t。这是因为耕地中的植被年内季节变化性较大，植被覆盖水平相对林、草地植被较低，当耕地的植被覆盖较高时，土壤中植被根系具有一定的固定土壤的能力；而未利用地以盐碱地、裸地为主，植被覆盖度极低，因此其土壤保持功能不如耕地，这导致了耕地和未利用地转成林、草地时的平均土壤保持量上存在较大差距(表5)。综上可以认为，北方农牧交错带近20 a来的土壤保持能力提升的主要原因是由于“退耕还林”的实施。其中，未利用地到林、草地的土地转入对土壤保持功能的提升贡献最大，转变区域内土壤保持量增加了近912.56万t，占研究区土壤保持总增量的89.95%；耕地转成林、草地的总面积占“退耕还林”区域面积的88.76%，转变区域内土壤保持增量仅占总增量的10.05%，约101.93万t。

表5 2000-2018年“退耕还林”区土壤保持变化

3 讨 论

3.1 误差分析

从本文结果的误差来源分析，土地利用数据的精度直接影响着对研究区土壤保持量模化的真实程度。不仅如此，InVEST模型对土壤保持过程的简化、输入数据的敏感性都会对结果产生影响[28-29]；从数据处理来看，土地利用数据经过了一定修正，如河北省张家口市张北县安固里淖，实际上2004年已经全部干涸[30]，在2018年土地利用数据中该区域被修正成未利用地；在气象数据处理上，本文采取了多年平均值，一方面是为了排除极端值给结果带来的不确定性，另一方面则是为了消除气象因素的干扰，从而保证土壤保持只受土地利用变化的单因素影响，这样的处理方法在同类生态系统服务的研究中[31-32]中也被使用。

3.2 土地利用变化对土壤保持影响

从本文结果来看，不同坡度的区域土壤保持能力变化(随坡度增大而增大)与参考文献[4]是相同的；对于研究区不同土地利用类型的土壤保持能力(平均土壤保持量)，林地>草地>耕地>未利用地的结果与已有研究[18]也是一致的，因此本研究结果具有可信度。从“退耕还林”对土壤保持能力的影响来分析，无论是土壤保持量总量还是均值，林、草地的相应水平都是最高的。林、草地作为植被覆盖度较高的土地利用类型，除了其地表的堆积的残物落叶可以抵挡降雨的直接冲刷，地底较长的根茎具有固定土壤以及良好的保水性，因此林、草地不易受侵蚀，土壤保持能力也相应较强；未利用地的土壤保持能力最低，这是因为盐碱地、沙地等土地土壤松散，不仅保水性差而且地表径流系数大，强降雨很容易造成沙化的土壤被雨水冲走而引起严重的水土流失。

4 结 论

(1) 在气候均质的条件下，中国北方农牧交错带2000年、2018年土壤保持量数量上分别为214.76亿 t，214.93亿 t，平均土壤保持量分别为459.033，459.395 t/hm2；在空间上呈现“东南强，西北弱”的分布。该区土地利用变化轻微增强了土壤保持能力，相比研究初期实现了0.356 t/hm2的增量。

(2) 2000—2018年，在省域尺度上，山西省土壤保持功能较其余6省最高(915 t/hm2)，河北省土壤保持能力提升最大(+1.83 t/hm2)；在不同土地利用类型上，各土地类型土壤保持功能依次为:林地>草地>耕地>未利用地，另外，这4种土地利用类型在2018年每1 hm2的土壤保持量分别为863.62，520.39，238.44，42.56 t；在坡度分布上，研究区土壤保持功能随着坡度梯级上升而显著增强，坡度最高区间内的单位土壤保持量是最低区间的26倍。

(3) 2000—2018年，“退耕还林”工程显著地提高了北方农牧交错带土壤保持功能，“退耕还林”区域内土壤保持量增加了1 014.50 t，占研究区总增量的近60%。本研究发现，在未利用地区域实施植被重建对土壤保持功能有显著的提升。但实际上，植被的恢复与重建往往会耗费大量人力、物力，且在盐碱地、荒地实施难度更大，这在改善区域土壤流失的生态议题中值得引起重视。