沂蒙山区土壤侵蚀时空变化及其驱动因素分析

黄 鑫,陈 红,徐小任,唐小倩,肖 新,陈世勇,梅大伟,王 梁① (.安徽科技学院资源与环境学院,安徽 滁州 3300;.山东省水土保持与环境保育重点实验室/临沂大学资源环境学院,山东 临沂 76005)

中国是世界上土壤侵蚀最严重的国家之一,2020年全国水蚀和风蚀土壤面积为269.27万km2,约占国土面积的28.0%[1]。北方土石山区是我国水力侵蚀最严重的区域之一,土壤侵蚀已成为制约区域发展的首要问题。沂蒙山区作为北方土石山区的典型代表,其特有的土壤地理环境、极不均匀且强度大的降水分布、复杂多样的地质构造以及不合理的人类活动成为土壤侵蚀的内外驱动力,浅薄的土层流失殆尽造成严重的水土流失问题,进而导致土地退化、自然灾害频发,生态系统濒临崩溃。因此,定量评价该区域的土壤侵蚀现状,科学认识和定量评估自然因素、人类活动及其耦合作用对土壤侵蚀的影响,对遏制该地区的水土流失及推进生态安全建设至关重要。

对于土壤侵蚀这一影响全球生态安全的重要问题,多年来国内外学者在土壤侵蚀机理、模型模拟和驱动因素等方面进行多领域研究,建立了具有不同特色的侵蚀预报模型,主要分为物理过程模型、分布式模型与经验统计模型3类。典型的物理过程模型主要有水蚀预报模型(WEEP)[2]、欧洲土壤侵蚀模型(EUROSEM)[3]等,分布式模型以分布式水文模型(SHE)[4]最为经典。但由于物理过程模型与分布式模型所需参数众多,难以推广应用,因此以通用土壤流失方程(USLE)/修正版通用土壤流失方程(RUSLE)为代表的经验统计模型应用最为广泛。RUSLE在大量实测数据基础上对USLE的局限进行修正,形式简单灵活,参数获取方便,便于与地理信息系统(GIS)和遥感(RS)技术结合,是目前发展较为成熟的土壤侵蚀模型[5-6]。在土壤侵蚀驱动因素研究方面,目前常用的主要方法,如线性与非线性拟合、主成分分析、偏最小二乘回归、结构方程等,对于主控因子识别、因子共线性处理以及剖析因子的直接或间接作用等方面起到了关键作用。但是由于驱动因素的时空变异性及其之间的交互作用致使土壤侵蚀具有时空异质性。地理探测器是王劲峰等[7]基于空间方差理论提出的一种分析地理现象空间异质性的新统计学方法,不仅可以揭示地理现象背后的驱动力,而且能够量化因子间的交互作用、有效识别高风险区,成为当前运用较为广泛的前沿方法之一[8]。

土壤侵蚀在年际之间存在一定变化,目前关于沂蒙山区土壤侵蚀的研究主要集中于单一年份土壤侵蚀,对于长时间序列土壤侵蚀与其驱动因素的交互作用研究关注较少。深入这方面的研究有助于进一步认识土壤侵蚀的机理与变化规律,特别是自然因素与人类活动及其交互作用下的土壤侵蚀定量归因。因此,笔者以北方土石山区典型代表沂蒙山区为对象,基于GIS和遥感RS技术,利用RUSLE探究沂蒙山区2000—2018年土壤侵蚀的时空变化特征,进而运用地理探测器模型量化不同驱动因素之间的交互作用及其对土壤侵蚀的贡献率,揭示不同时期影响土壤侵蚀的主导因素及其交互因素,结果可为沂蒙山区的水土保持规划及国土空间规划提供参考,也可以为该区域乃至北方土石山区生态环境建设提供科学依据。

1 研究区概况及数据来源

1.1 研究区概况

沂蒙山区位于山东省中南部(34°22～36°23′ N,116°34′～119°39′ E),范围涉及临沂市、淄博市、日照市、济宁市、枣庄市等27个县(市、区)(图1),总面积约3.41万km2。研究区的气候类型属于暖温带季风型大陆性气候,年平均气温12～14 ℃,年平均降水量700～900 mm,夏季(6—8月)降水占全年降水量的60%以上。土壤以棕壤和褐土为主,地势西北高东南低,地貌类型复杂,以山地丘陵为主。自然植被以落叶阔叶林为主,但受长期人类活动影响,自然植被只见于海拔较高的山区,以人工植被或次生林为主。

