基于CSLE模型的区域水土流失风险分析

苏新宇， 吴镇宇， 刘 霞†， 唐 俊， 李 想， 赵传普， 黎家作， 张春强

(1.南京林业大学南方现代林业协同创新中心，江苏省水土保持与生态修复重点实验室，210037，南京；2.淮河水利委员会淮河流域水土保持监测中心站，233001，安徽蚌埠)

水土流失受气象、地形、土壤、植被和人为活动等多种风险因素影响和制约[1]。区域中某处受一个或多个风险因素影响后，某处相比于其他位置水土流失可能发生的相对概率，称之为水土流失风险[2]，其结果为从最低到最高的水土流失风险等级，较好地反映区域内水土流失强度的空间差异性；因此，开展水土流失风险分析研究，对了解区域水土流失空间差异，合理配置水土保持措施以及制订水土保持政策具有重要的现实意义。

目前有关区域水土流失调查与评价方面的研究[3-6]积累较多，集中在定性与定量2方面。定性是以多因子综合评判和专家经验法为主，定性判断水土流失发生发展的可能[7-8]；定量则以经验模型和物理模型为主，通过模型计算土壤侵蚀模数，定量评价区域水土流失状况[9-11]。而开展区域水土流失风险研究的还比较少，虽有部分学者对水土流失引发的生态风险进行探索性研究[12-16]，但主要是采用多因子定性分析或美国通用土壤流失方程(universal soil loss equation,USLE；revised universal soil loss equation, RUSLE)，而基于中国土壤流失方程(Chinese soil loss equation, CSLE)开展区域水土流失风险分析的研究还未见报道。

蒙阴县属国家级重点治理区，该县经济林比例大，有“中国蜜桃之乡”的称号，山地和丘陵广泛分布，由于长期以来人类不合理的活动，自然植被几乎荡然无存，被人工防护林与用材林替代，坡耕地和坡地果园是水土流失的主要策源地，采矿用地水土流失也比较严重；因此，为了解蒙阴县水土流失发生可能性较高的区域，以及不同地形、植被覆盖和土地利用方式导致水土流失风险的差异性，探究未来水土流失治理和水土保持规划的重点区域。笔者以国家级水土流失重点治理区蒙阴县为研究区域，在遥感(remote sensing, RS)和地理信息系统(geography information systems, GIS)技术支撑下，基于CSLE模型定量计算土壤侵蚀模数，分析县域尺度水土流失风险因子空间分异特征，评估水土流失风险及其等级，探讨水土流失风险与地形地貌、植被覆盖以及土地利用的关系，以期为区域水土流失综合整治和水土保持管理决策提供依据。

1 研究区概况

蒙阴县地处北方土石山区(E 117°45′～118°15′，N 35°27′～36°02′之间)，土地面积1 602 km2。地势南北高，中间低，由西向东逐渐倾斜，地貌以低山丘陵为主，占97.77%。属暖温带季风型大陆性气候，多年平均气温12.8 ℃，多年平均降水量700 mm，降水时空分布不均匀，年际变化大。岩石类型以石灰岩和页岩为主，土壤包括粗骨土、棕壤、褐土和潮土等，以粗骨土为主。区内有梓河、东汶河、蒙河3条主要河流，属沂河水系。植被属暖温带落叶阔叶林区域，自然植被破坏严重，现多为人工植被，主要乔木树种有侧柏(Platycladusorientalis)、刺槐(Robiniapseudoacacia)、油松(Pinustableulaeformis)等；自然灌木与草本植物主要有黄荆(Vitexnegundo)、胡枝子(Lespedezabicolors)、三裂绣线菊(Spiraeatrilobata)等。

2 材料与方法

2.1 数据源及处理

基于1∶5万地形图，提取海拔、坡度与坡长，获取地形地貌特征数据；基于2017年2 m分辨率GF-1，通过野外调查和人机交互解译，获取土地利用与水土保持措施数据；基于MODIS遥感数据的NDVI产品，采用归一化植被指数计算林草植被覆盖度，并通过调查样地校核林草植被覆盖度计算值，获取林草植被覆盖特征数据；搜集沂蒙山区88个雨量站点30年日降雨数据计算降雨侵蚀力；搜集2018—2019年研究区及周边径流小区观测数据，为水土流失风险等级划分提供依据；搜集1∶5万土壤类型分布图、森林资源分布图，开展矢量化及配准，获取研究区土壤分布图。依据土地利用现状分类(GB/T 21010—2017)和区域实际，考虑人为活动及水土保持措施影响，将研究区土地利用类型划分为8个一级类、15个二级类(表1)。

表1 土地利用现状分类Tab.1 Land use classification

2.2 水土流失风险评估方法

