绥中县土壤侵蚀时空动态变化及其影响因素分析

张成哲

（营口市水利勘测建筑设计院，辽宁 营口 115000）

实践表明，土壤侵蚀是改变地貌景观的主要方式[1]。水土流失使土地退化、土质变差、农业减产以及山洪地质灾害加剧，给人类居住和生态环境安全构成潜在威胁。我国土壤侵蚀程度深、水土流失面积广，2019年动态监测资料显示，全国土壤侵蚀以风蚀和水蚀为主，水土流失面积达到271.08万km2，并以西北部最为严重[2]。目前，土壤侵蚀研究主要集中于区域、流域、坡面3个尺度上的定量分析，模型主要有EROSEM、WEPP等物理模型和USLE、RUSLE等经验模型，由于获取物理模型参数的难度较大，应用较多的是经验模型[3-6]。为了使得RUSLE模型本地化，许多学者考虑区域特点对其做了适当修正，并以刘保元的CLSE模型最为典型[7-9]。因具有所需数据易获取、模型简单和考虑影响因子全面等特点，RUSLE模型受到我国诸多学者的重点关注，并广泛应用于多个地区，如喀斯特地区、黄土高原地区、长江中下游和滇南山区等[10-12]。地理探测器是研究土壤侵蚀影响因素最常用的方法之一，可以识别地理现象的高风险区和揭示地理现象背后的驱动因子。

绥中县位于大兴安岭——太行山北北东向构造带东缘的交接部位，脆性特征非常明显，从辽宁省土壤侵蚀成果中提取的绥中县土壤侵蚀数据难以反映小区域实际特征。因此，本文以绥中县为例，基于2004年、2012年和2020年绥中县降雨和遥感影响数据，利用RUSLE模型揭示土壤侵蚀影响因子及其时空动态变化特征，并进一步识别出驱动力因素，为绥中县国土空间合理规划及其退耕还林、生态修复、土壤侵蚀防治等措施提供科学依据。

1 区域概况

绥中县隶属于葫芦岛市，介于东经119°34′－120°37′，北纬39°51′－40°32′之间，辖13个镇、11个乡，总面积2721.18km2。县内河流较多，在流域划分上属辽河流域和辽东湾西部沿渤海诸河水系，主要有王宝河、六股河、狗河、猫眼河、黑水河、金丝河、九江河、石河、黄家河秋皮河、小庄子河、塔子沟河、长滩河等。地貌特征呈明显的阶梯状，东南低、西北高，依次为沿海平原、丘陵、中低山，高程50m以下的平原区占34.3%，高程50～200m的丘陵占39.4%、高程超过200m的中低山区占26.3%，按其成因可以分为冲积海积平原、剥蚀堆积求和构造侵蚀中低山3种地貌类型。

绥中县属北温带大陆性气候区，其气候特征呈冬季寒冷干燥、夏季多雨炎热、秋季凉爽宜人、春季多风干旱，多年平均降水量617.4mm，其中汛期雨量471.1mm。绥中县土层疏松，地质灾害频发对道路、桥梁和耕地影响较大；矿产资源、采石、陡坡开垦使得自然灾害频发，土壤侵蚀范围扩大，自然环境变差，对当地经济发展产生极大阻碍。

2 研究方法

2.1 数据来源

2004年、2012年、2020年降雨数据和遥感数据来源于国家地球系统科学数据中心、Landsat TM/OIL影像(1～2月数据)；土壤有机碳以及粉砂含量等质地类型数据来源于HESD数据库；采用最大值合成法生成2004年、2012年和2020年1～12月NDVI数据，基础数据由MOD13Q1产品提供；应用ASTER GDEM V2全球数字高程数据提供分辨率30m的数字高程模型(DEM)；通过几何校正、辐射定标解译遥感影像以及野外验证等程序获取土地利用类型，随机从解译结果中选取50个样点到现场实地检验其准确度，经验证精准度达到85%，符合精度要求。其中，地理空间数据提供数字高程模型和Landsat遥感影像数据，统一坐标取WGS_1984_UTM_47N，数据栅格统一大小为30m。

