青藏高原典型样区2种土壤侵蚀评价与制图方法的对比

杜朝正,杨勤科,王春梅,庞国伟

(1 西北大学 城市与环境学院, 陕西 西安 710127;2 菏泽学院 城市建设学院, 山东 菏泽 274015)

土壤侵蚀是全球性生态环境问题之一[1]，被侵蚀的土壤有机碳是温室气体碳成分的重要来源之一[2]，土壤侵蚀速率大小直接影响水土流失区域的土壤侵蚀经济损失状况[3]。治理土壤侵蚀必须以土壤侵蚀定量评价和制图为基础[4-6]，即通常采用土壤侵蚀模型定量计算土壤侵蚀速率[7]。逐步建立的土壤侵蚀模型有美国的WEPP模型[8]、USLE[9]和RUSLE[10]，欧洲的EROSEM模型[11]，荷兰的LISEM模型[12]以及中国的CSLE[13]等。国外较多研究采用USLE/RUSLE计算区域土壤水蚀速率，如Panagos等[14]编制欧盟土壤流失图，Teng等[15]编制澳大利亚土壤侵蚀图。国内学者多采用CSLE计算区域土壤水蚀速率，如章文波等[16]采用CSLE预测坡面和田间土壤侵蚀，程琳等[17]采用CSLE对陕西省土壤侵蚀进行定量评价。

一般进行区域土壤侵蚀定量评价方法可归纳为2种：第一种是基于侵蚀因子图层整体运算的方法，可称为地图代数法，代表性研究有程琳等[17]、Borrelli等[18]、Yang等[19]、周来等[20]；第二种是基于抽样调查结果和插值的方法，在研究区布设基本抽样调查单元，结合土壤侵蚀模型计算抽样单元土壤侵蚀速率，使用空间插值模型来完成区域土壤侵蚀评价制图，简称为空间插值法，代表性研究有刘宝元等[21]、Yin等[22]。地图代数法和空间插值法均已被广泛应用，但这2种方法有何异同还鲜有报道。

本研究以青藏高原典型样区为例，以现有土壤侵蚀因子数据和抽样调查单元数据为数据源，分别采用基于侵蚀因子栅格图层的地图代数法和基于抽样调查单元的空间插值法，计算研究区土壤水蚀速率，并对计算结果进行对比分析，以期为土壤侵蚀定量评价与制图提供方法参考。

1 材料与方法

1.1 研究区域与基础数据

研究区位于青藏高原典型样区，总面积170万km2，海拔为331～7 182 m，平均海拔超过4 000 m[23]，年均降水量分布不均，由西北向东南呈逐渐递增的形式分布，是多条重要河流的水源地，也被称为“亚洲水塔”，为超过14亿人提供水源。目前该区域土壤侵蚀研究相对滞后[24]，迫切需要对土壤侵蚀进行试验观测和制图[25]。本研究所用青藏高原典型样区的基础数据见表1。

表1 青藏高原典型样区的基础数据Table 1 Basic data of typical sample areas of the Tibetan Plateau

1.2 空间插值法

1.2.1 抽样单元布设与土壤水蚀速率的计算 在研究区外围创建55 km缓冲区，在研究区和缓冲区内，按照0.25纬度×0.25经度布设样点，共计2 631个抽样单元，再根据《泛第三极土壤侵蚀抽样调查单元遥感解译工作大纲》，利用高分辨率遥感数据，解译抽样单元内土地利用类型和水保措施信息，用于支持抽样单元植被覆盖与生物措施因子(B)、工程措施因子(E)、耕作措施因子(T)的计算。再结合降雨侵蚀力因子(R)、土壤可蚀性因子(K)、坡度坡长因子(LS)等数据，根据水利部水土保持监测中心编写的《水土流失普查技术规定》和文献[26]，采用CSLE模型计算各抽样单元的土壤水蚀速率。

