基于地理探测器的岷江上游地区土壤侵蚀变化

2022-03-14刘敦龙李玉霞

水土保持研究 2022年2期

姚昆, 周兵, 何磊, 刘斌，罗涵，刘敦龙，李玉霞

(1.西昌学院资源与环境学院, 四川西昌 615000;2.国家气候中心气候服务室, 北京 100081; 3.成都信息工程大学软件工程学院, 成都 610225;4.四川省信息化应用支撑软件工程技术研究中心, 成都 610225; 5电子科技大学自动化工程学院, 成都 611731)

土壤侵蚀会造成河道淤塞、耕地质量下降和促进泥石流频发等环境问题[1-3]。区域的地形地貌、气候环境和土壤类型等条件是促进土壤侵蚀形成和改变的先天因素，而后天不合理的人类活动，又会驱动土壤侵蚀强度和分布面积的改变[4-6]。截至2020年，四川省仍有12.1万km2的地区存在水土流失现象，其中水力侵蚀就占据侵蚀总面积的94.21%[7]。因此，加强区域土壤侵蚀(特别是水力侵蚀)的实时监测，掌握其动态变化规律并分析背后的驱动因素是实现侵蚀综合治理重要基础。

为实现土壤侵蚀状况的定量评价，学者们将RS和GIS技术应用于该领域，并且提出了通用土壤流失方程(USLE)[8-10]、修正的土壤侵蚀模型(RUSLE)[11-12]、中国土壤流失方程(CSLE)[13-15]和四川省土壤流失方程(SCSLE)等[16]众多模型。田宇等[5]利用RUSLE模型完成1990—2015年三峡库区土壤侵蚀定量评价后，还对其时空变化规律和驱动因素进行了探讨；康琳琦等[12]以RUSLE模型为基础，分析了青藏高原近30 a来土壤侵蚀的时空变化规律；姜琳等[17]基于该模型完成了岷江上游流域2000—2010年土壤侵蚀时空格局的动态变化分析，并探讨了侵蚀状况与高程、坡度和土地利用类型的关系；龚雪梅等[18]基于RUSLE模型完成了1995—2014年岷江上游地区土壤侵蚀的敏感性评价，并分析了“汶川地震”和土地利用类型改变对其变化的影响。

虽然以上成果均为区域土壤侵蚀综合治理提供了重要的参考价值，但是也存在部分不足需要进一步完善。如以往成果在进行驱动力探索时，大都以定性分析或数理统计(相关性或回归分析)[5,17-18]为主要方法；前者存在较强的主观性未能实现驱动因素较准确的定量分析，而后者实施的前提条件是假设在整个时间序列中，土壤侵蚀与驱动因子间存在显著线性关系，而现实中此关系却不一定存在[19-20]；同时，以上两种方法都未将各因子交互产生的协同作用对侵蚀改变的驱动作用进行分析；而针对以上两点不足地理探测器却较好解决[3,20]。综上，研究以RUSLE模型为基础，完成岷江上游地区2000—2018年土壤侵蚀状况定量评价；借助斜率变化模型完成其时空变化规律分析；以地理探测器为工具对影响区域土壤侵蚀变化的驱动力进行探索。期望能为该地区土壤侵蚀综合治理提供科学理论参考。

1 研究区概况

岷江上游地区地处103°32′—104°15′E，30°45′—33°09′N，包括松潘、黑水和茂县等5个行政区，幅员面积约2.4万km2，海拔表现为东部比西部相对较高，地形起伏差异相对明显，降水主要集中于夏季(5—10月)，季节性变化明显，全域植被生态系统相对复杂，阔叶林、针叶阔叶林和灌丛等为主要植被类型，土壤则主要包括棕壤、粗骨土和草毡土等，区域地质构造复杂，地质灾害相对频发。

2 数据来源与研究方法

2.1 数据来源

基础数据如下：(1) 1980—2018年岷江上游地区及周边气象站逐日累计降水量数据，来源于国家气候中心；(2) 各年份1∶10万土地类型矢量成果，由中国科学院资源环境科学中心(http:∥www.resdc.cn)提供；(3) 90 m分辨率DEM数据，来源于地理空间数据云平台；(4) 研究区30 m土壤可侵蚀性(K)因子，由中科院山地所刘斌涛课题组提供；(5) 1∶100万中国土壤类型矢量数据，来自于中国土壤数据库；(6) 各年份NDVI为MOD13Q1成果，来自USGS。

