基于栅格像元和矢量地块的区域土壤侵蚀计算结果对比分析——以怀来县外井沟小流域为例

邹海天，李子轩，张彩云

（1.海河水利委员会海河流域水土保持监测中心站，天津 300170；2.怀来县水务局，河北 怀来 075400）

随着国家深入推进长江经济带发展、黄河流域生态保护和高质量发展以及京津冀协同发展等战略的实施，水土流失等生态环境问题进一步得到重视，水土流失综合治理力度逐年加大，迫切需要水土流失监测数据为水土流失综合治理、规划设计和效果评价提供数据支撑。近年来，全国采用遥感解译、定位观测与栅格模型计算等方法，开展水土流失动态监测，获取了区域的水土流失面积、强度和分布[1]。基于栅格像元计算，能够快速掌握区域内水土流失及分布情况，但在同一地块内存在不同土壤侵蚀强度等级的栅格像元，这就导致无法直接确定地块土壤侵蚀强度等级，结果难以用于水土流失综合治理图斑设计。本文选取河北省怀来县外井沟小流域，采用基于栅格像元和矢量地块2种方法分别计算获取小流域内水土流失状况，对比分析2种方法计算结果差异性和适用性，为水土流失监测成果应用提供参考。

1 研究区概况

外井沟小流域位于河北省怀来县南部，距县城37 km，面积56.93 km2，处于永定河上游，属燕山山地丘陵水源涵养生态维护区，以水力侵蚀为主。地势南高北低、东高西低，最高海拔1 653 m，最低海拔810 m，相对高差843 m，地貌形态主要为山地。属温带大陆性季风气候，四季分明，年均气温6℃，7月温度最高，1月温度最低，无霜期125 d，冻土深大于1 m。多年平均年降雨量550 mm，多集中在6—9月，多年平均径流深85 mm。土壤类型以褐土为主，岩石以白云岩为主，兼有少量石灰岩和石英砂岩。植物类型以次生类植物为主，多为森林植被以及灌丛。近年来，小流域内实施了京津风沙源治理工程等水土流失综合治理，具有梯田、水平阶等水土保持工程措施和人造林等水土保持生物措施。

2 材料与方法

2.1 数据及处理

（1）土地利用和水土保持措施数据。利用空间分辨率优于2 m的遥感影像开展解译工作，所有解译最小地块面积为400 m2。

（2）地形数据。采用空间分辨率为12.5 m的数字高程模型，通过坡度坡长因子计算工具[2]，提取坡度、坡长和CSLE模型中涉及的坡度因子、坡长因子数据。

（3）土壤数据。源于1∶50万土壤类型图。

（4）植被覆盖度数据。采用2016—2018年每年24期MOD13Q1 NDVI数据，依据章文波[3]等提出的计算方法计算植被覆盖度数据，确定CSLE模型中的植被覆盖与生物措施因子。

（5）CSLE模型中涉及的降雨侵蚀力因子、土壤可蚀性因子和耕作措施因子。采用全国水土流失动态监测项目数据。

2.2 土壤侵蚀模数计算方程

采用CSLE模型计算土壤侵蚀模数[4]。

式中：A为土壤侵蚀模数[t/（hm2·a）]；R为降雨侵蚀力因子[MJ·mm/（hm2·h·a）]；K为土壤可蚀性因子[t·hm2·h/（hm2·MJ·mm）]；L为坡长因子；S为坡度因子；B为植被覆盖与生物措施因子；E为工程措施因子；T为耕作措施因子。

2.3 基于栅格像元的土壤侵蚀计算和统计

运用基础数据计算各因子，经重采样，生成10 m×10 m分辨率栅格，运用ArcGIS软件，对7个图层进行乘积运算，得到每个栅格的土壤侵蚀模数，按照《土壤侵蚀分类分级标准》（SL 190-2007），判断每个栅格的土壤侵蚀强度，统计不同土壤侵蚀强度的面积和比例。

2.4 基于矢量地块的土壤侵蚀计算和统计

根据CSLE模型计算原理可知，计算结果受土地利用类型、水土保持措施类型、降水、土壤、地形、植被覆盖等因素影响。为尽可能保证同一地块内计算因子相对唯一，在遥感解译的土地利用和水土保持措施基础上，依据土壤类型、植被覆盖度、地形数据进一步细化确定矢量地块。依据影响因素确定矢量地块，获取每个地块的R、K、L、S、B、E和T因子值，在矢量文件属性表中新建属性列，运用ArcGIS软件对7个因子进行乘积运算，得到每个地块的土壤侵蚀模数，按照《土壤侵蚀分类分级标准》（SL 190-2007），判定每个地块土壤侵蚀强度，统计不同土壤侵蚀强度的面积和比例。