1.2.2 空间插值计算 以计算的抽样单元土壤水蚀速率平均值、GLC30中的土地利用类型数据为基础，通过空间插值法计算研究区土壤水蚀速率，本研究采用地统计协同克里金插值，该插值方法能够容纳多个解释变量，较好地反映多个环境因子对目标变量的影响[27]。与普通克里金方法相比，协同克里金插值可以输入多个协变量作为影响因素，使区域土壤水蚀速率的计算结果更加准确。具体步骤为：(1)统一土地利用分类系统，使解译土地利用类型与GLC30中的数据相对应；(2)按照一级土地利用类型，统计每个抽样单元的平均土壤水蚀速率；(3)按照一级土地利用类型进行空间插值，依据统计结果，并以降雨侵蚀力因子(R)、土壤可蚀性因子(K)、地形因子(LS)为协变量，按照不同土地利用类型依次进行地统计协同克里金插值，获得耕地、林地、草地、居住及建设用地、裸土等的插值面；(4)空间插值后处理：按照GLC30中的一级土地利用类型分布，将各土地利用类型插值结果图层进行合成，最终获得区域土壤水蚀速率栅格数据层，其中湖泊、坑塘、永久性冰川积雪、戈壁、裸岩等赋值为0。

为了确保计算精度及方便对比结果，本研究插值结果采用WGS84坐标系，区域水蚀速率专题层分辨率设定为1 km。

1.3 地图代数法

在1.2节中由于采用了CSLE模型计算抽样单元的土壤水蚀速率，为了便于对比，在地图代数法中也采用CSLE模型计算研究区土壤水蚀速率专题层，其基本形式如下：

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

(1)

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

1.3.1 降雨侵蚀力因子R、土壤可蚀性因子K、地形因子LS数据处理[26]根据项目组提供的大范围降雨侵蚀力R、中国范围土壤可蚀性因子K，通过裁剪处理，并重采样为1 km分辨率栅格数据，获得研究区的降雨侵蚀力因子R专题层和土壤可蚀性因子K专题层；将项目组提供的多幅地形因子基础数据进行拼接并裁剪出研究区域，再根据水利部水土保持监测中心编写的《水土流失普查技术规定》计算坡度因子上限阈值(坡度30°时坡度因子值，S上限=10.0)和坡长因子上限阈值(坡长100 m时坡长因子值，L上限=2.1)，进而计算得出LS上限阈值(LS上限= 21.2)，并重采样为1 km分辨率栅格数据，最终获得研究区LS数据专题层。

1.3.2 植被覆盖与生物措施因子B计算 根据水利部水土保持监测中心编写的《水土流失普查技术规定》和文献[26]中的方法计算B值。本研究中，取2014-2016年NPV均值作为乔木类林下地表植被盖度数据。草地、灌木林地、乔木林地B值采用下式计算：

草地：

(2)

灌木林地：

(3)

乔木林地：

SLRi=0.444 68×e(-3.200 96×GD3)-0.040 99×
e(YBD-GD3×YBD)+0.025，

(4)

(5)

式中：SLRi是第i个半月植被覆盖与生物措施因子值，无量纲；GD1、GD2、GD3分别为草地、灌木林地、乔木林地的地表植被盖度，取值0～1；YBD是郁闭度，指乔木林地植被覆盖度，取值0～1；WRi为第i个半月降雨侵蚀力占全年降雨侵蚀力的比例，取值0～1。

裸土地类B值：当24期半月植被盖度数据中的第14期半月植被盖度数据(the 14th fractional vegetation coverage，fvc14)>0.05时，B裸土=B草地；当fvc14≤0.05时，B裸土=0.5。耕地、居住及建设用地、湖泊、坑塘、永久性冰川积雪、戈壁、裸岩等土地类型B值按《水土流失普查技术规定》要求赋值计算。

1.3.3 工程措施因子E计算 依据抽样单元布设原则，按照0.25纬度×0.25经度间隔，绘制研究区内矢量网格，确保每个网格内只包括一个抽样调查单元，将抽样单元内E的平均值作为矢量网格的E值，再将矢量数据转换为1 km分辨率的栅格数据，获得研究区的E因子专题图。

1.3.4 耕作措施因子T计算 首先将刘巽浩等[28]编制的中国耕作制度区划数据转换为矢量图形数据，并与《水土流失普查技术规定》中的轮作分区的T因子值表进行数据关联，获得全国耕作区划T值图，对其进行裁剪，并将矢量数据转换为分辨率1 km的栅格数据，获得研究区的区域耕作区划T因子图；然后依据GLC30数据，将非耕地区域的T值设置为1，将耕地区域的T值设置为耕作区划T值，最终获得研究区的T因子专题图。

为了确保计算精度及方便对比结果，本研究计算采用WGS84坐标系，最终获得的R、K、LS、B、E、T、A专题层分辨率均设定为1 km。

