基于CSLE模型全覆盖和抽样调查法的天山北坡县域尺度土壤侵蚀计算方法的比较

李菊艳

（新疆维吾尔自治区水土保持生态环境监测总站，乌鲁木齐 830011）

土壤侵蚀造成的水土流失是一个世界性的严重生态环境问题。土壤侵蚀降低了土壤质量，导致区域生态经济损失和生态治理投入增大，及时把握区域水土流失状况、发展态势既是水土流失预防和治理的基础，更是水土保持规划制定与实施的基本前提。区域土壤侵蚀状况变化可作为判定水土流失治理成效的重要依据，对区域水土保持措施的优化配置、水土资源保护具有重大意义。

目前国内外对区域土壤侵蚀分析研究主要有模型法和抽样调查法，抽样调查法于1956 年始用于美国水保调查并沿用至今[1]。我国自20 世纪80 年代共开展了4 次土壤侵蚀普查[2]，并在第一次全国水利普查水土流失普查中采用了分层不等概率抽样法和模型法[3]。分层不等概率抽样法通过设定抽样调查单元密度实地调查采集单元数据，基于我国土壤流失方程CSLE采用不同推算方法定量估算土壤侵蚀量以获得区域土壤侵蚀强度面积及分布等[4-7]，该方法遵循统计学原理，样本代表性强同时结合实地调查，能够获取高精度的土地利用、植被覆盖度及水土保持措施等调查单元数据资料[8]。模型法通过遥感解译、GIS技术等方法获取覆盖全域土壤侵蚀因子的栅格数据，采用土壤侵蚀模型计算[7-8]以对调查区域土壤侵蚀进行评价。近年来，诸多学者对抽样调查密度设定、推算方法应用等方面进行研究和探讨[9-13]，认为栅格计算法受抽样密度影响较小。如赵维军等[9-10]通过对比分析吴起县不同抽样密度下土地利用、坡度及坡长等因子变化发现1%和4%抽样密度下计算的土壤侵蚀状况相差不大，认为1%抽样密度就能够很好地反映区域土壤侵蚀状况，适用性更高；邹丛荣等[11]、段倩等[12]在1%和4%的抽样密度下，通过采用单元直接外推法、单元插值外推法、栅格计算法等三种推算方法以及全覆盖计算法分析沂蒙山区土壤侵蚀状况，发现在这两个抽样密度下采用前两种方法估算的土壤侵蚀面积和强度相差较大，而栅格法计算的土壤侵蚀状况相差不大。齐斐等[13]以蒙阴县为对象采用单元插值外推法和栅格计算法对1%和4%抽样密度下土壤侵蚀因子进行精度评价，发现4%抽样密度下采用栅格计算法对土壤侵蚀结果影响较小。

CSLE 模型综合考虑了区域降雨、土壤、植被、地形及水土保持措施等，能较全面地预测区域土壤侵蚀过程，并客观、有效地反映水土流失治理成效[14]。基于CSLE模型的全域覆盖计算可获得区域总体情况并体现到图斑上，但工作量大，土壤侵蚀状况因子的区域本地化计算需要投入大量精力。新疆面积广大、生态系统脆弱，水土流失问题较为严重，当地政府予以高度重视。近年来经过大力治理，流失总面积显著减少、土壤侵蚀强度总体呈减弱趋势，其中强烈、极强烈侵蚀面积减少，轻度和中度侵蚀面积增加，水土流失状况显著改善，但局部地区问题仍然严重。因此水土流失预防和治理成效有待进一步巩固加强，治理任务仍然艰巨。在有限的技术力量与资金条件下，如何快速准确地掌握新疆地区土壤侵蚀状况、变化趋势及原因，探索适合当地的土壤侵蚀估算法已成为亟待解决的问题。

天山北坡经济和社会发展活跃，地形地貌和气候条件特殊，生态环境脆弱，受季节性强降水影响水蚀强烈，成为新疆水力侵蚀的典型区。本研究以天山北坡阜康市为例，基于CSLE 模型全域覆盖计算和密度抽样单元推算法估算土壤侵蚀量，通过对比分析不同估算方法结果的差异性及主要原因，明确适宜天山北坡县域尺度的抽样密度和土壤侵蚀推算方法，为区域水土流失动态监测野外抽样密度和方法的选取提供技术支撑和理论依据。

1 材料与方法

