土壤侵蚀时空分布规律及侵蚀热点分析——以陕西省王益区为例

陈四宾，霍艾迪，2，陈思名，赵志欣，陈 建

(1.长安大学 水利与环境学院，西安 710054； 2.长安大学 旱区地下水文与生态效应教育部重点实验室，西安 710054)

1 研究背景

土壤侵蚀是全球性的环境问题之一[1-2]，它会减少表层土壤的养分含量，破坏土壤肥力造成土地退化等问题[3-6]。严重的土壤侵蚀还会影响区域农业经济发展[7-8]。黄土台塬是宝贵的宜农、宜耕、宜居的土地资源，其塬面地势平坦，土地肥沃，易于耕作。然而受土壤侵蚀问题的困扰，黄土塬周围沟道丛生，沟头不断前进导致塬面萎缩[9]，严重影响当地的经济发展。自20世纪50年代以来，黄河流域开展的小流域治理、植树造林、退耕还林等治理工程成效显著，但是在治理程度方面仍然存在区域差异[10]。定量分析黄土高塬沟壑区的土壤侵蚀状况及其时空变化特征，对当地土壤侵蚀防治和生态建设具有重要意义。

一般而言，受研究区域客观条件限制，采用实地观测的土壤侵蚀评价方法往往存在较大困难[11]，实际应用过程中局限性较大，因此常采用土壤侵蚀预报模型进行土壤侵蚀的定量评价。从20世纪50年代开始，我国学者逐步结合中国土壤侵蚀的特点开展土壤侵蚀模型的研究，有关黄土高原地区侵蚀模型的研究也得到了快速发展[12]。刘宝元等[13]基于通用土壤流失方程(Universal Soil Loss Equation，USLE)利用黄土丘陵沟壑区实测资料建立了中国土壤流失预报方程(Chinese Soil Loss Equation，CSLE)；Zheng等[14]根据黄土高原不同小流域的实测资料建立了预测次洪产沙量的比例函数，流域模拟产沙结果比较准确。物理模型方面，王建勋等[15]研究发现水蚀预报模型(Water Erosion Prediction Project，WEPP)坡面版在黄土丘陵沟壑区具有良好的适用性；此外，Huo等[16]基于水土评价工具(Soil and Water Assessment Tool，SWAT)模型模拟了不同固沟保塬工程下陇东黄土高原典型流域的径流泥沙变化。随着现代技术的快速发展，采用遥感和地理信息系统(GIS)技术进行土壤侵蚀的评价分析，可以快速、准确地获取具有一定认可程度的研究结果[17]。

在生态保护与城市建设发展的双重影响背景下，如何科学防治土壤侵蚀，协调好人地关系成为当前城市建设和生态发展工作中亟须解决的问题。本文以典型的黄土高塬沟壑区——铜川市王益区为研究对象，分析了其2000—2019年的土地利用和植被覆盖度变化，依据水利部颁布的《土壤侵蚀分类分级标准》(SL 190—2007)[18](以下简称“标准”)估算分析王益区土壤侵蚀强度变化及其时空分布特征。基于土壤侵蚀强度风险指数(Soil Erosion Dangerous Index，SEDI)和探索性空间数据分析(Exploratory Spatial Data Analysis-GIS，ESDA-GIS)方法进一步分析了区域内的土壤侵蚀强度风险状况和土壤侵蚀冷热点分布(土壤侵蚀高强度集聚区和低强度集聚区)情况，以期为王益区水土流失重点防治区的划分、固沟保塬工程实施和生态环境保护提供一定的科学依据。

2 数据来源与研究方法

2.1 研究区概况

王益区位于陕西省铜川市中部，地处渭北黄土高原中部南缘地带(图1)，地理坐标为108°57′E—109°09′E，34°57′N—35°12′N，全区面积162.2 km2。经过多年的地壳运动，现整体地势西北高、东南低，区域内梁峁川塬山丘交错相间，沟壑纵横。西北部是以黄土梁为主的丘陵山地，中部的黄土川塬地带面积大，黄土厚度深，东南部为黄土丘陵，海拔1 000～1 400 m，多沟壑残塬，地形复杂。该地区属于暖温带大陆性季风气候，四季分明，年平均气温8.7～12.1 ℃，年均降水量589.2 mm，年蒸发量837.1 mm，年均无霜期182 d。