2 结果与分析

2.1 2种方法计算的土壤水蚀速率的统计特征对比

对2种方法计算的土壤水蚀速率的基本统计信息进行对比，以累积频率为99%时的土壤水蚀速率为极大值进行统计，结果如表2所示。

表2 基于2种方法的青藏高原典型样区土壤水蚀速率的基本统计结果比较Table 2 Comparison of basic statistical results of soil water erosion rates based on two methods in the typical sample areas of Tibet Plateau t/(km2·年)

从表2可以看出，从土壤水蚀速率平均值来看，空间插值法的结果(593.5 t/(km2·年))与地图代数法(549.4 t/(km2·年))较为接近，表明2种方法计算的研究区土壤侵蚀状况基本一致。标准差可以衡量计算结果的离散程度。从土壤水蚀速率标准差来看，地图代数法的结果是空间插值法的1.7倍，说明空间插值法的结果更平滑，更能反映研究区土壤侵蚀的整体状况。从土壤水蚀速率极大值来看，地图代数法的结果(5 000 t/(km2·年))是空间插值法(2 800 t/(km2·年))的1.8倍，所以地图代数法能更好地反映研究区土壤侵蚀的极大值状况。

为了进一步比较2种方法计算结果的差别，以土壤水蚀速率计算结果为基础，设定25为步长，统计各步长频率，绘制频率曲线和累积频率曲线，结果如图1所示。

图1 基于2种方法的青藏高原典型样区土壤水蚀速率频率曲线(A)和累积频率曲线(B)的对比Fig.1 Comparison of soil water erosion rate frequency curve (A) and cumulative frequency curve(B) based on two methods in the typical sample areas of Tibet Plateau

从图1可以看出，地图代数法与空间插值法的土壤水蚀速率频率曲线和累积频率曲线差异较大。从土壤水蚀速率频率曲线来看，地图代数法的结果为单峰曲线，空间插值法的结果为双峰曲线。从土壤水蚀速率累积频率曲线来看，在地图代数法的结果中，土壤水蚀速率小值累积频率百分比较大且土壤水蚀速率大值累积频率曲线更长，可知地图代数法结果更突出土壤水蚀速率小值和大值，空间插值法结果更突出均值附近土壤水蚀速率值。出现这种状况的原因如下：(1)在地图代数法的的工程措施因子E计算过程中，以一个抽样单元的E均值作为其周围0.25纬度×0.25经度范围的E因子值，由于0.25纬度×0.25经度范围内的E因子值不是均匀分布，这导致0.25纬度×0.25经度范围内的E因子值要么被缩小要么被放大，进而使该范围内土壤水蚀速率小值和土壤水蚀速率大值的百分比均增大。(2)在空间插值法的计算过程中，以抽样单元的土壤水蚀速率均值为输入数据，这使得空间插值法的结果主要反映每个抽样单元的均值状况，不能反映出每个抽样单元土壤水蚀速率最小值和最大值状况。所以地图代数法的结果中土壤水蚀速率小值和大值所占比例均大于空间插值法结果，空间插值法所得的结果中土壤水蚀速率在均值附近出现第二个波峰。由于地图代数法的结果中土壤水蚀速率小值占比较大，同时空间插值法的结果侧重反映抽样单元土壤水蚀速率均值状况，抽样单元的土壤水蚀速率最小值和最大值均未参与插值运算，从而导致当土壤水蚀速率较小时，地图代数法累积频率高于空间插值法；当土壤水蚀速率较大时，地图代数法的累积频率曲线更长。

以统计的频率数据为基础，计算直方图交(histogram intersection，HI)[29]，判断2种方法计算结果的直方图相似性，一般HI取值范围为0～1，HI值越大，说明对比的2种数据结果越相似。经计算，2种方法计算结果的HI值为0.715，说明这2种方法计算结果具有较高的相似性，即2种方法计算结果反映的研究区土壤侵蚀状况基本一致。

综上可知，在青藏高原样区，从统计特征方面来看，地图代数法与空间插值法的计算结果均能反映研究区土壤侵蚀基本状况，其中空间插值法的结果能更好地反映研究区土壤水蚀速率均值左右的状况，且更贴近真实值。

2.2 2种方法计算的土壤水蚀速率表面特征的对比