2.2 研究方法

2.2.1 RUSLE模型通用土壤流失方式表达式如下：

A=R×LS×K×C×P

(1)

式中:A为侵蚀模数[t/(hm2·a)];R为降雨侵蚀力[(MJ·mm)/(hm2·h·a)];K为土壤可侵蚀[(t·hm2·h)/(MJ·mm·hm2)];LS为坡度坡长；C,P分别为植被管理与水土保持因子,均无单位;侵蚀模数乘以100单位转换为t/(hm2·a)。

其中，R因子以章文波等[21]提出的日降雨量侵蚀力模型计算；K因子为基于EPIC方程[22]计算，结合张科利等[23]提出的修正模型进行校正，最后采用GIS空间插值与平滑优化等算法对数据进一步优化以提高精度的成果(图1)；L因子则以Liu等[24]提出模型计算；S因子则采用刘斌涛等[25]修正的方程(图1)；C因子以蔡崇法等[26]基于植被覆盖度的模型计算(图2)，而植被覆盖度则采用像元二分模型[27]；P因子参考已有成果[28-29]结合区域实际分别赋水田0.15，旱地0.35，草地0.8，林地和未利用为1及其他为0。

2.2.2 斜率变化斜率可以实现土壤侵蚀与时间的回归拟合，对其未来发展趋势进行预测[30]，本文以最小二乘法实现逐像元侵蚀模数随时间的变化分析。表达式如下：

(2)

式中:K为斜率;n为年份;Ai为第i年侵蚀模数。K<0代表土壤侵蚀状况有良好的发展趋势;相反,K>0则侵蚀状况有加重的发展倾向。同时,采用F检验,以p为0.05为分界点实现斜率变化的显著性检验。

2.2.3 地理探测器地理探测器是一种能揭示地理现象空间分异差异和分析自变量与因变量之间相互作用的数学模型。相比传统的相关系数模型，其不仅能实现定量数据分析，也能完成定性数据的处理，同时又能对各变量之间的交互作用进行分析[20,31]。

图1 坡度坡长和K因子

图2 研究区C因子

因子探测器主要用于探测各个自变量X对因变量Y空间分异的影响力度[20,31]，采用q值进行表征：

(3)

(4)

式中:q为因子解释力,数值在[0,1];h为变量X或Y的分类区间;Nh,N为各层或者全域的分区数量；SSW,SST为层内和全区总方差；q数值越大表明X对Y的影响力越明显。

交互探测器可用于分析当各不同自变量X因子存在交互作用时，其交互作用对因变量Y的驱动解释力是增强或者减弱[20,32-33]。其先分别单独计算自变量因子X1和X2分别对Y的解释力q(X1)和q(X2)，其次计算q(X1∩X2)，最后将这3个变量的q值进行比较以此完成交互影响作用的判定(表1)。

本文选取多年平均NDVI(X1)、高程(X2)、多年平均降雨量(X3)、土壤类型(X4)、土地利用类型(X5)和坡度(X6)共计6个因子为自变量，而将多年平均土壤侵蚀模数为因变量。

表1 交互作用判定

基于王劲峰等[20]提出的数据离散化和先验知识，在参考已有成果[4,19,30]的基础上利用“自然间断点分级”将高程、坡度和多年平均降雨量均分为9类，而多年平均NDVI分为13类；土壤类型参照《1∶100万中华人民共和国土壤图》将其分为17类，土地利用类型按照耕地、林地和草地等分为6个一级类。同时，以1 km×1 km格网完成该地区中心样点的规则取样。

3 结果与分析

3.1 土壤侵蚀评价

研究基于USLE模型完成2000—2018年各年份土壤侵蚀模数的计算，并按《土壤侵蚀分类分级标准(2008)》文件[34]完成各年份侵蚀程度的等级划分(图3)。

