摘 要: 为进一步研究基于 1%和 4%抽样密度，采用不同水土流失推算方法对县域水土流失监测成果的影响，为全国水土流动态监测提供技术支撑，本文以沂蒙山区沂水县为例，采用遥感与抽样单元调查相结合的方法，分析不同抽样密度下直接外推法、插值外推法和全覆盖计算法对水土流失状况的影响，并进行差异性分析。结果表明：（1）在 1%、4%抽样密度下，直接外推法、插值外推法水土流失面积比分别相差 6.49%和 7.69%，相对误差达 15.87%和16.22%，而全覆盖计算法水土流失面积比相差 0.88%，相对误差 1.97%；（2）1%抽样密度下，直接外推法与插值外 推法水土流失状况结果相近，但与全覆盖计算法结果有差异；4%抽样密度下，三种推算方法水土流失状况结果接近， 但是水土流失空间分布有明显区别。因此，基于抽样单元的直接外推法和插值外推法受抽样密度影响较大，全覆盖 计算法受抽样密度影响较小，建议采用 1%抽样密度下的全覆盖计算法进行水土流失动态监测。
关键词: CSLE 模型; 抽样密度; 推算方法; 县域尺度; 相对误差
        水土流失是全球普遍关注的生态环境问题。作为水土流失最严重的国家之一，我国已为治理水土流失投入了大量财力、物力和人力。为评价水土流失治理效果，开展水土流失动态监测具有重要的意义。目前，在国外水土流失调查计算过程中，美国主要采用抽样调查和通用土壤流失方程（USLE），澳大利亚主要采用基于 GIS 和 USLE 的水土流失网格估算，欧洲主要采用基于专家法和USLE 模型法[1-4]等。在国内，水土流失状况评价方法主要有定性评价法（综合评判法）以及 CSLE 模型法。前者在我国 80 年代、90 年代的水土流失遥感调查中广泛应用[5]；后者首次与抽样调查相结合大规模应用于 2010~2012 年开展的第一次全国水利普查水土保持情况普查中[6-9]，并在 2018 年开始的全国水土流失动态监测用于水力侵蚀计算。在 CSLE 模型应用中发现，抽样密度会影响小尺度尤其是县域尺度的水土流失状况推算结果，但是影响大小具有地域性差异；而且不同的水土流失推算方法也会影响最终的结果。赵维军、张岩等[10-12]以陕西吴起县为例，对比分析不同密度对水土流失因子精度影响，发现 1%密度和 4%密度均能很好地代表吴起县水土流失状况，并在 1%抽样密度基础上对比 CSLE 模型和综合评判法估算成果，认为前者更有优越性，但是二者以抽样单元内分析为主，均未涉及到不同抽样密度、不同推算方法对县域水土流失状况估算的影响；邹丛荣等[13]在山东蒙阴县研究 1%、4%抽样密度与三种水土流失推算方法对水土流失状况的影响，发现 4%抽样密度适用于直接外推法和插值外推法，采用 1%抽样密度适用于全覆盖计算法（栅格计算法）。故本研究以沂蒙山泰山国家级重点治理区沂水县为例，分别探索不同抽样密度、不同推算方法对县域水土流失状况的影响，对比分析差异和原因，进一步探讨适用于县域尺度水土流失调查与评价的抽样密度和推算方法，为县域水土流失动态监测方法提供数据支撑。
一、材料与方法
1、研究区概况
       沂水县地处沂山南麓，隶属山东省临沂市，面积 2420.08 km2，属北方土石山区。地貌主要以低山丘陵为主，地势自西北向东南倾斜，山丘区面积比 94.3%。土壤包括粗骨土、棕壤、褐土、红粘土和潮土五大类，以粗骨土为主。属暖温带季风型大陆性气候，多年平均气温 12.3 ℃，多年年均降水量 629 mm。水资源丰富，境内沂河为山东省第一大河，以及山东省第三大水库跋山水库。植被属暖温带落叶阔叶林区，主要乔木树种有刺槐（Robinia pseudoacacia L.）、侧柏（Platycladus orientalis （Linn.）Franco）、花椒（Zanthoxylum bungeanum Maxim.）、油松（Pinus tabulaeformis Carr.）等，自然灌木与草本植物主要有黄荆（Vitex negundo L.）、三裂绣线菊（Spiraea trilobata L.）、胡枝子（Lespedeza bicolors Turcz.）等。地理位置见图 1。
