2000-2014年徐州地区土壤侵蚀的变化趋势

王国斌, 陈报章, 孙少波, 李彦鹏, 梁 森, 渠俊峰

(1.中国矿业大学 环境与测绘学院, 江苏 徐州 221008; 2.天津大学表层地球系统科学研究院， 天津 300072; 3.徐州市水利科学研究所, 江苏 徐州 221003)

土壤侵蚀是指在水力、风力、重力及冻融等自然营力以及人类社会生产和建设活动下，水土资源和土地生产力的破坏和损失[1]。土壤侵蚀会消耗土壤肥力，改变土壤结构，降低有效的生根深度，扰乱所有自然过程的基础[2]，水土流失还会带走大量泥沙，不仅污染水环境，更导致河道淤堵，加大了洪涝灾害的风险。因此，在中国经济发展正处于转型期的背景下，掌握水土流失的变化规律对国计民生具有重要的现实意义。

随着3S技术的迅猛发展，基于遥感影像资料的定量计算土壤侵蚀程度的数学模型法得到了广泛的关注，《全国水土流失规划》明确指出在水土流失监测中应注重采用遥感和地理信息系统等技术结合地面观测和抽样调查方法对水土流失重点预防区和治理区进行监测，综合评价区域水土流失强度和分布状况、治理措施动态变化[3]。Wischmeier等[4]提出的USLE是全球常用的模型之一(包括热带地区)，用于估算水蚀。该模型形式上比一般的经验模型和物理模型更简单，公式仅由6个决定因子构成，许多学者对因子的确定方法进行了研究，以满足某一特定环境的需求。在降雨侵蚀因子(R)的经典算法中，需要高频的降雨资料，而在实际中很难获取如此充分的降雨数据， Beguería等[5]针对该问题，提出了基于中频数据或低频数据的最佳回归方法计算降雨侵蚀力，该方法被应用于西班牙地区，取得了不错的效果。Wu等[6]在实验室条件下，模拟了不同土壤类型在不同降雨强度下的侵蚀状况，得到了相同土壤不同降雨强度和不同土壤相同降雨强度的可侵蚀性规律，为以后模型的应用提供了宝贵的经验。Zhang等[7]开发了一种改进的方法(LS-TOOLMTFD)来评估LS因子，并在复杂变化斜率的流域集水区中得到了应用和验证，取得了较高的因子计算精度。Bambang[8]以印尼Marawi流域为试验区，将3种不同的植被覆盖因子(C)模型应用到侵蚀方程中，用统计学手段比较基于3种不同方法得到得侵蚀模数结果，从而分析了不同C模型的优劣，这在侵蚀计算本地化方面，是一个很有启发性的案例。Panos等[9]依据地球观测数据以及文献指导，估算欧洲范围内的P因子值，并评估农业保护措施(如石墙和草地边缘)对减少土壤流失的影响，对欧洲的耕地水土保持以及共同农业发展具有重大意义。迄今，USLE模型已经被广泛应用在中国各地区的侵蚀风险评估，周夏飞等[10]利用USLE模型，测算了江苏省2001—2015年土壤侵蚀模数，并进行了时空分析。张超等[11]基于USLE模型估算了渝东北生态涵养区的土壤侵蚀量，并在此基础上分析了不同降雨强度、坡度、土壤类型下的土壤强度分布特征。

徐州市位于江苏省西北部，华北平原的东南部，地跨东经116°22′—118°40′，北纬33°43′′—34°58′，东西长约210 km，南北宽约140 km，土地总面积11 259 km2，占江苏省总面积的11%。由于地处沂、沭、泗诸水的下游，境内河流纵横交错，湖沼、水库星罗棋布，加之土壤质地松散、季节性降水丰沛、人为扰动地表面广、量大等因素，使得辖区易发生土壤水蚀，近些年来，徐州市坚持不懈地抓好治水工作，以预防监督为主，防治并重，有效抑制了水土流失的发生和发展，但水土流失综合治理任务依然艰巨，其中水土保持新技术的引进、推广和应用力度不足是较为突出的问题[12]。迄今，徐州地区的土壤侵蚀时空分布特征鲜有研究，本文拟应用USLE模型计算得到土壤侵蚀模数，并以此数据为基础，以统计分析为手段旨在从时间尺度和空间尺度分析多年来徐州地区土壤侵蚀的变化规律。

1 数据与方法

1.1 数据来源