2.2 研究方法

1）计算土壤侵蚀模数A。对于土壤侵蚀量和侵蚀模数考虑利用RUSLE模型进行计算，其表达式为：

式中：A代表土壤侵蚀模数[t/(hm2·a)]；R、K代表降雨侵蚀力因子[(MJ·mm)/[hm2·h·a]和土壤可蚀性因子[(t·h)/(MJ·mm)]；LS代表地形因子(无量纲)，即坡度与坡长；C、P代表植被覆盖因子和水保措施因子(无量纲)。

a）降雨侵蚀力因子R反映了降雨诱发土壤侵蚀情况的潜力，本文选用杨子生等[13]修订的方法计算R因子值，即：

式中：M代表年降雨量(mm)。

b）土壤可蚀性因子K反映了土壤对侵蚀的敏感程度，认为土壤质地越细腻或越粗糙则越不易被侵蚀，土壤质地属于中间水平时则越易被侵蚀，本文利用Williams[14]在EPIC模型中提出的公式计算K因子值，即：

式 中：SAN、SIL代 表0.050～2.000mm的 沙 粒 含量(%)和0.002～0.050mm的粉粒含量(%)，其中SN1=1-SAN/100；C、CLA代表土壤有机碳含量(%)和＜0.002mm黏粒含量(%)。依据有关研究资料，必须对计算得到的美制单位K值进行单位系数换算，以国际制单位K国替换美制单位K美，转换公式为K国=0.1317K美。

c）地形因子LS用于表征地形加速土壤侵蚀的程度，其数值等于坡长因子L、坡度因子S的乘积值，本文借鉴刘保元[15]和Mc Cool等[16]的经验公式计算地形因子值，具体表达式为：

式中：λ、θ代代表坡长(m)和坡度(°)；m代表坡度指数(无量纲)，数值等于坡度百分比。

d）植被覆盖因子C反映了土地管理、耕作方式、林草植被等对土壤侵蚀的影响作用，本文利用蔡崇法[17]提出的修正公式求解C因子值，即：

依据遥感影像统计NDVI(归一化植被指数)的最小值与最大值，然后利用下述公式计算植被覆盖度F，即：

e）水保措施因子P是反映不同用地类型是否采取水保措施的主要依据，其取值区间0～1，地表裸露无水保措施时因子取1，水保措施越完善越趋近于0。本文结合实地调查资料和区域微观坡面形态特征、宏观地貌类型等因素，对不同用地类型合理赋予相应的P因子值，如表1。

表1 不同用地类型的P因子值

2）地理探测器。这是一种能够有效规避人为主观性影响，有效识别多因子之间的关系分析地理现象空间异质性的统计学方法，无假设条件是其最大的特征，理论依据为空间方差原理[18-20]。

本文利用交互作用探测器、风险探测器、因子探测器3个模块定量分析可能影响土壤侵蚀的土地利用类型、植被覆盖度、坡度、海拔、年均降雨量5类因子。属性M对N的解释程度(取值0～1)因子解释力(q)值来衡量，数值越小则代表影响因子的解释力越强；各因子的显著性检验用p值来判断人，若p＜0.05则代表通过显著性检验。两个因子之间的作用以交互作用探测器来判别，通过比较多因子与单一因子的q值判定其影响程度；对影响土壤侵蚀的因子取值范围主要利用风险探测器来确定。依据相关文献资料，离散化处理以上5类因子，采用自然间断法将海拔和年均降雨量划分成9类，植被覆盖度和坡度划分成8类，土地利用类型划分成未利用土地、水域、建设用地、耕地、草地和林地。对研究区利用Arc GIS软件划分网格，设定采样点间距500m，将异常数据剔除后最终确定12160个点数。

3 结果与分析

3.1 土壤侵蚀时空变化

通过计算处理各因子值确定绥中县土壤侵蚀模数平均值10.64t/(km2·a)，年均侵蚀量2.8953×104t，将绥中县土壤侵蚀强度按照土壤侵蚀标准分级，统计计算各时期的土壤侵蚀面积和侵蚀模数。

表2 土壤侵蚀强度

结果表明，绥中县土壤侵蚀强度以微度和轻度为主，侵蚀程度较低且侵蚀强度分布较广，3个时期微度和轻度侵蚀面积占总面积90%以上，中度侵蚀占1%～3%，土壤侵蚀整体处于较轻状态。2004年、2012年、2020年研究区土壤侵蚀模数表现出逐渐减小趋势，依次为14.25、9.10和8.57t/(km2·a)，2004～2020年不同强度的侵蚀面积变化显著。2004年、2012年、2020年微度侵蚀面积从1648.86km2不断增大到1702.65km2、1908.71km2，轻度侵蚀从836.02km2不断减小到811.00km2、614.52km2，轻度、中度和剧烈侵蚀面积表现出明显减小趋势，高侵蚀等级不断转移到低侵蚀等级，强烈和极强烈侵蚀面积呈现出先增加后减小的趋势。总体而言，绥中县土壤侵蚀程度以微度和轻度为主，土壤侵蚀模数呈不断减小趋势，近15a来轻度、中度和剧烈侵蚀面积明显减小，土壤侵蚀程好转趋势，这与国家生态文明建设和退耕还林草工程的实施密切相关。

为了更加直观地揭示绥中县土壤侵蚀时空动态变化特征，将土壤侵蚀强度数据利用Arc GIS软件的空间叠加功能生成时空变化转移矩阵，见表3。

表3 土壤侵蚀强度转移矩阵

结果表明，2004-2012年研究区土壤侵蚀强度转化主要集中于微度、轻度和中度之间，从轻度转化为微度的侵蚀面积有208.08km2，从中度转化为轻度的侵蚀面积有25.78km2；从低强度转化为高强度的侵蚀面积有275.2km2，从高强度转化为低强度的侵蚀面积有379.10km2。2012-2020年和2004-2012年的侵蚀强度转化基本一致，从轻度转化为微度的侵蚀面积有312.73km2；从低强度转化为高强度的侵蚀面积有233.67km2，从高强度转化为低强度的侵蚀面积有443.62km2。2004-2012年强烈和极强烈侵蚀面积增加41.63%，上升幅度达到59%；统计分析强烈和极强烈中各地类所占面积比重发现这两个等级中的耕地是所有用地类型面积占比最高的地类。

叠加统计强烈和极强烈侵蚀等级分布图、坡度图和土地利用类型图，如图3所示。结果显示，在强烈和极强烈侵蚀整体面积变化中坡度＞25°的耕地占据着主导作用，作物类型、耕作方式、耕作水平和区域降水是影响坡度＞25°的耕地中侵蚀等级差异性的主要因素。总体而言，2004-2020年绥中县土壤侵蚀强度未发生变化的区域占74.68%，侵蚀强度变化集中于中度、轻度和微度等级。

图3 坡耕地面积与侵蚀等级面积关系图

3.2 土壤侵蚀影响因子

将各因子对土壤侵蚀空间分异特征的显著性检验(p值)和解释力(q值)利用因子探测器进行计算，如表4。结果表明，各因子的p值均为0.00，均通过0.05的显著性检验，即土壤侵蚀分布受各因子的影响显著。从低到高个因子的解释力排序为坡度＜海拔＜年均降雨量＜植被覆盖度＜土地利用类型。q值最高的是土地利用类型的62.27%，即土地利用类型是主导因子，对土壤侵蚀的贡献率最高；q值最低的是坡度的0.18%，对研究区土壤侵蚀影响最小，以上人为和自然因素共同决定了绥中县土壤侵蚀分布特征。其中，年均降水量和海拔的q值比土地利用类型低，表明绥中县土壤侵蚀受人类工程活动的影响最为显著。

表4 各因子p、q值

地理探测器中的交互探测器计算结果见表5，单个因子的解释力小于交互后的解释力。土地利用类型及其与其它因子的叠加具有较强的解释力，均排在第一；其中，土地利用类型和坡度因子叠加的q值最高，表明土地利用类型与不同坡度之间的侵蚀程度存在显著差异。植被覆盖度与土地利用类型交互的q值为0.6701，而单个因子——植被覆盖度的q值为0.0271，表明退耕还林草工程的实施具有重要作用，因此为有效抑制土壤侵蚀的发生可通过生态修复未利用土地和耕地，植树种草提高林草地面积和植被覆盖率，增强下渗能力及减少地表径流，逐渐改善区域水土流失现状。

表5 交互作用下各因子的q值

对土壤侵蚀高风险区域利用风险探测器的计算结果进行识别，结果显示绥中县高风险区位于植被覆盖度＜30%、年均降雨量超过600mm、坡度＞35°、海拔200m以上未利用地土地的区域。通过叠加植被覆盖度和土壤侵蚀分布图发现，对抑制土壤侵蚀植被覆盖度发挥着明显作用，植被覆盖度＜30%属于高风险地区。建设用地相较于其它影响因子q最高，也是土地利用类型中最易引起侵蚀的地类，其侵蚀强度最高达到510.85t/(km2·a)；研究区域内的未利用土地以裸地为主，未采取任何水保措施且植被覆盖度低、环境差。坡度超过35°时极易发生土壤侵蚀，坡度越大则产生侵蚀的概率越高，受自身产生的重力作用树木在坡度达到一定阈值时无法存活，树木存活率和植被覆盖度较低，故未利用土地或建设用地的固土能力差，降雨冲刷地面会引起侵蚀量的增多[21]。在植被覆盖度较低且雨量比较丰富的区域，降雨径流冲刷坡面土壤比较强烈，极易产生土壤侵蚀。

绥中县土壤侵蚀影响因子主要是土地利用类型因子，不同用地类型的侵蚀强度具有明显差异，对比分析各用地类型土壤侵蚀面积和侵蚀模数，如表6所示。

表6 各用地类型面积及土壤侵蚀模数

从表6可以看出，从低到高不同用地类型的土壤侵蚀强度排序为水域＜草地＜林地＜耕地(含园地)＜未利用地＜建设用地，未利用土地面积0.12万hm2，主要分布于六股河沿岸，并且均未采取水保措施。草地和林地侵蚀强度较小，植被覆盖度较高，耕地(含园地)面积11.24万hm2，侵蚀模数为88.37t/(km2·a)，主要分布于人类活动密集区。根据绥中县土地利用类型分布特征、面积占比和地理探测器分析结果，坡耕地应作为研究区土壤侵蚀重点治理区，为抑制土壤侵蚀现状必须继续实施退耕还林草等工程措施。

4 结 论

1）2004-2020年绥中县土壤侵蚀强度以微度和轻度为主，微度和轻度侵蚀占总面积90%以上，中度侵蚀占1%～3%。研究区土壤侵蚀模数总体呈逐渐减小趋势，依次为14.25、9.10和8.57t/(km2·a)，近15a来轻度、中度和剧烈侵蚀面积明显减小，这与国家生态文明建设和退耕还林草工程的实施密切相关。

2）六股河两岸的土壤侵蚀比较严重，研究区土壤侵蚀强度转化主要集中于微度、轻度和中度之间，统计分析强烈和极强烈中各地类所占面积比重发现这两个等级中的耕地是所有用地类型面积占比最高的地类。

3）从低到高不同用地类型的土壤侵蚀强度排序为水域＜草地＜林地＜耕地(含园地)＜未利用地＜建设用地。建设用地相较于其它影响因子q最高，也是土地利用类型中最易引起侵蚀的地类，绥中县高风险区位于植被覆盖度＜30%、年均降雨量超过600mm、坡度＞35°、海拔200m以上未利用地土地的区域。根据绥中县土地利用类型分布特征、面积占比和地理探测器分析结果，认为坡耕地应作为研究区土壤侵蚀重点治理区。