利用CSLE 模型的东北漫川漫岗区土壤侵蚀评价

顾治家，谢 云，李 骜，刘 刚，史燕东

（1. 信阳师范学院地理科学学院，河南省水土环境污染协同防治重点实验室，信阳 464000；2. 北京师范大学地理科学学部，地表过程与资源生态国家重点实验室，北京 100875；3. 珠江水利科学研究院，广州 510611）

0 引 言

土壤侵蚀是指水力和风力等外营力对土壤的消损过程，它直接导致土地退化、江河库塘等淤积、水体富营养化等一系列的环境生态问题，已经成为了全球性环境问题[1-5]。防治土壤侵蚀首先要明确区域内土壤侵蚀的强度、面积及分布，才能有的放矢地布设水土保持措施。因此，区域土壤侵蚀评价受到世界上许多国家和地区的广泛关注。比如美国于1934 年组织的全国侵蚀勘探调查（National Erosion Reconnaissance Survey）[6-7]，这次调查初步确定了发生在农地上的土壤侵蚀强度分级面积[8]。澳大利亚运用土壤侵蚀模型进行了多次土壤侵蚀评价[9-10]。20 世纪90 年代至21 世纪初，欧盟实施的土壤侵蚀危险性评价项目，先后在整个欧洲、欧洲内不同区域或国家采用不同方法进行了侵蚀危险性评价[11]。中国开展的4 次大规模的“全国土壤侵蚀调查”，为今后的区域土壤侵蚀动态监测与评价打下了坚实的基础[12-13]。在土壤侵蚀评价中，通用土壤流失方程（Universal Soil Loss Equation，USLE）[14]和修订版的通用土壤流失方程（Revised Universal Soil Loss Equation，RUSLE）[15]得到了广泛的应用[16-18]。Liu 等[19]在USLE模型的基础上，充分考虑中国坡面土壤侵蚀的实际特点和防治土壤侵蚀所采取的植被覆盖与生物措施、工程措施和耕作措施，建立了中国土壤流失方程CSLE。该方程在中国得到了众多科研工作者的认可，被应用于中国许多地区的土壤侵蚀评价中[20-26]。在中国开展的“第一次全国水利普查”（2010—2012 年）中，CSLE 模型用来评价全国范围的土壤水蚀状况。尽管如此，该模型在东北黑土区应用仍有待深入研究。东北黑土区是中国重要的粮食生产基地，土壤侵蚀使得该区黑土层厚度逐渐变薄，土地生产力下降，严重破坏宝贵的黑土资源，威胁到中国粮食生产安全[27-28]。东北黑土区土壤侵蚀主要发生在漫川漫岗区，是土壤侵蚀防治的重点区。土壤侵蚀防治与水土保持规划需要科学合理的土壤侵蚀评价。但土壤侵蚀具有显著的多尺度特征，评价的标准不一，在坡面、沟道、小流域和区域尺度上，土壤侵蚀都有不同的表现方式与特征。在坡面上大家更关注土壤侵蚀模数，在沟道上更关注土壤侵蚀量，而在区域上土壤侵蚀面积往往是关注的对象。土壤侵蚀评价“落地”是指评价结果能“落”在具体的地块或小尺度范围上，给决策者提供明确具体的土壤侵蚀防治区域。纵观已有的研究，区域土壤侵蚀评价的结果在空间上通常以土壤侵蚀模数、侵蚀量或者侵蚀强度方式呈现，但受土地利用、地形和水保措施的影响，这些指标的空间变异性很大，并不能很好地实现侵蚀评价的“落地”，对水土保持规划的指向性不明确。在地块和小流域范围内，土壤侵蚀评价能够很好地实现“落地”，但在区域尺度上，需要探索更为合适的空间呈现方式。因此，本文选择位于侵蚀较为严重的东北漫川漫岗区拜泉县为例，运用CSLE模型，评价其土壤侵蚀状况，进而确定土壤侵蚀的重要源区。尝试提出土壤侵蚀评价“落地”的合理方案，以期为区域土壤侵蚀防治和水土保持规划提供科学依据与合理参考。

1 研究区概况

