基于CSLE和抽样调查的泰国土壤水蚀评价

杜朝正, 杨勤科, 张加琼, 马 波, 王春梅,4, 庞国伟,4

(1.西北大学 城市与环境学院, 西安 710127; 2.菏泽学院 城市建设学院, 山东 菏泽 274015; 3.中国科学院水利部 水土保持研究所, 陕西 杨陵 712100; 4.旱区生态水文与灾害防治国家林业和草原局重点实验室, 西安 710127)

土壤是人类生存与发展的重要自然资源之一[1]，土壤可以为人类提供食物、能源、水源，在维持生态系统服务、生物多样性方面也具有重要作用[2]。但是，全球大部分地区的土壤受到水蚀的威胁[3]。土壤水蚀是粮食生产和可持续农业面临的最严重威胁之一[4-5]。土壤水蚀不仅会导致土地退化，生态服务水平降低，还会导致下游河流淤积，从而制约区域的可持续发展[6]。

治理土壤水蚀必须以掌握区域土壤水蚀状况为前提[7-8]，一般可归纳为2种方法：第1种是基于区域各侵蚀因子图层相乘，计算区域土壤水蚀速率专题图，可称为地图代数法。该方法计算结果具有较好的空间分布，得到广泛应用，代表性研究有Borrelli[9]、Panagos[10]、查良松[11]等。第2种是抽样调查的方法，在研究区布设基本抽样调查单元，结合土壤水蚀模型，计算各抽样单元水蚀速率，进而推算出研究区土壤水蚀状况，可称为抽样调查法，代表性研究有Nusser[12]、USDA[13]、刘宝元等[7]、Liu等[14]、Yin等[15]。但如何将地图代数制图结果与抽样调查结果相结合，实现对区域土壤侵蚀的定量评价则鲜见报道。

泰国属于发展中国家，农作区占全国的40%以上，大力发展农业导致农地不断扩张，而森林不断萎缩，加之粗放的耕作和管理方式，土壤侵蚀日益加剧，造成农业地区土地严重退化[16]。在农作区，农民习惯将坡面裸露，这使土壤完全暴露于暴雨及径流中，土壤流失现象普遍存在[17]。在山区和高原，大规模砍伐森林、陡坡耕种(45°～60°)加之游移耕作，进一步加剧了土壤侵蚀[18-19]。目前，虽然已有学者对泰国进行土壤水蚀评价，如Rangsiwanichpong等[16]，但数据分辨率较粗，地形因子和水保措施因子计算有待改进。需进一步提高泰国土壤水蚀定量评价的准确性，为泰国土壤水蚀治理提供更为准确的数据支撑。

本研究以泰国为例，在前期工作的基础上[20]，采用中国土壤流失方程(Chinese Soil Loss Equation，CSLE)模型，基于各侵蚀因子基础数据和抽样调查数据，定量计算泰国土壤水蚀速率专题图，以期全面准确掌握泰国土壤水蚀基本状况，为该国土壤水蚀治理提供数据支撑。

1 研究区域与基础数据

泰国(5°30′—21°N，97°30′—105°30′E)位于东南亚的中南半岛中部，总面积513 115 km2。泰国地势西北高东南低。根据地形特征，可将泰国分为北部(山地为主)、东北部(高原为主)、中部(平原为主)、南部(丘陵地带)4个自然区域带。泰国地处热带地区，主要受热带季风气候影响，年均降水量约1 150 mm。由于降雨量大，盲目开垦耕地，并且存在不合理的耕作方式，水土流失问题较为严重[21]。泰国是重要的农产品出口国，土壤侵蚀问题已严重制约该国的可持续农业发展[22-23]。

本研究所用基础数据主要包括2个部分，一是用于地图代数法制图整体计算的各水蚀因子数据(表1)；二是土壤水蚀抽样调查数据，研究区共计175个抽样调查单元(图1)，所有抽样调查单元数据均来自野外调查解译并通过了质量检验[8,24]。

表1 基础数据一览表

图1 抽样单元分布图

2 研究方法

本研究首先采用地图代数法计算研究区土壤水蚀速率专题图，再以175个抽样调查单元计算结果为参照，对地图代数结果进行直方图匹配，将匹配之后的制图结果作为研究区最终的土壤水蚀速率专题图。

2.1 土壤水蚀评价模型

