基于InVEST模型的太行山区土壤侵蚀与土壤保持

朱建佳，彭晓伟，刘耀亮，刘金铜

(1 河北科技师范学院a园艺科技学院，b农学与生物科技学院，c工商管理学院，河北 秦皇岛，066600；2 中国科学院遗传与发育生物学研究所农业资源研究中心)

早在20世纪60年代，美国水土保持学家Wischmeier率先提出了通用土壤流失方程(The Universal Soil Loss Equation，USLE)，对区域土壤侵蚀和生态系统服务的研究具有重大意义[1]。该方程经过2次修正后得到修正的通用土壤流失方程RUSLE，与USLE的结构相同，但引入了土壤侵蚀过程的概念，改进了各因子的含义和算法，因而使用范围更加广泛[2]。然而，RUSLE只有在美国境内参数才不用调整，且忽略了上游地块自身的泥沙截留能力。由美国斯坦福大学、大自然保护协会和世界自然基金会共同开发的生态系统服务价值与权衡综合评价工具——InVEST模型(the Integrate Valuation of Ecosystem Services and Tradeoffs Tool)，旨在通过模拟不同土地利用下生态系统服务物质量和价值量的变化，为决策者权衡人类活动的影响和效益提供科学依据。InVEST模型的出现，弥补了USLE和RUSLE的上述缺陷。

近年来，InVEST模型已经广泛应用到我国不同地区土壤保持模拟的研究中。例如，李苒[3]研究表明，榆林市土壤保持效益与土地利用方式、流域内部地形地貌、生态系统类型、植被覆盖等因素有很大关系；王敏等[4]发现，福建宁德水土保持功能与社会经济发展密切相关；智长贵等[5]对河南淅川县土壤保持功能进行评估，得出有林地土壤保持量最大、未成林地土壤保持量最小的结论；陈童尧等[6]分析了不同土地利用、海拔高度、坡度以及土壤类型下祁连山国家级自然保护区的土壤保持状况，得出草地和栗钙土的土壤保持量最高。以上研究表明，土壤侵蚀与土壤保持的空间分布特征受多种人为和自然因素的影响，然而，这种空间异质性以及影响因素在山区表现尤为明显和复杂。因此，笔者拟在以上研究的基础上，探究山区土壤保持的垂直分布特征，试图明晰其与土地利用、地形等因素的相关关系。

太行山是我国华北地区主要山系之一，由于太行山的特殊地理位置，人口急骤增加，工农业的发展和战争的连年爆发，原始森林破坏殆尽，造成该区生态环境恶化，水土流失严重，自然灾害频繁。源自太行山的诸河流，雨季水势猛，泥沙含量高，酿成泥沙泛滥；旱季又干旱无雨，河流干涸，造成该区农业生产长期低而不稳，对华北平原和京津地区的经济建设和安全带来极大的威胁[7]。本次研究以太行山区为研究区域，利用InVEST模型探究太行山区的土壤保持空间分布特征及其影响因素，并分析空间上的土壤侵蚀强度，有助于决策者明晰土壤侵蚀的主要影响因素，从而更好的防止水土流失，有利于区域生态系统平衡。

1 研究区概况

太行山位于北纬34°35′～40°19′，东经110°15′～116°27′，包括山西、河北、北京、河南等省市共计101个县(区)(表1)，海河和黄河2个一级流域，4个二级流域和11个三级流域(表2)，总面积12.78万km2(图1)。太行山地处我国中纬度的温带和暖温带，气候温和，降雨量集中，光照充足，属于大陆性气候[8]。太行山是一个以山地为主，拥有黄土丘陵、山间盆地分布的复合山地地貌，为综合发展农林牧副生产提供了有利条件，同时也为林业发展提供了良好基地。太行山区土壤主要为褐土带，一般较低地带为山地褐土，海拔较高山地分布有淋溶褐土、山地棕壤及粗骨性土等，海拔2 000 m以上还分布有山地草甸土和亚高山草甸土。由于气候条件的分异，太行山植被从南到北也存在一定的差异，全区主要为暖温带落叶阔叶林带[7]。

表1 太行山区行政区划

图1 太行山区位置及地形

表2 太行山区一、二、三级流域

2 材料与方法

2.1 数据来源