拜泉县位于黑龙江省的中西部、齐齐哈尔市的东北部，地处小兴安岭余脉向松嫩平原过渡的地带，属于东北漫川漫岗区。地理位置为125°30′E～126°31′E，47°20′N～47°55′N，总面积约为3 600 km2，耕地面积比例为57.4%，总人口56.85 万。拜泉县在地理位置上处于中高纬度、欧亚大陆的东部地区，属于典型的中温带大陆性季风气候，多年的平均气温为1.2 ℃、年均降水为 490.1 mm[29]。全县除沿河地带水量较丰富外，东部、中部和西南部地区的地表水资源十分匮乏，地下水埋藏很深，开采困难，利用率不高。土壤类型主要包括黑土、黑钙土、草甸土、沼泽土和盐土等，黑土为主要的耕作土壤。拜泉县的原始植被早已破坏殆尽，整体上具有从森林和草甸草原或草原化草甸交错分布的特点。拜泉县曾经是黑龙江省土壤侵蚀最为严重的县之一，早在20 世纪80 年代，拜泉县就被列入了省水土保持重点县，土壤侵蚀以水蚀为主，风蚀次之。

2 资料与方法

2.1 CSLE 模型及各因子计算

土壤侵蚀模数的计算采用中国土壤流失方程CSLE，其基本结构形式如下

式中，A 为土壤水蚀模数，t/(hm2·a)；R 为降雨侵蚀力因子，MJ·mm/(hm2·h·a)；K 为土壤可蚀性因子，t·hm2·h/ (hm2·MJ·mm)；L、S 为坡长坡度因子，无量纲；B 为植被覆盖与生物措施因子，无量纲；E 为工程措施因子，无量纲；T 为耕作措施因子，无量纲。

1）降雨侵蚀力因子R

R 因子的计算选用干雨季日降雨量估算模型[30]。搜集黑龙江省及其周边地区共46 个气象站点的逐日雨量资料，降雨数据年限为1986—2015 年。通过地统计插值得到黑龙江省24 个半月降雨侵蚀力比例的栅格图层和年降雨侵蚀力的栅格图层，空间分辨率为30 m。在此基础上，用拜泉县的矢量图层裁切得到研究区的年降雨侵蚀力栅格图层和24 个半月降雨侵蚀比例图层。

2）土壤可蚀性因子K

土壤可蚀性因子K 是直接采用了第一次全国水利普查中土壤可蚀性因子成果。所用的土壤可蚀性因子的计算采用Wischmeier 提出的K 值估算公式[31]和Williams 提出的K 值估算公式[32]。

3）坡长因子L 和坡度因子S

基于DEM 的坡长因子L 选用坡长法来进行计算。依据栅格所在的位置分别采用Wischmeier 等[14]、Foster 等[33]提出的坡长因子计算公式。

10°以下的区域的坡度因子S 选用McCool 等[34]的公式（2）和公式（3）；10°以上的坡度区域选用Liu 等[35]的公式（4）。

式中S 为坡度因子；θ 为坡度(°)。

覆盖拜泉县的DEM 数据由1∶5 万等高线转换而来，并生成空间分辨率为30 m 的栅格图层。基于“土壤侵蚀模型地形因子计算工具”[36]求得坡长因子与坡度因子。

4）植被覆盖与生物措施因子B

植被覆盖与生物措施因子赋值与计算需要依据不同土地利用来进行相应的赋值与计算。林地和草地的B 值计算公式如下

式中WRi为第i个半月降雨侵蚀力占全年降雨侵蚀力的比例，取值范围为0～1。SLRi为第i 个半月的林地和草地的土壤流失比例，取值范围为0～1。第i 个半月的有林地和其他林地的SLRi计算公式如下

灌木林地的 iSLR 计算公式如下

草地的 iSLR 计算公式如下

式中FVC 为基于遥感影像的NDVI 计算的植被覆盖度，取值范围为0～1；GD 为乔木林的林下盖度，取值范围为0～1，包括除乔木林冠层以外的所有植被（灌木、草本和枯落物）构成的林下盖度，按实地调查或经验取值。

植被覆盖度FVC 由250 m 分辨率的MODIS NDVI与基于30 m 分辨率的TM 影像计算的植被覆盖度进行融合得到。其中MODIS NDVI 影像和TM 影像时间段都来自于2015—2017 年。

