基于遥感与抽样单元调查的县域尺度水土流失推算方法

2021-01-09罗梦琦张春强吴镇宇姚孝友

山东农业大学学报（自然科学版） 2020年6期

段倩,齐斐,罗梦琦,刘霞*,唐俊,张春强,吴镇宇,李想,姚孝友

1.南京林业大学南方现代林业协同创新中心,江苏省水土保持与生态修复重点实验室,江苏南京210037

2.淮河水利委员会淮河流域水土保持监测中心站,安徽蚌埠233001

水土流失是全球普遍关注的生态环境问题。作为水土流失最严重的国家之一，我国已为治理水土流失投入了大量财力、物力和人力。为评价水土流失治理效果，开展水土流失动态监测具有重要的意义。目前，在国外水土流失调查计算过程中，美国主要采用抽样调查和通用土壤流失方程（USLE），澳大利亚主要采用基于GIS 和USLE 的水土流失网格估算，欧洲主要采用基于专家法和USLE 模型法[1-4]等。在国内，水土流失状况评价方法主要有定性评价法（综合评判法）以及中国土壤流失方程（Chinese Soil Loss Equation,CSLE）模型法。前者在我国80 年代、90 年代的水土流失遥感调查中广泛应用[5]；后者首次与抽样调查相结合大规模应用于2010～2012 年开展的第一次全国水利普查水土保持情况普查中[6-9]，并在2018 年开始的全国水土流失动态监测用于水力侵蚀计算。在CSLE 模型应用中发现，抽样密度会影响小尺度尤其是县域尺度的水土流失状况推算结果，但是影响大小具有地域性差异；而且不同的水土流失推算方法也会影响最终的结果。赵维军、张岩等[10-12]以陕西吴起县为例，对比分析不同密度对水土流失因子精度影响，发现1%密度和4%密度均能很好地代表吴起县水土流失状况，并在1%抽样密度基础上对比CSLE 模型和综合评判法估算成果，认为前者更有优越性，但是二者以抽样单元内分析为主，均未涉及到不同抽样密度、不同推算方法对县域水土流失状况估算的影响；邹丛荣等[13]在山东蒙阴县研究1%、4%抽样密度与三种水土流失推算方法对水土流失状况的影响，发现4%抽样密度适用于直接外推法和插值外推法，采用1%抽样密度适用于全覆盖计算法（栅格计算法）。故本研究以沂蒙山泰山国家级重点治理区沂水县为例，分别探索不同抽样密度、不同推算方法对县域水土流失状况的影响，对比分析差异和原因，进一步探讨适用于县域尺度水土流失调查与评价的抽样密度和推算方法，为县域水土流失动态监测方法提供数据支撑。

1 材料与方法

1.1 研究区概况

沂水县地处沂山南麓，隶属山东省临沂市，面积2420.08 km2，属北方土石山区。地貌主要以低山丘陵为主，地势自西北向东南倾斜，山丘区面积比94.3%。土壤包括粗骨土、棕壤、褐土、红粘土和潮土五大类，以粗骨土为主。属暖温带季风型大陆性气候，多年平均气温12.3 ℃，多年年均降水量629 mm。水资源丰富，境内沂河为山东省第一大河，以及山东省第三大水库跋山水库。植被属暖温带落叶阔叶林区，主要乔木树种有刺槐（Robinia pseudoacaciaL.）、侧柏（Platycladus orientalis（Linn.）Franco）、花椒（Zanthoxylum bungeanumMaxim.）、油松（Pinus tabulaeformisCarr.）等，自然灌木与草本植物主要有黄荆（Vitex negundoL.）、三裂绣线菊（Spiraea trilobataL.）、胡枝子（Lespedeza bicolorsTurcz.）等。地理位置见图1。

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

1.2.1 野外抽样单元布设采用分层不等概抽样法[7,9,14,15]，以国家第一次水利普查布设的野外抽样单元为基础，在山丘区增加抽样单元，勾绘0.2～3 km2小流域，使其抽样密度增加至4%，平原区维持不变。全县1%抽样密度调查单元24 个，单元控制区域100 km2；4%调查单元96 个，单元控制区域25 km2。野外调查点分布见图2。

图1 沂水县地理位置Fig.1 The location of Yishui County

图2 沂水县1%和4%抽样单元分布Fig.2 The location of field units in Yishui County

1.2.2 野外抽样单元数据采集与处理基于野外调查底图和现场实际状况进行地块边界勾绘，填写野外调查信息表，详细记录单元内各地块土地利用、水土保持措施等信息[7]。调查结束后，将野外调查地图配准矢量化，逐地块填写地块信息。

1.3 水土流失推算方法

1.3.1 水土流失计算模型本研究采用中国水土流失通用模型（CSLE）[11]，根据《北方土石山区水土流失综合治理技术标准》（SL 665-2014）进行水土流失状况统计分析。CSLE 模型基本形式为：

式中：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值，并基于全县土壤类型图，获取全县土壤可蚀性因子图层。