结合图3分析，该地区土壤侵蚀强度空间分布变化规律与姜琳等[17]基本吻合，间接证明研究成果基本可靠。以2018年岷江上游地区侵蚀强度分级结果为例，该地区侵蚀以微度和轻度为主，它们占据了全域总面积的77.59%，这些土壤侵蚀现象并不明显的地区在整个研究区的分布范围最广，分布在全域的大部分地区。结合土地利用类型和地形地貌数据分析可知，微度和轻度区地貌大多以中—高海拔的山地为主，而这些地区土地类型则又主要为有林地、灌木林地和高度盖度草地等，这些区域植被生态环境相对良好且覆盖度相对较高，良好的植被生态环境对调节气候变化和生态环境好转又起到了有利的促进作用，微度侵蚀全域分布范围广可能就是得益于此。

图3 2000-2018年土壤侵蚀强度分级

相反，中度及以上程度在全域分布范围仅占据全域总面积的22.41%，它们绝大部分分布于高海拔且地形起伏差异明显的高山和极高山地区，这些地区呈现出植被覆盖度相对较低，地形陡峭，特殊的气候易引发岩石风化和土壤剥离搬运等现象，因此土壤侵蚀程度相对明显。此次，它们还有极小部分分布于岷江、黑水河和杂谷脑河等江河水系两侧，这些地区虽然地势相对平缓，但它们是社会经济发展的核心地带，也是受人类活动影响相对明显的区域，全域主要的农牧生产均主要集中于此，促成了这些地区土壤侵蚀现象相对明显现象的产生。

3.2 趋势变化分析

结合岷江上游地区土壤侵蚀斜率变化(图4)计算结果分析可知：2000—2018年内，土壤侵蚀呈现好转(K<0)的地区全域面积占比为75.60%；而整个岷江上游地区仅有25.40%的侵蚀呈现出严重化发展趋势。结合显著性检验结果又可知：整个时段内，有21.14%的地区土壤侵蚀呈现明显好转的发展状态，它们主要分布于松潘县东部地区、岷江、黑水河以及杂谷脑河等江河水系两侧地形相对平缓的区域；侵蚀状况未发生明显改变(好转或严重化)的地区全域面积占比最大为76.99%，在全域的大部分范围均有分布；仅有1.87%的地区呈现出土壤侵蚀明显严重化的发展倾向，结合地形地貌数据分析可以发现，它们主要分布在海拔4 000 m以上的高山地区。

图4 斜率及显著检验

3.3 侵蚀驱动因子探测

3.3.1 侵蚀驱动因子探测结合(表2)对土壤侵蚀驱动因子的关系进行分析可知：各个因子对区域土壤侵蚀空间分布格局变化的驱动作用存在不同差异；其影响力大小关系为多年平均NDVI>多年平均降水>高程>土壤类型>土地利用类型>坡度；6个因子中多年平均NDVI的q值最大为0.331 4，说明其对土壤侵蚀空间分布格局形成和改变的影响力最大，坡度的q值最小仅为0.017 2，和其他5个因子相比坡度对侵蚀强度空间分布格局变化的驱动作用最不明显。

表2 驱动因子探测结果

结合q值变化与实际分析可知：岷江上游地区森林和中高覆盖度草地分布范围相对较广，全域95%以上为林草区(其中，中高植被覆盖区又占据林草区总面积的80%以上)，全域植被覆盖状况整体相对良好，植被对土壤侵蚀的形成与变化起到相对明显的抑制作用；区域海拔差异明显，地形起伏相对较大，降雨又主要集中于下半年(5—10月)和海拔2 500 m以下的河谷地带，夏季暴雨相对频发，降雨强度也相对较大，是促进岷江上游地区土壤侵蚀面积扩张和程度加深的主要因素；此外，多年平均植被覆盖度、高程和多年平均降雨量这3个因子的解释力均明显远大于土地利用类型；综上，客观程度可以判定，自然因子是影响区域土壤侵蚀空间分布格局变化的主要驱动因素。

3.3.2 交互探测交互探测主要用于判定，若两个驱动因子产生交互关系时，其交互作用对土壤侵蚀变化的影响作用为增强、减弱或者相互独立。通过交互探测器完成岷江上游地区土壤侵蚀与各驱动因子的交互探测，结果表明各驱动因子间不存在相互独立作用，因子交互作用表现为非线性增强和双因子增强两种。

