土壤侵蚀已严重威胁世界许多地区的社会经济发展^[1]，对其治理受到广泛关注^[2]，土壤侵蚀的治理需以区域土壤侵蚀调查制图为数据基础。自20世纪90年代以来，国内外学者对区域土壤侵蚀评价进行一系列研究^[2]。Lu等^[3]和Teng等^[4]先后用RUSLE对澳洲土壤侵蚀进行评价，Borrelli等^[5]在全球尺度上进行土壤侵蚀定量分析，刘宝元等^[6]基于抽样调查方法进行全国土壤水力侵蚀普查制图，杨志成等^[7]在抽样调查基础上进行县域和省域土壤侵蚀评价。目前已初步形成对坡面和小流域尺度土壤侵蚀速率的精度评价方法和对区域尺度侵蚀危险性评价方法，但区域尺度上实现土壤侵蚀速率的调查和制图方法有待研究。

珠江是中国第三大河流，珠江三角洲是我国乃至世界最具活力的经济区之一。由于地处南亚热带，降水量大、土层薄，加上人类活动干扰，珠江流域土壤侵蚀愈发严峻^[8-9]。自1990年以来，研究者对珠江流域土壤侵蚀及石漠化等进行一些研究^[9-13]，但目前的研究对流域土壤侵蚀空间分异特征和主控因子的认识不足，与珠江三角洲繁荣的经济不相适应，不能有效支撑珠江流域水土保持规划和决策。为此，迫切需要利用现代科学技术手段，对该区土壤侵蚀进行系统研究。

笔者基于抽样调查单元获取的相关信息和区域土壤侵蚀因子数据，结合CSLE模型^[14]与GIS技术，采用地图代数和空间插值2种方法对研究区土壤侵蚀速率制图并进行土壤侵蚀分析，以期为该地区生态环境治理及水土保持工作实施提供可靠的科学依据，同时促进我国亚热带地区土壤侵蚀的研究。

研究区包括珠江流域及海南省等相邻区(图 1)，总面积约59.14万km²。珠江流域地貌类型复杂，山地、丘陵占总面积94.4%，平原面积小且分散；西江中上游地区以石灰岩山地为主，东江流域以花岗岩丘陵为主，北江中上游地区干流西侧为石灰岩山地、东侧为红岩盆地、下游以花岗岩丘陵为主。流域地处南亚热带，年平均气温在14~22 ℃之间；流域内年平均降水量1 470 mm，自然植被以常绿阔叶林为主，主要土地利用类型为林地、草地和耕地。主要土壤类型有赤红壤、红壤、黄壤和水稻土等。

图 1(Fig. 1)

图 1 研究区位置图 Fig. 1 Location map of the study area

本研究的数据主要包括2部分，一是土壤侵蚀抽样调查数据，研究区(包括陆地部分55 km缓冲区)共布设抽样调查单元1 695个(图 2)，该数据基于高分辨率影像解译而来，主要包含土地利用类型和水土保持措施2部分解译内容，精度优于1∶5 000^[15]；二是用于计算区域土壤侵蚀因子专题数据。数据及其来源详见表 1。

图 2(Fig. 2)

图 2 抽样单元分布图 Fig. 2 Distribution of sampling units

表 1(Tab. 1)

