沱江流域土壤侵蚀动态变化及驱动力分析

钟旭珍, 张 素, 吴瑞娟, 敬远兵, 门雷雷, 周 婷

(1.内江师范学院 地理与资源科学学院, 四川 内江 641100; 2.中煤科工生态环境科技有限公司, 北京 100013)

土壤侵蚀是指发生在外营力作用条件下的土体迁移过程[1]。土壤侵蚀严重影响了工农业的正常发展，造成土壤肥力下降、致使生态环境恶化、破坏了水利交通工程设施等[2]。目前全球土壤侵蚀形势严峻，我国土壤侵蚀现状也不容乐观，有研究表明我国土壤侵蚀面积高达32%，四川省约有24.90%的地区存在土壤侵蚀现象，研究土壤侵蚀的原因和演变规律，有效地进行防治，保护生态环境已经成为全球关注的焦点[3]。

对土壤侵蚀的研究主要是将RS和GIS技术与数学模型相结合进行，常用方法如通用侵蚀方程(USLE)[4-7]、修正通用方程(RUSLE)[8-10]、水蚀预报模型(WEPP)[11-13]和中国土壤流失方程(CSLE)等[14-15]众多计算模型，其中修正通用方程(RUSLE)是具有统一规范性的国际上应用较为广泛的模型，具有结构简明、计算因子含义明显的特点。然而，由于不同区域间的复杂性和特殊性以及研究方法的多样性，特别是针对属于丘陵地形的沱江流域，仍然存在一些问题需要进一步探索。例如对土壤侵蚀较为典型的地区研究较多，而对于丘陵、坡耕地等地区的土壤侵蚀研究较少；对土壤侵蚀的现状研究较多而对于其影响因素归因的研究较少，而由王劲峰等[16]提出的地理探测器是用于探测空间分异性，以及揭示其背后驱动力的一种统计学方法，近年来广泛应用于自然科学、社会科学、环境科学和人类健康等方面。目前已有少数学者将其用于水土流失驱动力的探测，并取得了良好的效果，得到了学者们的认可，如陈锐银等[15]基于CSLE模型及地理探测器对四川省水土流失重点防治区土壤侵蚀进行研究；黄硕文等[17]运用RULSE模型及地理探测器分析河南省近10 a来土壤侵蚀时空变化及归因；王欢等[18]基于地理探测器研究不同喀斯特地貌形态类型区的土壤侵蚀定量归因等。

针对沱江流域，《四川水土流失综合治理蓝图绘就》中将其划定为国家级水土流失重点治理区，同时也是长江上游重要生态屏障最大的环境风险带。2018年，《中共四川省委关于全面推动高质量发展的决定》将沱江流域综合环境治理作为四川省“十三五”时期全面推动高质量发展的重要组成部分，但目前，极少有人研究沱江流域的水土流失，更没有人基于地理探测器研究其土壤侵蚀的驱动因子。综合以上分析，研究基于RUSLE模型完成沱江流域2000—2018年土壤侵蚀模数定量计算与分级，对其时空动态演变规律进行分析，并借助地理探测器对土壤侵蚀进行定量归因研究，有利于了解流域内土壤侵蚀状况，为改善流域生态环境及对该国家级水土流失重点治理区的防治工作提供参考依据。

1 研究区概况及数据来源

1.1 研究区概况

沱江流域位于长江上游，四川省中部，坐标为(103°40′—105°50′E,28°53′—31°41′N)，范围涉及德阳、成都、资阳、内江、自贡、泸州、重庆等11个地市，其东南部为重庆市的荣昌县、大足县、双桥区等区县，总面积约2.79万km2，是四川省人口密度最大、城市分布最密集、经济社会最发达的地区，经济总量占全省的30.8%，人口占全省的26.2%，有全省经济发展“金腰带”之称。流域气候属于亚热带湿润气候，多年平均气温在17℃左右。流域形状瘦长，地形以丘陵平原为主，地势由西北向东南倾斜，起伏较小，植被覆盖的空间差异大，流域内森林覆被率仅6.1%，为四川各河中最低者。近年来，由于以城市化为主的人类经济社会活动影响的加剧、气候变化以及人口增长，引发了植被退化、水土流失等一系列生态问题[19]，沱江流域长期是四川省环境治理的重点区域，开展该流域土壤侵蚀时空演变及定量归因研究，有利于了解流域生态环境状况，并进行有针对性地生态整治。