王益区地处半干旱地带，属于典型的黄土高塬沟壑区，土壤侵蚀以水力侵蚀[19]为主，黄土塬塬面受土壤侵蚀影响严重，塬面萎缩，农耕地肥力下降。区域内土壤侵蚀风险可能进一步加剧。

2.2 数据来源

本文采用王益区数字高程模型(Digital Elevation Models，DEM)数据、土地利用数据以及植被覆盖数据。DEM数据来自地理空间数据云网站(www.gscloud.cn)，分辨率为30 m。基于地理空间数据云网站(www.gscloud.cn)下载2000年的Landsat ETM+影像和2019年的Landsat OLI影像，提取王益区植被覆盖度变化，两期影像成像时间分别为9月份和8月份，轨道号127/36，影像成像时间接近且质量较好。在1.96 m分辨率遥感影像支持下，结合Landsat遥感影像的光谱特征，采用人机交互解译的方法，判别地物的几何形状，颜色纹理特征和空间分布情况，目视勾绘出地类边界，解译得到所需土地利用数据。

2.3 指标因子提取

本文基于遥感影像，采用像元二分模型[20]进行植被覆盖度的估算，李苗苗等[21]基于归一化植被指数(NDVI)构建的像元二分模型在植被覆盖度提取方面可操作性强，精度较高。相比较于直接利用NDVI反映植被覆盖度，该方法在削减植被类型、土壤和大气的影响方面更具有优势[22]。其计算公式为

(1)

式中：FVC表示植被覆盖度(%)；NDVIsoil表示裸土或无植被覆盖区域的NDVI值，其变化范围在-0.1～0.2之间[21，23]；NDVIveg表示完全由植被覆盖区域的NDVI值。参照NDVI的频率统计表，分别取累计频率5%和95%的NDVI值作为NDVIsoil和NDVIveg的取值。

基于DEM数据，利用ArcGIS软件进行研究区坡度的提取，并根据相应标准进行重分类，由于研究区坡度>35°的区域面积较小，故将(25°，35°]和>35°区间统一合并为>25°一个区间。依据相应标准划分为(0，5°]、(5°，8°]、(8°，15°]、(15°，25°]、>25°五个坡度区间。

2.4 土壤侵蚀量估算

采用水利部颁布的《土壤侵蚀分类分级标准》(SL 190—2007)，可以基于植被覆盖度、坡度和土地利用较为快速准确地获取研究区土壤侵蚀强度分布状况[24]。

王益区土壤侵蚀以水力侵蚀为主[19]，根据面蚀分级指标(表1)，将土壤侵蚀强度划分为“微度”“轻度”“中度”“强烈”“极强烈”5类，侵蚀模数分别为<1 000、[1 000，2 500)、[2 500，5 000)、[5 000，8 000)、[8 000，15 000) t/(km2·a)。对于王益区城乡居民用地和工矿用地，参考城市化土壤侵蚀等级分级评分标准中以斜坡高差、植被覆盖度等为评价项目进行分级的方法[25-27]，以及侵蚀强度调查中，对非耕地和未利用地室内测算采用坡度和地面物质为主要指标的分类方法[28]。对大面积地面硬化条件较好的建成城区，土壤侵蚀强度以微度侵蚀为主[29]；其余工矿裸露用地主要考虑坡度、下垫面条件及区域植被覆盖情况进行侵蚀强度划分。此外，基于2019年的高分辨率影像及野外调查发现，王益区梯田面积较大，在使用“标准”进行土壤侵蚀强度划分时，为区别于坡耕地，参考通用土壤流失方程模型中P因子(采取水土保持措施与不采取任何措施的地块土壤流失量的比率)的取值[30]，判断不同坡度下梯田水土保持作用的折减效果，并对区域内梯田的土壤侵蚀强度进行划分。而对于坡度>24°的梯田，水土保持效果不明显[31]，仍将其划分为极强烈侵蚀。

