基于InVEST模型的低山丘陵区土壤侵蚀变化与驱动因素分析

韩 晶, 崔金芳, 杨 威, 徐阳吉 哲, 秦冬晖, 高凤杰

(东北农业大学 公共管理与法学院, 哈尔滨 150030)

土壤侵蚀是全球最严重的环境问题之一[1-2]，在造成土地退化与土壤养分流失的同时还会导致洪水频发、河流淤积和水污染等次生环境问题[3-5]。迁西县地处燕山南麓，低山丘陵景观占比36.1%，是典型北方土石山区，富有铁矿，且经济作物板栗驰名中外。20世纪90年代初期，奉行经济发展优先战略，采矿业无序开采导致林地被毁，山体滑坡、碎石崩落、尾矿库堆积等重力侵蚀灾害的危险性长期存在[6]；经济林板栗规模的盲目扩张侵占大量天然次生林，为生产方便，林下植被被人为清除，加之山区坡度大、土层薄的天然劣势，在降雨与地形双重作用下，水土流失问题十分严峻。大量泥沙随地表径流进入河流，造成滦河水系河道淤积与水质下降，生态环境质量不断下降。根据《中国水土保持区划》[7]，迁西县属于燕山山地丘陵水源涵养生态维护区，国家级水土流失重点预防区。2020年，国家级生态文明示范县建设对区域内水土流失等生态问题治理与修复要求进一步提升。因此，科学认识并定量评价区域人类不合理土地开发活动导致的土壤侵蚀变化对水土保持治理与生态恢复具有重要意义。

GIS与通用水土流失方程(Universal Soil Loss Equation)的结合是土壤侵蚀研究应用最广泛的方法之一[8-10]，其考虑了自然和人为扰动的影响，能够科学地反映区域土壤侵蚀的空间格局[11]，但该模型忽略了地块自身对泥沙的拦截能力使评估结果存在误差[12]。InVEST模型弥补了通用水土流失方程的不足，近年来成为生态服务功能研究的重要方法之一[13]。该模型综合考虑流域水库及地块泥沙持流能力的影响，其评价结果更为精准并具针对性[14]，因此基于InVEST模型研究土壤侵蚀在国内日益广泛[12,15-17]。以往在土壤侵蚀驱动因素的研究中，多采用相关分析或回归分析等方法，这些常规手段无法判断影响因子对土壤侵蚀的作用程度[18]且无法统计各因子间的交互作用，地理探测器弥补了传统分析方法的不足，可定量探究因子及因子间交互作用对结果的影响程度，更适用于土壤侵蚀这一复杂生态问题的研究[19-20]。本文采用InVEST模型的泥沙输移比例模块及地理探测器，定量评估1990—2020年迁西县土壤侵蚀时空变化规律及其驱动因子，研究结果为后续生态治理与区域可持续发展提供决策支撑。

1 研究区概况

迁西县位于唐山市北部燕山南麓，滦河流域中下游，地理坐标为118°10′—118°30′E，40°00′—40°20′N(图1)。地势由北部燕山脉系山地丘陵向滦河三角洲平原过渡，山地丘陵占比36.1%，是典型的低山丘陵地区，坡耕地比重高。滦河在迁西境内全长67.5 km，流域面积1158 km2，占全县面积的80%以上。境内有潘家口和大黑汀2座国家级大型水库，在上游山洪来水调蓄与下游平原农田灌溉发挥重要作用，同时是“引滦入津”工程基地。属温带大陆性半湿润的季风气候，四季分明，雨热同季，多年平均气温10.1℃，多年平均降雨804.2 mm。迁西号称“板栗之乡”，板栗种植已有2 000多年历史，创造了享誉世界的“围山转”工程，森林覆盖率达63%，位居全省第二位。富有铁、金、铜等矿产资源，其中铁矿储量4.7亿t，境内津西集团成为全国最大的H型钢生产基地。20世纪90年代起，随着市场经济的发展和技术水平的提升，板栗经济林的盲目扩张及铁矿的无序开采，对山地丘陵地表植被造成大面积破坏，在降雨侵蚀作用下，境内水土流失问题十分突出，由此衍生的次生灾害，如河库淤积及水体污染等问题也日益严峻。

图1 迁西县地理位置及地形

2 数据来源与研究方法

2.1 数据来源与处理