2.2.1 土壤侵蚀模数计算 笔者基于CSLE模型计算土壤侵蚀模数，CSLE模型如下：

A=RKLSBET。

(1)

式中：A为土壤流失量，t/(hm2·a)；R为降雨侵蚀力因子，MJ·mm/(hm2·h·a)；K为土壤可蚀性因子，t·hm2·h/(hm2·MJ·mm)；L为坡长因子，量纲为1；S为坡度因子，量纲为1；B为生物措施因子，量纲为1；E为工程措施因子，量纲为1；T为耕作措施因子，量纲为1。降雨侵蚀力因子(R) 基于沂蒙山区88个雨量站点1986—2017年日降雨数据，通过逐日雨量公式法计算[17]，采用普通Kriging插值法，生成10 m空间分辨率降雨侵蚀力空间分布图。

土壤可蚀性因子(K) 基于1∶5万土壤类型分布图，采集典型土壤样品，测定土壤粒径和有机碳，根据Williams模型[18-19]计算土壤可蚀性因子K值，并采用径流小区观测数据修正。经重采样生成10 m空间分辨率的K因子数据图层。

坡度、坡长因子(LS) 基于1∶5万地形图，采用符素华等[20]和刘宝元等[21]的修正计算公式，提取坡度坡长因子，结合野外调查单元修正，获取10 m空间分辨率LS因子图层。

水土保持措施因子(B、E、T) 基于MODIS标准化植被指数产品和TM遥感数据，采用参数修订法[22]，获取24个半月植被覆盖度，依据不同土地利用类型确定各地类生物措施因子B，镶嵌生产10 m空间分辨率B因子图层。依据水土保持工程措施、全国轮作区名称及代码，赋值法确定工程措施因子E、耕作措施因子T，经重采样，生成10 m空间分辨率的ET因子图层。

2.2.2 水土流失风险评估 基于计算的土壤侵蚀模数，结合野外调查结果和研究区径流小区资料，探索多阶自然断点法与区域实际情况的结合，合理确定风险断点，研究土壤侵蚀风险等级及其空间分异性，探讨与风险影响因素之间的耦合关系。

3 结果与分析

3.1 水土流失风险因子空间分异性

3.1.1 降雨侵蚀力因子 降雨侵蚀力反映降雨能引起土壤侵蚀的潜在能力。研究区降雨侵蚀力因子(R)分布呈现由北到南依次升高趋势，见图1a。降雨侵蚀力介于3 333.05～4 033.19 MJ·mm/(hm2·h·a)之间，最小值位于北部岱崮镇，最大值在南部垛庄镇。在保持其他因子不变的情况下，仅考虑降雨单一因素影响，联成镇南部、桃墟镇、垛庄镇存在较高的水土流失风险。

图1 水土流失风险因子及水土流失风险等级分布Fig.1 Risk factors of soil and water loss and risk level distribution of soil and water loss

3.1.2 土壤可蚀性因子 土壤可蚀性反映土壤受到侵蚀的难易程度。研究区土壤可蚀性因子(K)值介于0.001～0.013 t·hm2·h/(hm2·MJ·mm)之间，≥0.012 t·hm2·h/(hm2·MJ·mm)的K因子在空间上呈两带两片状分布，见图1b。岱崮镇、坦埠镇、垛庄镇、以及旧寨镇东部和桃墟镇北部土壤类型主要是粗骨土，土壤可蚀性因子K值较高，若考虑土壤单一因素影响，上述区域面临较高的水土流失风险。研究区中部及南部K因子值平均较低，水土流失风险小，主要涉及野店镇、高都镇、旧寨镇、联成-桃墟-垛庄南部山区。

3.1.3 地形因子 地形是影响水土流失的重要因素之一，地形因子是坡度(S)与坡长(L)因子综合。总体来看，研究区地形因子分布呈现南北高中间低的特点，LS值介于0～31.56之间，见图1c。从空间上看，岱崮镇、野店镇、坦埠镇、桃墟镇以及垛庄镇的南部LS值较大，而中部的蒙阴镇、常路镇、高都镇和联城镇LS值较小，只考虑地形因素，南部与北部较中部区域更易引发水土流失。

3.1.4 生物措施因子 生物措施因子(B)能够定量衡量植被覆盖对土壤侵蚀的抑制作用，植被覆盖度高的区域，生物措施因子较低。研究区B因子介于0～1之间，从空间上看，南部蒙山山脉，中部蒙阴镇附近，B因子较小；其次为北部野店镇、高都镇、旧寨镇和岱崮镇；B因子值最大的区域为常路镇、联城镇、桃墟镇和垛庄镇(图1d)，若考虑植被单一因素影响，这些区域则面临较高的水土流失风险。