采用中国土壤流失方程(Chinese Soil Loss Equation，CSLE)，结合GIS技术，计算研究区土壤水蚀速率。CSLE基本形式为：

A=R×K×L×S×B×E×T

(1)

式中：A为土壤水蚀速率[t/(hm2·a)]；R降雨侵蚀力因子[MJ·mm/(hm2·a·h)]；K为土壤可蚀性因子[t·hm2·h/(hm2·MJ·mm)]；L为坡长因子，无量纲；S为坡度因子，无量纲；B为植被覆盖与生物措施因子，无量纲；E为工程措施因子，无量纲；T为耕作措施因子，无量纲。

2.2 水蚀因子数据处理

(1)R，K，LS因子数据处理：根据项目组提供的全球降雨侵蚀力因子R、泛第三极土壤可蚀性因子K，通过裁剪、重采样处理，获得研究区30 m分辨率的R，K因子栅格数据。将30 m分辨率的多幅分块LS因子基础数据进行拼接，并裁剪出研究区范围。参考中国水利部的关于水土流失动态监测的有关技术文件和Brychta等[25]的研究，计算LS因子上限阈值(坡度30°、坡长100 m时，LS上限=21.2)，最终获得研究区LS因子专题层数据。

(2)B因子计算：植被覆盖与生物措施是一种重要的水保措施[26]，参考Borrelli等[9]的研究方法，计算非农业土地类型B因子值。

下载MOD44 B数据，并通过MRT(MODIS Reprojection Tool，MRT)工具提取TC(Tree Cover，TC)、NTV(NonTree Vegetation，NTV)图层。计算各土地利用类型B因子值采用如下公式。

林地：

BNA=MINB+(MAXB-MINB)×(1-TC)

(2)

非林地：

BP=MINB+(MAXB-MINB)×(1-NTV)

(3)

式中:BNA为林地B因子值；TC为林地覆盖度,取值为0～1；BP为非林地(即:草地、灌木、裸地等)的B因子值,NTV为植被(除林地外)覆盖度,取值为0～1。maxB,minB为对应土地利用类型B因子值的最大值和最小值,取值如表2所示。

根据土地利用类型数据，合并各地类B因子值，最终获得研究区B因子图层。

表2 各土地利用类型B因子取值范围

(3) E因子计算：工程措施是用人工或机械挖掘和修筑完成[26]，在前期工作的基础上，依据抽样单元布设原则，将研究区划分为规则矢量网格，确保每个矢量网格内只包括1个抽样调查单元，将抽样单元的E因子平均值作为矢量网格的E因子值，再将矢量网格数据转换为栅格数据，最终获得研究区的E因子专题图。

(4) T因子计算：根据水利部水土保持监测中心编写的《水土流失动态监测技术指南》(2020年8月版)，参考中国各农业熟制区T值表的规则，为Fischer等[27]制作的全球熟制区赋予T因子值。再依据研究区土地利用类型数据，将非耕地区域的T值设置为1，将耕地区域的T值取值为全球T值图的T值，最终获得研究区的T因子专题图。

2.3 土壤水蚀制图方法

2.3.1 地图代数法制图 地图代数法制图是将事先准备好的每个土壤侵蚀因子栅格数据相乘，获得研究区的土壤水蚀速率专题图。为了确保计算精度及便于显示，本研究采用WGS84坐标系，最终获得的区域R，K，LS，B，E，T，A专题数据分辨率均设定为30m。

2.3.2 抽样单元水蚀速率计算 研究区内共解译了175个均匀分布的抽样调查单元，根据解译单元点位置的地形特征，确定调查单元的形状为小流域(在山区)或矩形(在平原)，面积大小为0.2～3km2。解译每个抽样单元内的土地利用和详细水土保持措施，其中水土保持措施包括生物措施、工程措施、耕作措施。在研究区范围内，CSLE模型适用于计算抽样单元水蚀速率[28]。

对区域R因子、区域K因子、区域LS因子进行多重采样和裁剪，获得每个抽样单元R，K，LS因子专题数据。对各土地利用类型的区域B值图(250m分辨率)进行多重采样和裁剪，再根据解译的土地利用类型和生物措施，计算每个抽样单元的B因子专题数据。根据解译的工程措施，为每个图斑赋E因子值，再将矢量数据转换为栅格数据，获得每个抽样单元的E因子专题数据。根据解译的土地利用类型和水保措施信息，计算每个抽样单元的T因子专题数据。再采用CSLE模型计算每个抽样单元的土壤水蚀速率(针对抽样调查单元的每个图斑)，借助python语言实现数据批量处理，最终获得175个以各抽样单元为边界的土壤水蚀速率专题图(2.5m分辨率)。