图1 沂蒙山区行政区域及数字高程模型Fig.1 Administrative division and digital elevation model of Yimeng Mountainous Area

沂蒙山区是我国水土流失治理的重点区域,在全国水土保持区划中被定为鲁中南低山丘陵土壤保持功能区,在全国水土流失重点防治区划分中被列为国家级重点治理区,在国家水土保持重点工程建设中被纳入第4期(2008—2012年)和第5期(2013—2017年)规划项目。截至2018年底,国家和地方在沂蒙山区开展了以修建梯田、淤地坝为代表的水土保持综合治理工程和以退耕还林还草为代表的一系列生态修复工程,小流域水土流失治理度达70%以上,林草覆盖度达30%以上,水土流失情况得到了一定改善。

1.2 数据来源

研究所使用的主要数据有:沂蒙山区矢量范围,来源于全国地理信息资源目录服务系统(http:∥www.webmap.cn);1∶100万中国土壤数据集,来源于中国科学院南京土壤研究所(http:∥vdb3.soil.csdb.cn/),包含土壤类型、质地、有机质等数据;数字高程数据ASTER GDEM,来源于地理空间数据云平台(http:∥www.gscloud.cn);沂蒙山区及周边29个气象站点2000—2018年的逐日降雨数据,来源于中国气象局气象数据中心(http:∥data.cma.cn);土地利用数据、地貌形态数据,来源于中国科学院资源环境科学数据中心(http:∥www.resdc.cn);植被覆盖指数NDVI数据,采用NASA发布的MOD13Q1产品(http:∥lpdaac.usgs.gov)。

2 研究方法

2.1 RUSLE模型

采用RUSLE模型评估沂蒙山区土壤侵蚀,计算公式为

A=R×K×L×S×C×P。

(1)

式(1)中,A为年土壤侵蚀模数,t·hm-2·a-1;R为降雨侵蚀力因子,MJ·mm·hm-2·h-1·a-1;K为土壤可蚀性因子,t·hm2·h·MJ-1·mm-1·hm-2;L、S为坡长、坡度因子;C为植被覆盖与管理因子;P为水土保持措施因子。

2.1.1降雨侵蚀力因子(R)

选用WISCHMEIER提出并由ARNOLDUS等[9]修正的基于月平均降雨量和年平均降雨量的经验公式来估算降雨侵蚀力。

式(2)中,pi为月平均降雨量,mm;p为年平均降雨量,mm。

2.1.2土壤可蚀性因子(K)

选用WILLIAMS等[10]在侵蚀/生产力影响模型(EPIC)中土壤可蚀性因子的计算方法进行估算。

(3)

(4)

式(3)～(4)中,SAN、SIL、CLA、Co分别为砂粒、粉粒、黏粒和有机质的含量,%;SNI为计算参数。

2.1.3坡长、坡度因子(L、S)

坡长因子的计算选用WISCHMERIER等[11]提出的计算公式。当坡度≤18%,选用MCCOOL等[12]提出的计算公式;当坡度>18%,选用LIU等[13]改进后的计算公式进行估算。

(5)

(6)

式(5)～(6)中,λ为坡长,m;α为坡长指数;β为坡度修正值;θ为DEM提取的坡度值,(°);标准小区坡长为22.13 m。

2.1.4植被覆盖和管理因子(C)

选用蔡崇法等[14]建立的c值和植被覆盖度之间的回归方程计算C值,核心算法如下

c=(INDV-INDV,s)/(INDV,v-INDV,s),

(7)

(8)

式(7)～(8)中,INDV为归一化植被指数;INDV,s为像元在裸地覆盖状态下的NDVI值;INDV,v为在植被全覆盖状态下的NDVI值。

2.1.5水土保持措施因子(P)

