基于InVEST模型的矿区土壤保持功能评估——以银顶格矿区为例

武翼飞,陈国光

(1.中国地质科学院，北京 100037；2.中国地质调查局南京地质调查中心,江苏 南京 210016;3.中国地质大学(北京),北京 100083)

土壤保持功能是生态系统服务功能的一个重要方面,为土壤形成、植被固着、水源涵养等提供了重要基础[1]。生态系统服务和权衡综合评估模型(InVEST模型),是目前应用最广泛的生态系统服务功能评估模型[2]。InVEST模型中的土壤保持模块(SDR)在通用土壤流失方程(RUSLE)的基础上进行改进,考虑到地块自身拦截上游沉积物的能力,弥补了RUSLE方程中忽略这一水文过程的缺陷,使得土壤保持功能评估结果更加直观、合理[3]。近些年来,InVEST模型在我国不同地区的生态系统服务功能评估和土壤保持研究领域得到了应用。周彬等[4]应用InVEST模型对北京山区土壤侵蚀进行模拟,探究不同森林类型土壤保持功能差异;刘英等[5]基于InVEST模型构建综合生态系统服务评估模型,对神东矿区生态系统服务功能进行分析。随着InVEST模型的广泛应用,唐尧等[6]、肖微[7]等发现模型需根据地区具体情况不断调整参数等不足之处。

福建省大田县银顶格矿区曾经是福建省主要铁矿产区,矿区开采历史久,尤其是近30年运用大型机械深度开采,导致矿山开采区域山体植被破坏严重、生物多样性降低;矿山边坡土层结构松散、透水强、土壤肥力差,且土壤和水体污染严重,水源涵养能力不足[8]。近些年，随着福建省闽江流域山水林田湖草生态保护项目的推进,银顶格废弃矿山成为重点生态修复对象。本文以银顶格矿区所处小流域为研究对象,以土壤保持功能作为评估指标,基于InVEST模型应用银顶格矿区的降雨侵蚀力因子、土壤可蚀性因子、植被覆盖度、DEM和生物物理系数(如土壤保持措施因子)等数据,对比分析生态修复前后土壤保持功能的变化,探讨土壤保持功能的影响因素,为矿区生态修复效果评估和矿区生态系统服务功能评估提供借鉴。

1 研究区概况

研究区位于福建省大田县前坪乡、太华镇、均溪镇交界处(117°48′38″～117°50′34″E、25°50′0″～25°51′28″N)的银顶格矿区及所处的上华村流域,总面积为31.12 km2(图1)。区内地质构造复杂,岩石以钙质砂岩、石英砂岩为主,其次为花岗岩、凝灰岩和变粒岩[9]。大田县银顶格矿区生态修复项目于2019年10月启动,2021年4月完成修复。生态修复工程主要采取工程措施与生物措施相结合的方式。工程措施以水土保持、土地恢复为主,采取削坡减载措施,将采矿区挖成多级平台,共挖土方约729 062 m3;对原来剥离挖损的土地采取回填推平措施,共回填推平土方约686 155 m3;各平台设置钢筋混凝土截排洪沟合计长约39 552 m; 对强酸性土壤采取“播撒石灰石+覆盖专用防水毯”治理措施。生物措施以改良土壤、景观恢复为主,采取挖穴种植、播种草籽、栽植草坪等措施,恢复植被绿化面积约2.29 km2。在修复工程完成后,矿山景观得到改善,边坡稳定性提高,坡面流水通道畅通,存在的地质灾害隐患基本消除,酸性土壤得到有效治理,共恢复耕地面积0.21 km2,增加建设用地0.32 km2和农业设施用地0.94 km2。

图1 福建省大田县(a)及研究区(b)地理位置图Fig. 1 Geographical location map of Datian County (a) and the study area (b) in Fujian Province

2 研究方法

2.1 InVEST模型

本文利用InVEST模型中的土壤保持功能模块(SDR)评估银顶格矿区的土壤保持能力,通过对比修复工程前后的土壤保持状况,分析矿区土壤保持功能恢复效果。SDR模块又称为泥沙输移比模块,计算公式为

RKLSx=R·Kx·LSx，

(1)

USLEx=R·Kx·LSx·Cx·Px，

(2)

(3)

SDRx=SEDRx+RKLSx-USLEx，

(4)