坡长、坡度因子（L、S）：基于1:1 万地形图，经投影变换后，生成10 m 栅格大小DEM，采用刘宝元[19]修正算法，提取坡度坡长因子。

水土保持措施因子（BET）：基于径流小区数据和水利普查中提供的措施因子参考值，对地块进行赋值[14,16]。

1.3.2 水土流失推算方法主要采用直接外推法、插值外推法和全覆盖计算法[13,14]。

直接外推法和插值外推法均是以抽样单元内水土流失状况推算县域水土流失状况。前者分别将1%和4%密度抽样单元内水土流失状况按照各级土壤侵蚀强度面积百分比推算至其控制区域，并汇总，此法无法直接成图；后者则对轻度及其以上土壤侵蚀强度面积比和水土流失面积比，采用普通克里金插值，经平衡计算后获取县域水土流失状况，采用累积百分比分级法获取县域水土流失状况空间分布示意图。

全覆盖计算法是基于详细的土地利用矢量数据（根据SPOT 6 获取），进行BET 因子赋值，叠加县域R、K、L、S、BET 因子图层，获取县域水土流失状况及空间分布。其中，B因子是根据县域林草植被覆盖度结合野外调查数据进行赋值；E、T因子分别根据野外调查数据采用面积加权平均法计算并赋值。

1.3.3 相对误差分析采用测量学相对误差的概念，以某一密度或某一推算方法结果为基准，来反映不同密度不同推算方法的差异程度。

式中：δi为第i级土壤侵蚀强度相对误差（%）；Δi为第i级土壤侵蚀强度绝对误差，此文中为两种对应结果之差，在此为方便后续计算取其绝对值；Li为真值，此文中同种密度下不同推算方法结果比较分析时以全覆盖计算法结果为真值，同种推算方法下不同密度结果比较分析时以4%抽样密度结果为真值；δ¯为平均相对误差（%）；Pi为第i级土壤侵蚀强度面积百分比（无量纲）；i=1,2,…,6，土壤侵蚀强度等级。

2 结果与分析

2.1 同一抽样密度不同推算方法水土流失特征

2.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%抽样密度下三种推算方法各级土壤侵蚀强度面积对比Fig.3 The results of three estimation methods based on 1%sampling ratio

图4 1%密度下插值外推法和全覆盖计算法土壤侵蚀分布图Fig.4 Soil erosion intensity by different calculation methods based on 1%sampling density

2.1.2 4%抽样密度下水土流失数量及空间差异4%密度下，直接外推法水土流失面积990.03 km2，占40.91%；插值外推法水土流失面积998.04 km2，占41.24%；全覆盖计算法水土流失面积1074.17 km2，占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%抽样密度下三种推算方法各级土壤侵蚀强度面积对比Fig.5 The results of three estimation methods based on 4%sampling ratio in Yishui

从空间分布来看，4%抽样密度下，插值外推法和全覆盖计算法水土流失空间分布差异较大，前者水土流失区域片状分布，中度及其以上侵蚀主要分布在西部，东部各乡镇几乎均为微度侵蚀，而后者呈均匀异质性零散分布，在全县各个乡镇广泛存在。见图6。

图6 4%密度下插值外推法和全覆盖计算法土壤侵蚀分布图Fig.6 Soil erosion intensity by different calculation methods of 4%sampling density

2.2 同一推算方法下不同抽样密度水土流失特征与分布

2.2.1 1%和4%密度下直接外推法结果差异对比1%、4%抽样密度，沂水县水土流失面积随抽样密度增大而升高，相对误差15.87%；随着侵蚀强度的增加，各级土壤侵蚀强度相对误差变大，其中剧烈侵蚀相对误差达96.00%。直接外推法只能根据抽样单元水土流失状况，大体估算水土流失分布区域。

2.2.2 1%和4%密度下插值外推法结果差异对比1%、4%抽样密度，水土流失面积同样随抽样密度增大而升高，相对误差达18.65%；随着侵蚀强度的增加，各级土壤侵蚀强度相对误差变大，其中剧烈侵蚀相对误差达达129.99%。

插值外推法可插值形成水土流失分布示意图，但空间分布受抽样单元影响较大。两种密度下水土流失区域空间分布基本一致，但是4%抽样密度东部和北部范围扩大，且中度侵蚀面积增加。见图4 和图6。

2.2.3 1%和4%密度下全覆盖计算法结果差异对比1%、4%抽样密度，水土流失面积差异较小，仅相差21.11 km2，相对误差1.97%；从各级土壤侵蚀强度面积来看，强烈及以上侵蚀面积差异较大，尤其是剧烈侵蚀，由于其面积基数小，相对误差仍达54.71%。

全覆盖计算法可成详细的水土流失图，两种密度下，土壤侵蚀强度和水土流失分布规律相似，这与全覆盖计算法是基于全县因子图层，以栅格为单元进行计算有关。见图4 和图6。

3 讨论

（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%密度抽样单元，其余因子差异不明显。