3 结果与分析

3.1 基于栅格像元计算的外井沟小流域土壤侵蚀结果

基于栅格像元计算的外井沟小流域共有水土流失面积31.13 km2，占小流域面积的55.02%，其中轻度侵蚀31.05 km2，占小流域面积的54.88%；中度侵蚀0.03 km2，占小流域面积的0.05%；强烈侵蚀0.04 km2，占小流域面积的0.07%；极强烈侵蚀0.01 km2，占小流域面积的0.02%。

按土地利用类型划分，小流域内林地水土流失面积占总水土流失面积的66.65%，草地28.22%，耕地2.70%，建设用地、园地和交通用地分别为1.35%、1.04%和0.04%。

3.2 基于矢量地块计算的外井沟小流域土壤侵蚀

（1）基于土壤侵蚀影响因素确定的外井沟小流域矢量地块解译结果。外井沟小流域解译土地利用和水土保持措施矢量地块912个，平均16个/km2，最大地块为有林地，面积2.05 km2；最小地块为灌木林地，面积429 m2。叠加土壤类别数据后，矢量地块数量为1 117个；叠加植被覆盖度分级数据并消除面积小于400 m2地块后，矢量地块数量为2 512个；叠加地形数据并消除面积小于400 m2地块后，矢量地块数量为7 833个。

结果显示，植被覆盖和地形对地块提取数量影响较大，小流域内提取的园地、林地和草地内不同等级植被覆盖地块有1 545个，提取的不同坡面和坡度地块有1 808个，在叠加植被覆盖和地形因素后，矢量地块较未叠加该因素前分别增加了2.25倍和3.12倍。不同影响因素确定矢量地块情况，详见表1。

表1 不同影响因素确定矢量地块情况

（2）基于矢量地块计算的外井沟小流域土壤侵蚀结果。基于各土壤侵蚀影响因素确定后的矢量地块进行土壤侵蚀强度计算，结果显示：外井沟小流域共有水土流失面积33.77 km2，占小流域面积的59.68%，其中轻度侵蚀33.71 km2，占小流域面积的59.58%；中度侵蚀0.02 km2，占小流域面积的0.04%；强烈侵蚀0.04 km2，占小流域面积的0.06%；无极强烈侵蚀。

按土地利用类型划分，小流域内林地水土流失面积占总水土流失面积的67.98%，草地28.94%，建设用地1.88%，耕地、园地和交通用地分别为0.57%、0.52%和0.11%。

3.3 基于栅格像元和基于矢量地块计算方法结果分析

对比2种方法计算结果发现，外井沟小流域基于栅格计算的水土流失面积比基于矢量地块计算的结果小2.64 km2，二者相对差异率为8.47%。外井沟小流域2种计算方法的土壤侵蚀结果，详见表2。从土壤侵蚀强度等级上看，基于栅格像元计算的结果在中度及以上侵蚀等级面积大于基于矢量地块计算的结果。从空间分布上看，基于栅格像元计算的结果与基于矢量地块计算的结果在小流域内的水土流失位置整体没有明显变化，但有2点区别较大：一是基于矢量地块计算的结果呈空间连续状态，基于栅格像元计算的结果多呈离散分布状态；二是中度以上侵蚀的位置变化较大，如图1所示。从土地利用类型上看，二者计算结果中林地、草地水土流失面积占比较大，但基于栅格计算结果的耕地水土流失面积大于基于矢量计算结果。

表2 外井沟小流域2种计算方法的土壤侵蚀结果

图1 基于栅格像元和基于矢量地块计算的土壤侵蚀

4 讨论

4.1 土壤侵蚀因子数据对计算结果的影响

基于栅格像元和基于矢量地块2种计算方法均采用CSLE模型。模型中R、K因子均采用全国统一数据，比例尺小于1∶50万；L、S因子采用空间分辨率为12.5 m的数字高程模型数据计算获得；B因子采用空间分辨率为250 m的MODIS数据计算获得；E、T因子基于遥感解译的水土保持措施和土地利用数据，依据《区域水土流失动态监测技术规定（试行）》直接赋值。基于栅格像元计算可直接、精确地反映每个栅格像元的水土流失情况，但也受因子数据的空间位置精度、准确性和时效性影响，在同一地块上因子值有多种情况。基于矢量地块计算采用土地利用类型、水土保持措施类型、降水、土壤、植被覆盖度、地形等因素确定矢量地块，能有效地保证同一地块上各因子值的相对唯一。此处理虽然会概化地块内水土流失情况，但能规避因栅格像元值不一致导致的同一地块计算出多种土壤侵蚀强度的问题，实现土壤侵蚀计算结果分布相对集中，更好地体现地块总体水土流失情况。