降雨数据来源于中国气象数据国家气象科学数据中心(http:∥data.cma.cn/)，本文选取研究区境内及周围31个站点的逐月降雨数据计算降雨侵蚀力因子；土壤数据来源于中国科学院南京土壤研究所(http:∥www.issas.ac.cn/)和土壤属性数据库(http:∥westdc.westgis.ac.cn)，用于计算土壤可蚀性因子；DEM(Digital Elevation Model)来源于NASA地球科学数据网站(https:∥nasadaacs.eos.nasa.gov/)的12.5 m数据，用于研究区小流域提取和坡度坡长因子计算；基于GEE平台对landsat系列数据提取1990年、2000年、2010年、2020年土地利用类型与植被覆盖度，在此基础上通过赋值获取水土保持措施因子。GDP与人口数据来源于资源环境科学与数据中心(https:∥www.resdc.cn/)的空间分布公里网格数据。将以上数据重采样至12.5 m栅格并统一到WGS_1984坐标系统。

2.2 研究方法

2.2.1 InVEST模型 InVEST模型3.9.0版本中的泥沙输移比例模块(Sediment Delivery Ratio)是基于通用水土流失方程(USLE)通过像元尺度来描述坡面土壤侵蚀和流域输沙空间过程。在模型运行前需通过GIS将所有数据处理为模型运行所需的格式，主要参数包括降雨侵蚀力R、土壤可蚀性K、地形因子DEM、土地利用/覆被、生物物理系数表、小流域矢量、汇水面积阈值以及模型运行所需的其他参数。通用水土流失方程为：

A=R×K×LS×C×P

(1)

式中:A为年土壤侵蚀量〔t/(hm2·a)〕;R为降雨侵蚀力因子〔MJ·mm/(hm2·h·a)〕;K为土壤可蚀性因子〔t·hm2·h/(MJ·mm·hm2)〕;LS为坡度坡长因子;C植被覆盖与管理因子;P为水土保持措施因子。

(1) 降雨侵蚀力因子R。降雨是泥沙输移的主要因素，本研究采用Wischmeier在1978年提出的简化算法[21](公式2)，结果乘以17.02转换成国际单位制。

(2)

式中:R为年降水侵蚀力〔MJ·mm/(hm2·h·a)〕；Pi为月降水量(mm);P为年降水量(mm)。

(2) 土壤可蚀性因子K。本研究采用Williams等提出的EPIC模型[22](公式3)，土壤颗粒和有机碳含量来源于土壤属性数据库[23]，结果乘0.131 7转为国际单位制。

(3)

式中:K为土壤可蚀性〔t·hm2·h/(MJ·mm·hm2)〕；SAN是砂粒的含量(%)；SIL是粉粒的含量(%)；CLA是黏粒的含量(%)；SN1=1-SAN/100；C为有机碳含量(%)。

(3) 地形因子LS。地形是导致泥沙输移的主要因素，InVEST模型根据DEM自动提取LS因子。通常情况下，坡度越大、坡长越长土壤侵蚀发生的可能性越大，计算原理如公式(4)。

(4)

式中:Si表示栅格单元坡度因子,当坡度θ<9%,Si=10.8·sinθ+0.03当坡度θ≥9%,Si=16.8·sinθ-0.50；Ai-in表示栅格径流入口以上产沙区域面积(m2),D表示栅格尺寸(m)；xi=|sinαi|+|cosαi|；αi表示栅格单元i的输沙方向；m表示USLE长度指数因子,当θ≤1%,m=0.2;当1%<θ≤3%,m=0.3;当3%<θ≤5%,m=0.4;当θ>5%,m=0.5。

(4) 植被覆盖与管理因子C。植被叶片通过削弱降雨对土壤的冲击来减少土壤流失，此外，植被可以有效抑制地表泥沙的输移，对水土流失起到很好的抑制作用，用植被覆盖度来度量(公式5)。本研究采用王万忠对中国土壤侵蚀因子的赋值方法[24](表1)计算植被覆盖与管理因子。

VFC=(NDVI-NDVIsoil)/(NDVIveg-NDVIsoil)

(5)

式中:VFC为植被覆盖度;NDVI为像元植被指数;NDVIsoil为裸土或无植被覆盖指数;NDVIveg为完全被植被覆盖区域指数。

表1 不同土地利用类型C，P因子值

