基于CSLE模型的珠江流域土壤侵蚀强度评价

陈羽璇， 杨勤科,2†， 刘宝元， 黄晨璐， 王春梅,2， 庞国伟,2

(1.西北大学城市与环境学院，710127，西安；2.旱区生态水文与灾害防治国家林业局重点实验室，710127，西安；3.中国科学院 水利部 水土保持研究所，712100，陕西杨凌；4.北京师范大学地理学院，100875，北京)

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

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

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

1 研究区概况

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

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

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

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

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

2 研究方法

2.1 土壤侵蚀因子计算

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

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

CNA=MINCF+(MAXCF-MINCF)(1-TC)

(1)

式中：CNA为林地、灌木覆盖C因子值；MINCF=0.000 1、MAXCF=0.003；TC为林地覆盖的比例。

CP=MINC+(MAXC-MINC)·(1-NVS)

(2)

式中：Cp为非林地覆盖(草地覆盖等)C因子值；MINC=0.01、MAXC=0.05；NVS为非林地(草地等)覆盖的比例。

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

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

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

2.2 流域土壤侵蚀速率图制作

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

2.3 主控性因子分析方法

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

(3)

(4)

3 结果与分析

3.1 土壤侵蚀影响因子分布特征

影响土壤侵蚀速率各因子中(图3)，R因子总体上呈从东至西递减的规律，其值域为1 786.81～16 324.18 MJ·mm/(hm2·h·a)，平均值为6 813.56 MJ·mm/(hm2·h·a)。K因子值域为0～0.012 5 t·hm2·h/(MJ·mm·hm2)，均值为0.004 8 t·hm2·h/(MJ·mm·hm2)，有从东到西递增趋势。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 CSLE各因子空间分布图Fig.3 Spatial distribution of CSLE factors

3.2 2种计算方法下土壤侵蚀的对比分析

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

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

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

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

表2 土壤侵蚀速率对比参照表Tab.2 Soil erosion rates of this study and extracted from references t/(km2·a)

表3 主要土地利用土壤侵蚀速率对比表Tab.3 Average soil erosion rates of land use types t/(km2·a)

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

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

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

3.3 土壤侵蚀主控因子分析

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

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

4 结论

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

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

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

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

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