4.2 矢量地块确定方法

确定矢量地块实际上就是确定计算土壤侵蚀计算的最小单元，且保证每个单元各属性相对唯一。确定矢量地块的方法有多种，比如网格法、地块法、叠置法等[5]。为尽可能避免计算误差，本文直接依据土壤侵蚀模型计算因子来分割矢量地块，直接获取矢量地块上模型计算因子值。外井沟小流域范围较小，降水空间分布差异不大，矢量地块确定时暂不考虑降水因素影响。首先，基于高分辨率卫星遥感影像目视识别地块土地利用和水土保持措施类型，判定其属性。其次，利用土壤类型数据，确定属于同一土种的地块。再次，在园地、林地和草地等地类上，依据《土壤侵蚀分类分级标准》（SL190-2007），确定植被高、中高、中、中低和低覆盖等级，进一步细分地块，减少地块内因植被覆盖度数据差异对水土流失计算结果产生影响。最后，由于地面坡度的大小直接影响水土流失情况[6]，因此，本文提取小流域内坡面和坡度，并采用坡面和坡度分级数据，确定处于同一地形的矢量地块，减少因坡度差异较大、地形突变等对水土流失计算结果的影响。

4.3 矢量地块确定影响因素

矢量地块确定方法直接影响计算效率和结果准确性。在数据处理过程中发现，土地利用、水土保持措施、土壤类型和降水较容易获取矢量地块计算所需的数据。但植被覆盖度和地形数据因存在像元值空间上的连续性，在分级后植被覆盖度、坡度和坡面数据进一步细化矢量地块时，叠加后的矢量地块特别破碎，较难实现地块单元化。本文应用了2种方法优化植被覆盖度和地形数据：①将分级后植被覆盖度、坡度和坡面数据由栅格数据转换为矢量数据，再与之前确定的矢量地块叠加计算，软件自动优化叠加后图斑边界和消除碎屑多边形；②将前序获取的矢量地块边界叠加分级后植被覆盖度、坡度和坡面栅格数据，人工提取植被覆盖度、坡度和坡面对矢量地块的分割线，逐个确定新的矢量地块边界。这2种方法前者用时较短，但需处理大量碎屑多边形等问题；后者精准度较高，但用时较长，在进行大尺度范围的土壤侵蚀计算时不易实现。

5 结论

（1）基于栅格像元和基于矢量地块计算均可实现小流域土壤侵蚀评价且精度相对一致，整体空间分布大致相同，结果差异为8.47%。产生差异原因与二者的最小计算单元、土壤侵蚀各因子计算方法等因素有关。基于栅格像元计算时每个计算单元都是10 m×10 m的栅格，但呈现在地块上的计算结果存在多样性；基于矢量地块计算时计算单元存在不确定性，依据计算单元所处位置的土地利用、水土保持措施、降水、植被覆盖度和地形综合情况而定，但呈现在地块上的计算结果是唯一的。

（2）在外井沟小流域，基于栅格像元和基于矢量地块计算结果在高等级侵蚀强度、耕地上水土流失面积差异大，主要原因是在确定矢量地块时对模型计算因子进行均值处理，从而对地块上土壤侵蚀情况刻画不如基于栅格像元计算结果细致，类别面积越小差别幅度越大。

（3）从计算过程和结果来看，基于栅格像元计算更快捷、简单，可快速掌握区域水土流失及分布情况，但因模型计算因子数据精度不一，计算结果空间上呈现插花状态，连续性不强；基于矢量地块计算结果空间分布集中连片，水土流失地块边界清晰，有助于实现土壤侵蚀地块单元化管理，但受计算数据精度、质量和计算单元划分影响，数据处理分析工作量较大。

（4）在今后的水土流失监测工作中，需进一步优化模型计算因子，提高数据的精度和质量，确保不同数据层在空间上能有效匹配，使计算结果既有数值计算上的准确性，又有空间分布上的合理性，从而保证水土流失监测结果在小流域水土流失综合治理规划设计和效益评价中得到有效利用。