(5) 水土保持措施因子P。该因子表征人为因素对土壤侵蚀的抑制作用，介于0—1之间。其中，不发生土壤侵蚀的区域P值为0，未采取任何保护措施的地区P值为1，其他用地类型P值赋值方法沿用前人的研究成果[11]，见表1。

(6) 小流域提取。InVEST模型主要通过小流域的泥沙输移过程与坡面的侵蚀状况评估区域的土壤流失与保持情况。通常情况下，汇水面积阈值越小，河网越稠密，划分的小流域越多。通过反复尝试，确定汇水面积阈值为5 000最为合适，最终形成172个小流域。

(7) 其他参数。kb(特定地块与径流的空间连接程度)和泥沙输移比IC0(进入河谷泥沙量与坡面侵蚀量之比)决定小流域水文过程空间联系与泥沙输移比关系形态的校准参数。IC0和k值为定义SDR(泥沙输移比)与IC关系的校准参数(递增函数)，采用模型默认值：kb=2，IC0=0.5。SDRmax栅格最大泥沙输移比由土壤质地决定，本文将其设为0.8。

2.2.2 地理探测器 地理探测器用以探索空间分异性，揭示其驱动力的一组统计学方法[25]。其以离散化的空间数据为基础探索采样点自变量与因变量之间空间分布的一致性，并以q值度量自变量对因变量的解释度，主要分为：因子探测器、风险探测器、生态探测器和交互作用探测器4大模块。因子探测器用于探测因变量的空间分异性，用q值大小描述自变量(影响因子)与因变量(土壤侵蚀量)空间一致性强弱，见公式(6—7)；风险探测器用于判断两个子区域间的属性均值是否有显著的差别，用t统计量来检验；交互作用探测器探测影响因子X1与X2共同作用时是否会增强对因变量Y(土壤侵蚀量)的作用程度，或这些因子共同作用时q值变小，即这些因子对Y的影响是相互独立的。

(6)

(7)

本文利用地理探测器判定年降雨量、土壤类型、坡度、土地利用类型、植被覆盖度、GDP和人口密度7个影响因子对土壤侵蚀的作用程度。模型运行前需对连续的因子进行离散化处理，通过将每层的连续型因子数据的属性值设置断点，利用断点划分出不同区间并依次编号，这样每个区间内所有的单元都具有相同数值。常用的离散化方法较多，经过反复试验，发现在本研究中采用等距离划分方式的探测结果中自变量与因变量的关系更符合实际，因此，将年降雨量、GDP和人口密度数据等距离划分成6份(表2)。坡度根据水利部关于水土保持坡度等级划分标准分为0°～5°，5°～8°，8°～15°，15°～25°，25°～35°和>35°共6个等级；土地利用类型沿用表1；土壤类型分为：棕壤、褐土、褐土性土、石质土、棕壤性土、土娄土、珊瑚砂土、钙质石质土、酸性粗骨土、钙质粗骨土、潮土共11类。最后，将以上数据在GIS中重采样到1 km×1 km网格点，剔除异常值后得到1 400个样本点。以影响因子为自变量，对应的土壤侵蚀量为因变量，运行地理探测器。

表2 土壤侵蚀模数及影响因素分级标准

3 结果与分析

3.1 土壤侵蚀变化

InVEST模型运行结果表明：1990年、2000年、2010年与2020年土壤潜在侵蚀量分别为：2.98亿t，3.08亿t，2.78亿t，2.80亿t，实际侵蚀量依次为：1.25×107t，1.41×107t，1.77×107t与2.00×107t。土壤保持量分别为：2.86亿t，2.94亿t，2.60亿t，2.60亿t。根据《土壤侵蚀分类分级标》(SL190-2007)，得到各期土壤侵蚀情况(表3)。1990—2020年研究区土壤侵蚀以微度和轻度侵蚀为主，二者面积占比之和均达到90%以上，4期分别为94.49%，93.47%，92.63%和91.68%；随着时间推移，二者侵蚀量占比呈下降趋势，4期分别为50.84%，47.17%，36.71%和33.10%。中度、强度、极强和剧烈侵蚀面积占比之和虽然不足10%，但逐期增加，其中剧烈侵蚀等级侵蚀量占比最高，并呈逐年增加趋势。总体而言，研究区尽管在空间上以微度和轻度侵蚀为主，但中度以上侵蚀等级土壤流失量却占据主要地位，特别是剧烈侵蚀，近20年呈愈演愈烈态势，水土流失生态治理需求迫切。