2、野外抽样单元数据采集与处理
（1）野外抽样单元布设 采用分层不等概抽样法[7,9,14,15]，以国家第一次水利普查布设的野外抽样单元为基础，在山丘区增加抽样单元，勾绘 0.2~3 km2 小流域，使其抽样密度增加至 4%，平原区维持不变。全县 1%抽样密度调查单元 24 个，单元控制区域 100 km2；4%调查单元 96 个，单元控制区域 25 km2。野外调查点分布见图 2。
图 1 沂水县地理位置    图 2 沂水县 1%和 4%抽样单元分布

（2）野外抽样单元数据采集与处理 基于野外调查底图和现场实际状况进行地块边界勾绘，填写野外调查信息表，详细记录单元内各地块土地利用、水土保持措施等信息[7]。调查结束后，将野外调查地图配准矢量化，逐地块填写地块信息。
3、水土流失推算方法
（1）水土流失计算模型 本研究采用中国水土流失通用模型（CSLE）[11]，根据《北方土石山区水土流失综合治理技术标准》（SL 665-2014）进行水土流失状况统计分析。CSLE 模型基本形式为：
A = R·K·L·S·B·E·T      （1）
式中：A 为土壤流失量（t×hm-2a-1）；R 为降雨侵蚀力因子（MJ×mm × hm-2  ×h-1×a-1）；K 为土壤可蚀性因子（t×hm2×h ×hm-2×MJ-1×mm-1）；L、S 为坡长、坡度因子（无量纲）；B、E、T 为水土保持措施因子（无量纲）。
降雨侵蚀力因子（R）：基于沂蒙山区 30 年 88 个雨量站点的日降雨数据，采用章文波[16]公式计算 R，利用普通克里金插值法获取 30 m 区域栅格图层[17]，通过裁剪获取沂水县 R 栅格图层。
土壤可蚀性因子（K）：根据 Williams 模型[18]和径流小区观测资料计算 K 值，并基于全县土壤类型图，获取全县土壤可蚀性因子图层。
坡长、坡度因子（LS）：基于 1:1 万地形图，经投影变换后，生成 10m 栅格大小 DEM，采用刘宝元[19]修正算法，提取坡度坡长因子。
水土保持措施因子（BET）：基于径流小区数据和水利普查中提供的措施因子参考值，对地块进行赋值[14,16]。
（2）水土流失推算方法 主要采用直接外推法、插值外推法和全覆盖计算法[13,14]。
     直接外推法和插值外推法均是以抽样单元内水土流失状况推算县域水土流失状况。前者分别将1%和 4%密度抽样单元内水土流失状况按照各级土壤侵蚀强度面积百分比推算至其控制区域，并汇总，此法无法直接成图；后者则对轻度及其以上土壤侵蚀强度面积比和水土流失面积比，采用普通克里金插值，经平衡计算后获取县域水土流失状况，采用累积百分比分级法获取县域水土流失状况空间分布示意图。
       全覆盖计算法是基于详细的土地利用矢量数据（根据 SPOT 6 获取），进行 BET 因子赋值，叠加县域 R、K、L、S、BET 因子图层，获取县域水土流失状况及空间分布。其中，B 因子是根据县域林草植被覆盖度结合野外调查数据进行赋值；E、T 因子分别根据野外调查数据采用面积加权平均法计算并赋值。
（3）相对误差分析 采用测量学相对误差的概念，以某一密度或某一推算方法结果为基准，来反映不同密度不同推算方法的差异程度。

    （2）
     （3）