表 1 基础数据一览表 Tab. 1 List of basic data 编号No. 基础数据Basic data 数据来源Sources of the data 分辨率Resolution 1 降雨侵蚀力Rainfall erosivity 项目组提供From project team 1 km 2 土壤可蚀性Soil erodibility 项目组提供From project team 30 m 3 坡度坡长因子Slope steepness and slope length factor 项目组提供From project team 30 m 4 抽样调查单元土地利用和水土保持措施数据Land use and soil conservation measures data in sampling survey unit 项目组提供From project team 矢量数据 Vector data 5 土地利用数据(GLC30)Land use type (http://data.ess.tsinghua.edu.cn/)网站下载Website download 30 m 6 MODIS数据集(percent_Tree_cover、percent_NonTree_vegetate) MODIS dataset (https://ladsweb.modaps.eosdis.nasa.gov/search/)网站下载Website download 250 m 7 研究区耕作区划图Cultivation area map of study area 《区域水土流失动态监测技术规定》(2018年8月版)① “Technical Rules for Dynamic Monitoring of Regional Water and Soil Erosion ”(August 2018 edition) 矢量及表格数据 Vector and tabular data 8 水土保持分区Regional divisions by soil and water conservation 根据赵岩等[16]附图数字化Digitizing based on the figure by Zhao Yan et al[16] 矢量数据 Vector data 注：1~4均系泛第三极土壤侵蚀定量评价与分区防控对策(XDA20040202)项目组的阶段成果数据。Notes: 1-4 are all the stage result data of the project team for “Quantitative Evaluation of Pan-third Pole Soil Erosion and Regional Prevention and Control Countermeasures” (XDA20040202).

表 1 基础数据一览表 Tab. 1 List of basic data

① 水利部水土保持监测中心《区域水土流失动态监测技术规定(试行)》2018年8月The Center of Soil and Water Conservation Monitoring, Ministry of Water Resources. Technical rules for dynamic monitoring of regional water and soil erosion (trial). (2018.8)

笔者使用CSLE模型^[14]计算土壤侵蚀速率。所用参数为全球降雨侵蚀力因子(R因子)、中国土壤可蚀性因子(K因子)、全球L坡长因子与S坡度因子(LS因子)。这3个因子的计算方法详见第一次全国水利普查^[17]和相关文献，本文主要介绍水土保持措施因子的计算。

植被覆盖与生物措施因子(B因子)：基于Borrelli等的^[5]C因子计算方法，根据研究区实际情况改进。对农业用地C因子赋值为1；对非农业用地C因子，结合GLC30土地利用分类，定义相应土地利用C因子值的约束范围。

式中：C_NA为林地、灌木覆盖C因子值；MIN_CF=0.000 1、MAX_CF=0.003；T_C为林地覆盖的比例。

式中：C_p为非林地覆盖(草地覆盖等)C因子值；MIN_C=0.01、MAX_C=0.05；N_VS为非林地(草地等)覆盖的比例。

裸土按照式(2)计算，MIN_C、MAX_C取值为0.1、0.5，其计算结果分辨率为250 m，通过重采样得到30 m分辨率数据。最后基于GLC30土地利用分类合并C_NA、C_p、裸土C和耕地C，得到B因子图层。

工程措施因子(E因子)：利用抽样单元土地利用数据，根据《区域水土流失动态监测技术规定》中计算方法编写程序计算得到抽样单元E因子均值，再对E因子均值的点数据进行反距离权重空间插值得到E因子图层。

耕作措施因子(T因子)：根据《区域水土流失动态监测技术规定》，对抽样单元内耕地的轮作措施属性字段赋值，非耕地赋值1，并基于Borrelli等^[5]对农业中水田的处理方法，将水田赋值为0.15，最终得到T因子图层。

笔者采用2种方法：空间插值法和地图代数法。空间插值法是基于抽样单元的侵蚀速率图^[15]，用地统计协同克里金插值方法，以R、K、B因子作为协变量，对耕、林、草、灌木及裸地分别插值生成不同地类图，并将水体、湿地、不透水面设置为0，最终得到土壤侵蚀速率专题图。地图代数法是根据准备好的各因子图层数据，直接进行各因子相乘得到土壤侵蚀图，考虑到LS、K因子分辨率均为30 m，本研究的运算将分辨率设置为30 m。

基于地理探测器进行土壤侵蚀主控因子分析。地理探测器是一种通过探测要素的空间分层异质性及揭示背后驱动力的一种统计学方法^[18]，其大小由q来衡量，q表达式如下：

式中：h=1，2…；L为自变量X的分层数；N_h和N分别为对应层为h和区域内单元数；δ_h²和δ分别为对应层h和全部区域Y值的方差；SSW为层内方差之和；SST为区域总方差；最终计算q的值域范围为[0, 1]，q值越大，自变量X对土壤侵蚀速率Y的分异性解释力越强。本研究以R、K、LS、B、E、T、植被覆盖度、坡度、土地利用作为影响因子X。根据研究区范围及软件可处理容量，将研究区离散化为5 km×5 km的网格，采样其中间点(共2万9 050个)作为地理探测器的输入数据。

影响土壤侵蚀速率各因子中(图 3)，R因子总体上呈从东至西递减的规律，其值域为1 786.81~16 324.18 MJ·mm/(hm²·h·a)，平均值为6 813.56 MJ·mm/(hm²·h·a)。K因子值域为0~0.012 5 t·hm²·h/(MJ·mm·hm²)，均值为0.004 8 t·hm²·h/(MJ·mm·hm²)，有从东到西递增趋势。LS因子值域为0.04~21.24，平均值为8.70，高值见于云贵高原、东部山地和海南五指山等地。R、LS因子值普遍较高，是土壤侵蚀的主要诱发因子。B因子、E因子和T因子，值域分布范围依次为0~1、0.15~1、0.15~1，其均值依次为0.21、0.91、0.87。B因子值较低，而E和T因子值较高，研究区植被覆盖度较好，是土壤侵蚀主要抑制因子。

图 3(Fig. 3)

图 3 CSLE各因子空间分布图 Fig. 3 Spatial distribution of CSLE factors

在水利部批准SL190—2007《土壤侵蚀分类分级标准》^[19]的基础上进行细分，编制了土壤侵蚀分级图(图 4和图 5)，研究区土壤侵蚀均主要集中分布在研究内的贵州省及云南省、广西中部和广东省沿海区，其中强烈和极强烈侵蚀区域范围较小且分布在较零散的坡耕地上。2张图上土壤侵蚀速率较高的地区与土地利用图的耕地分布较吻合，而林、草地的侵蚀则相对较弱。从微观角度看，2图在细节纹理结构上略有不同。基于地图代数法的结果局地变化比较明显，空间插值结果整体较为平滑，主要原因在于插值所用的侵蚀速率为抽样单元均值，且插值过程中有平滑作用。

图 4(Fig. 4)

图 4 地图代数法土壤侵蚀速率图 Fig. 4 Soil erosion rate map based on map algebra method

图 5(Fig. 5)

图 5 空间插值法土壤侵蚀速率图 Fig. 5 Soil erosion rate map based on spatial interpolation method

采用地图代数和空间插值法制图的土壤侵蚀速率均值分别为791.78 t/(km²·a)和615.37 t/(km²·a)，两者值较接近。地图代数法计算的土壤侵蚀速率均值较高，可能与水保措施难以全面纳入计算有关。笔者计算的土壤侵蚀速率值与前人研究^[20-22]对比(表 2)，比较接近，可以认为本研究的结果较为可信。对比研究区2种方法的主要土地利用类型的土壤侵蚀速率结果(表 3)发现，土壤侵蚀速率值较为接近。空间插值法制图结果，耕地、林地的土壤侵蚀速率均值略小，这是由于地图代数法是从整个研究区的计算，分辨率相对较粗，难以获取详细的水土保持措施，空间插值法从抽样单元小流域可计算得到准确的水土保持措施值，计算的土壤侵蚀速率值较小。

表 2(Tab. 2)

表 2 土壤侵蚀速率对比参照表 Tab. 2 Soil erosion rates of this study and extracted from references   t/(km2·a) 编号 No. 研究区域 Study area 分析方法 Analysis method 分辨率 Resolution 侵蚀速率 Erosion rate 本研究侵蚀速率 Erosion rate in this study 地图代数法 Map algebra 空间插值 Spatial interpolation 1 珠江流域 Pearl River Basin 泥沙输[19] Sediment transport ratio — 498.56 736.86 655.04 2 贵州省 Guizhou province RUSLE[20] 50 m 880.81 784.02 813.44 3 广西西北部岩洞区 Cave area of northweast Guangxi RUSLE[21] 30 m 716.6 763.97 683.80 注：—表示无数据。Notes: — refers to no data.

表 2 土壤侵蚀速率对比参照表 Tab. 2 Soil erosion rates of this study and extracted from references

表 3(Tab. 3)

表 3 主要土地利用土壤侵蚀速率对比表 Tab. 3 Average soil erosion rates of land use types   t/(km2·a) 计算方法 Calculation method 耕地 Cropland 林地 Forest 草地 Grassland 地图代数法Map algebra 2 767.11 166.16 1 235.63 空间插值法Spatial interpolation 2 090.35 79.52 1 620.54

表 3 主要土地利用土壤侵蚀速率对比表 Tab. 3 Average soil erosion rates of land use types

根据SL190—2007《土壤侵蚀分类分级标准》^[19]对土壤侵蚀强度进行统计(表 4)，2种方法计算的水土流失面积(轻度及其以上侵蚀)差值为10.17%，空间插值法略大；土壤侵蚀均主要集中于微、轻度侵蚀；但其在微、轻度侵蚀频率分布略有差异，这可能与空间插值时耕地土壤侵蚀速率值较集中有关；在中度及以上的侵蚀强度上，土壤侵蚀速率结果频率相差不大。

表 4(Tab. 4)

表 4 空间插值法与地图代数法土壤侵蚀强度统计 Tab. 4 Soil erosion intensity statistics by spatial interpolation and map algebra   % 侵蚀强度 Soil erosion intensity 空间插值法 Spatial interpolation 地图代数法 Map algebra 二者差值 Difference value 微度侵蚀Micro erosion 71.87 82.04 -10.17 轻度侵蚀Mild erosion 20.95 11.34 9.61 中度侵蚀Moderate erosion 5.91 2.99 2.92 强度侵蚀Serious erosion 0.98 1.33 -0.35 极强度侵蚀Extreme erosion 0.29 1.26 -0.97 剧烈侵蚀Severe erosion 0.00 1.04 -1.04

表 4 空间插值法与地图代数法土壤侵蚀强度统计 Tab. 4 Soil erosion intensity statistics by spatial interpolation and map algebra

总体上看，2种方法计算的土壤侵蚀速率结果均可表达研究区土壤侵蚀状况，考虑到土壤侵蚀防控的需要，优先选用空间插值法。

基于地理探测器的分析结果(表 5)表明，整个研究区和水土保持区划的各个区，各因子影响程度从大到小的排序基本为土地利用、水保措施(特别是B因子)、地形、土壤和降水。除浙闵山地丘陵区外，土地利用方式为土壤侵蚀的主控因子，影响因子q值均达到47%以上；其次，B因子对土壤侵蚀速率及其空间分布的影响较大；T因子的影响排第三，华南沿海丘陵台地区和海南及南海诸岛丘陵区的q值达到40%以上，可能与该区复杂的耕作方式有关，其原因有待进一步探究分析。K因子和R因子的影响较小，这与研究区以林草植被为主、植被覆盖普遍较高有关。

表 5(Tab. 5)

表 5 各影响因子q值统计表 Tab. 5 Statistical table of q values of each impact factor 水土保持分区 Regional division by soil and water conservation R K LS B E T FVC SLP 土地利用类型 Landuse type 滇黔桂山地丘陵区Hilly and mountainous area of Yunnan-Guizhou-Guangxi 0.004 0.008 0.109 0.458 0.002 0.260 0.062 0.091 0.570 南岭山地丘陵区Hilly and mountainous area of Nanling 0.002 0.036 0.208 0.483 0.001 0.172 0.137 0.209 0.579 华南沿海丘陵台地区Hilly platform area of South China coastal 0.028 0.021 0.096 0.355 0.011 0.426 0.110 0.101 0.532 浙闵山地丘陵区Hilly and mountainous area of Zhejiang-Fujian 0.006 0.028 0.113 0.549 0.011 0.319 0.086 0.118 0.383 海南及南海诸岛丘陵区Hilly areas of Hainan and South China Sea islands 0.048 0.039 0.161 0.515 0.010 0.404 0.235 0.152 0.565 研究区Study area 0.018 0.023 0.082 0.370 0.002 0.206 0.069 0.075 0.470 注：SLP为坡度。Notes: SLP trefer to slope.

表 5 各影响因子q值统计表 Tab. 5 Statistical table of q values of each impact factor

1) 研究区土壤侵蚀特征环境是降雨侵蚀力较大、坡度较陡及植被覆盖度较高，前2项是土壤侵蚀的主要诱发因子，后者则是土壤侵蚀的主要抑制因子。

2) 本文通过地图代数和空间插值2种方法制图得到土壤侵蚀速率均值分别为791.78 t/(km²·a)、615.37 t/(km²·a)，研究区土壤侵蚀广泛分布，坡耕地是土壤侵蚀主要的发生部位。

3) 土地利用类型为影响土壤侵蚀速率的主控因素，B因子对土壤侵蚀速率及其空间分布影响较大，总体上讲，由于研究区植被覆盖度普遍较高，因而土壤和降雨的影响较小。调整土地利用结构、不断优化植被的水土保持生态功能是值得注意的问题。

4) 地图代数和空间插值两种方法制图结果相比，空间插值法的结果更接近土壤侵蚀的实际情况，可作为区域土壤侵蚀制图首选方法，但两者之间的差异和相互结合有待讨论。

中国科学院水利部水土保持研究所张晓萍研究员、水利部水土保持监测中心王爱娟高工、珠江水利委员会水土保持监测中心金平伟主任对本研究给予多方支持，一并致谢。