式中VH(ti)为高分辨率像元的NDVI 融合值；VM(ti)为此高分辨率像元对应地类MODIS 多年平均值序列；VT(tj)为此像元对应某时段TM 影像的高分辨率NDVI 的数据，总计有n 景；ti为MODIS-NDVI 数据获取时的儒略日（DOY=16×时相-7）；tj为高分辨率NDVI 数据获取时对应的儒略日；ω(ti,tj)为tj时的高分辨率NDVI 的权重，表达为

将融合后的NDVI 转换为植被覆盖度。将上式中融合的24 个半月30 m 空间分辨率的NDVI 转换为相对应的植被覆盖度FVC

式中NDVI 是像元的NDVI 值；NDVImax与NDVImin为像元所在地类的转换系数；k 为非线性系数，在同一气候类型确定MODIS影像不同植被的NDVI最大值、裸土NDVI最小值所在像元，取该像元内TM 的NDVI 的平均值为转换系数。

最后，计算3 a 平均24 个半月植被覆盖度。将3 a对应的24 个半月植被覆盖度取平均值参与到B 值的计算当中，即为公式（10）中的FVC。

除林地和草地外的土地利用类型的B 值依据下表进行直接赋值（表1）。

表1 不同地类植被覆盖与生物措施因子B 赋值表 Table 1 Vegetation cover and biological practices factor value (B) of different land use types

5）工程措施因子E 和耕作措施因子T

在研究区的实地调查发现，工程措施主要包括梯田、地埂和极少量的水平阶等。工程措施因子的赋值借鉴第一次全国水利普查的水土保持工程措施因子赋值表（表 2）。对E 值进行重采样生成空间分辨率为30 m 的栅格大小（对应1∶5 万DEM）。

耕作措施因子反映2 种影响：一是轮作制度形成的作物覆盖影响，二是水土保持耕作措施。依据全国轮作区划，研究区属于东北平原丘陵半湿润喜温作物一熟区、三江平原长白山地温凉作物一熟区，反映作物覆盖影响的轮作措施T 因子赋值为0.331。研究区水土保持耕作措施主要是等高沟垄种植，对应T 值赋值为0.251。对T 值进行重采样生成空间分辨率为30 m 栅格大小（对应1∶5万DEM）。

表2 工程措施因子E 赋值表 Table 2 Value of engineering practices (E)

研究区范围的土地利用数据是基于2016 年覆盖研究区域的资源三号卫星（ZY-3）影像解译获得全区的土地利用数据，其空间分辨率为2.1 m，属高分影像，可一定程度上保障土地利用和工程措施解译的精度。

2.2 土壤侵蚀模数计算及强度判断

依据CSLE 模型，计算的土壤侵蚀模型进行侵蚀强度的判断，判断的依据标准为水利部《黑土区水土流失综合防治技术标准（SL446—2009）》中的分级标准[37]。东北黑土区土壤侵蚀的容许流失量为2 t/(hm2·a)，所对应的土壤侵蚀分级标准如表3。

表3 土壤侵蚀强度分级标准表 Table 3 Classification standard of soil erosion intensity

3 结果与分析

3.1 土壤侵蚀空间分布

将土壤侵蚀模型因子各图层乘积后得到拜泉县的土壤侵蚀模数（图1 a），基于土壤侵蚀模数和侵蚀强度分级标准，分别得到土壤侵蚀量分布图（图1 b）和土壤侵蚀强度分级图（图1 c）。结果表明，土壤侵蚀模数变 化 于 0 ～557 t/(hm2·a)，平 均 土 壤 侵 蚀 模 数 为6.77 t/(hm2·a)，土壤侵蚀面积比例为39.35%。模数小于2 t/(hm2·a)的面积比例最大，为60.65%，主要分布于拜泉县西南部平原区和靠近河流的低洼区。其次为轻度侵蚀和中度侵蚀，所占面积比例分别为16.56%和15.06%。强烈及其以上侵蚀的面积比例较小为7.73%。极强烈和剧烈侵蚀的高侵蚀主要分布于拜泉县的东南部丘陵区和西北部的漫川漫岗区，虽然所占的面积比例较小，仅为3.75%，但对侵蚀量的贡献达1/4，是侵蚀的主要源区。