表1 土壤侵蚀分类分级标准

王益区属于西北黄土高原区，容许土壤流失量为1 000 t/(km2·a)，即该地区允许存在微度侵蚀。研究过程中，取各个土壤侵蚀模数区间的中值作为平均侵蚀模数进行王益区土壤侵蚀量的估算。

2.5 土壤侵蚀强度风险分析

为进一步分析王益区土壤侵蚀风险状况，采取SEDI分析研究区内的土壤侵蚀严重程度。SEDI可以综合反映出不同土壤侵蚀强度对区域内侵蚀状况的贡献程度大小，评判土壤侵蚀严重情况，计算公式[32-33]为

(2)

式中；Fi表示第i级侵蚀强度的量化分级值；Ai表示评价区域第i级侵蚀强度的面积；A表示评价区域的总面积；n为土壤侵蚀强度级别数目；100为扩大系数。微度侵蚀、轻度侵蚀、中度侵蚀、强烈侵蚀、极强烈侵蚀对应的量化分级值分别为0、2、4、6、8。相应的土壤侵蚀强度风险等级分为轻微(SEDI<50)、一般(50200)。

2.6 空间数据分析

全局空间自相关分析是利用全局Moran’ I指数统计不同土壤侵蚀强度之间的平均关联程度和显著性，描述土壤侵蚀强度空间上整体的关联程度和差异性[34-35]，计算公式为

(3)

(4)

式中：Z为正态检验统计量；E(I)和VAR(I)分别为期望值和方差，若其通过5%的显著性水平检验且Z>1.96，则说明具有显著的空间自相关。

局部空间自相关可以反映区域内不同土壤侵蚀强度之间的空间相关程度，弥补Moran’ I不能反映局部相关类型和聚集区的不足[34-35]，计算公式为

(5)

式中Ii表示空间关联局域指标，取值在-1～1之间，其指标值判定与Moran’ I指数一致。

(6)

(7)

3 结果与分析

3.1 王益区土地利用时空变化特征

3.1.1 土地利用时间变化特征

王益区现有土地利用类型中，耕地(坡耕地+梯田)面积最大(图2)，2000年和2019年分别为107.12 km2和79.85 km2，其中2019年耕地面积中梯田占比约47.63%，为38.03 km2。研究期间，水体面积占比极小且变化不明显，除坡耕地外，其他土地利用类型面积均不同幅度的增加。坡耕地作为其他土地利用类型转入量的最大来源，2019年其面积占比减少至25.78%。与此同时，林草地和城乡居民用地呈现较大幅度增长，前者增幅最为显著(39%)，高达17.05 km2。城乡居民用地增加了10.15 km2。耕地、林草地、城乡居民用地的变化，侧面反映出研究区内退耕还林的成效以及城市化的进程。

研究期间土地利用变化主要集中在林草地、耕地和城乡居民用地、梯田的变化上。因此土地利用转移的空间分布主要考虑6种土地利用变化类型(图2)。

3.1.2 土地利用空间变化特征

王益区沟壑纵横，地势起伏变化大，(8°，25°]坡度区间面积为94.24 km2，占全区面积的58.10%。空间分布上，随着坡度增加，各土地利用类型的面积总体呈现出先增后减的趋势(图3)，不同时期林草地多集中在坡度>8°区域内，在(15°，25°]坡度区间内达到峰值，该坡度区间内的林草地面积分别占林草地总面积的38.11%和38.54%。2000年(0°，25°]坡度区间内的优势土地利用类型均为耕地，到2019年，(15°，25°]区间内的优势土地利用类型则转化为林草地。城乡居民工矿用地的总体面积占比较小且主要均匀分布在(0°，15°]坡度区间内。

总体而言，各土地利用类型虽呈现出不同幅度的变化，但是除(15°，25°]坡度区间外，其他坡度区间内的优势土地利用类型未发生变化。

3.2 林草地植被覆盖度变化