式中：o_i为第 i 级土壤侵蚀强度相对误差（%）；11_i为第 i 级土壤侵蚀强度绝对误差，此文中为两种对应结果之差，在此为方便后续计算取其绝对值；L_i为真值，此文中同种密度下不同推算方法结果比较分析时以全覆盖计算法结果为真值，同种推算方法下不同密度结果比较分析时以 4%抽样密度结果为真值；o扩为平均相对误差（%）；P_i 为第 i 级土壤侵蚀强度面积百分比（无量纲）；i=1,2，……，6，土壤侵蚀强度等级。
二、结果与分析
1、同一抽样密度不同推算方法水土流失特征
（1）1%抽样密度下水土流失数量及空间差异 在 1%抽样密度下，沂水县直接外推法水土流失面积832.93 km2，占 34.42%；插值外推法水土流失面积 811.94 km2，占 33.55%；全覆盖计算法水土流失面积 1053.06 km2，占 43.51%。土壤侵蚀强度均是以微度侵蚀为主，其次为轻度侵蚀。见图 3。
      通过不同推算方法结果对比发现，直接外推法、插值外推法各级土壤侵蚀强度面积和水土流失面积差异较小，水土流失面积仅相差 20.99 km2，占全县 0.87%。但二者与全覆盖计算法的结果存在差异， 水土流失面积分别偏低 220.13 km2 、241.12 km2 ； 轻度及其以上侵蚀相对误差较高为12.18~250.33%。在 1%密度下，以全覆盖计算法水土流失面积为基准，直接外推法和插值外推法相对误差分别为 20.90%和 22.90%，各级土壤侵蚀强度平均相对误差分别为 22.69%和 24.84%。见图 3。
       从空间分布来看，1%抽样密度下插值外推法和全覆盖计算法存在明显差异。插值外推法水土  流失区域集中分布在西部和西北部，而全覆盖计算法分布相对零散。见图 4。







图3  1%抽样密度下三种推算方法各级土壤侵蚀强度面积对比
 
图4  1%密度下插值外推法和全覆盖计算法土壤侵蚀分布图
（2） 4%抽样密度下水土流失数量及空间差异 4%密度下，直接外推法水土流失面积 990.03 km2，占40.91%；插值外推法水土流失面积 998.04 km2，占 41.24%；全覆盖计算法水土流失面积 1074.17km2，占 44.39%。土壤侵蚀强度均是以微度侵蚀为主，其次为轻度侵蚀。见图 5。
       通过对比发现，直接外推法和插值外推法计算结果差异极小，水土流失面积仅相差 8.01 km2， 占 0.33%，各级土壤侵蚀强度面积相差 15 km2 以下。且二者与全覆盖计算法的计算结果差异减小， 水土流失面积分别低 84.14 km2 和 76.13 km2，各占 3.48%和 3.15%，但中度及其以上侵蚀仍有较大差异，相对误差为 17.47~31.71%。4%密度下，以全覆盖计算法水土流失面积为基准，直接外推法和插值外推法相对误差分别为 7.09%和 7.83%，各级土壤侵蚀强度平均相对误差分别为 10.04%和9.66%。见图 5。







图5  4%抽样密度下三种推算方法各级土壤侵蚀强度面积对比
 