水土保持措施因子P是指采用一定的水土保持措施后的土壤流失量与不采取任何措施的土壤流失量之比。P值范围在0～1之间,0代表不发生土壤侵蚀,1代表未采取水土保持措施或措施完全失效。对于P值,目前通常采用结合实地考察资料给不同土地利用类型赋值的方法确定。结合文献[15-17]和研究区土地利用情况,将水田赋值为0.15,旱地赋值为0.35,林地、疏林地、草地、未利用地赋值为1,水域、居民点、建设用地、裸岩赋值为0。

2.2 地理探测器

研究主要利用地理探测器中的因子探测器、交互探测器和风险探测器,定量分析沂蒙山区土壤侵蚀的驱动因子以及因子间的相互影响,并识别土壤侵蚀高风险区域。

(9)

(10)

ST=Nσ2。

(11)

式(9)～(11)中,q为土壤侵蚀因子的影响强度;L为因变量Y或自变量X的分层;Nh和N分别为层h和全区单元数;σh2和σ2分别为层h和全区的Y值方差;SW为层内方差和;ST为全区总方差。q的范围在0～1之间,q值越大说明土壤侵蚀Y的空间分异性越明显。如果分类是由因子X造成的,q值越大表示因子X对土壤侵蚀Y的解释力越强,反之则越弱。

土壤侵蚀的发生与发展受众多因素的影响和控制,是自然因素和人类活动综合作用的结果。土壤侵蚀的驱动因素可以分为内在因素和外在因素2类,内在因素在较长时间尺度上具有高度的稳定性,其空间异质性显著影响土壤侵蚀的空间分布,而外在因素在较小尺度的时空变化是驱动土壤侵蚀过程的直接诱因。基于此,选取地形地貌、坡度和海拔等内在因素,降雨、土地利用类型和植被覆盖度等外在因素作为自变量X,2000、2005、2010、2015和2018年土壤侵蚀强度作为因变量Y代入地理探测器中进行运算。由于地理探测器的输入变量要求为类别数据,首先需对连续型变量做离散化处理,结合数据离散化方法及先验知识[9,18-19],地形地貌使用《中华人民共和国地貌图集1∶100万》中的类别编号,分为26类,土地利用类型分为耕地、林地、草地、水域、建设用地、未利用地6类,将植被覆盖度分为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.0),海拔分为3类(≤500、>500～1 000、>1 000～1 500 m),年平均降雨量等间距分为9类,坡度分为6类(≤5°、>5°～10°、>10°～15°、>15°～20°、>20°～25°、>25°),并将土壤侵蚀强度及驱动因子分类值赋予1 km×1 km格网点上,作为地理探测器的运行数据。

3 结果与分析

3.1 土壤侵蚀强度的时间变化

沂蒙山区2000、2005、2010、2015与2018年的平均土壤侵蚀模数分别为5 129.01、6 661.76、6 344.00、3 493.23和5 669.03 t·km-2·a-1,2005年侵蚀最为严重,2015年侵蚀程度最低。将研究区5期土壤侵蚀强度按照SL 190—2007《土壤侵蚀强度分类分级标准》划分为微度、轻度、中度、强烈、极强烈、剧烈6个级别(图2),各侵蚀级别面积占比见表1。由表1可知,2000—2018年各侵蚀级别整体变化不大,微度、强烈、极强烈和剧烈侵蚀面积分别增加2.95%、0.81%、0.72%和0.71%,轻度和中度侵蚀面积分别下降4.18%和1.01%。

表1 2000—2018年沂蒙山区不同土壤侵蚀等级面积比例Table 1 Percentage change in the area of different soil erosion levels in the Yimeng Mountainous Area from 2000 to 2018

图2 2000—2018年沂蒙山区土壤侵蚀强度分级Fig.2 Soil erosion intensity level in the Yimeng Mountainous Area from 2000 to 2018

为进一步探讨不同土壤侵蚀强度之间的时间变化特征,分别对5个时期内不同侵蚀强度面积数据进行统计,制成土壤侵蚀强度转移矩阵(表2)。

表2 2000—2018年沂蒙山区土壤侵蚀强度转移矩阵Table 2 Soil erosion intensity transfer matrix for Yimeng Mountainous Area from 2000 to 2018