采用地图代数法和空间插值法2种方法获得的研究区土壤水蚀速率表面特征分布结果如图2所示。图2显示，从宏观方面来看，2种方法计算结果均显示由北向南研究区土壤水蚀速率呈逐渐增强的趋势，分布趋势与实际相符。这可能是因为北部降雨量较少、地势相对平坦，南部降雨较多、且地形起伏比较大。土壤水蚀速率较高的区域大多分布于河流两侧，水动力的搬运能力较强，较易产生水力侵蚀。

图2显示，从微观来看，空间插值法的结果较为平滑，土壤水蚀速率较大值(>3 000 t/(km2·年))的分布较为集中；而地图代数法的结果较为粗糙，土壤水蚀速率较大值(>3 000 t/(km2·年))的分布较为分散。以2种结果为基础，计算局部方差(local variance，LV)，判断2种结果的局部变异情况[30]。LV平均值越大，表明数据表面的局部变化越大、表面越粗糙。经计算，地图代数法的LV均值(696 138.3)大于空间插值法(107 957.6)，说明地图代数法的结果能较好地反映局部区域变化情况。这是因为在地图代数法中，获得准确的面状水保措施较为困难，难以表达真实的土壤水蚀速率；但该方法采用的是逐栅格单元运算，每个栅格单元都具有侵蚀因子值，所以地图代数法的计算结果能够更好地反映局部区域水蚀速率变化情况。

图2 基于2种方法的青藏高原典型样区土壤水蚀速率的表面特征分布Fig.2 Surface characteristic distribution of soil water erosion rate based on two methods in the typical sample areas of Tibet Plateau

2.3 2种方法计算的土壤侵蚀强度面积所占比例的对比

对2种方法计算得土壤水蚀速率进行分级，并统计各等级侵蚀强度面积所占比例，结果如图3所示。图3显示，2种方法计算结果中，当土壤水蚀速率为0～25 t/(km2·年)时，其面积所占比例最大，并且不同土壤水蚀速率等级面积所占比例的变化趋势基本一致，说明2种方法计算结果反映的研究区土壤侵蚀强度级别状况基本一致。当土壤水蚀速率为>250～3 000 t/(km2·年)时，空间插值法的面积比例均高于地图代数法；当土壤水蚀速率≤250 t/(km2·年)或者>3 000 t/(km2·年)时，地图代数法的面积比例均明显高于空间插值法，说明空间插值法的结果更突出中间土壤侵蚀强度级别，而地图代数法的结果更突出较小和较大土壤侵蚀强度级别。表2中，采用空间插值法和地图代数法的所得的土壤水蚀速率平均值分别为593.5和549.4 t/(km2·年)，均位于500～750 t/(km2·年)侵蚀强度级别，而在此级别下空间插值法的面积比例高于地图代数法，可知空间插值法的结果相对更合理。

图3 基于2种方法的青藏高原典型样区不同土壤水蚀速率等级面积所占比例的对比 Fig.3 Comparison of area ratios of different soil water erosion rate grades based on two methods in the typical sample areas of Tibet Plateau

综上可知，空间插值法和地图代数法的结果均能反映研究区的土壤侵蚀强度级别状况，其中空间插值法的结果相对更合理。

2.4 本研究结果与其他研究结果的对比

Wang等[31]基于USLE和数学模型，综合评价西藏地区土壤侵蚀敏感性，判断西藏地区土壤侵蚀发生的可能程度，评价的趋势结果与本研究计算结果相吻合。此外，已有部分学者对研究区的平均土壤水蚀速率进行了定量计算，故将本研究采用空间插值法和地图代数法对不同研究区的平均土壤水蚀速率进行了计算，并将所得结果与前人研究结果(文献值)进行了对比，结果如表3所示。

表3 基于本研究2种方法所得不同研究区土壤水蚀速率与前人结果的对比Table 3 Comparison of soil water erosion rate in different study areas between the two methods and results from previous studies