1.1 研究区概况

阜康市位于天山山脉东段北坡，准噶尔盆地南缘，总面积8 528.78 km2，东临吉木萨尔县，西与乌鲁木齐市米东区接壤，南以天山分水岭为界，北与古尔班通古特沙漠和阿勒泰地区富蕴县相接。地势自东南向西北倾斜，地貌由南向北形成高山、平原、沙漠的垂直景观地带，海拔高程450～5 445 m（图1）。研究区属于温带大陆性干旱气候，年均气温6.7 ℃，最高气温40.3 ℃，年均降水量205 mm，年均蒸发量2 064 mm；自西向东分布有水磨河、甘河子河、三工河等七条河流；土壤主要以风沙土、灰漠土为主；植被分布具有典型的垂直地带性，分布有高山亚高山草甸、中山森林植被、平原人工和荒漠植被，平原农区以人工植被为主，植被覆盖度以中覆盖度为主。

图1 研究区地理位置Figure 1 Location of the study area

1.2 研究区水土流失类型区划分

依据区域水土流失、地貌、气象及农业生产特征等，运用GIS 空间叠置分析工具采用先定量后定性分析的研究方法，综合考虑研究区地貌特征、气候特征、土壤类型、植被、土地利用、社会经济等定量指标及主要侵蚀外营力和集中连片等综合定性指标，开展阜康市水土流失类型区划分，并采用“地理位置+地貌类型+主导水土流失类型”方式对类型区进行命名。

1.3 数据获取及处理

遥感影像采用2015 年和2016 年2.0 m GF1 号影像，通过人机交互解译获取土地利用数据[15]、地形数据资料；结合区域土地利用现状和《土地利用现状分类》（GB/T 21010—2017），基于30 m 分辨率Landsat 8多光谱影像，提取归一化植被指数（NDVI）；采用1∶1万地形图，获得30 m DEM数据；降雨数据采用阜康市及其周边20 个站点1981—2018 年日降雨资料；土壤数据利用阜康市行政区划图和阜康市1∶1 万土壤类型图，基于ArcGIS软件，获得研究区土壤类型分布图。

1.4 野外抽样单元数据采集与处理

1.4.1 野外抽样单元布设

以第一次全国水利普查水土保持普查分层系统抽样原则，采用分层不等概率抽样方法，抽样密度设定为0.25%、1%和4%。综合考虑人为扰动和可到达性，剔除难以到达的沙漠区和深山区抽样单元，最终确定0.25%抽样密度调查单元20 个，1%抽样密度调查单元48个，4%抽样密度调查单元108个（图2）。

图2 阜康市抽样调查单元分布Figure 2 Distribution of sampling survey units in Fukang City

1.4.2 野外抽样单元数据采集与处理

以高分遥感影像和1∶1 万地形图为基础，叠加地形要素、单元边界，制作野外调查底图，勾绘边界。于2016 年和2017 年7—8 月，开展野外实地调查并详细记录单元地块土地利用、植被郁闭度和盖度、水土保持措施等数据信息[6，16]，获得单元调查成果。

1.5 水土流失推算方法

1.5.1 水土流失计算模型

采用我国土壤流失方程CSLE（Chinese soil loss equation）计算土壤侵蚀模数[17]：

式中：A为土壤侵蚀模数，t·hm-2·a-1；R为降雨侵蚀力因子，MJ·mm·hm-2·h-1·a-1，基于1980—2018 年阜康市及周边20 个站点日降雨资料，通过数据插补，获得逐日降雨序列数据，采用降雨侵蚀力因子计算法[18-20]获取24 个半月降雨侵蚀力占年降雨侵蚀比例的栅格图层；K为土壤可蚀性因子，t·h·MJ-1·mm-1，通过野外采集土壤样品，采用Williams 模型[21]计算土壤可蚀性因子K值，叠加土壤类型空间分布图，获得研究区土壤可蚀性K值栅格图层；L和S分别为坡长因子（无量纲）和坡度因子（无量纲），基于1∶1 万地形图生成30 m 栅格大小DEM，采用刘宝元等[6]修正算法，提取坡度、坡长因子，形成坡度、坡长因子栅格图层；B、E、T为水土保持措施因子（无量纲），分别指生物措施、工程措施、耕作措施，基于水利普查中提供的措施因子参考值，对地块进行赋值。