结果表明,2000—2005年微度、轻度、中度、强烈、极强烈和剧烈侵蚀的稳定率(2005与2000年土壤侵蚀强度等级相同的面积之比)分别为99.82%、0.38%、0.43%、0.00%、59.37%和0.00%。除微度侵蚀外,其他级别侵蚀强度均发生极大变化,整体趋势为低等级向高等级转移。2005—2010年微度、轻度、中度、强烈、极强烈和剧烈侵蚀的稳定率分别为67.46%、3.80%、9.14%、15.20%、11.39%和27.66%。除微度侵蚀外,其他级别侵蚀强度均发生较大变化,整体趋势为高等级向低等级转移。2010—2015年土壤侵蚀强度整体变化趋势仍为高等级向低等级转移,轻度、中度、强烈和极强烈侵蚀转化为微度侵蚀的面积比例分别为77.14%、47.27%、45.45%和38.46%,剧烈侵蚀全部转化为极强烈侵蚀。2015—2018年微度、轻度、中度、强烈、极强烈和剧烈侵蚀的稳定率分别为70.17%、18.55%、10.37%、12.15%、10.18%和30.94%,整体趋势依旧表现为高等级向低等级转移,但也存在低等级向高等级转移的现象:轻度侵蚀中27.62%转化为极强烈侵蚀,中度侵蚀中15.49%转化为极强烈侵蚀,强烈侵蚀中12.55%转化为极强烈侵蚀,极强烈侵蚀中21.25%转化为剧烈侵蚀。2000—2018年微度、轻度、中度、强烈、极强烈和剧烈侵蚀的稳定率分别为98.74%、28.27%、12.50%、0.00%、21.43%和50.00%,轻度、中度、强烈、极强烈和剧烈侵蚀向低等级转移的面积占比分别为67.70%、60.00%、80.00%、78.57%和50.00%。可以看出,2000—2018年土壤侵蚀强度向低等级转化趋势明显,说明研究区土壤侵蚀整体上呈现改善趋势。

3.2 土壤侵蚀强度的空间变化

叠加分析沂蒙山区不同时期土壤侵蚀强度分布图与沂蒙山区行政区划图,获得沂蒙山区2000—2018年各县区土壤侵蚀强度的面积比例(图3),提取各县域平均土壤侵蚀模数,根据SL 190—2007划分土壤侵蚀等级(图4)。微度、轻度、中度、强烈、极强烈和剧烈侵蚀在不同年份各县区均有分布,其中微度和轻度侵蚀分布面积广泛,其他4个等级在不同年份的分布特征差异明显。2000年各县区强烈侵蚀、极强烈侵蚀和剧烈侵蚀面积占比分别为1.05%～8.33%、1.12%～6.75%和1.05%～7.49%,其中山亭区占比最大,微山县最小;中度侵蚀在沂源县和山亭区面积占比较大,分别为18.36%和21.41%。2005年各县区强烈侵蚀和极强烈侵蚀面积占比分别为2.18%～10.86%和1.58%～10.23%;中度侵蚀在沂源县和山亭区的面积占比最大,分别为15.80%和15.40%。2010年各县区强烈侵蚀和极强烈侵蚀面积占比分别为2.33%～14.81%和1.23%～11.45%;中度侵蚀在蒙阴县、费县、沂源县和岚山区的面积占比较大,分别为16.02%、14.59%、14.02%和14.29%。2005和2010年剧烈侵蚀面积占比在山亭区(13.08%和13.97%)和邹城市(10.28%和11.00%)最大。2015年各县区中度侵蚀和强烈侵蚀面积占比分别为1.81%～13.45%和0.78%～6.06%;极强烈侵蚀在莒南县面积占比最大(4.86%);剧烈侵蚀除莒县(5.72%)外其他26个县区面积占比均较小。2018年各县区极强烈侵蚀和剧烈侵蚀的面积占比分别为0.76%～10.02%和1.88%～11.92%;强烈侵蚀在沂源县(10.69%)和山亭区(10.55%)的面积占比较大;中度侵蚀在蒙阴县、平邑县、沂源县和山亭区的面积占比较大,为14.98%～15.34%。