根据相应标准将植被覆盖度等级划分为5类，分别为(0，30%]，(30%，45%]，(45%，60%]，(60%，75%]，>75%。研究期间王益区植被覆盖度变化如图4所示。

植被覆盖度FVC>60%的林草地面积大幅增长，增加了18.93 km2，约占2000年林草地总面积的42.94%；其中植被覆盖度>75%的林草地面积增幅最大，增加了约0.7倍。此外，植被覆盖度<60%的林草地面积小幅下降，其面积占比由25.5%减少到15.4%。林草地面积增加的同时伴随着高植被覆盖度面积的增长，研究区内的植被覆盖整体上明显好转。

结合研究区DEM(图1)，土地利用及其空间分布(图2、图3)，王益区林草地多分布于>8°的高坡度区域，高植被覆盖度比例的增长，说明王益区封山禁牧取得了一定的成效，对退耕还林成果的巩固较好。

3.3 2000—2019年王益区土壤侵蚀时空变化特征

3.3.1 土壤侵蚀时间变化特征

根据相应标准，2000—2019年王益区土壤侵蚀面积变化如图5所示：2000—2019年王益区土壤侵蚀面积由115.04 km2减少至90.58km2，侵蚀比例由70.92%减少至55.84%。侵蚀强度结构由中度侵蚀>强烈侵蚀>轻度侵蚀>极强烈侵蚀变化为轻度侵蚀>中度侵蚀>强烈侵蚀>极强烈侵蚀。极强烈侵蚀面积减少虽小(5.27 km2)，但减幅较大(48.75%)，强烈侵蚀面积减少最大(22.45 km2)，减幅也最大(74.26%)。一方面由于陡坡耕地的减少，另一方面也得益于梯田的修建以及高植被覆盖度林草地面积的增加。

研究期内，王益区土壤侵蚀量减少幅度明显(39.36%)，由53.81×104t减少至32.63×104t。年均土壤侵蚀模数由0.332×104t/(km2·a)变化为0.201×104t/(km2·a)。与此同时，王益区土壤侵蚀强度风险指数(SEDI)由2000年的310减少至198，区域内的土壤侵蚀状况整体呈好转趋势。

3.3.2 土壤侵蚀空间变化特征

通过全局空间自相关(Global Moran’ I)分析不同土壤侵蚀强度之间的平均关联程度和显著性，结果显示2000年和2019年王益区的全局Moran’ I指数分别为0.53，0.51，相应的Z得分分别为20.78、22.43，且通过了5%的显著性水平检验，说明王益区土壤侵蚀存在明显的空间正相关性。

基于ArcGIS的聚类和异常值分析(Anselin Local Moran I)计算得到王益区土壤侵蚀强度高值HH(High-High)、低值LL(Low-Low)的空间聚类和空间异常值HL(High-Low)，LH(Low-High)的分布图(LISA图)。其能够直观反映不同土壤侵蚀强度的空间聚集及相异性特征，弥补全局空间自相关在空间聚集位置表达上的不足。结果(图6)显示，2000年和2019年的土壤侵蚀强度空间聚集状态明显，2000年的HH聚集区面积为32.85 km2，占全部聚集区面积的54.32%，在王益区内呈连片集中分布。到2019年LL聚集区面积增幅明显(49.60%)，其增幅与研究期间微度侵蚀面积增幅(51.86%)相近。与此同时HH聚集区面积减少了8.4 km2，减幅高达25.57%。区域内LL聚集程度的增强和HH聚集程度的降低与全区侵蚀强度变化趋势相符。

总体上，研究期间土壤侵蚀强度的空间演化较为明显，LL聚集区面积的增加与HH聚集区面积减少，也反映土壤侵蚀状况得到了有效缓解。

叠加分析研究区侵蚀热点分布与坡度分级。由表2可知，2000年和2019年土壤侵蚀热点主要集中在(15°，25°]区间内，分别约占热点区总面积的39.95%和40.34%；其次为(8°，15°]和>25°区域。此外，(8°，15°]和(15°，25°]区域内的热点区面积减少较明显且极为接近。这与前述林草地主要分布在>8°区域且(15°，25°]坡度区域内的优势土地利用类型转化为林草地关系密切。