图 6 4%密度下插值外推法和全覆盖计算法土壤侵蚀分布图
2、同一推算方法下不同抽样密度水土流失特征与分布
（1）1%和 4%密度下直接外推法结果差异 对比 1%、4%抽样密度，沂水县水土流失面积随抽样密度增大而升高，相对误差 15.87%；随着侵蚀强度的增加，各级土壤侵蚀强度相对误差变大，其中剧烈侵蚀相对误差达 96.00%。直接外推法只能根据抽样单元水土流失状况，大体估算水土流失分布区域。
（2）1%和 4%密度下插值外推法结果差异 对比 1%、4%抽样密度，水土流失面积同样随抽样密度增大而升高，相对误差达 18.65%；随着侵蚀强度的增加，各级土壤侵蚀强度相对误差变大，其中剧烈侵蚀相对误差达达 129.99%。插值外推法可插值形成水土流失分布示意图，但空间分布受抽样单元影响较大。两种密度下水土流失区域空间分布基本一致，但是 4%抽样密度东部和北部范围扩大，且中度侵蚀面积增加。见图 4 和图 6。
（3）1%和 4%密度下全覆盖计算法结果差异 对比 1%、4%抽样密度，水土流失面积差异较小，仅相差 21.11 km2，相对误差 1.97%；从各级土壤侵蚀强度面积来看，强烈及以上侵蚀面积差异较大， 尤其是剧烈侵蚀，由于其面积基数小，相对误差仍达 54.71%。全覆盖计算法可成详细的水土流失图，两种密度下，土壤侵蚀强度和水土流失分布规律相似， 这与全覆盖计算法是基于全县因子图层，以栅格为单元进行计算有关。见图 4 和图 6。
三、讨 论
1、主要受坡度影响，直接推算法和插值推算法结果在两种抽样密度下差异较大。在 CSLE 模型各因子中，4%抽样单元 S 因子均值明显高于 1%抽样密度结果。根据统计数据分析，沂水县山丘区面积百分比 94.33%，平缓坡占 51.80%，而 1%密度平缓坡占 68.06%，4%抽样密度平缓坡占52.45%，4%密度平缓坡比例接近全县比例，并明显低于 1%抽样密度，从而导致了 1%抽样单元的S 因子明显偏高。
2、与邹从荣等[13]结果发现一致，直接外推法和插值外推法结果受抽样密度影响大，全覆盖计算法结果受抽样密度影响小。邹丛荣等[13]研究发现蒙阴县三种推算方法下水土流失面积随抽样密度增加而降低，直接外推法和插值外推法结果受抽样密度影响大，全覆盖计算法结果受抽样密度影响小；而本文研究发现沂水县采用三种推算方法时水土流失面积随抽样密度增加而升高，在 4%密度下，三种推算方法结果数量特征基本一致。蒙阴县和沂水县同属于沂蒙山区，但是二者土地利用方式和石灰岩山区分布有明显不同，蒙阴县土地利用以林地和园地为主，全县内石灰岩山区成条带状分布，沂水县以林地和耕地为主，石灰岩山区在西部成片状分布。此外，蒙阴县 1%密度下抽样单元土壤可蚀性因子、坡度因子和工程措施因子均高于 4%密度抽样单元，而沂水县 1%密度下抽样单元坡度因子低于 4%密度抽样单元，其余因子差异不明显。
四、结论
1、1%抽样密度下，直接外推法与插值外推法计算数据结果相近，但与全覆盖计算法结果有差异；4%抽样密度下，三种推算方法计算数据结果接近，但是水土流失空间分布差异明显。
2、直接外推法与插值外推法受抽样密度影响较大，各级土壤侵蚀强度相对误差一般大于10%，水土流失面积相对误差和平均相对误差均高于 15%；全覆盖计算法受抽样密度影响较小，水土流失面积相对误差和平均相对误差均小于 10%，但各级土壤侵蚀强度级别分布存在一定差异。
3、受坡度影响，1%抽样密度下直接外推法和插值外推法水土流失面积低于 4%抽样密度结果。
4、需要根据工作需要选择合适的抽样密度和推算方法。若需要县域详细的水土流失分布图，可采用 1%抽样密度下的全覆盖计算法；若采用直接外推法或插值外推法，宜采用 4%抽样密度布设调查单元。
参考文献：
1、Jeffery GJ. The national resources inventory and its role in U S agriculture[C]//Agricultural Statistics 2000: Proceedings of the conference on agricultural statistics organized by the National Agricultural Statistics Service of the US Department of Agriculture, under the auspices of the International Statistical Institute, 1998:181-192
2、US Department of Agriculture. Summary report: 2007 national resources inventory [R/OL]. (2009-12-01) [2016-03-27]. http://www.nrcs.usda.gov/Internet/FSE_DOCUMENTS/stelprdb1041379 .pdf.
3、Lu H, Gallant J, Prosser IP, et al. Prediction of sheet and rill erosion over the Australian continent, incorporating monthly soil loss distribution[R/OL]. (2001-01) [2016-03-27]. https://www.researchgate.net/publication/236770265_Prediction_of_Sheet_and