3.1.5 工程措施与耕作措施因子 水土保持工程措施可以改变地形条件并控制水土流失，耕作措施则通过增加地表覆盖和土壤入渗等提高土壤抗蚀性能。结合研究区实际，本文重点考虑梯田措施和耕作制度影响。研究区内梯田措施有土坎和石坎2类，土坎梯田主要分布在常路镇、高都镇、联城镇以及野店镇的北部，石坎梯田则在全县均匀分布。在中国耕作制度区划中，研究区属山东丘陵水浇地二熟旱坡地花生棉花一熟区(图1e)。

3.2 区域水土流失风险评估

3.2.1 水土流失风险等级

1)特殊土地利用土壤侵蚀风险界定：结合野外调查结果与专家经验，将城镇村建设用地、交通运输用地及河湖库塘直接归为无风险；对应生产建设活动扰动图斑(包括正在修建的城镇村建设用地、交通运输用地以及采矿用地)按照其类型，将采矿类项目归为高风险、非采矿类归为低风险；其他土地归为低风险。

2)多阶自然断点法：结合径流小区观测数据和水土流失较严重区域的侵蚀模数，先在0.2、0.3和0.5中确定0.5 t/(hm2·a)为“无”风险断点，又在25和50的侵蚀模数中确定了25 t/(hm2·a)为高风险断点，再根据自然断点法得到0.5～25之间的2个断点，最终得到0.5、2、10和25共4个断点侵蚀模数，直方图如图2，将水土流失风险等级划分为无风险、低风险、中等风险、较高风险和高风险5级。

图2 土壤侵蚀模数直方图Fig.2 Histogram of soil erosion modulus

3.2.2 水土流失风险分布 通过对不同区域水土流失风险等级统计分析，研究区水土流失风险等级的侵蚀模数阈值及不同等级面积见表2。总体看，随着风险等级升高，风险分布面积呈减少趋势。无风险等级区域占比最大，达45.75%，其次为中等风险、低风险等级，较高和高风险等级仅占8.78%。

表2 研究区水土流失风险等级分布统计Tab.2 Distribution statistics of soil and water loss risk levels in the study area

从空间分布看(图1f)，研究区水土流失风险空间分异性明显，中部以“无”风险和“低”风险分布为主，两边以“较高”风险和“高”风险等级分布为主，其中“高”风险与“较高”风险分布面积最大的乡镇为岱崮镇、野店镇、垛庄镇和坦埠镇，占各乡镇面积比例均>20%；而“无”风险分布面积最大的乡镇为蒙阴镇、高都镇和常路镇，占各乡镇面积比分别为64.34%、55.11%和55.16%。

3.3 水土流失风险性关系分析

3.3.1 水土流失风险与地形地貌 为分析不同地形地貌水土流失的风险情况，将研究区地貌与坡度，与水土流失风险进行叠置分析，得到不同地形地貌水土流失风险分布特征(图3和图4)。从地貌类型来看，“无”风险等级主要集中在平原区，其次为丘陵区，分别占其面积的80.72%和47.97%；“低”风险和“中等”风险在各地貌类型中均有分布，主要集中在丘陵与山地，但随着海拔与地表起伏高度增加，在各地貌类型中分布面积占比呈增加趋势；“较高”风险和“高”风险则主要分布在小起伏山地与中起伏山地，其次为丘陵区，平原区没有分布。

图3 不同地貌类型水土流失风险面积分布Fig.3 Area distribution of soil and water loss risk by different geomorphic types

图4 不同坡度等级水土流失风险面积分布Fig.4 Area distribution of soil and water loss risk at different slope levels

从不同坡度等级来看，“无”风险集中分布在0～2°，比例达67.99%，随坡度等级增加面积比例显著下降；“低”风险在>2°～5°、>5°～8°范围内面积比例最高，分别为27.40%和26.82%；“中等”风险随坡度等级增加面积比例增大，其中>8°～15°、>15°～25°、>25° 3个坡度带内的面积比例均>50%；“较高”风险主要分布在>5°～8°，面积比例为13.90%，8°以后的3个坡度等级范围内面积比例相近；“高”风险主要分布在>8°～15°坡度带，面积比例为10.82%。

3.3.2 水土流失风险与土地利用 为分析水土流失风险与土地利用关系，叠加分析水土流失风险图层与土地利用图层，获取不同土地利用类型水土流失风险分布特征(图5)。可见，城镇村建设用地、交通运输用地、水域及水利设施用地均为“无”风险，说明这3种土地利用类型不易产生水土流失，旱平地中“无”风险等级分布达62.35%，有林地36.26%；“低”风险主要分布在裸土地、裸岩石砾地，2种土地利用类型面积比例均>90%，平缓地果园中“低”风险面积比例达52.11%；“中等”风险主要分布在灌木林地、其他林地和其他草地，面积比例均达50%左右；“较高”和“高”风险主要分布在坡耕地、坡地果园与采矿用地中，其中采矿用地“高”风险面积比例达82.91%，这3种土地利用类型是区域水土流失治理的重点。