1.2 数据来源

研究数据主要有：(1) 研究区1∶10万土地利用类型矢量数据、1∶100万中华人民共和国地貌图集、中国年度植被指数(NDVI)空间分布数据集、中国人口、GDP空间分布km格网数据集等，来源于中国科学院资源环境数据云平台(http:∥www.resdc.cn)；(2) 研究区数字高程模型(DEM)，空间分辨率30 m，来源于地理空间数据云平台；(3) 2000—2018年沱江流域及周边地区28个基础气象站的各月累计降雨量，来源于中国气象数据网(http:∥data.cma.cn/)；(4) 1∶100万中国土壤数据库，来源于中国科学院南京土壤所；(5) 其他数据。研究区河流矢量图、行政区划界线、社会经济要素等，来源于中国科学院资源环境科学数据中心及《四川省统计年鉴》。

2 研究方法

2.1 RUSLE模型

土壤侵蚀是指土壤及成土母质被水力、风力、冻融或重力等外力作用下，土壤物理结构或化学成分发生破坏、搬运等的过程[20]。研究采用修正通用土壤侵蚀方程RUSLE，该模型的表达式如下：

A=R·K·L·S·C·P

(1)

式中：A为土壤侵蚀量[t/(hm2·a)],乘以100后单位转换为t/(km2·a);R为降雨侵蚀力因子[(MJ·mm)/(hm2·h·a)];K为土壤可蚀性因子[(t·hm2·h)/(MJ·mm·hm2)];LS为坡长坡度因子;C为植被覆盖因子;P为水土保持因子;LS,C,P为无量纲。

为保证各因子在GIS中叠加重合，研究将各因子图层统一设置为Albers投影系统，空间分辨率均为30 m×30 m。

2.1.1 降雨力侵蚀因子 降雨侵蚀力能较客观地描述出降雨对地区土壤侵蚀的影响作用。研究利用Wischmeier模型完成沱江流域降雨侵蚀力的计算，并采用反距离权重插值模型完成其空间插值。Wischmeier表达式如下：

(2)

式中:R为降雨侵蚀力;Pi为月降雨量;P为年降雨量。

2.1.2 土壤可侵蚀因子 土壤可侵蚀性是驱动土壤发生侵蚀变化的内因，不同类型的土壤因为粒径、黏性和化学成分等存在差异，其可侵蚀性也存在明显的不同。土壤类型不同，K值不同。K值越大，土壤越容易受到侵蚀，K值越小，则土壤受侵蚀的可能性越小。研究使用1∶100万中国土壤数据，采用EPIC模型[21]完成K因子的计算。

K=[0.2+0.3e-0.0256SAN(1-SIL/100)]×

(3)

式中：SAN为砂粒含量%；SIL为粉粒含量%；CLA为黏粒含量%；C为有机质含量%；SNI=1-SAN/100。

2.1.3 坡度坡长因子 地形是影响区域土壤侵蚀类型和强度形成与变化的主要因素之一。通常情况下，地势平缓的地区坡度和坡长值相对较小，土壤侵蚀程度也相对较轻。坡长因子L和坡度因子S一般基于DEM进行提取。本研究采用刘宝元等[21]提出的坡长模型完成该地区坡长因子的计算，表达式如下：

(4)

坡度因子的计算采用刘斌涛等[23]针对西南山地区提出的10°～25°范围的S修正模型完成，模型如下：