(1) 降雨数据。本文采用的降雨数据由徐州市水文局提供，监测站点包括双沟、新店、睢宁、高流、邳城闸、贾汪、运河、宋楼、新楼、五段、安国、徐州、正集。

(2) 土壤数据。土壤数据由寒区旱区科学数据中心网站(http:∥westdc.westgis.ac.cn/)下载获取，其中国境内的数据源为南京土壤所提供的中国第二次土壤调查1∶100万数据。

(3) 地形数据。地形数据为ASTGTM2 DEM数据，空间分辨率为30 m。

(4) 植被覆盖数据。来源于比利时VITO网站(http:∥www.vito-eodata.be/PDF/portal/Application.html)的SPOT_VEGETATION数据集中的NDVI数据，空间分辨率为1 km。辅助的遥感数据为Landsat TM数据。

(5) 土地利用类型数据。从中国科学院资源环境科学数据中心(http:∥www.resdc.cn/)获得。

1.2 土壤侵蚀量计算

通用土壤流失方程(USLE)形式为：

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

(1)

式中：A——土壤年流失量(t/hm2·a)，R——降雨侵蚀因子〔MJ·mm/(hm2·h·a)〕，K——土壤可蚀性因子〔(t·hm2·h)/(MJ·mm·hm2)〕，S和L——地形坡度和坡长因子，C和P——地表植被覆盖度因子和水土保持措施因子。

该模型考虑的各个因子，基本包容了土壤侵蚀力与抵抗力这一矛盾统一体的宏观轮廓，为土壤侵蚀量估算的因子选择提供了成熟的基本骨架[13]。方程中各因子值的确定均要求对监测区的相关地理要素进行详尽的分析,所以应用USLE的关键在于对各相关因子进行科学的确定[14]，本研究中各侵蚀因子的确定方法如下。

(1) 降雨侵蚀因子R的确定。由于R因子的经典计算方法以次降雨量为基础资料[15]，但数据采集难度较大且无法保证其连续性，故在本研究中采用了国内广泛应用的章文波等[16]提出的基于日降雨量的中国降雨侵蚀力模型。

(2) 土壤可侵蚀因子K的确定。采用Sharply等[17]发展的EPIC模型计算K值，并根据张科利等[18]的研究，将结果修订，使其更符合中国土壤特性。

(3) 地形因子LS的确定。本研究中坡度因子S和坡长因子L据区域数字高程模型DEM进行估计[19]。

(4) 植被覆盖因子C的确定。该因子是植被防止土壤侵蚀程度的定量指标，C的值取决于控制地表侵蚀的植被覆盖的规模、表面积，以及植物根部、表面粗糙度、含水量的状态[20]。本研究选用了蔡崇法等[21]提出的基于归一化植被指数(NDVI)的C值估算公式。

(5) 水土保持因子P的确定。P值的确定结合了地形坡度数据以及土地利用类型数据，具体参照了方广玲等[22]提出的选取标准。

1.3 土壤侵蚀时空变化分析

本文中计算年最大、平均侵蚀模数变化率以及整体变化分布，应用的数学手段为一元线性回归。回归分析，尤其是一元线性回归分析，做为统计学中的一个重要分支，以其成熟的原理，高可信度的预测结果，运用的简洁性等优点，在工程技术经济领域有着十分广泛的应用[23]。具体表达式为：

(2)

式中：变量i——1～15的年序号；Ei——第i年的侵蚀模数值。变化趋势分布图则反映了在15 a的时间序列中，徐州土壤侵蚀的变化趋势，某像点的像元值是这个像点15 a间的侵蚀模数用一元线性回归模拟出的变化直线的斜率。

2 结果与分析

2.1 土壤侵蚀强度分布

《全国水土保持区划》(2015—2030)采用三级分区体系，一级区为总体格局区，确定全国水土保持工作战略部署与水土流失防治方略，是水土资源保护、开发和合理利用的总体格局；二级区为区域协调区，协调跨流域、跨省区的重大区域性规划目标、任务及重点；三级区为基本功能区，确定水土流失防治途径及技术体系，作为重点项目布局与规划的基础[3]。徐州市位于北方土石山区(一级区)、华北平原区(二级区)、黄泛平原防沙农田防护区和淮北平原岗地农田防护保土区(三级区)[12]。本研究以全国水土保持区划为指导，按照三级区划将徐州区县划分为黄泛平原和丘陵地区(图1)。