图3 2000—2018年沂蒙山区各县区土壤侵蚀强度的面积比例Fig.3 Percentage change in the area of different soil erosion intensity in each county in the Yimeng Mountainous Area from 2000 to 2018

图4 2000—2018年沂蒙山区各县区土壤侵蚀分级Fig.4 Soil erosion classification in each country in the Yimeng Mountainous Area from 2000 to 2018

综上可知,2000—2018年强烈及以上等级的侵蚀主要分布在西部和北部,集中于山亭区和沂南县一带,2个县区的平均土壤侵蚀模数分别为9 217.55 和8 005.60 t·km-2·a-1,属于极强烈侵蚀,强烈以下等级的侵蚀主要分布在西南部和中东部。随时间变化,土壤侵蚀呈2000—2005年向北和西北方向扩散、2005—2010年向西北方向扩散、2010—2015年向中部扩散、2015—2018年向西和东北方向扩散的趋势。

3.3 土壤侵蚀驱动因素变化分析

运用地理探测器的因子探测器获取2000、2005、2010、2015和2018年地理探测器的q值(表3)。2000年土壤侵蚀驱动因子的贡献力由大到小依次为q(X5)、q(X1)、q(X2)、q(X6)、q(X4)、q(X3),其中土地利用类型起到主导作用,其次具有较强解释力的是地形地貌。2005年土壤侵蚀驱动因子的贡献力由大到小依次为q(X5)、q(X1)、q(X4)、q(X2)、q(X6)、q(X3),其中起主导作用的仍然是土地利用类型,地形地貌、植被覆盖度都具有较强的解释力。2010年土壤侵蚀驱动因子的贡献力由大到小依次为q(X5)、q(X1)、q(X6)、q(X4)、q(X2)、q(X3),其中起主导作用的因子与上一时期相同,其他有较强解释力的因子为地形地貌、年平均降雨量与植被覆盖度。2015年土壤侵蚀驱动因子的贡献力由大到小依次为q(X5)、q(X1)、q(X4)、q(X2)、q(X6)、q(X3),该时期主导因子是土地利用类型,另外具有较高解释力的是地形地貌。2018年土壤侵蚀驱动因子的贡献力由大到小依次为q(X5)、q(X1)、q(X6)、q(X2)、q(X4)、q(X3),该时期主导因子依然是土地利用类型,地形地貌、年平均降雨量都具有较强的解释力,这一时期驱动因子的解释力较上一时期均有明显提高,土壤侵蚀强度的变化也最剧烈。

运用地理探测器的交互探测器获取2000、2005、2010、2015和2018年各个驱动因子对土壤侵蚀的交互作用(表4～5)。2000年土壤侵蚀的交互作用均为双因子增强与非线性增强,没有单独起作用的因子,土地利用类型、地形地貌可以增强各个驱动因子对土壤侵蚀的解释。2005年各个因子对土壤侵蚀的交互作用大多为非线性增强,双因子增强较少,在这一时期解释力较强的交互因子与2000年趋同。2010年各个因子对土壤侵蚀的交互作用大多为双因子增强,非线性增强较少,其中与土地利用类型交互的因子解释力增强,土地利用类型、地形地貌、年平均降雨量和植被覆盖度在这一时期对土壤侵蚀有显著影响。2015年各驱动因子对土壤侵蚀的交互作用与前3个时期相同,仍然还是非线性增强和双因子增强,土地利用类型和地形地貌在这一时期还是对土壤侵蚀有显著影响。2018年各因子对土壤侵蚀的交互作用仍为非线性增强和双因子增强,在这一时期土地利用类型、地形地貌、年平均降雨量可以增强各驱动因子的解释力,坡度对土壤侵蚀有一定的引导作用。

表4 2000—2018年驱动因子交互探测结果Table 4 Results of the interaction of driving factors from 2000 to 2018

表5 2000—2018年沂蒙山区因子交互贡献率Table 5 Interaction contribution of factors in the Yimeng Mountainous Area from 2000 to 2018