1.5.2 水土流失估算方法

（1）全覆盖计算法是基于土地利用矢量数据，对B、E、T因子进行赋值。B因子基于NDVI计算得到的24 个半月植被覆盖度，查表得园地、林地和草地土壤流失比例，以24个半月降雨侵蚀力比例为权重加权平均计算得植被覆盖度因子[20]；耕地、居民点及工矿用地、交通运输用地、水域及水利设施用地或其他土地，直接赋值为1或0，获得植被覆盖与生物措施因子栅格数据；E因子基于遥感影像解译和统计调查获得的水土保持措施和土地利用数据，进行工程措施因子赋值[20]，阜康市工程措施因子统一赋值为1；依据第一次全国水利普查水土保持普查中的轮作制度分区轮作措施赋值表[20]，研究区属西北干旱灌溉一熟兼二熟区，北疆灌溉一熟填闲区，T因子统一赋值0.281。

（2）抽样调查法是以土地利用图层为基础，对调查单元数据采用栅格计算法进行加权平均[11]，进行B、E、T因子赋值，得到各调查单元的B、E、T因子。

运用ArcGIS，将R、K、L、S、B、E、T等7 个因子经重采样生成10 m 分辨率栅格图层，采用CSLE 模型（公式1）进行所有图层的空间乘积运算，得到阜康市10 m×10 m 的基于栅格土壤侵蚀模数分布图。根据《土壤侵蚀分类分级标准》（SL 190—2007）[22]中的土壤侵蚀强度分级标准，判断每个栅格土壤侵蚀强度，获得区域各土地利用类型土壤侵蚀强度面积及比例。

2 结果与分析

2.1 水土流失类型区划分

根据划分原则、依据、方法等，将阜康市划分为4个水土流失类型区（图3），即北部沙漠风力侵蚀主导区、中部山前冲积扇及倾斜平原风水复合侵蚀主导区、南部中低山丘陵水力侵蚀主导区和南部高山冻融侵蚀主导区。据调查，阜康市水力侵蚀主要影响范围为南部中低山丘陵区和中部山前冲积扇及倾斜平原区，侵蚀面积总计3 758.06 km2，本研究仅在上述区域范围内开展水力侵蚀的定量计算与分析。

图3 阜康市水土流失类型区Figure 3 Soil and water conservation regionalization of Fukang City

2.2 基于CSLE模型的全覆盖计算结果

研究结果表明，基于CSLE模型的全覆盖计算，阜康市水土流失面积为2 626.62 km2，占土地总面积30.80%。在水土流失区域中，阜康市水力侵蚀以轻度侵蚀为主，其次为中度，两者总占比达81.73%；而强烈及以上等级占比不足20%，高强度侵蚀随侵蚀强度增大呈降低趋势，剧烈侵蚀面积占比最小（表1）。

从图4 可知，阜康市水力侵蚀强度空间分异性明显。轻度侵蚀主要分布在中部耕作区、植被盖度较好的冲积扇区域及海拔1 200～2 800 m间郁闭度、盖度较高的中山森林区域，且呈集中连片分布；中度侵蚀主要发生在植被盖度中等的区域，分布较分散；强烈及以上等级侵蚀主要分布在海拔700～1 200 m间植被盖度低于30%、采矿或放牧等人为活动频繁的区域，以及海拔2 800 m 左右坡度大于35°且植被盖度低于45%的区域。

图4 基于CSLE模型法估算的阜康市土壤侵蚀图Figure 4 Soil erosion distribution in Fukang City based on CSLE model calculation

2.3 基于CSLE 模型的不同密度抽样单元土壤侵蚀计算结果

研究结果表明，基于0.25%、1%和4%三种抽样单元密度采用栅格计算得到阜康市水土流失面积分别为2 771.89、2 740.94 km2和2 635.26 km2，占土地总面积比例分别为32.50%、32.14%和30.90%。随着抽样密度增大，水土流失总面积和轻度侵蚀面积呈逐渐减小趋势，而微度和强烈侵蚀面积则呈逐渐增大趋势，但均表现为以微度和轻度侵蚀为主（表1）。在水土流失区域中，三种抽样密度下轻度与中度侵蚀总面积占比分别为78.65%、71.75%和77.42%，在强烈及以上侵蚀强度中，除了1%抽样密度表现为极强烈侵蚀面积占比大于强烈侵蚀外，0.25%和1%抽样密度下均表现为强烈侵蚀>极强烈侵蚀>剧烈侵蚀。