图1 徐州市县级行政区划及其所属水土保持三级区

本研究利用上述的USLE模型计算出徐州土壤侵蚀模数影像，裁剪得到徐州黄泛平原及丘陵地区的土壤侵蚀数据。根据水利部颁发的《土壤侵蚀分类分级标准》(SL190-2007)，将徐州地区土壤侵蚀强度分为微度、轻度、中度、强烈和极强烈5个不同等级[24]，具体划分标准如表1所示。

表1 北方土石山区土壤侵蚀强度分级标准

黄泛平原包括丰县、沛县、睢宁、新沂和邳州的部分区域，由计算出的侵蚀模数分布得到，黄泛平原绝大多数地区的侵蚀程度为微量侵蚀，其余地区以轻度侵蚀为主，夹杂着小区域的其他程度的侵蚀，按照行政区域划分，轻度及以上侵蚀等级的区域主要为新沂市的棋盘镇、马陵山镇，邳州市的铁富镇、邳城镇，睢宁县的姚集镇、王集镇、岚山镇以及丰县的华山镇、宋楼镇、大沙河镇等。

丘陵地区的绝大多数区域以微度侵蚀为主，其余地区主要为轻度侵蚀，呈条状或带状分布，小部分区域为中度及以上侵蚀。从行政区划来看，轻度侵蚀主要发生在贾汪区的贾汪镇、利国镇、江庄镇、柳泉镇、青山泉镇、茅村镇，铜山区的大彭镇、汉王镇、铜山镇、三堡镇、棠张镇、徐庄镇、张集镇、伊庄镇，邳州市的八义集镇、占城镇，鼓楼区的庙山镇，以及泉山区、云龙区、鼓楼区的部分地区，中度侵蚀、强度侵蚀、极强度侵蚀、剧烈侵蚀范围逐级递减，分布在轻度侵蚀范围之中。比较丘陵地区和平原地区的侵蚀状况可知，丘陵地区各等级(除微度)侵蚀面积明显比黄泛平原大，并且侵蚀等级跨度更大。坡度是坡面土壤侵蚀中影响最大的因素，丘陵地区坡度更大，降雨冲刷更加剧烈，易造成更大程度的水土流失[25]。

2.2 年最大与平均侵蚀模数年际变化

逐年统计2000—2014年土壤侵蚀模数最大值、平均值，根据统计数据，分别绘制黄泛平原和丘陵地区最大侵蚀模数和平均侵蚀模数变化特征如图2—3所示。

图2 黄泛平原、丘陵地区年最大侵蚀模数变化

图3 黄泛平原、丘陵地区年平均侵蚀模数变化

由2—图3可知，从2000—2014年，黄泛平原和丘陵地区年平均侵蚀模数和年最大侵蚀模数的变化都处于下降的趋势，侧面说明多年来黄泛平原和丘陵地区的土壤侵蚀程度在逐渐减弱，并且丘陵地区年最大侵蚀模数变化的回归线通过了显著性水平为0.05的假设检验。

综合两组结果可以发现，丘陵地区侵蚀模数下降程度相较黄泛平原更大、更明显。从徐州市治理水土流失的历史资料可以发现，近些年来徐州有关部门修建了一大批蓄水保土工程，显著降低了徐州市的水土流失总量，尤其在主要的水土流失区——丘陵地区，水土流失治理率达到57.6%[3]，土壤侵蚀强度显著降低。

2.3 侵蚀面积变化

将黄泛平原和丘陵地区对应等级的侵蚀面积相加，统计多年来徐州市各等级土壤侵蚀面积变化过程，由于徐州属于北方土石山区，而北方土石山区的容许流失强度为200 t/km2，即微度侵蚀[19]。所以在统计面积时没有将微度侵蚀考虑在内，结果如图4所示。

图4 2000-2014年徐州市各等级侵蚀面积变化

可以发现除微度侵蚀外徐州市轻度侵蚀占主要比重，与其他侵蚀类型差距明显，各侵蚀等级面积呈波浪式变化，总体呈下降态势，分别将徐州黄泛平原及丘陵地区各年不同等级的侵蚀面积以所占总侵蚀面积比例的形式统计，结果如图5所示。

图5 2000-2014年黄泛平原和丘陵地区各等级侵蚀面积所占比例/%