图1 拜泉县土壤侵蚀空间分布 Fig.1 Spatial distribution of soil erosion in Baiquan County

土壤侵蚀模数越高，单位面积土地上发生侵蚀量越大，在流域、小流域或者坡面尺度上，模数可以用来计算土壤流失量，并与泥沙、径流、面源污染等密切相关。侵蚀量的分布图表示土壤被侵蚀的多少，侵蚀量越多土壤面临的危险性越高，河道、湖泊、库塘等面临的泥沙压力越大。依据土壤侵蚀模数，不同的区域有不同的分级标准，将结果以土壤侵蚀强度图的形式呈现。但这3种结果呈现方式大多呈现“斑点状”，不能确定侵蚀严重区的确切位置，并不能很好地实现侵蚀评价的“落地”，对水土保持规划的空间指向性不明确。因此，需要在此基础上探索更为合适的呈现方式。目前，中国正在进行水土流失重点预防区与重点治理区的划分[38-39]。水土流失重点防治区的“落地”则需要在县级尺度上进行科学的划定。在以往的水土保持实践中，以小流域为单元的水保规划与治理已经取得显著的成效，在拜泉县尤为如此。因此，首先本文在1∶5 万等高线生成的DEM基础上进行拜泉县小流域的划分，小流域平均面积约为11.7 km2，共划分306 个小流域，其中小流域面积最小的为2.1 km2，最大的为45.5 km2。依据土壤侵蚀强度分级标准将土壤侵蚀模数分为6 级，即分别为微度侵蚀、轻度侵蚀、中度侵蚀、强烈侵蚀、极强烈侵蚀和剧烈侵蚀。通常将轻度及其以上侵蚀面积之和视为土壤侵蚀面积，该面积与总面积的比值即为水蚀比例。其次，在拜泉县306 个小流域的范围内统计计算水蚀面积比例，求得的结果见图1 d。结果表明，约有1/5 的小流域土壤侵蚀面积比例大于 50%，占总土地面积的23.4%。因此，该1/5 的小流域就是土壤侵蚀防治的重点区域，这些侵蚀小流域都进行了编号，在空间位置上也相对固定，这就实现了区域土壤侵蚀评价的“落地”，可更好地为区域水土保持规划服务。拜泉县的北部和西北地区都是侵蚀较为严重的地区，是水土保持规划与治理的重点区域。

3.2 不同土地利用的土壤侵蚀

基于高分影像解译得到的研究区的土地利用类型，采用中科院土地利用覆盖分类体系，并对个别地类进行了合并，结果表明旱地所占的比例最大，为75.0%，其次为居民点及工矿交通用地和有林地，所占面积比例分别为10.5%和5.2%。有林地主要分布于拜泉县东部和西北部的丘陵地区，灌木林地和草地所占的面积比例仅分别为2.9%和2.5%，但分布面积相对较广（图2）。水田在拜泉县中也有一定的比例，但较小，为1.5%，主要分布于水源丰富、地势低洼的河流沿岸低地，集中于东部地区的通肯河流域和中西部地区的双阳河流域。水域及水利设施的面积比例略大于水田，为1.8%。其他用地的面积比例非常少，包括其他未利用地、裸岩或裸地和沙地。剩余土地利用类型都占很小的面积比例，不足0.5%。

图2 拜泉县土地利用类型分布图 Fig.2 Distribution map of land use in Baiquan County

拜泉县各主要土地利用类型的土壤侵蚀模数进行对比分析（图 3），旱地的侵蚀模数变化于 0 ～557 t/(hm2·a)，平均侵蚀模数最大，为8.9 t/(hm2·a)，发生于旱地上的土壤侵蚀面积比例高达62.7%，旱地又占据研究区绝大多数土地面积，因此，旱地为土壤侵蚀治理的重点区。灌木林地和草地的平均土壤侵蚀模数仅次于旱地，分别为1.4、1.3 t/(hm2·a)，变化范围分别为0～44 t/(hm2·a)和0～176 t/(hm2·a)，土壤侵蚀面积比例分别为20.8%和25.0%，由于其在研究区的分布面积较广且比较分散，应防止人为开发破坏。有林地和其他林地的土壤侵蚀强度最小，这主要得益于林地较高的植被覆盖度。