3.2 空间分布特征

为确定土壤侵蚀空间分布特征，将InVEST模型输出的172个小流域的平均土壤侵蚀模数分为0～50，50～100，100～200，200～300，300～400和大于400 t/(hm2·a)6个等级，依据《土壤侵蚀分类分级标》(SL190-2007)的分级标准，分别命名为微度、轻度、中度、强度、极强度和剧烈(图2)。结果表明：整体上土壤侵蚀北部比南部严峻，与研究区北部山地丘陵南部平原地形地貌直接相关。1990年与2000年不存在剧烈等级，极强度和强度侵蚀等级占比极低，2010年和2020年土壤侵蚀等级明显加重，剧烈、极强度和强度等级面积占比不断增大。

3.3 土壤侵蚀驱动因素分析

3.3.1 影响因子显著性分析 根据地理探测器输出结果(表4)，7个影响因素对土壤侵蚀解释力由强至弱依次为用地类型、坡度、土壤类型、植被覆盖度、年降雨量，GDP和人口密度作用微乎其微。用地类型的解释能力呈“倒U”形变化，2010年土地利用对土壤侵蚀的作用强度达到峰值，降雨、GDP和人口密度对土壤侵蚀量解释力均较弱，降雨在1990年与2000年未通过显著性检验，GDP与人口密度四期均未通过显著性检验，可能是三者在县域内无显著差异所导致的。

3.3.2 风险区域识别 风险探测器输出见表5，4期的土地利用类型风险高低顺序相似，依次为：工矿用地>未利用地>建设用地>经济林地>耕地>林地>草地>水域。工矿用地与未利用地由于无植被覆盖且没有水土保护措施，土壤侵蚀风险最高。由图3可知，研究区土壤侵蚀风险随坡度变化呈先增后减再回弹的趋势，两个拐点所在坡度范围分别为15°～25°和25°～35°，最大值普遍出现在35°以上地区，坡度为25°～35°时土壤侵蚀呈波谷低值，与该坡度范围土地利用类型大多为林地和草地有关，林地和草地可以有效抑制水土流失。由图4可以看出，土壤侵蚀风险随着植被覆盖度的增加先增大后减小，波峰位置各期不同，但拐点过后土壤侵蚀风险与植被覆盖度大体上呈负相关关系。土壤类型的高风险区域为潮土、棕壤和钙质石质土。降雨、GDP和人口密度不具有显著性差异。

表3 1900-2020年土壤侵蚀分级结果

图2 1990-2020年小流域土壤侵蚀等级

表4 1990-2020年土壤侵蚀驱动因素q值统计结果

表5 1990-2020年影响因子高风险区域

图3 1990-2020年不同坡度风险变化

图4 1990-2020年不同植被覆盖度风险变化

3.3.3 驱动因素交互作用分析 从表6可以得出各因子的交互作用均为双向促进关系，即因子间的交互作用均大于独立因子对土壤侵蚀的作用。用地类型与其他因子协同作用对土壤侵蚀解释力最强，这一结果与前文3.3.1中地利用类型对土壤侵蚀影响最大相吻合。其中，坡度与用地类型交互作用除2020年，其余年份均最大，说明坡度与土地利用类型交互作用对土壤侵蚀影响最大。降雨单独作用q值较小，但与其他因素协同作用q值大幅提高，特别是与用地类型及坡度，因此，加强陡坡区域退耕还林、植树造林等水土流失治理措施十分必要。GDP与人口密度两经济社会因素单独作用对土壤侵蚀作用微乎其微，但二者与其他自然因素共同作用时对土壤侵蚀的解释力大幅提高，说明人为干扰下自然因素对土壤侵蚀的决定作用更为强烈，因此，土壤侵蚀治理要兼顾自然因素与经济社会因素。

表6 1990-2020年影响因子交互作用q值

4 讨论与结论

4.1 讨 论