结合(表3)对各因子之间的交互影响作用进行分析：多年平均NDVI分别与高程、多年平均降雨量或土地利用类型；高程分别与多年平均降雨量、土壤类型、土地利用类型、多年平均降雨量或坡度。以上几类驱动因子产生交互关系时，它们的交互协同作用对土壤侵蚀分布格局变化产生双因子增强影响；其他几类因子产生交互关系时，交互协同作用对土壤侵蚀分布格局变化产生的影响作用为非线性增强；当多年平均NDVI分别与多年平均降雨量、高程或土壤类型产生交互关系时，它们的协同作用对土壤侵蚀分布格局变化的影响力均达到40%以上，客观程度表征了3类组合产生的交互作用对土壤侵蚀分布格局变化影响的贡献作用最大；相反，土壤类型和坡度产生交互关系时，其协同作用对土壤侵蚀分布格局变化的影响力最小，仅产生5%左右的驱动作用；多年平均NDVI与其他驱动因子产生交互关系时，它们产生的协同作用对土壤侵蚀分布格局变化的影响力都在30%以上，均比其他因子间产生交互的影响力大，这也意味着植被覆盖度差异较大的地区土壤侵蚀程度通常也呈现出较明显的差异。

表3 因子交互探测结果

4 讨论与结论

4.1 讨论

岷江上游地区既是四川省西部生态高原建设的重要功能区，也是土壤侵蚀问题相对典型的地区之一。因此，较全面掌握其空间分布特征、时间变化规律和变化驱动因素对科学实现该地区土壤侵蚀综合治理有着重要意义。研究以RUSLE为基础，完成该地区2000—2018年土壤侵蚀状况的定量评价与侵蚀动态变化规律的探索，结合地理探测器对影响岷江上游地区土壤侵蚀空间分布格局和变化的驱动因素进行分析。对研究结果分析可知，本文采用的方法可行，预期目标也基本实现。

相比以往成果，本研究引入用于分析长时间序列数据空间变化规律分析的线性回归模型，其能较好避免由于某时期数据突变而造成整体趋势变化规律出现异常现象的发生，这对较全面掌握区域土壤侵蚀真实的变化规律更具价值；采用地理探测器对影响区域土壤侵蚀空间分布格局变化因素进行探索，其不仅能较有效探测地理现象空间差异性，揭示其背后驱动力，还能分析因子交互协同作用对侵蚀强度空间分布格局变化的影响作用；其较有效弥补了以往成果仅从单因子对驱动力进行分析的不足，进一步促进了对驱动力的较全面掌握。

植被覆盖度(NDVI)、海拔和降水是影响岷江上游地区土壤侵蚀分布及程度改变的主要驱动因素。特别是植被覆盖度，单因子探测结果显示其影响力几乎是海拔和降水的两倍。双因子交互探测结果也显示，植被覆盖度(NDVI)与其他因子产生交互协同作用的驱动力也比其他因子两两交互产生的影响力大。以上两点均表明，植被覆盖度(NDVI)是驱动该地区土壤侵蚀形成与变化的最主要因素，其影响力远远大于其他因子的驱动作用。在整个研究时段内，岷江上游地区土壤侵蚀状况整体得到有效遏制，对祝聪等[30]研究成果分析可知，整个阶段内岷江上游地区植被覆盖度整体处于良好发展的趋势，植被生态环境的改善对区域土壤侵蚀的改变起到了有效的抑制作用。因此，今后对该地区土壤侵蚀进行综合治理时，加强区域植被的保护与努力提高植被覆盖度仍然是一项重要的工作。

由于影响土壤侵蚀变化的因素相对较多，同时受数据收集的限制，本研究仍然存在一定的不足尚需完善。如在计算土壤侵蚀模数时，本文以250 m的栅格尺寸为尺度，虽然利用遥感技术进行时空变化规律分析时大都更侧重宏观角度，该地区范围也相对较大，该栅格尺度未能对动态变化规律分析产生较明显干扰，但若能采用更高分辨率数据进行分析，那研究结果的可靠性会进一步提高；同时，研究在进行驱动力因子选取时，主要选择了部分典型因子，今后的研究中若加入地层岩性、人口和GDP等更多因子进行探讨，那么驱动力探讨的结果也更具价值。