由图5可以发现，总体来说，徐州黄泛平原轻度侵蚀面积占总侵蚀面积的80%以上，个别年份有所上升或下降，各等级侵蚀面积逐级递减，从多年来看，轻度以上等级的侵蚀比例在逐渐缩小。多年来丘陵地区各等级侵蚀所占比重的变化规律相较于黄泛平原更加明显，1 a中轻度侵蚀占总侵蚀面积的70%左右，自2007年之后，轻度侵蚀所占比重呈阶梯状上升趋势，轻度以上侵蚀比重呈阶梯状下降趋势，2009年以后极强度侵蚀控制在5%左右，剧烈侵蚀基本得到了遏制。

2.4 像素尺度侵蚀模数年际变化

对徐州市15 a间侵蚀模数影像逐像素计算回归系数，若回归系数大于0，则代表侵蚀模数呈上升趋势，即侵蚀程度恶化，反之，若回归系数小于0，则代表侵蚀模数呈下降趋势，即侵蚀程度减弱。本研究计算回归系数应用了Python编程语言，结合ArcGIS自带的Arcpy站点包以及处理数组的Numpy模块和进行数学分析的Scipy模块，分别计算了黄泛平原和丘陵地区的回归系数，并进行了显著性检验。

由年际变化率分布得到徐州市绝大多数地区的回归系数小于0，表明徐州市的土壤侵蚀模数整体呈下降趋势。对于丘陵地区而言，贾汪区大部、铜山区西南部以及东南部等轻度及以上侵蚀程度的地区下降趋势尤为明显，按照镇级行政区划分，这些地区主要包括贾汪区的利国镇、江庄镇、贾汪镇、青山泉镇、柳泉镇、茅村镇，铜山区的铜山镇、汉王镇、徐庄镇、张集镇，以及鼓楼区和泉山区的部分地区。对于黄泛平原来说，下降趋势较大的区域主要位于睢宁县北部，邳州市北部、东部，丰县南部，沛县东部，按镇级行政区划分包括邳州市的邳城镇、铁富镇、戴庄镇，睢宁县的姚集镇、古邳镇、岚山镇，沛县的五段镇、魏庙镇，丰县的大沙河镇、宋楼镇。

结合侵蚀变化率分布和侵蚀显著性分布得知，侵蚀程度降幅较大的区域多分布于各河流流域及湖泊周围，并且变化显著性强的区域也多分布于各河流附近，对于黄泛平原，侵蚀模数降低较大的区域分布于黄河故道北部，正集河与大沙河交叉处以及京杭运河部分段沿岸，通过了显著性水平为0.01或0.05的假设检验的地区为邳州北部京杭大运河(中运河)上游，京杭运河(中运河段)东部，徐洪河流域中部以及复兴河中段区域。对于丘陵地区，降幅较大的区域分布于京杭运河(梁济运河、湖西航道、不牢河)以北，奎河以西以及黄河故道以北，丘陵地区的某些部分通过了显著性水平为0.01或0.05的显著性检验，它们零星散布于黄河故道和京杭运河周围。对于丘陵地区，废黄河、奎河、京杭运河流域的部分地区仍存在着侵蚀程度加重的现象，对于黄泛平原，各流域均存在着侵蚀加重的情况，尤其在骆马湖沿岸存在着大面积侵蚀加重的状况。

2.5 成因讨论与分析

从以上的分析结果来看，多年间徐州黄泛平原及丘陵地区的土壤侵蚀程度整体呈减弱态势，产生此结果的原因主要包括以下3个方面。

(1) 徐州地区作为黄河故道区，历史上多次遭遇决口泛滥，使得该地区生态环境恶劣，森林覆盖率极低，20世纪80年代以后，徐州开始了全面的生态环境建设，其中封山造林、退耕还林是主要的手段，据统计，徐州地区的林地面积在1987—2016年增加了1 448.44 km2，增加量约是1987年林地总面积的2倍[26]。尤其在2013年，徐州市被水利部确认为首批生态文明城市建设试点后，有关部门加强了林业建设、修复和保护力度，截至2013年底，徐州的森林覆盖率达到了32%，位列全省第一名。森林覆盖度的增加直接减小了USLE模型中P因子的取值，进而促使徐州地区土壤侵蚀强度逐年减弱。