表2 不同坡度区间土壤侵蚀热点面积分布

研究期间侵蚀热点区面积大幅减少，土壤侵蚀强度风险指数降低，进一步说明王益区土壤侵蚀状况明显好转。究其原因，一方面由于王益区封山禁牧政策的实施；另一方面，是由于区域塬面涝池的修建以及部分坡改梯工程的实施，对区域内土壤侵蚀防治起到一定的促进作用。

3.3 讨 论

(1)采用《土壤侵蚀分级分类标准》(SL 190—2007)进行土壤侵蚀量的估算，可以基于土地利用、坡度、植被覆盖度快速地获取具有一定认可程度的区域土壤侵蚀强度变化，在土壤侵蚀相关研究[4，30，36-37]中应用较多。探索性空间数据分析方法在土壤侵蚀空间聚集的分布演化方面更具优势，但不能反映土壤侵蚀驱动因素，后续研究中应综合考虑降雨、坡长、水土保持措施等多个影响因子在内的土壤侵蚀变化情况，进一步分析各个影响因子与研究区域土壤侵蚀的相关性，或结合地理探测器等方法进一步揭示研究区土壤侵蚀可能的驱动因素。

(2)已有研究发现[38]，黄土塬区耕地面积减少主要由于建设用地的增加以及其他诸多因素的影响；这与研究期间王益区土地利用变化趋势相符，此外还得益于退耕还林工作取得的进展。本文估算的2000年、2019年王益区侵蚀模数分别为0.332×104t/(km2·a)和0.201×104t/(km2·a)，与郭力宇等[4]估算的渭北高原沟壑区年平均土壤侵蚀模数在0.3×104t/(km2·a)上下浮动较为接近。王益区沟壑面积大，沟头延伸，沟岸扩张问题严重，一方面可以进一步推进塬面保护工程，沟坡排水工程的实施，栽植水保林和经济林，防治沟道侵蚀、沟头前进；另一方面，结合当地封山禁牧政策，保护好现有耕地及林草地，避免出现过度开垦和滥牧行为，完善水土保持治理体系，实现蓄水保土，固沟保塬和区域生态良性循环。

(3)针对土壤侵蚀较为严重的侵蚀沟，积极开展区域内小尺度的土壤侵蚀模拟研究，一方面可评判固沟保塬等水土保持措施的效益，为其提供科学合理的改善建议；另一方面推进沟头前进主要致灾因子的研究，制订更加科学合理的土壤侵蚀防治体系。

4 结 论

本文基于植被覆盖度及土地利用等研究了王益区土壤侵蚀变化，结论表明：

(1)王益区耕地面积大幅减少，耕地面积减少了27.75 km2，林草地和城乡居民用地增长明显，林草地面积增长最为显著(17.06 km2)。空间上，林草地面积大幅增加，取代耕地成为(15°，25°]区域内的优势土地利用类型。

(2)王益区植被覆盖度明显好转。植被覆盖度>75%的林草地面积增加了0.7倍，>60%植被覆盖度的林草地面积增幅高达42.94%，封山禁牧的实施取得了一定的效果，对早期退耕还林工作成果保护较好。

(3)王益区土壤侵蚀状况明显好转。侵蚀强度结构由中度侵蚀和强烈侵蚀为主转化为中度侵蚀和轻度侵蚀为主，2000年和2019年王益区土壤侵蚀面积分别为115.04 km2和90.58 km2，占全区总面积的70.92%和55.84%。SEDI由310减少到198，土壤侵蚀量减少了21.18×104t，年均土壤侵蚀模数由0.332×104t/(km2·a)变为0.201×104t/(km2·a)。

(4)王益区土壤侵蚀强度分布呈明显的空间聚集效应。全局Moran’ I指数由0.53变为0.51，侵蚀热点由大面积连续分布转化为局部集中的零散分布，其空间变化与陡坡林草地变化关系密切，这也表明退耕还林工作取得了很好的效果。