表3显示，康琳琦等[32-33]采用USLE模型计算的青藏高原1980-2013年土壤水蚀速率为400～800 t/(km2·年)，与本研究区2种方法所得结果相近。Shao等[34]、李俊杰等[37]、宫奎方等[38]计算结果与本研究中2种方法计算结果也均相近，说明2种方法计算结果具有一定的可信度。Borrelli等[18]采用地图代数法计算的青藏高原典型样区潜在土壤水蚀速率均值为660.9 t/(km2·年)，本研究采用空间插值法和地图代数法计算的土壤水蚀速率均值分别为593.5和549.4 t/(km2·年)，比Borrelli等的结果稍小。分析其主要原因是，虽然本研究采用的CSLE与Borrelli等[18]采用的USLE/RUSLE没有本质差别，但Borrelli等的计算结果未考虑详细的水保措施，计算结果仅仅为土壤水蚀潜在研究，并不是真实的土壤水蚀速率，本研究2种方法均考虑了详细的水保措施，计算结果更接近真实值。本研究与文献[35-36]关于纳通河小流域计算结果差异较大，可能原因为文献[35-36]中的结果为137Cs采样点的平均值，未计算非采样点的区域土壤水蚀速率；而本研究结果为区域的平均值，这两种平均值之间可能存在一定差异。本研究与文献[39]关于拉萨河流域计算结果差异较大，可能原因是文献[39]中的水保措施较粗放，仅根据土地利用信息进行了主观赋值，且赋值偏高，其中文献[39]将盐碱地、沙地、冰川/永久积雪等土地类型赋值为1，这与水利部水土保持监测中心编写的《水土流失普查技术规定》赋值为0差异较大。通过对比可知，本研究的2种方法计算结果较为合理，具有一定的可信度。

3 讨 论

空间插值法与地图代数法2种方法的差异在于对基础数据的掌握程度以及对地面真实侵蚀状况的估算程度不同。

土壤侵蚀发生在微观尺度上，为了掌握土壤侵蚀真实值，必须在微观尺度上计算土壤水蚀速率，为此美国农业部的自然资源清查和中国水利部的水利普查(水土流失状况普查)等工作发展了土壤侵蚀抽样调查方法。为了准确获取微观尺度上的各侵蚀因子信息，特别是小流域尺度上的水土保持措施信息，必须进行抽样单元布设并对高分辨率遥感影像解译；在编制区域土壤侵蚀图时，为了计算无采样地区的土壤侵蚀速率，必须借助空间插值方法来实现[21-22]。受到方法自身的限制，空间插值法不一定很好地表达区域土壤侵蚀的局部变异；但该方法以抽样单元真实水蚀速率为基础，以降雨侵蚀力因子R、土壤可蚀性因子K、坡度坡长因子LS为协变量，计算的区域水蚀速率结果稳定且贴近真实值。

地图代数法是基于各侵蚀因子专题层进行整体运算，对于大范围研究区而言，各侵蚀因子专题层的分辨率与土壤侵蚀发生的尺度可能不完全匹配，水保措施信息难以用专题层的方式精确表达。因此，在缺少足够水保措施信息的情况下，地图代数法计算的结果是潜在土壤侵蚀，难以完全反映区域土壤水蚀速率真实情况；但该方法是逐栅格单元运算，每个栅格单元都具有侵蚀因子值，所以，地图代数法计算的区域水蚀速率结果能够更好地反映局部变化情况。

虽然Borrelli等[18]已在全球范围内研究数据的分辨率可达250 m，但本研究侧重计算方法的对比，考虑到本研究区大部分区域环境条件下土壤水蚀速率空间变异较小，暂时将计算结果分辨率设定为1 km。虽然本研究也得到了青藏高原典型样区不同土壤水蚀速率计算结果，但需要更多的实测数据来验证。以NPV数据作为林下盖度计算林地水蚀速率，可能使林地水蚀速率有些偏差，如何提升林下盖度计算精度，并结合更多的实测数据进行验证，还需要进一步研究。

4 结 论

在青藏高原典型样区，空间插值法和地图代数法均可用于区域土壤侵蚀定量评价与制图，2种方法计算结果均能反映研究区土壤水蚀速率总体分布状况，2种方法计算结果的宏观分布趋势相同，平均值相近，直方图相交值HI达0.715。地图代数结果能较好地反映局部土壤侵蚀变异情况，但面域范围内获得详细的水保措施信息较为困难；空间插值法以抽样单元真实水蚀速率为基础，其结果能更好地反映研究区的平均土壤水蚀状况，结果稳定且接近真实值，故在青藏高原典型样区可优先选择空间插值法。

志谢：本研究在刘宝元教授指导下完成，章文波老师和殷兵博士参与了计算过程的讨论，一并表示感谢。