运用地理探测器的风险探测器识别2000、2005、2010、2015和2018年土壤侵蚀高风险区,结果显示,不同时期土壤侵蚀高风险区分布较为一致。地形地貌为中起伏的区域土壤侵蚀风险较高,而海拔分类中易发生土壤侵蚀的级别为500～1 000 m。叠加这2个因子分布图发现,2个高风险区有较大重叠,包括沂源县以及蒙山地区(蒙阴、平邑、费县交界处),这部分区域地表起伏大,陡坡较多,侵蚀较为严重。通过对比不同坡度分区对应的土壤侵蚀量可以发现,随着坡度的增加,土壤侵蚀量呈先增加后减少的趋势,在15°～ 20°区间出现拐点,即为坡度的高风险区域。植被覆盖度处于0.4～0.6之间的区域土壤侵蚀风险较高。比较发现,植被覆盖度介于0.4～0.6之间时为临界值,低于这个临界值土壤侵蚀随植被覆盖度的增加而增加,高于这个临界值土壤侵蚀随植被覆盖度的增加而降低,说明植被覆盖度对土壤侵蚀有一定的抑制作用。各土地利用类型中,草地发生土壤侵蚀的风险最高,这是因为研究区草地主要分布在坡度较大的山区,环境条件较差,同时科学有效的水土管理措施实施并不到位,造成土壤侵蚀现象较为严峻。叠加因子分布图发现,年平均降雨量的高风险区在空间上也对应了土地利用类型的高风险区,主要分布在沂源县、沂水县、蒙阴县、平邑县以及山亭区等地,这部分区域山多坡陡,降雨强度大,加之植被覆盖度低,造成侵蚀作用强烈,水土流失严重。

土地利用类型是该研究区土壤侵蚀的主导因子,为进一步分析不同土地利用类型的土壤侵蚀状况,对沂蒙山区不同土地利用类型的土壤侵蚀强度平均值和面积占比进行统计(表6)。土地利用类型侵蚀强度排序为草地、林地、耕地、水域、建设用地和未利用地,不同土地利用类型侵蚀强度差异较大,2005、2010和2018年草地的平均侵蚀强度达到极强烈侵蚀,2010、2015与2018年未利用地与2015年建设用地的平均土壤侵蚀强度处于轻度侵蚀。

表6 2000—2018年不同土地利用类型土壤侵蚀平均模数及面积占比Table 6 Average modulus of soil erosion and area share of different land use types in the Yimeng Mountainous Area from 2000 to 2018

在不同土地利用类型中,侵蚀面积最大的为耕地,占侵蚀总面积的60%,其次是草地、建设用地和林地,未利用地的侵蚀面积较小,侵蚀面积占比不足1%。结合地理探测器分析结果和研究区土地利用类型面积占比和分布情况,草地是研究区亟需治理的土地利用类型,这与钟旭珍等[20]、王猛等[21]进行的土壤侵蚀定量归因研究相符合,耕地不仅侵蚀面积占比最大,而且侵蚀强度突出,应成为研究区应加强预防和治理的重点区域,而对于林地、水域、建设用地和未利用地则需要重点关注并预防进一步扩展。

4 讨论

4.1 土壤侵蚀时空变化

RUSLE模型模拟结果表明,沂蒙山区2000、2005、2010、2015与2018年的平均土壤侵蚀模数分别为5 129.01、6 661.76、6 344.00、3 493.23和5 669.03 t·km-2·a-1,2005年侵蚀最为严重,2015年侵蚀程度最低。2005年土壤侵蚀强度最高,主要原因在于2005年是研究区工业化、城镇化进程发展最快的时期,农业及城乡建设等人类活动不断扩大,导致人为水土流失加重。2015年土壤侵蚀强度最低,主要是由于2015年沂蒙山区降雨出现异常情况,侵蚀性降雨量大幅减少且空间分布较为均匀。2015年平均降雨侵蚀力为1 323.93 MJ·mm·hm-2·h-1·a-1,较2000、2005、2010和2018年平均降雨侵蚀力(1 805.02、2 474.22、2 519.69和2 466.41 MJ·mm·hm-2·h-1·a-1)分别下降26.7%、46.5%、47.5%和46.3%。采用模拟降雨法[22]计算2000、2005、2010与2018年在2015年降雨下的土壤侵蚀情况,得出沂蒙山区2000、2005、2010、2015与2018年的平均土壤侵蚀模数为3 435.29、3 837.27、3 541.73、3 493.23和3 434.30 t·km-2·a-1,计算2000、2005、2010与2015年在2018年降雨下的土壤侵蚀情况,得出沂蒙山区2000、2005、2010、2015与2018年的平均土壤侵蚀模数为6 306.29、6 344.00、6 189.78、6 115.75和5 669.03 t·km-2·a-1,可以看出降雨之间的年际差别是2015年土壤侵蚀强度大幅度下降的原因。近20 a来强烈及以上等级的侵蚀主要分布在西部和北部,集中在山亭区、沂南县一带。这除了与不可避免的人类活动不断扩大有关外,这部分地区长期以来开采矿产资源,形成大量废弃矿山,也在一定程度上加剧了该区域的土壤侵蚀。