表1 阜康市不同土壤侵蚀强度面积及比例Table 1 Area and ratio of different erosion intensities in Fukang City

三种抽样密度计算结果均显示阜康市水力侵蚀空间分布具有明显的差异（图5、图6、图7），且总体较一致：轻度侵蚀主要分布在中部耕作区，呈集中连片分布；中度侵蚀主要在植被盖度30%～60%、坡度8°～25°的区域，分布较分散，且0.25%抽样密度下轻度和中度侵蚀区域比较大；而强烈及以上侵蚀强度在1%抽样密度下主要分布在三工河至白杨河之间的区域，4%抽样密度下主要分布在四工河至白杨河之间的区域。整体上随着海拔的逐渐升高，土壤侵蚀强度等级逐渐增大。

图5 基于0.25%抽样密度估算阜康土壤侵蚀图Figure 5 Soil erosion distribution in Fukang City based on 0.25%sampling density estimation

图6 基于1%抽样密度估算土壤侵蚀图Figure 6 Soil erosion distribution in Fukang City based on 1%sampling density estimation

图7 基于4%抽样密度估算土壤侵蚀图Figure 7 Soil erosion distribution in Fukang City based on 4%sampling density estimation

2.4 分层系统抽样下不同抽样密度结果比较

本研究以0.25%、1%和4%为抽样梯度，以4%抽样密度数据为基准，对比分析三种抽样密度下水土流失面积、强度及空间分布的差异性，结果表明0.25%密度下水土流失面积结果偏高136.63 km2，相对差异5.18%；1%密度下水土流失面积结果偏高105.68 km2，相对差异4.01%（表2）。0.25%抽样密度下各级土壤侵蚀面积平均相对差异8.76%，相差较大的侵蚀强度等级是微度和轻度，分别为偏低136.63 km2、偏高140.30 km2，相对于4% 抽样密度差异分别是12.17%和9.06%，差异主要位于中部山前冲积扇及倾斜平原区域。1%抽样密度下各级土壤侵蚀面积平均相对差异28.99%，相差较大的土壤侵蚀强度等级是中度、极强烈和剧烈，分别偏低70.77 km2、偏高98.08 km2和90.16 km2，相对差异分别是14.39%、42.42%和104.32%，差异主要位于南部中低山丘陵区域。

表2 0.25%、1%与4%抽样密度土壤侵蚀强度面积对比分析Table 2 Comparative analysis of soil erosion intensity area at 0.25%，1%and 4%sampling density

各抽样密度土壤侵蚀空间分布规律基本一致。栅格计算法采用的R、K、L、S因子区域统一，仅B、E、T因子根据野外调查单元数据面积加权平均，略有差异。总体上看，基于高分遥感数据进行水土流失定量分析，0.25%、1%和4%三种抽样密度侵蚀面积总体相差不大。因此，对县域尺度来说，基于CSLE模型采用栅格计算，0.25%的抽样密度即可基本满足样本推算需要，同时还能减轻外业工作量。

2.5 CSLE模型全覆盖计算与抽样调查推算结果比较

研究结果（表1）表明，基于CSLE 模型全覆盖计算阜康市水土流失面积与0.25%、1%和4%抽样密度相比结果均偏小。全覆盖计算与0.25%和1%抽样密度计算结果流失面积相差较大，分别偏小145.27 km2和114.32 km2，差异为5.53%和4.35%；与4%抽样密度计算结果相比仅偏小8.64 km2，两者相对差异仅为0.33%。全覆盖计算结果和不同密度抽样调查结果均以轻度侵蚀为主，且仅0.25%密度下轻度侵蚀面积高于全覆盖计算。0.25%密度下除中度和剧烈侵蚀面积与全覆盖结果相差不大外，其余侵蚀强度面积均显著高于全覆盖计算结果，而1%和4%密度下强烈及以上侵蚀强度面积均与全覆盖计算面积差异显著。全覆盖和不同密度抽样调查水土流失在空间分布上有一定差异，主要表现在前者强烈及以上侵蚀强度分布相对较分散，而抽样调查结果相对集中于四工河和白杨河之间区域。

3 讨论