图3 拜泉县主要土地利用类型的土壤侵蚀模数 Fig.3 Soil erosion modulus of main land use types in Baiquan County

3.3 地形对土壤侵蚀的影响

地形是影响土壤侵蚀的基础因素之一，控制着土壤侵蚀的面积、强度及其空间分布。拜泉县的平均坡度较小为1.1°，最大坡度为19.8°。≤0.25°的面积比例最大为34.5%，其次为1～2°和0.25～1°坡度范围，所占的面积比例分别为23.1%和18.2%（图4 a）。3°以上的坡度共占比约11.1%，大于15°的坡度面积很小，不足0.1%。土壤侵蚀模数随坡度的增加而增大。坡度为8°以上的区域土壤侵蚀模数高达50 t/(hm2·a)（图4 b），为剧烈侵蚀，但一方面由于其面积比例较小（为0.3%），侵蚀多发生于旱地上，另一方面该坡度区域林地分布较广，因此并不是侵蚀的主要源区。从对侵蚀量的贡献来说，经计算，来源于1～2°的坡度范围侵蚀量最大，占比为27.9%，其次为3～5°和2～3°，贡献比例分别为26.6%和24.3%。拜泉县的平均坡长超过100 m，最大坡长超过3 000 m，达3 254 m。50 m 以上的坡长面积比例为43.8%，其中100～200 m 的坡长所占比例最大，约为16.2%，其次为50～100m 和200～400 m，所占比例分别为12.6%和9.8%。1 000 m 以上的坡长面积比例较小，约为0.4%（图 4 c）。侵蚀模数最大的坡长为200～400 m，其次为100～200 m 和50～100 m（图 4 d）。经计算，100～200 m的坡长对侵蚀的贡献最大，达35%，其次为200～400 m和50～100 m 的坡长，对侵蚀的贡献分别为22.0%和20.7%。东北黑土区的九三水土保持试验站（48°59′N～49°03′N，125°16′E～125°21′E）的野外径流小区土壤侵蚀监测数据表明（表4），坡长与土壤侵蚀模数的关系并非简单的线性关系，随着坡长的变长，土壤侵蚀模数大致呈现出先增大后减小的变化趋势。这与本文的研究结果相对一致。

图4 不同坡度坡长下的土壤侵蚀面积比例及其土壤侵蚀模数 Fig.4 Area proportion and soil erosion modulus at different slopes and slope lengths

表4 2011—2015 年野外径流小区的土壤侵蚀模数 Table 4 Soil erosion modulus of field runoff plots from 2011 to 2015

4 结 论

本文以地处东北漫川漫岗区拜泉县为例，针对土壤侵蚀问题开展研究，通过定量评价区域的土壤侵蚀状况，并提出了解决土壤侵蚀评价的“落地”方案。土壤侵蚀面积比例作为土壤侵蚀评价的结果呈现方式便于理解，在空间上能较好地实现“落地”，对水土保持规划和土壤侵蚀防治的指向性确定。本文得到的主要结论如下：

1）拜泉县平均土壤侵蚀模数为6.77 t/(hm2·a)，土壤侵蚀面积比例为39.35%，强烈及其以上侵蚀的面积比例较小为7.73%。极强烈和剧烈侵蚀的高侵蚀主要分布于拜泉县的东南部丘陵区和西北部的漫川漫岗区，是侵蚀的主要源地。

2）旱地的侵蚀模数最大，发生于旱地上的土壤侵蚀面积比例高达62.7%，旱地又占据研究区绝大多数土地面积，因此，旱地为土壤侵蚀治理的重点区。

3）土壤侵蚀模数与坡度成正比，坡度越大侵蚀强度越大。但随着坡长的变长，土壤侵蚀模数呈现出先增大后减小的变化趋势。对侵蚀量的贡献来说，来源于1～2°的坡度范围侵蚀量最大，其次为3～5°和2～3°。100～200 m 的坡长对侵蚀的贡献最大，其次为200～400 m 和50～100 m 的坡长。

4）在小流域的尺度上统计其土壤侵蚀面积比例，是实现区域土壤侵蚀评价“落地”的有效途径，可以更好地服务于区域水土保持规划。