4.2 驱动因素时空变化

不同时期各驱动因子对土壤侵蚀的影响程度显著不同。土地利用类型在近20 a对土壤侵蚀始终具有最高的贡献率,表明土地利用类型是影响沂蒙山区土壤侵蚀的最主要因素。研究区地貌复杂,以山地、丘陵为主,这样的地形地貌从宏观上控制地表过程的发生和发展,对土壤侵蚀也具有主导性作用。虽然坡度在近20 a对土壤侵蚀的贡献率不足20%,但坡度与其他因子交互作用可大大增强对土壤侵蚀的控制作用,由此也表明坡度是影响土壤侵蚀的主要因素之一。植被覆盖度受人类活动的影响强烈,1995—2000年间通过实施退耕还林政策,整体的植被覆盖度得到了一定程度的增加,导致了植被覆盖度空间分布趋于均匀的变化。受此影响,植被覆盖度对土壤侵蚀的影响程度在2000年减弱。年均降雨量对土壤侵蚀的影响在2000、2005和2015年较小,这可能与这3个时期降雨量的空间分布较为均匀有关。随着2010和2018年降雨量分布空间差异的增大,年均降雨量成为影响土壤侵蚀的重要因素之一。海拔对土壤侵蚀的影响程度在5个时期均未达到显著水平,其原因可能在于海拔差异小且空间离散化程度较低。

由于影响土壤侵蚀的因素较多,同时受数据收集的限制,研究仍存在一些不足:RUSLE模型因子计算过程中各个因子的数据精度不同,在栅格计算过程中会产生一定的误差。在地理探测器的应用中主要选择了部分典型的驱动因子,今后的研究中若加入地层岩性、人口等更多因子进行探讨,将更有利于沂蒙山区土壤侵蚀重点治理区的判别及因地制宜的进行治理。

5 结论

(1)沂蒙山区2000、2005、2010、2015与2018年的平均土壤侵蚀模数分别为5 129.01、6 661.76、6 344.00、3 493.23和5 669.03 t·km-2·a-1。2000—2018年轻度、中度、强烈、极强烈和剧烈侵蚀向低等级转移的面积占比分别为67.70%、60.00%、80.00%、78.57%和50.00%,土壤侵蚀强度向低等级转化趋势明显。从研究区的总体特征来看,土壤侵蚀先升高后降低,土壤侵蚀呈改善趋势。

(2)近20 a来沂蒙山区土壤侵蚀较为严重的区域主要分布在西部和北部,主要发生强烈及以上等级的侵蚀(≥5 000 t·km-2·a-1),强烈以下等级的侵蚀(<5 000 t·km-2·a-1)主要分布在西南部和中东部。随时间变化,土壤侵蚀呈现2000—2005年向北和西北方向扩散,2005—2010年向西北方向扩散,2010—2015年向中部扩散,2015—2018年向西和东北方向扩散的趋势。

(3)驱动因子在不同时期的解释能力不同,土地利用类型是影响沂蒙山区土壤侵蚀的主导因子,对土壤侵蚀的解释力在42.52%～66.84%之间。不同时期驱动因子的交互作用均会增强其对土壤侵蚀的解释能力,土地利用类型和地形地貌与其他因子交互后对侵蚀强度的贡献率明显提高。土地利用类型为草地、地形地貌为中起伏中山、坡度位于15°～20°区间以及植被覆盖度处于0.4～0.6区间的地区是水土防护治理的重点区域。