本次研究使用的TM卫星图像来源于美国地质调查局网站(https://earthexplorer.usgs.gov/)，采用目视解译得到研究区土地利用图，数据精度为87.96%。30 m分辨率DEM来源于地理空间数据云(www.gscloud.cn/)，用于研究区地形分析和流域划分。2005～2014年日降水量来源于中国气象数据网(http://www.nmic.cn/)，用于计算降水侵蚀力因子R。土壤机械组成、有机碳含量来源于国家地球系统科学数据中心(http://www.geodata.cn/)，用于计算土壤可蚀性因子K。典型流域径流输沙量来源于《中华人民共和国水文年鉴》，用于土壤保持模型结果的验证。

2.2 InVEST模型原理

应用InVEST模型以及ArcGIS 10.5平台对太行山区的土壤侵蚀与土壤保持进行了模拟，土壤保持为栅格单元土壤潜在侵蚀减去实际侵蚀。土壤潜在侵蚀为裸地的土壤侵蚀量RKLS，其计算公式为：

RKLS=R·K·LS

(1)

土壤实际侵蚀为自然状态下的侵蚀量USLE，其计算公式为：

USLE=R·K·LS·P·C

(2)

SR=RKLS-USLE

(3)

式中：R为降雨侵蚀力因子[MJ·mm/(hm2·h)]；K为土壤可蚀性因子[t·h/(MJ·mm)]；LS为地形因子；C为植被覆盖因子；P为水土保持措施因子；SR为土壤保持量(t/hm2)；h为时间。

2.3 模型所需参数

2.3.1降雨侵蚀力因子R采用基于日降水量的计算方法[9]。

(4)

式中：Mi为第i个半月时段的侵蚀力值；k为该半月时段内的时间，Dj为半月时段内第j天的侵蚀性日雨量，要求日雨量大于等于12 mm，否则以0计算；α和β为模型待定参数。分别以每个气象站资料进行回归分析，由于不同站点得到的α和β不同，α和β的公式为：

α=21.586×β-7.189 1

(5)

(6)

式中：Pd12为日降雨量≥12 mm的日平均雨量；Py12为日降雨量≥12 mm的年平均雨量。

2.3.2土壤可蚀性因子K采用EPIC模型[10]进行K值的计算。

(7)

(8)

其中：SAN，SIL，CLA分别为沙粒、粉粒、粘粒的质量分数(%)；C为有机碳质量分数(%)。

2.3.3地形因子LS在地形要素中，坡度、坡长是形成有侵蚀能力径流的最主要因素，反映了地形地貌对土壤侵蚀的影响。LS指单位面积坡面流失与标准侵蚀小区(坡度为9%，坡长为22.12 m)的比值。在小流域和区域尺度上，往往通过目标区域DEM来提取该因子[11]。

对于坡度较缓的区域(<25%)：

(9)

式中：flowacc为栅格的汇流量，用栅格数量表示；cellsize为栅格大小，本次研究为1 km2；slope为坡度。

对于坡度较陡的地区(≥25%)：

LS=0.08λ0.35×prctslope0.6

(10)

式中：flowacc为栅格的汇流量，用栅格数量表示；cellsize为栅格大小，本次研究为1 km2；m为坡长指数；slope为坡度；prctslope为栅格百分坡度。

2.3.4C，P因子植被覆盖因子C是在某种条件下有作物种植地块上的土壤流失量与对应的无作物种植地块上流失量的比值，一般介于0与1之间；水土保持措施因子P指的是某种水土保持措施支持下的土壤流失量与对应的顺坡耕作条件下流失量的比值，一般介于0与1之间[1]。C，P因子反映的是土壤流失与地块种植的植被、作物以及地块管理措施等因素之间的关系。本次研究通过查阅文献[11～14]，得到太行山区C，P的取值。

2.4 模型验证

选取太行山区域内8个水文站对2005～2014年平均输沙量进行了验证。结果表明，8个子流域的模拟输沙量总和为19.61×104t，实际观测的输沙量总和为19.67×104t，仅相差0.06×104t，模拟效果较好(表3)。而且，模拟的输沙量数值与实际输沙量数值具有显著的线性关系(y=0.84x+0.38，R2=0.98，P<0.001)，故本次研究对于土壤侵蚀和土壤保持的模拟结果具有很好的参考价值。

表3 太行山区域内8个水文站控制子流域的集水面积、实际输沙量与模拟输沙量

3 研究结果

3.1 USLE各因子空间分布特征

太行山区USLE各因子均表现出明显的空间分异特征，其中降雨侵蚀力因子R受气候条件影响，自东南向西北逐渐降低(图2)；其余4个因子表现出随海拔发生变化的垂直地带性，并与土地利用类型密切相关。土壤可蚀性因子K，植被覆盖因子C随海拔升高而降低；地形因子LS，水土保持措施因子P随海拔升高而升高。USLE五大因子的空间分异规律共同决定了太行山区土壤侵蚀以及土壤保持的空间分布特征。

图2 太行山区USLE各因子空间分布

3.2 土壤侵蚀空间分布特征

2005～2014年，太行山区多年平均土壤侵蚀总量为0.37×108t，平均土壤侵蚀强度为4.3 t/hm2。参照《土壤侵蚀分类分级标准(SL 190—2007)》可将土壤侵蚀强度分为微度、轻度、中度、强度、极强度、剧烈共6个等级[15]。太行山区土壤侵蚀以微度和轻度为主，各等级土壤侵蚀面积从微度到剧烈在全区所占的比例分别为71.08%，27.01%，1.22%，0.36%，0.23%，0.10%。太行山区土壤侵蚀强度自南向北逐渐增高，在县域尺度上，河北省邢台、行唐，河南安阳、孟州土壤侵蚀强度较低，均低于0.12 t/hm2；山西省繁峙、平顺、五台、左权等县土壤侵蚀强度较高，均超过9 t/hm2。在一级流域中，海河流域土壤侵蚀强度较高，达到了3.84 t/hm2；黄河流域较低，为2.06 t/hm2。二级流域土壤侵蚀强度由高到低顺序依次为：海河南系(4.16 t/hm2)，海河北系(2.79 t/hm2)，龙门至三门峡(2.28 t/hm2)，三门峡至花园口(1.77 t/hm2)。三级流域中，位于太行山区东北部的子牙河、北三河、大清河土壤侵蚀强度较高，均高于4.42 t/hm2；太行山区东南部的沁丹河、小浪底至花园口土壤侵蚀强度较低，约为1 t/hm2(图3)。

图3 太行山区土壤侵蚀空间分布

3.3 土壤保持空间分布特征

太行山区2005～2014年10年平均土壤保持量为4.58×108t，单位面积土壤保持量为40.22 t/hm2。土壤保持与土壤侵蚀之间具有显著的正相关关系(R=0.51，P<0.01)，即对于太行山区来说，土壤侵蚀严重的地区，土壤保持功能也较好。

在县域尺度上，山西省黎城、左权土壤保持量较高，分别为104.73 t/hm2和101.78 t/hm2；河北省邢台、行唐，河南省安阳、孟州土壤保持量均低于2 t/hm2。在一级流域中，海河流域的土壤保持量较高，为46.04 t/hm2；黄河流域较低，为28.84 t/hm2。二级流域土壤保持量由高到低顺序依次为：海河南系(49.26 t/hm2)，三门峡至花园口(36.88 t/hm2)，海河北系(35.15 t/hm2)，龙门至三门峡(22.69 t/hm2)。三级流域中，北三河、三门峡至小浪底土壤保持量较高，分别为78.44 t/hm2和71.59 t/hm2；小浪底至花园口土壤保持量最低，仅为12.83 t/hm2(图4)。

图4 太行山区土壤保持空间分布

3.4 土壤保持的影响因素

3.4.1与土地利用类型之间的关系太行山区土壤保持总量随土地利用类型的变化从大到小依次为：林地，草地，水域，未利用土地，城乡、工矿、居民用地，耕地。而林地和草地提供的土壤保持量占总土壤保持量的54.1%。说明人为干预对土壤保持有很大的影响。

3.4.2与海拔之间的关系太行山区单位面积土壤保持量与海拔之间的关系为：海拔0～700 m时，土壤保持量随着海拔的升高而线性上升；700～800 m时，土壤保持量随海拔升高而下降；在800～2 200 m时，土壤保持量随海拔升高而缓慢上升；2 200～2 850 m时，土壤保持量随海拔升高波动性下降(图5)。

图5 太行山区土壤保持量随海拔的变化

3.4.3与坡度的关系按照《土地利用现状调查规程》规定，凡坡度小于或等于2°的称之为平地，坡地分为2°～6°，6°～15°，15°～25°，>25°等4个不同的坡级。太行山区土壤保持量随着坡度的增加而逐渐增加，呈明显的线性关系。4个坡级对应的单位面积土壤保持量依次为：20.74，32.32，49.54，101.12 t/hm2。

图6 太行山区土壤保持量随坡度的变化

4 讨 论

土壤侵蚀和土壤保持不仅受生态系统结构影响，更与人口、地形、土地利用类型等因素密切相关。本次研究利用InVEST模型对太行山区土壤侵蚀和土壤保持进行模拟，得出近10年太行山区土壤侵蚀强度为4.30 t/hm2，单位面积土壤保持量为40.22 t/hm2，土壤侵蚀强度以微度和轻度为主。其主要原因为，虽然太行山区地形陡峭，坡度较大，导致土壤潜在侵蚀较大。然而，植被覆盖类型主要以林地和草地为主，植被覆盖度广，生态系统结构较复杂，而且太行山地处我国中纬度的温带和暖温带，气候温和，降雨量集中，侵蚀性降雨较少。因此，土壤保持服务在该区作用显著，导致实际的土壤侵蚀强度较低。

太行山区不同土地利用类型土壤保持能力由高到低依次是：林地，草地，水域，未利用土地，城乡、工矿、居民用地，耕地。与前人得出的林地和草地对缩减侵蚀面积贡献率最明显的研究结果相符合[5,16]。太行山区由于森林和草地等植被覆盖率较高，提供的土壤保持量占总保持量的54.1%，因此太行山区具有较强的土壤保持能力。然而，土壤保持功能与区域自然条件及人类活动密切相关，不当的人类活动会导致水土流失的增加[17，18]，故城乡居民用地和耕地的土壤保持能力远低于林地和草地。以上结论表明，如果人类对于土地的扰动强烈，会增加土壤侵蚀，导致土壤保持能力下降，而对于一些侵蚀严重的地区退耕还林可以很好的改善土壤保持能力。

太行山区土壤保持能力与海拔和坡度因子之间关系密切。海拔800 m以下，土壤保持随海拔升高呈现先升高后降低的变化趋势；800 m以上，土壤保持随海拔升高呈现先缓慢升高再波动下降的变化趋势。分析原因为，海拔700 m以下主要为太行山东部与华北平原过渡地带，随着海拔升高，人口和耕地逐渐减少，由于不当的人类活动减少，故土壤保持量随海拔升高线性上升；在700～800 m主要为太行山西部的山西地区，是与黄土高原的过渡地带，人口和耕地多集中在这个海拔高度，由于人类活动影响增大，且黄土高原植被覆盖较少，因此土壤保持量有减少的趋势；随着海拔的升高，人口分布和耕地逐渐减少，土地利用类型以林地为主，所以土壤保持量有增加的趋势；2 200 m以上，植被逐渐由森林过渡到草甸，土壤保持量开始下降。

太行山区土壤保持量随坡度增加而显著增加，这与陆传豪[19]、郜红娟[20]等的研究结果一致。张彪等[21]研究表明，在植被覆盖、土壤属性等条件相近的情况下，土壤侵蚀量随坡度增加而显著增加[21]。本次研究结果表明，太行山区NDVI(归一化植被指数)随坡度的变化呈显著的二次函数关系，在坡度35°以下时，NDVI随坡度增加而增加，即植被覆盖率逐渐增大；坡度35°以上，NDVI随坡度变化平稳，植被覆盖率基本保持不变(图7)。由于植被覆盖率直接决定了土壤保持能力的高低，因此太行山区土壤保持量随坡度增大而增大，这也在一定程度上解释了太行山区土壤保持与土壤侵蚀呈正相关的原因。

图7 太行山区NDVI随坡度的变化

综合以上研究可以得出，太行山区土壤保持服务高值区主要集中在海拔和坡度较大高山区，土地利用类型以林地和草地为主，该区域同时也是土壤侵蚀相对严重区。坡面水与泥沙的空间重分布过程主要受控于植被的空间镶嵌格局，植丛斑块与裸地斑块和谐的镶嵌格局，能够最大程度地防止水土流失[22～25]。因此，在合理保护高山区土壤保持服务的基础上，提高中低山区土地利用的空间异质性，也有利于水土流失的缓解。

5 结 论

(1)2005～2014年，太行山区多年平均土壤侵蚀总量为0.37×108t，平均土壤侵蚀强度为4.3 t/hm2，全区侵蚀强度以微度和轻度为主；空间上土壤侵蚀强度自南向北逐渐增高。

(2)2005～2014年，太行山区土壤保持总量平均为4.58×108t，单位面积土壤保持量为40.22 t/hm2；土壤保持与土壤侵蚀之间显著正相关，土壤侵蚀严重的地区，土壤保持功能也较好。

(3)太行山区土壤侵蚀与土地利用类型、海拔、坡度等因子关系密切，林地和草地对土壤保持的贡献率最大，占总量的54.1%，城乡、工矿、居民用地和耕地贡献率最小；土壤保持服务高值区主要集中在海拔和坡度较大的高山区。保护高山区较高的土壤保持能力，同时提高中低山区土地利用类型的空间变异性，有利于进一步缓解区域水土流失和保证山区生态安全。