2.3.3 直方图匹配制图 以均匀分布的175个抽样调查单元为基础，抽样单元计算结果可全面表达水土保持措施的影响，可真实反映研究区的统计特征；地图代数法制图结果具有较好的空间分布特征。为了结合两者的优点，尝试以175个抽样单元水蚀速率专题图结果(针对每个抽样单元的每个栅格值统计)的累积频率为参照，对地图代数制图结果进行累积频率直方图匹配，最终获得研究区新的土壤水蚀速率专题图。

3 结果与分析

3.1 土壤水蚀因子分析

按照上述方法获得的研究区各水蚀因子专题层如图2所示。研究区R因子空间分布上呈现出由东南向西北逐渐减弱的趋势，主要原因为泰国属于热带季风气候，降雨量由东南沿海向西北内陆逐渐减小。K因子空间分布呈现出与土地利用类型相关的空间格局，K因子值高的区域与耕地分布区域相对应，K因子值低的区域主要对应林地区域。LS因子呈现由西北向东南逐渐减小的分布格局，主要原因为泰国西北部为山地，地形起伏较大，坡度较陡，中部为平原区，地形较为平缓。B因子值的空间分布与土地利用类型的分布有关，按照林地、草地、耕地依次增大。E因子的分布与国家政策和当地耕作习惯有关，在泰国范围内工程措施较少。T因子只针对农地地块，T因子值的大小与轮作方式、管理措施相关，轮作方式、管理措施又存在区域间差异，在南部丘陵地带T因子值较大，在东北部T因子值较小。

3.2 泰国水蚀速率特征

3.2.1 空间分布特征分析 泰国水蚀速率专题图显示水蚀速率高值区主要集中于泰国北部和东北部，水蚀速率低值区主要分布于中部平原区(图3)。这与泰国的地形和耕地分布密切相关，泰国北部属于山地，地形起伏较大，降雨容易引起土壤流失；东北部地区存在大量耕地区域，耕地区域B因子值较大，并且该区域的K，R因子值也较大。中部水蚀速率较小，主要原因为中部地形以平原为主，并且R因子相对较小。

3.2.2 统计特征分析 参考Bahadur等[29]、Borrelli等[9]的分级方案，结合泰国区域的实际情况，对泰国土壤侵蚀速率的侵蚀等级划分结果显示：泰国侵蚀速率以<500t/(km2·a)为主(74.2%)，随着侵蚀等级升高，侵蚀面积比率呈减少趋势。土壤侵蚀速率>15 000t/(km2·a)的面积仅占研究区总面积的0.6%，但年侵蚀量约占研究区侵蚀总量的21.5%，可见研究区局部侵蚀剧烈。泰国平均水蚀速率为687.9t/(km2·a)，是全球平均土壤水蚀速率[287.2t/(km2·a)]的2.4倍，个别地区达到1 000t/(km2·a)以上(占面积13.2%，占侵蚀总量72.0%)(表3)。可见，与全球平均水蚀速率相对比，研究区土壤水蚀较为严重。

3.2.3 热点地区分析 参考Borrelli等[9]的研究，结合本研究分级方案，根据水蚀治理的紧迫程度，将水蚀速率以500，2 500t/(km2·a)为阈值设置为热点地区。>500t/(km2·a)的热点地区沿泰国周边分布，该热点地区地形LS因子值(均值为7.5)是泰国区域LS因子值(均值为3.2)的2.3倍，降雨侵蚀力也较大[均值为6 419.2MJ·mm/(hm2·a·h)]。>2 500t/(km2·a)的热点地区主要分布于泰国北部和东北部，如此分布主要是地形LS因子、植被覆盖与生物措施B因子、土地利用方式等因素共同作用的结果。该热点地区LS因子值(均值为7.3)、B因子值(均值为0.849)是泰国区域的LS因子值(均值为3.2)、B因子值(均值为0.507)的2.3倍和1.7倍，降雨侵蚀力也较大[均值为6 587.6MJ·mm/(hm2·a·h)]，这导致该热点地区水蚀较为严重(图4)。对该热点地区进一步分析发现，热点地区84.1%为耕地，平均坡度为13.4°，因此，坡耕地进一步加剧了局部土壤水蚀。>500t/(km2·a)的热点地区面积占研究区面积的25.8%，年侵蚀量占研究区侵蚀总量的84.3%，该热点地区需制定长期治理规划，逐步减少土壤水蚀的影响；>2 500t/(km2·a)的热点地区面积占研究区总面积的4.7%，但年侵蚀量占研究区侵蚀总量的53.1%，该热点地区迫切需要尽快治理，并且治理效果显著。