图5 不同土地利用类型水土流失风险面积分布Fig.5 Area distribution of soil and water loss risk by land use type

3.3.3 水土流失风险与林草植被覆盖 为分析不同植被的水土流失风险状况，将研究区植被类型分为林草地和园地，覆盖度划分为低覆盖(0～30%)、中低覆盖(>30%～45%)、中覆盖(>45%～60%)、中高覆盖(>60%～75%)、高覆盖(>75%)5个等级。叠置分析水土流失风险图层与植被覆盖度图层，得到不同植被类型及覆盖度下水土流失风险分布情况(图6)。

由图6a可看出，林草地“无”风险主要分布在中覆盖及其以上等级；“低”风险在低覆盖和中低覆盖面积比例较高，分别为36.1%和30.83%；“中等”风险在各覆盖度等级中面积比例相近；“较高”和“高”风险在各覆盖度等级中分布较少，均<10%，“较高”风险主要分布在中覆盖及其以下等级，“高”风险则集中在低覆盖和中低覆盖等级。

由图6b可看出，园地中“无”风险在各覆盖度等级中分布均最多，但随着覆盖度降低面积比例呈下降趋势；“低”风险与“中等”风险在各覆盖度等级中面积比例相近；“较高”风险集中分布在低覆盖和中低覆盖，面积比例分别为14.33%和14.86%，“高”风险等级则在低覆盖中面积比例达到最大为10.14%，这2种风险等级面积比例均随覆盖度增加而降低，表明随着植被覆盖度增加，水土流失相比于其他区域发生的相对概率也随之降低。因此，中覆盖及其以下等级的园地和林草地水土流失风险较高，是区域水土流失综合治理和封育保护的重点。

图6 不同植被覆盖度等级水土流失风险面积分布Fig.6 Area distribution of soil and water loss risk at different coverage levels

4 讨论

自然断点法是依据数据固有的特征来进行断点分级，实现各类别之间的差异最大化，断点选择的个数直接影响分级体系。本研究在分析不同土地利用对水土流失风险响应的基础上，探索自然断点和多阶段自然断点分析法在断点及其个数选择上的应用，合理确定研究区0.5、2、10和25 4个断点侵蚀模数，水土流失风险等级相应划分为无风险、低风险、中等风险、较高风险和高风险5级。结合周边径流小区定位观测数据，依据多阶段自然断点法确定的风险等级适合于研究区域，但是否适用于其他区域仍有待进一步研究。

本研究中采矿用地不单单指采矿、采石、采(砂)沙场、砖瓦窑等地面生产用地、排土(石)及尾矿堆放地等生产建设项目，还包括正在扰动的在建生产建设项目。CSLE模型不能完全适应这类项目的土壤侵蚀计算，目前，采矿用地主要采用遥感和实地调查手段，依据坡度和水土保持措施实施程度判定，但措施实施程度判定受人为因素影响较大，同时因数据源限制，解译的措施数量往往与实地调查措施状态有一定差距。如何开展采矿用地土壤侵蚀的定量计算是CSLE模型应用中面临的问题，仍需进一步研究探索。

5 结论

1)根据多阶自然断点法确定断点模数为0.5、2、10、25 t/(hm2·a)，分为无、低、中等、较高和高风险5个等级。蒙阴县水土流失风险以“无”风险为主，占45.75%，集中在该县中部；“较高”和“高”风险等级面积比例最少，主要分布在该县南部和北部，是区域水土流失防治重点。

2)山地丘陵水土流失风险明显高于平原，>8°～15°及其以上坡度等级水土流失风险高。“较高”风险和“高”风险等级主要集中于>8°～15°、>15°～25°坡度等级，分别占10.82%和9.00%；随着海拔与地表起伏高度增加，“低”风险以上等级在各地貌类型中面积比例呈增加趋势。

3)土地利用方式导致水土流失风险差异，“较高”和“高”风险等级集中分布在坡耕地、坡地果园与采矿用地，“高”风险在这3种土地利用中分别占41.59%、30.52%和82.91%；“低”风险主要集中在裸岩石砾地、裸土地和平地果园；“无”风险集中在城镇村建设用地、交通运输用地、水域及水利设施用地。

4)中覆盖及其以下等级的园地和林草地是区域水土流失综合治理和封育保护的重点。“较高”和“高”风险在林草地中主要分布在中覆盖及其以下等级，在园地集中在低覆盖和中低覆盖等级，“中等”和“低”风险在各植被覆盖度等级中面积比例相近。