(5)

式中:L和S分别为坡长和坡度因子;λ为坡长；θ为坡度;m为坡度计算指数。

2.1.4 植被与水土保持因子 植被覆盖管理因子是一项用于描述植被覆盖状况对土壤侵蚀综合影响作用强弱的指标，水土保持因子是表示实施专门水土保持措施后的土壤流失量与未实施任何专门措施而种植时的土壤流失量之比，C和P的值大小变化均在0～1。

研究基于区域相似性，根据相关研究成果[24-27]，结合区域特征，确定沱江流域地表植被覆盖因子C和水土保持因子P的经验值(表1)。将不同土地利用类型的C值和P值在ArcGIS 10.2中赋值于土地利用类型矢量属性表中，然后利用要素转栅格功能转为栅格数据，得到流域植被覆盖因子和水土保持因子栅格图。

表1 沱江流域各土地利用类型的CP因子赋值

2.2 地理探测器

地理探测器即Geographical Detector，简称Geodetector，是一种统计工具，可以探测地理现象的空间分异性并揭示其背后驱动力，既适用于点数据，也适用于面数据，对于面数据一般要先进行重分类和离散化处理。其包括4个探测器[16]。

因子探测器，用于探测因变量Y的空间分异性，及某个因素X对因变量Y空间分布的解释力，其大小可用q值来度量(0≤q≤1)，本研究中，若某影响因子的q值越接近于1，说明该因子对土壤侵蚀的解释力越强，反之越弱[16-18]。

交互探测器用于识别不同因子之间的交互作用，即评估两个因子共同作用时是否会增加或减弱对因变量的解释力，或这些因子对因变量Y的影响是相互独立的[16]，其交互方式见表2[17]。通过交互探测器，可探究不同影响因子交互作用时对土壤侵蚀空间分布的影响。

表2 地理探测器交互作用方式

风险探测器用于判断2个子区域间的属性均值是否有显著的差别，可用于识别土壤侵蚀高风险区域[17]。生态探测器则是用于比较2个因子对属性Y的空间分布的影响是否有显著的差异，以F统计量来衡量，可比较影响因子对土壤侵蚀空间分布的影响是否有显著差异[15]。

研究选取了海拔、坡度、植被覆盖度、地形地貌、GDP、人口、土地、降水8个因子作为土壤侵蚀的影响因子，分别为X1—X8，土壤侵蚀为Y，作为地理探测器的数据输入。由于这些因子及土壤侵蚀都是连续变量，而地理探测器的输入数据要求为类别数据，根据王劲峰等[16]提出的数据离散化方法，首先对这些连续变量进行重分类处理，将海拔按照<500 m，500～1 000 m，1 000～1 500 m，>1 500 m分为 4类；坡度按照<5°，5°～10°，10°～15°，15°～20°，20°～25°，25°～30°，30°～35°、>35°分为8类；植被覆盖度按照<0.3，0.3～0.4，0.4～0.5，0.5～0.6，0.6～0.7，0.7～0.8，0.8～0.9，0.9～1分为7类；地形地貌使用《中华人民共和国地貌图集(1∶100万)》中的类别编号，分为24类；GDP、人口数据是来源于中国科学院资源环境数据云平台的由徐新良[28]处理的空间分布km网格数据，为栅格数据类型，利用GIS裁剪出沱江流域范围的GDP、人口km格网栅格数据，再运用GIS重分类功能，将其等间距分为8类作为地理探测器的输入变量；土地利用类型分为耕地、林地、草地、水域、建设用地、未利用地6类；降雨量等间距分为9类；土壤侵蚀强度按照土壤强度侵蚀分级标准分为6类。接下来通过ArcGIS创建渔网，设置其间隔为1 000 m，得到研究区共26 530个点，运用Spatial Analyst—提取分析—采样工具，输入栅格为所有重分类的X和Y，采样点为研究区渔网点，将自变量和因变量的值提取到点，对数据进行筛选，去掉小于0的值，最后得到22 487个点，作为地理探测器输入数据。