(2) 徐州境内水系众多，流域内土壤受河流的冲刷作用，易被河流分离、搬运和沉积，造成水土流失，并且在河流的治理过程中，人为的清淤、挖土和堆土，导致地貌、土壤和植被破坏，在水力、重力及风力等外营力作用下易引发强烈的水土流失[27]，针对这些问题，徐州市以小流域为主要突破口，发展了适用于徐州地区的小流域生态治理管理模式，取得了显著的水土保持效果。以新沂市山徐小流域为例，该流域治理前坡耕地面积为1 400 hm2，多为坡岗耕地，林地面积为386 hm2，其他用地(居民点)面积为105 hm2，水面、未利用土地面积为179 hm2。进行生态治理后，该流域内农地面积为1 000 hm2，优质果树、花木林面积300 hm2，林业用地面积达到486 hm2，其他土地利用类型面积保持不变，治理后该流域土地利用结构更加合理，充分利用现有土地资源，保护流域的生态环境。辖区内小流域生态治理模式的各种措施使得某些流域的土壤侵蚀大幅减弱，尽管如此，某些流域的侵蚀仍在加重，这是由流域侵蚀易发生并且易加重的自然本质决定的。

(3) 随着徐州市社会经济的快速发展，工业化、城市化进程的加快，生产建设项目造成的水土流失已成为本地水土流失的主要成因，《江苏省水土保持条例》规定：在易发生水土流失区域开办易导致水土流失的生产建设项目，生产建设单位应当编制水土保持方案，并通过有关行政主管部门审批。为此，徐州市近些年来进行了多项扰动土试验，为生产建设项目水土保持设计、施工、管理提供基础数据，为建设项目水土流失防治提供了技术支撑。随着试验科学性、完备性的逐渐提高，以及前人研究成果的不断积累，方案编写的合理性也逐步提高，使其应用于实际施工时，愈能减小水土流失的风险，减弱水土流失的强度。

3 结 论

(1) 15 a间徐州地区的土壤侵蚀模数总体呈下降趋势，即水土流失程度在逐渐减弱，其中丘陵地区的下降程度和显著性强于黄泛平原。近十几年来，徐州市完善了治理水土流失的方案，以监督预防为主，防治并重，同时加强水土保持监督、执法建设，一步步推进治理的科学化和合理化，使得水土流失工作取得了巨大的成功，尤其在徐州市区及其周围，即地理上的徐州市丘陵地区，效果更为明显，这与数据分析所得的结论相一致。

(2) 徐州黄泛平原和丘陵地区的土壤侵蚀以轻度(排除微度)侵蚀为主，从侵蚀的范围和等级跨度来看，丘陵地区的侵蚀程度更大。从15 a时间尺度的平均值来看，黄泛平原的轻度侵蚀面积占总侵蚀面积的90%左右，丘陵地区占总侵蚀面积的70%左右。各等级侵蚀面积随侵蚀程度的加重逐级递减，相较于黄泛平原，丘陵地区各等级侵蚀面积的变化规律更加明显，轻度以上等级土壤侵蚀面积随年份增加呈阶梯状下降态势。多年来，黄泛平原极强度侵蚀面积控制在5%以内，丘陵地区极强度侵蚀面积控制在10%以内，两区域剧烈侵蚀均得到了遏制。

(3) 对于黄泛平原，侵蚀模数下降趋势较大的区域主要位于睢宁县北部，邳州市北部、东部，丰县南部，沛县东部，对于丘陵地区而言，贾汪区、铜山区西南部以及东南部侵蚀模数下降趋势尤为明显。从侵蚀变化分布可以发现侵蚀减弱程度较大的地区多分布于各河流及湖泊周围，对于黄泛平原，这些区域分布在黄河故道、郑集河、大沙河、京杭大运河、徐洪河、复兴河周围，对于丘陵地区，这些区域分布于京杭运河、奎河、黄河故道附近。虽然流域治理取得了部分成果，但仍然存在着大面积的流域侵蚀加重情况，水土保持工作依旧任重道远，相信随着水土流失治理机制的日益完善、治理力度的日益加大、治理手段的日益科学化，徐州市的水土流失治理会收获更大的成效。

本研究存在的不足主要集中在数据获取上，尤其在计算植被覆盖因子(C)时，分别使用了SPOT_VEGETATION数据产品和Landast TM遥感影像两种数据源，从某种程度上影响了计算结果的统一性，若能提高获取数据的连续性、一致性和分辨率等属性，尽可能降低原始数据的误差，相信会得到更加精确的结果。后续工作的重点是提高数据精度，并利用实测数据验证结果的可靠性。