在以往的研究中USLE与地理探测器结合探索土壤侵蚀驱动因素的较多，本研究充分考虑InVEST模型的优势，确保低山丘陵区土壤侵蚀结果的准确性，尝试将地理探测器与该模型结合，分析土壤侵蚀变化规律与驱动因素。InVEST模型输出的子流域土壤侵蚀情况打破了行政边界的束缚，使土壤侵蚀情况在空间上反映的更为直观，胡胜等[26]也提出基于水文意义的子流域边界反映了沉积物质沿水文路径迁移的完整性，相较于行政边界其评估结果更为科学、准确。本研究发现土地利用类型对土壤侵蚀作用最强烈，其次为坡度，潘美慧等[11]曾提出土地利用类型和流域地形、坡度的变化是造成土壤侵蚀空间差异的主要原因，与本研究结果一致。迁西县土壤侵蚀发生的高风险区域为工矿用地与未利用地，与该两种用地类型地类植被覆盖较低且无水土保持措施有关，此结论与邹雅婧等[19]关于渭北矿区土壤侵蚀评估及驱动因素分析的研究结果一致。何莎莎等[16]的研究中土壤侵蚀主要发生在林地、草地和耕地，结果与本研究结论背道而驰，原因在于其研究区海拔为31～1 864 m，其中林地和草地多分布在高海拔地区，而本研究中海拔最高区域仅为829 m，进一步说明地形与土地利用类型共同作用结果对土壤侵蚀的作用变强。有学者认为土地利用类型是土壤侵蚀产生的主要因素[27-29]，因此，应严格控制用地类型调整方向。

陈思旭等[30]通过分析不同坡度的土壤侵蚀结构得出随坡度变化土壤侵蚀先增大后减小，15°～25°土壤侵蚀最剧烈，而何莎莎等[16]提出土壤侵蚀模数随坡度升高逐渐增大，本文基于地理探测器风险识别模块得出，土壤侵蚀发生的风险性随坡度增大呈“N”型变化，15°～25°与>35°是土壤侵蚀发生的高风险区域，且四期土壤侵蚀模数随坡度变化趋势一致，与该坡度上林地和草地分布相关。植被覆盖度范围为0～50%是发生土壤侵蚀的高风险区域，土壤侵蚀量与植被覆盖度大体上呈负相关关系，研究结果与李娜等[31]的探究结论一致。由于县级及以下降雨逐年数据获取较为困难，降雨因子插值结果在空间上无显著差异，从而导致地理探测器结果显示降雨与土壤侵蚀的相关性较小，已有研究[32-33]表明降雨强度决定着土壤侵蚀的剧烈程度，是土壤侵蚀治理的不可忽视的重要因素。本研究交互探测结果显示因子间交互作用对土壤侵蚀的解释力均大于单一因子的解释力，降雨单一因子与土壤侵蚀相关性不强，但与其他因子共同作用后相关性大幅提高，李桂芳[34]研究表明坡面土壤侵蚀过程是降雨强度、坡度与坡长共同作用结果，且三者存在互相促进的作用。此外，土地利用类型与坡度的交互作用结果对土壤侵蚀的作用最强烈，王欢等[18]针对喀斯特不同地貌形态对土壤侵蚀进行定量归因，结果显示土地利用类型与坡度的协同作用对土壤侵蚀的解释力最强。未来应增强不同坡度上土地利用类型与耕作方式土壤侵蚀的研究。

4.2 结 论

本文基于InVEST模型对迁西县1990年、2000年、2010年与2020年土壤侵蚀模数进行估算，结果显示：迁西县土壤侵蚀总量逐年增加，其中微度与轻度侵蚀量与侵蚀面积占比之和逐年减少，中度侵蚀以上等级的侵蚀量与侵蚀面积占比之和逐年增加，充分说明了区域土壤侵蚀逐年加剧。土壤侵蚀在空间上呈北高南低格局，与区域北部地势高坡耕地比重大直接相关。地理探测器对区域土壤侵蚀变化的驱动因素探索表明：(1) 土壤侵蚀影响因子的解释力由强至弱依次为：土地利用类型、坡度、土壤类型、植被覆盖度、降雨、GDP、人口密度。因子交互作用对土壤侵蚀的解释力均大于单一因子的解释力。(2) 土地利用类型与坡度与的交互结果对土壤侵蚀解释力最强。(3) 迁西县易产生土壤侵蚀区域分别为：工矿用地和未利用土地、>35°和15°～25°坡度、棕壤、潮土和钙质石灰土及植被覆盖度集中在17%～50%的区域。因此，区域水土流失治理的重点方向为：加强工矿用地生态修复整治与未利用土地生态保护，植树造林增强其地表植被覆盖度；控制板栗经济林扩张规模，恢复林下草皮，易发生土壤侵蚀的坡耕地退耕还林草。此外，注重小流域水土流失综合治理。