3 结果与分析

3.1 土壤侵蚀计算结果

根据RUSLE模型各因子计算方法，得到各因子计算结果，运用ArcGIS 10.2栅格计算器，将各因子相乘，得到沱江流域土壤侵蚀模数，并将单位转换为t/(km2·a)，同时参照水利部颁布的《土壤侵蚀分类分级标准》(SL190—96)规定的分级指标[29]，将该地区2000—2018年土壤侵蚀分为微度、轻度、中度、强度、极强烈和剧烈6个等级(图1)

图1 沱江流域2000-2018年土壤侵蚀强度等级

3.2 土壤侵蚀时空变化特征

表3为沱江流域2000—2018年土壤侵蚀各等级栅格数比例，可以看出流域主要以微度侵蚀为主，占比均大于50%，主要分布在低矮的平原和丘陵等耕地地区，且微度侵蚀比例随时间在逐渐增大，2018年相比2000年增加了7.03%，其中2010年微度侵蚀比例相比2005年有所减少，主要转化为了轻度侵蚀，原因主要是2010年降雨侵蚀力相对较大，另外考虑到2008年汶川地震，研究区北部岷山一带植被覆盖、地形地貌等受到影响，从研究结果可以看出到2010年土壤侵蚀改善程度不明显，微度侵蚀面积减少，但随着震后重建和恢复，研究区土壤侵蚀得到逐步改善，特别从图1可以看出，岷山一带改善较为明显，说明植被恢复及水土保持治理措施具有成效；剧烈侵蚀的比例随时间在逐渐减小，2018年相比2000年减小了2.00%，该侵蚀等级主要分布在西北部的岷山、中部的龙门山等海拔较高的山地地区；中度及以上侵蚀等级整体上都在随时间减小，说明沱江流域土壤侵蚀逐渐得到改善，尤其是岷山、龙门山地区土壤侵蚀改善明显，这主要是因为水土保持措施的实施。

借助IDRISI软件的CA-Markov模型，以2010年为时间分界点，得到2000—2010年和2010—2018年土壤侵蚀面积转移矩阵(表4—5)，并利用ArcGIS空间分析Combine功能，分析各阶段土地的空间分布变化(图2)。

表3 沱江流域2000-2018年不同土壤侵蚀等级面积比例 %

表4 2000-2010年沱江流域土壤侵蚀强度转移矩阵

表5 2010-2018年沱江流域土壤侵蚀强度转移矩阵

注：0代表土壤侵蚀强度等级不变，1代表土壤侵蚀强度等级升高，-1代表土壤侵蚀强度等级降低。

2000—2010年沱江流域微度、轻度、中度、强度、极强度、剧烈侵蚀的稳定率分别为78.53%，55.53%，24.93%，15.24%，12.73%，15.52%，2010—2018年稳定率分别为81.73%，53.72%，30.00%，16.94%，18.83%，42.84%。从表中可以分析到大部分土壤侵蚀强度向更低等级转变，2000—2010年有36.84%的轻度侵蚀转变为微度侵蚀，57.08%的中度侵蚀转变为轻度侵蚀，39.89%的强度侵蚀转变为轻度侵蚀，52.70%的极强度侵蚀转变为中度侵蚀，26.40%的剧烈侵蚀转变为强度侵蚀。2010—2018年有43.66%的轻度侵蚀转换为微度，62.32%的中度转为轻度，60.31%的强度转为中度，51.72%的极强度转为强度，33.91%的剧烈转为极强度。