注：耕地区域的B因子值为1；无工程措施区域的E因子值为1；非耕地区域的T因子值为1。

3.3 各土地利用类型的土壤水蚀分析

通过比较不同土地利用类型的土壤水蚀状况(表4)，我们发现从耕地到草地、林地的平均水蚀速率明显降低。耕地平均水蚀速率是草地、灌丛、林地平均水蚀速率的2.0倍、2.2倍、2.8倍，说明相对于耕地，林草地具有更好的水土保持作用。在泰国区域，耕地面积占研究区面积的比率最大(49.6%)，平均水蚀速率最高[1 020.2t/(km2·a)]，进而耕地侵蚀量占研究区总侵蚀量的比率最大(73.6%)，这是研究区平均水蚀速率较高的重要原因。

图3 泰国土壤水蚀速率

表3 不同侵蚀等级的土壤水蚀面积及水蚀量

进一步统计后我们发现，在耕地区域中，>5 000t/(km2·a)的区域占耕地总面积的比率高达4.6%。草地、灌丛、林地在>5 000t/(km2·a)的区域占该地类面积的1.6%，0.6%，0.1%，明显比耕地小。林地、灌丛、草地具有较好的储水固土能力，因此，在有条件的区域，将高侵蚀级别的耕地改变为林地、灌丛、草地，对于减少水土流失具有重要意义。

3.4 主控因子分析

分类决策树状图显示了水蚀因子的重要性，终端节点映射出每个像元值归因于分类决策树中最频繁的节点类(图5—7)。结合因子节点可视化图以及泰国土地利用类型图可知，节点11，12，14，15主要受R因子(降雨侵蚀力因子)影响，主要分布于林地区域，在该区域由于降雨量较大而容易发生侵蚀。节点7，8主要受K因子(土壤可蚀性因子)影响，主要分布于东北部高原区的耕作区域，该区域常年耕作(包括一年多熟)，并且不合理的耕作方式和管理措施而引起土壤可蚀性较强，进而容易引发土壤侵蚀。节点4，5主要受LS因子(坡度坡长因子)影响，主要分布于中部平原区的耕作区和南部丘陵区的耕作区，在中部平原区，耕地的坡长较长；在南部丘陵区的地面起伏较大，容易引发土壤侵蚀。

4 讨 论

4.1 计算结果合理性评估

本研究在UpperNamWa流域、NamChun流域、KhlongYai流域计算的土壤水蚀速率均值[2 154.6，706.5，223.1t/(km2·a)]与Bahadur[29]、Shrestha等[30]、Shrestha等[31]计算的均值[2 127，600，200t/(km2·a)]相近。本研究计算泰国平均土壤水蚀速率的结果为687.9t/(km2·a)，与Borrelli等[9]的结果类似[713.4t/(km2·a)]而与Rangsiwanichpong等[16]采用RUSLE模型的结果[812.5t/(km2·a)]存在差异，主要与采用不同分辨率数据导致LS因子准确度不同有关。Rangsiwanichpong等在计算水蚀速率过程中，采用1km分辨率的DEM计算LS因子，依据1km的DEM提取坡度，并对坡耕地赋水保措施值，由于分辨率较粗，这使计算的LS因子、水保措施因子误差较大。而本研究采用30m分辨率的DEM数据，计算结果更贴近真实值。通过对比可知，本研究计算结果较为合理，具有一定的可信度。

图5 土壤侵蚀因子重要性分类决策树

与抽样单元结果对比，泰国区域内175个抽样单元土壤侵蚀速率的总体均值为706.8t/(km2·a)，研究区计算结果的175个相应区域的总体均值为642.7t/(km2·a)，两者较为接近，总体计算精度为90.9%，进一步验证了计算结果具有一定的可信度。