式中：SDRx是栅格x的土壤保持量，SEDRx为泥沙持留量，RKLSx为潜在土壤侵蚀量，USLEx和USLEy分别为栅格x和其上坡栅格y的实际土壤侵蚀量，指植被覆盖下土壤流失量，以上单位均为t；SE为泥沙持留效率；R是降雨侵蚀力因子，MJ·mm·hm-2·h-1·a-1；K是土壤可蚀性因子，t·hm2·h·hm-2·MJ-1·mm-1；LS代表坡长与坡度因子；C为植被覆盖管理因子；P为土壤保持措施因子。LS因子、C因子和P因子均是无量纲参数。

2.2 数据来源

为评估银顶格矿区的土壤保持功能,本文采用遥感调查和资料收集相结合的方法开展工作。土地利用数据和植被覆盖数据是利用研究区2019年11月7日和2021年12月7日两期空间分辨率为2 m的高分六号遥感影像数据,经过预处理、土地利用类型解译[10-11]、实地验证后获得。InVEST模型的坡度与坡长因子分别选取两期DEM数据:第一期选用从NASA网站获取的2019年12.5 m空间分辨率DEM数据;第二期是2021年11月在研究区利用大疆无人机获得的矿山重点修复地区数字高程模型。土壤类型数据来源于福建省主要土壤类型专题图,由福建省1∶250 000土壤类型分布图[12]裁剪而成;土壤粒径、有机碳数据来源于中国土种数据库[13]。气象数据来源于中国气象数据网(http://data.cma.cn),本文选取福建省19个气象站点在2019年6月—2021年5月逐月降雨量数据,通过ArcGIS克里金插值获得研究区像元大小为12.5 m×12.5 m的月均降雨量栅格数据,用于计算降雨侵蚀力因子。

2.3 各因子计算方法

2.3.1 降雨侵蚀力因子

降雨侵蚀力因子(R)是降雨引起土壤侵蚀的潜在能力,是导致土壤侵蚀的关键性因素之一[14],主要受区域降雨量和降雨强度的影响。周伏健等[15]提出一套适用于福建省的基于月降雨量的降雨侵蚀力估算模型,该模型在降雨较充沛地区得到了广泛应用。本文采用此模型,利用站点月降雨量数据,估算银顶格矿区的降雨侵蚀力因子，其公式为

(5)

式中:R为降雨侵蚀力因子,MJ·mm·km-2·h-1·a-1;Pi为第i个月降雨量,mm。

经计算，研究区内降雨侵蚀力因子最大值为343 822 MJ·mm·km-2·h-1·a-1,最小值为331 897 MJ·mm·km-2·h-1·a-1,整体相差不大。利用ArcGIS软件克里金空间插值生成2019年和2021年研究区降雨侵蚀力因子分布图(图2),与陈晓瑜等[16]关于福建省降雨侵蚀力的研究结论相符。

图2 研究区2019年(a)和2021年(b)降雨侵蚀力因子分布图Fig. 2 Rainfall erosivity factor map of the study area in 2019 (a) and 2021 (b)

2.3.2 土壤可蚀性因子

土壤可蚀性因子(Kx)反映土壤对剥蚀与搬运的敏感性,不同土壤类型导致土壤侵蚀差异,其数值越高表示土壤越容易被侵蚀[17]。目前计算土壤可蚀性因子最为普遍是运用EPIC模型,方纲等[18]根据此模型计算出福建省主要土壤类型的土壤可蚀性因子,并应用到福建省水土流失研究中。EPIC模型中Kx的计算公式为

(6)

式中:Sa代表砂粒(粒径0.05～2 mm)含量,%;Si代表粉粒(粒径0.002 ～0.05 mm)含量,%;Sc代表黏粒(粒径<0.002 mm)含量,%;C代表有机碳含量,%。

研究区内的土壤类型主要包括黄壤、红壤、水稻土[19],其中红壤以黄红壤、硅钙质红壤、硅铝铁质红壤为主,水稻土以渗育水稻土和潴育水稻土为主(图3)。在矿区生态修复工程中主要采用了削坡平整的工程措施,这对研究区土壤类型的影响有限,因此，修复前后使用同一组土壤可蚀性因子。本文从中国土种数据库中选取邻近地区成土母质的土壤典型剖面的土壤粒径数据和有机碳数据计算土壤可蚀性因子,如黄红壤选择位于武夷山市星村的石英砂岩成土母质的典型剖面,红壤选择位于明溪县花岗岩成土母质典型剖面。经公式(6)计算后,得到研究区土壤类型典型剖面的土壤可蚀性因子(表1),在土壤类型图中赋值生成像元12.5 m×12.5 m的栅格数据,最终得到研究区土壤可蚀性因子分布图(图3)。

表1 研究区土壤类型可蚀性因子Table 1 Soil erodibility factor of main soil types in the study area

2.3.3 坡长与坡度因子

坡长与坡度因子(LSx)表示地形地貌因素对土壤侵蚀的影响程度,是影响土壤侵蚀的重要因素,InVEST模型一般通过DEM数据计算坡度与坡长因子[20]。

运用ArcGIS软件对修复前后的DEM数据(图4)进行填洼处理和地形反转处理,分别对正反地形DEM提取水流方向、计算汇流累积量、设置汇流累积量阈值(Fa)、生成河网、生成集水区。把集水区与反向地形集水区融合,得到汇水线与分水岭所组成的区域即为斜坡单元。应用坡度工具提取研究区坡度,研究区内坡度范围为0°～66°。依据《SL277—2002水土保持监测技术规程》[21]划分方案,将研究区坡度划分为6个等级:微坡(0°～5°)、较缓坡(5°～8°)、缓坡(8°～15°)、较陡坡(15°～25°)、陡坡(25°～35°)、急陡坡(>35°)。

图4 2019年(a)和2021年(b)重点修复地区高程图Fig. 4 DEM map of key restoration areas in 2019 (a) and 2021 (b)

InVSET模型可根据不同坡度自动选用不同的公式进行计算,由于模型是以美国坡度较小的试验区为原型,因此需调整边坡阈值,即大于此坡度的地区应当在采取水土保持措施条件下进行农业活动。根据王敏等[22]相关研究,坡度受土壤侵蚀影响的临界坡度为25°,因此本文将边坡阈值设为25°。

当坡度<边坡阈值时,坡度坡长因子的计算方法公式为

(7)

式中:LSx为坡度坡长因子;Fa为汇流累积量阈值;Cs为栅格大小,本文DEM数据栅格大小为12.5 m;s为百分比坡度,%;n为坡长指数,取值参考公式(8),

(8)

当坡度>边坡阈值时,坡度坡长因子的计算方法为

LS= 0.08β0.35s0.6，

(9)

(10)

式中:β为根据流向确定的栅格大小,m;s为百分比坡度,%。

2.3.4 植被覆盖管理因子

植被覆盖管理因子(C)是植被覆盖对土壤侵蚀的抑制能力,与植被覆盖度、土地利用类型相关。C值介于0～1,数值越大说明该区域表面的植被覆盖度越低,土壤侵蚀越严重[23]。应用EVNI软件提取出修复前后的归一化植被指数(NDVI)[24],再通过NDVI估算出研究区植被覆盖度。植被覆盖管理因子的计算采用目前适用较广的蔡崇法模型[25],公式为

(11)

式中:C为植被覆盖管理因子,取值范围为[0,1];FVC为植被覆盖度,%。

在ArcGIS软件运用上述公式计算得到研究区修复前后的植被覆盖管理因子栅格图。为符合InVEST模型数据输入格式,需要以土地利用类型为

单位确定植被覆盖管理因子的平均值。利用ArcGIS软件的空间分析工具统计出各个土地利用类型的植被覆盖管理因子平均值,其中将水面和交通运输用地类型分别修正为0和1,最后将因子赋值到研究区土地利用类型图中,得到研究区2019年和2021年植被覆盖管理因子分布图(图5)。

图5 研究区2019年(a)和2021年(b)植被管理覆盖因子分布图Fig. 5 Vegetation cover-management factor maps of the study area in 2019 (a) and 2021 (b)

2.3.5 土壤保持措施因子

土壤保持措施因子(Px)是反映采取水土保持措施后削弱土壤侵蚀能力的因子,取值范围为[0,1]。取值越小则代表该地区采取的水土保持措施对土壤侵蚀的抑制越有效,数值为0表示该地区无侵蚀;取值越高则代表水土保持措施越差,数值为1表示该地区未采取任何土壤保持措施。现今常采用土地利用类型赋值的方法[26],同一种土地利用类型采用的土壤保持措施大致相同。因此，本文生态修复前后同一种土地利用类型的P值也相同。借鉴周来等[27]、杨冉冉等[28]相关研究成果,确定研究区的土壤保持措施因子(表2):河流、坑塘、沟渠等水面赋值为0,该土地利用类型不会发生土壤侵蚀;交通用地、农村宅基地、工业用地等建设类用地,由于地面硬化几乎不发生土壤侵蚀,因此，赋予有效位数的最小值0.001;生态修复前的采矿用地P值赋值为1;其余植被覆盖的用地类型根据植物的抗侵蚀能力赋予0～1的值。

表2 不同土地利用类型的土壤保持措施因子表Table 2 The support practice factors of different land use types

2.3.6 模型参数及运行

植被覆盖管理因子和土壤保持措施因子需要制作成生物物理参数表格(csv格式),表中对相应的土地利用类型进行编码。矫正参数K、IC0和SDRmax值分别采用InVEST模型系统[22]默认值2、0.5、0.8。打开InVEST模型中的土壤保持模块(SDR),将所需数据导入运行模型,得到运算结果、栅格数据和矢量数据。

3 结果与分析

3.1 土壤保持功能分析

运用InVEST模型的土壤保持模块获得研究区修复前后的潜在土壤侵蚀量、实际土壤侵蚀量、泥沙持留量和土壤保持量等数据。

3.1.1 土壤侵蚀量空间分布及时空变化

通过统计得出,2019年银顶格矿区的实际土壤侵蚀总量为149 599.7 t,平均土壤侵蚀模数为4 806.4 t·km-2·a-1;2021年银顶格矿区的实际土壤侵蚀总量为81 373.3 t,平均土壤侵蚀模数为2 614.4 t·km-2·a-1。总量上来看,2021年的实际土壤侵蚀量比2019年减少了68 226 t,平均土壤侵蚀模数降低了2 192.1 t·km-2·a-1。根据中华人民共和国水利部颁布的《SL190—2007土壤侵蚀分类标准》[29]划分方案,土壤侵蚀强度划分为微度侵蚀、轻度侵蚀、中度侵蚀、强烈侵蚀、极强烈侵蚀、剧烈侵蚀6种类型(表3),依据此方案得到2019年和2021年银顶格矿区土壤侵蚀强度分级图(图6)。对比发现，2021年中度及以上侵蚀强度地区面积明显减少,强烈及以上侵蚀强度地区面积占比从2019年的21.0%降低至2021年的8.7%,表明生态修复工程对土壤侵蚀治理效果明显。

表3 不同土壤侵蚀强度统计结果Table 3 Statistical results of different soil erosion intensity

图6 研究区2019(a)和2021年(b)土壤侵蚀强度图Fig. 6 Soil erosion intensity map of the study area in 2019 (a) and 2021 (b)

经计算,修复后仍有部分地区土壤侵蚀较严重。实地调查后发现,目前仍有部分修复地区尚未被利用或裸露,以及部分强酸土壤治理区的植被恢复进度较慢,这些地区植被覆盖情况尚未改善,加之2021年降雨侵蚀力因子相对降低,都对修复结果产生一定影响。因此，生态修复后还应继续加强植物管护,推进矿区土地资源高效、合理利用。

3.1.2 土壤保持量空间分布及时空变化

应用ArcGIS软件得到2019年和2021年银顶格矿区土壤保持强度图(图7)。对运算结果进行统计分析表明,2019年银顶格矿区的土壤保持总量为17 546 658.6 t,土壤保持强度为563 748.2 t·km-2·a-1;2021年银顶格矿区的土壤保持总量为17 822 085.7 t,土壤保持强度为572 597.1 t·km-2·a-1。相较于2019年,2021年土壤保持总量和土壤保持强度都明显增加,土壤保持总量增加了275 427.1 t,土壤保持强度增加了8 848.9 t·km-2·a-1。尤其是研究区北部的矿山修复重点区在矿山修复工程的作用下,植被覆盖增加、土壤泥沙截留能力提升、土壤保持强度明显增强,表明生态修复工程对提升整体土壤保持能力效果较好。

图7 研究区2019年(a)和2021年(b)土壤保持强度图Fig. 7 Soil conservation intensity map of the study area in 2019 (a) and 2021 (b)

3.2 土壤保持功能的影响因素

3.2.1 土壤保持与土地利用类型的关系

由于不同土地利用类型的土壤保持能力存在差异,本文将2021年研究区的土壤保持量依据不同土地利用类型进行分区统计(表4)。统计结果显示:土壤保持总量从大到小依次是乔木林地、其他林地、水田、其他草地、采矿用地、旱地、竹林地、农村宅基地、交通运输用地、园地、设施用地;土壤保持强度从高到低依次是水田、其他林地、乔木林地、其他草地、采矿用地、竹林地、旱地、交通运输用地、农村宅基地、园地、设施用地。其中,乔木林地、其他林地的土壤保持总量最大,共占研究区土壤保持总量的77.1%,表明植被因子对土壤保持具有关键性作用。在今后的土壤保持防治工作中,应加大对林草地的保护力度,采取适当的水土保持措施来增加其土壤保持能力。由于水田采取了较好的土壤保持措施,且在地形较平缓地区,有效减弱了土壤侵蚀现象。研究区园地的植被覆盖相对较低,抗侵蚀能力较低,故土壤保持强度差。设施用地的土壤保持强度最差,主要因为研究区大部分的设施用地是由原来的采矿用地平整而来,相关设施还未完成建设,大片土地裸露,水土保持能力差。

表4 不同土地利用类型的土壤保持量统计结果Table 4 Statistical results of soil conservation in different land use types

3.2.2 土壤保持与坡度的关系

将研究区土壤保持量分布图与DEM图进行对比,发现土壤保持能力与坡度具有相关关系。将2021年坡度分级数据与土壤保持数据、植被覆盖度进行分区统计,得出不同坡度下的植被覆盖度、土壤保持强度和土壤保持总量(表5)。结果表明,坡度为0°～5°的地区土壤保持强度最弱,植被覆盖度也最低;从较缓坡至陡坡,植被覆盖度提高5%,土壤保持强度增加59 269 t·km-2·a-1。进一步统计发现,随着坡度的升高,植被覆盖度也相应提高,土壤保持强度也相应提高。经过分析和实地验证发现,矿区内此前受到严重破坏的高坡度地区已基本得到治理,现有的高坡度地区植被未受破坏,从而保持了较高的土壤保持强度。这也反映出,在使用土壤保持模块(SDR)时考虑到地块自身拦截上游沉积物能力的情况下,可以计算出更为真实的土壤保持强度和保持总量。

表5 不同坡度等级的土壤保持总量Table 5 Soil conservation for different slope gradients

不同坡度等级的土壤保持强度结果表明,在植被覆盖保持较高的情况下,泥沙的截留能力和植物根系对土壤的固结能力也会逐渐增强,从而抵消坡度对土壤保持的不利影响。因此，在银顶格地区水土保持工作中,要加强植被保护,提高土壤泥沙截留能力,以增强土壤保持功能。

4 结论

(1)银顶格矿区修复后,研究区土壤保持强度增加了8 848.9 t·km-2·a-1,土壤保持总量增加了275 427.1 t,表明研究区土壤保持强度明显提高，土壤保持总量有所提升。土壤侵蚀强烈及以上强度地区面积明显减少,占比从21.0%降低至8.7%,表明生态修复工程提升了研究区土壤保持能力。

(2)不同的土地利用类型中,水田、林地的土壤保持强度较强,有植被覆盖的土地类型的土壤保持能力远大于人类开发地区。从较缓坡至陡坡,植被覆盖度提高5%,土壤保持强度增加59 269 t·km-2·a-1,随着坡度的升高及植被覆盖度提高,泥沙截留能力和植物根系对土壤的固结能力也会逐渐增强,表明在地表植被不受破坏、保持较高植被覆盖的情况下,坡度对土壤保持强度的影响逐渐减弱。在土壤侵蚀防治工作中可以通过改善土地利用类型、加大对林草地的保护力度、提高植被覆盖度、加强陡坡生态植被保护等措施来提升土壤保持能力。

(3)InVEST模型的土壤保持模块(SDR)考虑到地块自身拦截上游沉积物的能力,弥补了RUSLE方程的缺陷,使得计算结果更加真实、准确。模型评估结果为栅格格式,实现了结果的可视化,可以更直观、便利地反映研究区土壤侵蚀和土壤保持情况。