结合图2，2000—2010年土壤侵蚀强度等级不变的区域主要分布在沱江流域东南部及北部岷山与龙门山之间，面积占比约为81.00%；土壤侵蚀强度等级降低的地区主要分布在龙门山一带，面积占比约为12.73%；土壤侵蚀强度等级升高的地区主要分布在龙泉山一带，面积占比约为6.7%。2010—2018年土壤侵蚀强度等级不变的区域主要分布在沱江流域东南部，面积占比约为84.69%；土壤侵蚀强度等级降低的地区主要分布在岷山及中部等植被覆盖较高的地区，面积占比约为12.34%；土壤侵蚀强度等级升高的地区面积占比约为2.96%。可以看出，两个时间段内土壤侵蚀强度等级降低的范围均大于侵蚀等级升高的范围，说明研究区土壤侵蚀整体上趋于改善。

3.3 土壤侵蚀影响因素分析

从因子探测器结果(表6)可以看出，不同土壤侵蚀影响因子对土壤侵蚀的解释力具有差异性，其中解释力最强的为坡度，达到48.32%，说明坡度对沱江流域土壤侵蚀具有主导性的影响，这与黄硕文等[17]对河南省土壤侵蚀定量归因的研究相似。接下来依次是地形地貌、海拔、土地利用、年降雨量、人口、植被覆盖度，解释力最小的影响因子为GDP。沱江流域大部分地区属于低矮的平原丘陵，北部的岷山、中部龙门山和龙泉山将流域划分为了山地丘陵平原相间分布几部分，这样的地形地貌从宏观上控制着地表过程的发生和发展，对土壤侵蚀也具有主导性作用，这与王欢等[18]进行的基于地理探测器的喀斯特不同地貌形态类型区土壤侵蚀定量归因研究相符合。说明本研究的结论相对可靠。由于沱江流域土地利用主要以耕地为主，2018年耕地面积占比达到78.68%，林地仅占12.65%，从而决定整个流域的植被覆盖度较低，因此土地利用和植被覆盖对土壤侵蚀的解释力明显低于了坡度、地形地貌和海拔。降水对土壤侵蚀也有着重要影响，对于以平原、丘陵地形为主的沱江流域，极端气候较少，降雨对土壤侵蚀的影响较为均匀，其解释力低于由地貌形态控制的坡度、海拔等。降水是否会引起水土流失，植被覆盖的拦截有重要作用，但沱江流域土地利用又以耕地为主，流域内林草地分布较少，林地和草地分布面积仅占15%左右，要提高植被覆盖度更多的靠人类活动，但在这样宏观地形地貌和土地利用结构的控制下，即使国家实行退耕还林还草政策，对土地利用结构也主要起微调作用，所以土地利用的影响大于降水和植被覆盖。此外，土壤侵蚀的影响因素错综复杂，主要以自然因素为主，但人类社会经济活动也会对土壤侵蚀产生加剧或减缓作用，因此研究选用了人口和GDP两个社会经济因子来进行定量归因探讨，结果表明人口对土壤侵蚀的解释力大于GDP，对于植被覆盖较低的沱江流域，人口活动对土壤侵蚀影响甚至大于植被覆盖的影响，这表明，在进行水土流失保护工作的时候，针对不同的地区要同时注重自然因素的影响和社会经济因素的影响，沱江流域土壤侵蚀的控制，要特别注意人类活动的影响。

表6 土壤侵蚀因子的p值和q值

从交互作用探测器结果(表7)可以看出，任意两个影响因子之间的交互作用对土壤侵蚀的解释力都大于单因子作用。这表明，要综合考虑各种影响因素对土壤侵蚀的影响。其中坡度与土地利用的交互作用解释力最大，达到61.58%，意味着不同坡度下土地利用的土壤侵蚀相差悬殊。坡度与年降雨量的交互作用解释力达到52.32%，这主要是由地貌形态决定的，沱江流域地貌形态主要由几大山脉主导，降雨量在低矮的丘陵对土壤侵蚀的解释力较小，当与坡度交互，特别是与坡度大的地区交互时，其对土壤侵蚀的解释力则猛增。坡度与植被覆盖和年降雨量的交互作用，地貌形态与植被覆盖和年降雨量的交互作用，都明显大于植被覆盖和年降雨量的单独作用，这表明在海拔、坡度较高的地方，降雨量和植被覆盖对土壤侵蚀的解释力更大，因此，要控制沱江流域的土壤侵蚀，在地形起伏较大的地区进行植被建设是关键。