与实地考察对比，2018年11月，张加琼等[20]对泰国北部的土壤侵蚀状况进行了实地考察，并沿考察路线布设了18个抽样单元。其18个抽样单元水蚀速率的总体均值为1 767.7t/(km2·a)，本研究在相应18个区域计算的水蚀速率总体均值为1 267.3t/(km2·a)，准确度为71.7%，差异较大，主要原因为：计算方法不同。张加琼等在计算考察的18个抽样单元的B因子值时，采用NPV数据作为林下植被盖度，这将低估生物措施的作用，使计算的B因子值偏高，进而使计算的水蚀速率偏高。本研究采用Borrelli等方法计算B因子值，结果更加合理。

图6 主控因子节点可视化图

图7 土地利用类型

4.2 治理建议

(1) 进一步加强土壤水蚀的基础理论研究。虽然泰国政府意识到土壤侵蚀的危害，但目前的土壤水蚀观测和调查资料匮乏[20]，难以满足土壤侵蚀治理的需求，现有资料也难以支撑准确评价水土保持效益，因此，建议进一步加强土壤侵蚀的基础理论研究，加大人才、科技、经费支持力度，构建土壤侵蚀的学科理论体系，为制定合理的水保措施和防治土壤流失提供理论支撑。

(2) 因地制宜制定土壤侵蚀治理措施。土地利用方式及变化对土壤侵蚀有重要影响[32]，经过多年耕作，土壤抗侵蚀能力明显下降[33]，合理的水保措施能够提高土壤抗蚀性[34]。离散分布于林地内部的耕地大多属于热点地区，侵蚀较为强烈，建议在林地区域内部禁止毁林开垦，并进行退耕还林草。东北部高原区的耕作区域主控因子为土壤可蚀性因子，针对该区域内的耕地，建议制定合理的耕作措施，选择合适的耕作措施类型和轮作制度，减少该区域的土壤流失。南部丘陵区的耕作区主控因子为坡度坡长因子，建议将坡耕地改为林草地或免耕。

(3) 提高作物单位面积产量。耕地的水蚀速率最大，而泰国49.6%区域属于耕地，建议控制耕地面积数量，提高单位面积作物产量。泰国作物单产具有较大提升潜力。如清莱府的稻米单产3 500kg/hm2，而土壤、气候条件类似的云南稻米单产为12 500kg/hm2[35]。

(4) 加强国际交流，引进切实可行的水土保持措施。从其他国家引进水保措施，并因地制宜加以改进，提高水土保持效益。如在坡面上覆盖的植被及其残留物，能明显减少地表径流和土壤流失，杂草覆盖是一种保护梯田的最便宜、最有效的方法[19]，这些方法也适用于泰国局部区域。

5 结 论

(1) 在泰国区域进行土壤水蚀定量评价，应优先选择直方图匹配方法制图。地图代数法制图结果具有较好的空间分布特征，但由于受到空间数据分辨率的限制，不能精确表达详细水保措施；在抽样单元尺度上，由于包括详细水保措施，计算的抽样单元水蚀速率较为准确，具有较为准确的统计特征，但难以精确表达空间分布。直方图匹配制图结果兼具两者的优点，因此，优先选择直方图匹配方法在泰国区域进行土壤水蚀定量评价制图。

(2) 虽然已有文献(如文献[9]、文献[16])对泰国进行土壤水蚀定量评价，但采用分辨率较粗，而本研究采用30m分辨率数据计算泰国土壤水蚀速率，并且通过与175个抽样调查结果对比、与实地考察结果相对比可知，本研究计算结果较为理想。文献[16]未对研究区进行主控因子分析，而本研究通过主控因子分析发现，东北部高原区的耕作区、中部平原区和南部丘陵区的耕作区、林地区域的主控因子分别为土壤可蚀性因子、坡度坡长因子、降雨侵蚀力因子，针对不同区域的主控因子，可采取更改土地利用方式、制定合理的耕作措施等不同的治理方案，治理研究区土壤流失问题。

(3) 泰国土壤水蚀较为严重，局部侵蚀较为剧烈。在各土地利用类型中，耕地水蚀最为严重，林草地水蚀速率较小，应因地制宜地制定不同区域的土壤流失治理措施，建议控制耕地面积数量，提高耕地单位面积产量。水蚀速率>2 500t/(km2·a)的热点地区主要分布于泰国北部和东北部，该热点地区84.1%区域为耕地，在有条件的区域，将该热点地区的耕地改变为林草地，对于减少水土流失具有显著效果。