基于CSLE 模型的全覆盖计算，土壤侵蚀措施因子为遥感图斑解译结果，解译判读的准确性对各因子尤其是土地利用、植被覆盖及水土保持措施等的判断和赋值有直接影响，若数据分辨率相对较粗，会影响土壤侵蚀结果精度[23]；而抽样调查通过单元实地调查土地利用类型、水土保持措施等，能够提高地形因子量化精度，准确反映郁闭度和植被覆盖度对水土流失的影响[8]。本研究基于CSLE 模型的全覆盖计算和4%抽样调查估算阜康市土壤侵蚀，结果表明侵蚀强度和空间分布上均差异不显著。全覆盖计算土地利用类型采用2.0 分辨率遥感影像获得，植被指数基于30 m 分辨率数据获得，R因子通过数据插补获取，坡度、坡长基于30 m分辨率的DEM数据获取，因子数据源精度存在差异、匹配性差，对计算结果精度会造成一定影响，而对于大尺度区域，因子的率定和本地化则需要大量的精力；抽样调查采用的10 m×10 m 栅格计算法，数据精度较高，受抽样密度影响较小[13]，植被通过实地调查，考虑了林下植被覆盖问题，更接近实际，但抽样调查受地形等影响，致使有些深山区域无法开展调查。本研究全域覆盖和抽样密度在土壤侵蚀强度和空间上差异不大，但都因计算方法受限，计算结果可能都不能精准地反映实际情况。对于天山北坡县域尺度土壤侵蚀调查，若采用全覆盖计算，应提高侵蚀因子数据源的精度和匹配性；抽样调查中难以到达区域，可借助高分影像解译完成土壤侵蚀抽样。

基于高分辨率遥感影像、利用分层不等概率系统空间抽样方法，可快速提取土地利用和水土保持措施信息，完成区域土壤侵蚀抽样调查[7]，而抽样密度大小会对推算结果产生一定影响，且影响程度具有一定的地域差异性。本研究结果表明，0.25%、1%、4%三种抽样密度下采用栅格计算法估算的阜康市土壤侵蚀强度和空间分布规律基本一致，总体相差不大，这与邹丛荣等[11]、齐斐等[13]研究结果相似。本研究野外调查单元布设时，除难以到达的深山区域，其他区域各密度调查单元均匀分布，不同抽样密度调查B、E、T因子依据野外调查单元数据面积加权平均，基础数据差异不大，因此0.25%的抽样密度土壤侵蚀量计算结果与4%抽样密度下差异不大。

人类活动对原有地表的破坏会加剧水土流失，土地利用方式变化与土壤侵蚀具有密切联系[24]，区域的土地利用破碎程度和人为活动干扰强度与抽样密度精度损失呈正相关，土地利用破碎化程度越强，人为活动越强烈，其地形、土地利用、土壤等侵蚀因子空间异质性越大，抽样密度精度损失就越大[13]。阜康市地貌空间垂直分异明显，但土地利用类型和人类活动相对集中、相对成片，各抽样密度下土壤侵蚀强度空间分布也相对集中。不同抽样密度下，轻度和中度侵蚀面积在0.25%抽样密度下所占比例相对较大，而强烈侵蚀以上等级面积在1%和4%抽样密度下相对较大，这可能与1%和4%的抽样单元布设数目分布相对较多有关。通过对比分析，0.25%抽样密度既能满足土壤侵蚀精度计算要求，又可使工作效益最大化。

4 结论

（1）0.25%、1%和4%三种抽样密度下采用栅格计算法推算的基础数据相差不大，三种密度下推算的土壤侵蚀结果差异不显著，0.25%抽样密度既能满足土壤侵蚀精度计算要求，又可使工作效益最大化；基于CSLE模型的全域覆盖计算和抽样调查估算天山北坡阜康市土壤侵蚀状况，均可较好地反映阜康市水土流失空间分布特点。

（2）在天山北坡进行县域尺度水土流失动态监测时，若对区域土壤侵蚀状况进行基础调查，在有限的资源、资金和技术条件下，0.25%抽样密度下采用栅格计算法即可满足工作要求；若开展土壤侵蚀定量分析，既要掌握土壤侵蚀模数，也要把握土壤侵蚀强度分级，还需与全国水土流失动态监测工作衔接，建议采用基于CSLE 模型全覆盖计算方法，其中植被因子可通过遥感影像提取并结合0.25%的分层系统抽样获取。