表7 土壤侵蚀因子交互作用下的q值

因子探测和交互探测表明单个因子及不同因子之间交互作用都会对土壤侵蚀产生不同影响，而通过风险探测器可探测同一影响因子不同分层间的土壤侵蚀平均值是否有显著差异，识别土壤侵蚀的高风险区(表8)。探测结果表明坡度>35°的区域为土壤侵蚀风险最高的地区，主要分布在北部岷山地区，对于坡度的8个分层，土壤侵蚀强度随坡度等级的增加逐渐增强。土地利用的高风险区为草地，也主要分布于北部岷山地区。年降雨量和植被覆盖高风险区均主要分布于北部海拔较高起伏较大的山区，极易发生水土流失。海拔的高风险区为>1 500 m的地区。对于人口，高风险区主要分布于人口密度小的地区，人口少的地方主要在崎岖的山地和未利用地，这些地方都容易发生水土流失。

表8 各影响因子侵蚀高风险区域

生态探测结果表明，坡度与土地利用、地形地貌与土地利用、地形地貌与海拔等对土壤侵蚀的影响具有显著差异，说明这几个因子主导着沱江流域是土壤侵蚀，应该充分根据流域的地貌形态，在不毁坏原始生态结构的基础上，加强对坡度大、海拔高、地形起伏大的地区的植被防护，对于坡度小的耕地、草地，应该合理调整农业产业结构和土地利用方式，种植水土保持功能更强的农作物，注重耕作区的造林防护。另外，沱江流域水域面积较广，沱江为长江的一级支流，极易发生洪涝灾害引起河流两岸水土流失，因此，应该加强河流两岸的水土保持工程，将两岸的耕地改为固坡固土更强的作物。

4 结 论

(1) 沱江流域土壤侵蚀以微度为主，土壤侵蚀强度较低的等级主要分布于低矮的坡耕地，强度及以上侵蚀区主要分布于坡度、海拔较高的地形地貌区，如岷山、龙门山、龙泉山等。微度侵蚀比例随时间在逐渐增大，2018年相比2000年增加了7.03%，剧烈侵蚀的比例随时间在逐渐减小，2018年相比2000年减小了2.00%，2000—2010年和2010—2018年土壤侵蚀等级中微度的变化稳定率都大于75%，侵蚀强度等级降低的范围均大于侵蚀等级升高的范围，说明研究区土壤侵蚀整体上趋于改善。

(2) 研究选取的8个影响因子对土壤侵蚀的解释力具有差异性，解释力最强的为坡度，这主要是受沱江流域特殊的地形地貌和土地利用方式影响，人口、GDP等社会经济因素也会对土壤侵蚀产生影响，因此，在进行土壤侵蚀归因分析时，应综合考虑自然和人类社会经济活动。

(3) 因子间交互作用均能增强对土壤侵蚀的解释力，坡度与土地利用和降雨量的交互最为显著，因此应该禁止陡坡开耕及坡耕地的管理，在降雨量较大的陡坡地区应该采取相应的防治措施，减少因降雨引起的自然破坏和经济损失。在土壤侵蚀的治理工作中，应考虑地貌类型的宏观控制作用，针对不同影响因子内部重点治理高风险区，主要为坡度大于35°，海拔大于1 500 m等起伏较大的山区。

运用地理探测器进行土壤侵蚀的定量归因还处于初步探索阶段，在今后的研究中需进一步探讨各个影响因子内部对土壤侵蚀影响的差异性，为土壤侵蚀的治理提供更加准确的依据。