20世纪初以来，水土流失日益严重，侵蚀问题已成为全球性的生态与环境问题^[1-3]，严重威胁着社会-生态的安全与稳定，关系着粮食安全、气候变化弹性以及地理社会的稳定性^[4]。土壤侵蚀危害在黄土高原地区表现的尤为突出^[5]，其水土流失面积高达45万km²，占总面积的71%，多年平均输沙量达16亿t，其水土流失面积之广，强度之大，流失量之多堪称世界之最^[6-7]。为了针对性地防治土壤侵蚀，必须加强土壤侵蚀的定量评价研究^[8]。

由降雨导致的土壤侵蚀受气候、物理、水文、化学、矿物质及生物因素的综合影响，如降雨强度、径流速率、土壤可蚀性等^[9]。在这些因素中，土壤可蚀性(soil erodibility)被认为是定量计算土壤流失的关键指标^[10]，是反映土壤性能和土壤侵蚀预报的必要参数，通常用K值来衡量其大小^[11-12]。土壤可蚀性研究对认识土壤侵蚀原理、定量估算土壤侵蚀量、合理进行水土流失防治有重要意义^[13]。

确定土壤可蚀性大小的最好方法是通过自然径流小区直接测定^[14]；然而由于径流小区研究必须具有足够大且观测历时较长的数据库资料^[15]，因此，国内外学者通常采用数学模型法估算K值。其中，应用广泛的代表性模型有美国学者Wishchmeier等^[16-17]建立的诺谟图方程和修正诺谟图方程、Williams等^[18]开发的侵蚀-生产力影响模型(EPIC)、Torri等^[19]建立的非线性最佳拟合公式、Shirazi等^[20]基于土壤几何平均粒径建立的估算模型等。上述估算方法自建立起，在我国不同区域土壤侵蚀敏感性评价、土壤流失量预测等方面均得到了广泛应用^[21-24]；但是由于土壤实际情况有所不同，上述经验公式在应用前，应对其预测结果进行必要的不确定性评价。对此，我国学者在不同地区进行了相应的研究。如史东梅等^[23]和张文太等^[24]先后对我国亚热带和紫色丘陵区不同K值估算方法的模拟精度进行研究，然而，黄土丘陵沟壑区不同K值估算方法的区域适用性仍尚待讨论。

安塞集水区地处水土流失十分严重的黄土丘陵沟壑区安塞县境内。安塞水文站输沙量监测值是指站点出口所控制的整个集水区的实测输沙量；泥沙输移比是实测输沙量与流域总侵蚀量的比值，是研究流域侵蚀产沙关系的重要依据，可用于推求流域的土壤侵蚀量^[25]。龚时旸等^[26]、牟金泽等^[27]和景可^[28]在黄土丘陵沟壑区的研究均证实黄土地区输沙量与流域产沙量基本一致，泥沙输移比约为1，即可直接用输沙量代替土壤侵蚀量。本研究选择安塞集水区作为黄土丘陵沟壑区典型集水区，基于5种K值估算方法，结合CSLE模型计算集水区2006—2014年平均土壤侵蚀量，对比相应年份输沙量监测值，探讨适用于黄土丘陵沟壑区集水区尺度的K值最优估算方法，以期服务于未来黄土丘陵沟壑区土壤侵蚀预报和水土流失的综合治理。

1 研究区概况

安塞集水区(E 108°5′44″-109°26′18″，N 36°30′45″-37°19′3″)地处延河流域上游，西北内陆黄土高原腹地，属典型的黄土丘陵沟壑区。集水区总面积1 334 km²，属中温带大陆性半干旱季风气候区。年均降水505.3 mm，且74%发生在6—9月份。年均温8.8 ℃，> 10 ℃积温2 876~3 270 ℃，年均日照时间2 397.3 h。地貌类型主要包括黄土梁涧、梁峁状黄土丘陵和沟谷阶地，地形起伏、沟壑纵横，地表支离破碎，地势西北高、东南低，海拔在997~1 731 m之间。土壤肥力低，以黄土母质上发育而来的黄绵土为主，土质疏松，土壤抗冲抗蚀性较差，水土流失严重，水土流失面积约占土地总面积的97%，属于强度水土流失区，也是西北典型的生态环境脆弱区^[29]。

2 材料与方法 2.1 数据来源

笔者采用的基础数据包括：空间分辨率为25 m×25 m的研究区数字高程模型(DEM)，来源于国家基础地理信息中心1:50 000数据库；2005和2015年2期土地覆被矢量数据，来源于中国科学院资源环境科学数据中心；降水数据来源于《中华人民共和国水文年鉴》黄河流域水文资料黄河中游区上段(河口镇至龙门)；研究区2006—2014年生长季遥感影像，来源于地理空间数据云网站；梯田、淤地坝等工程措施数据来源于《安塞县统计年鉴》；土壤理化性质来源于2014年7—8月土壤调查获取的151个典型样点数据集。

本次野外调查涵盖研究区主要分布的8种典型植被类型(天然草地、人工苜蓿草地、农田、苹果园、沙棘、柠条、刺槐、野山桃)。各植被类型的取样间距＞2 km，样地选择具有代表性(每种植被类型典型植被面积大于30 m×30 m，具有一定的主导作用)，且在研究区均匀分布。其中，农田和草地的样方大小为2 m×2 m，灌木样方大小为5 m×5 m，林地样方大小为10 m×10 m，手持GPS精确定位样地位置。

2.2 研究方法

中国土壤流失方程(China Soil Loss Equation)是在USLE模型和RUSLE模型的基础上，结合中国土壤侵蚀实际提出来的适用于中国土壤侵蚀特征的土壤流失预报模型^[30]，模型表达式为

$ A = R \cdot K \cdot L \cdot S \cdot B \cdot E \cdot T。$

(1)

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

2.2.1 土壤可蚀性因子K

本研究采用诺谟方程^[16](NOMO)、修正诺谟方程^[17](M-NOMO)、EPIC公式^[18]、Torri模型^[19]和Shirazi模型^[20]5种方法估算K值^[16-20]。考虑到不同方法间的可比性，对于K值笔者统一采用国际单位制单位。

2.2.2 降雨侵蚀力因子R

降雨侵蚀力因子反映降雨对土壤侵蚀的影响，是侵蚀预报的重要因子^[31]。鉴于日降雨与次降雨不是一一对应的关系，笔者将日降雨资料以半月为时段进行合并，利用章文波等^[32]提出的半月侵蚀力简易算法模型计算降雨侵蚀力。

2.2.3 坡度坡长因子LS

坡长因子和坡度因子反映土壤侵蚀对地形的响应特征。小区和坡面尺度土壤侵蚀评价中，LS因子一般通过野外实测获得，而在流域和区域尺度上，通常是基于DEM提取的；但是由于黄土丘陵沟壑区地形破碎，基于DEM提取的LS因子不能反映当地的真实情况^[33-34]，且相关研究发现，相较于坡长而言，坡度是影响黄土区LS因子的主导因素^[35]，因此，本研究以25 m DEM为基础数据，采用刘宝元等^[36]在黄土高原地区建立的方法计算坡度和坡长。在此基础上，采用汤国安等^[33]提出的坡度转化图谱法对所提取的坡度进行修正。该转换图谱是基于典型的黄土丘陵沟壑区地貌类型区获得的，平均纠正率在90%以上，适用于研究区坡度数据的修正^[33]。

2.2.4 植被覆盖与生物措施因子B

植被覆盖与生物措施因子反映地表覆盖对土壤侵蚀的作用，是指一定条件下有植被覆盖或实施田间管理的土地土壤流失总量与同等条件下实施清耕的连续休闲地土壤流失总量的比值^[37]，介于0~1之间。笔者基于研究区生长季遥感影像，提取NDVI值，计算植被盖度。根据B因子与土地利用类型和植被盖度的关系，对B因子赋值^[38]。

2.2.5 工程措施因子E

工程措施因子代表实施某种水土保持工程措施的土壤流失量与同等条件下无工程措施的土壤流失量之比^[36]。考虑到工程措施资料收集的困难性，本研究参考谢红霞等^[38]计算延河流域工程措施因子的方法计算E值。

2.2.6 耕作措施因子T

耕作措施因子是指采取专门措施后的土壤流失量与顺坡种植时的土壤流失量之比^{[31, 39]}，介于0~1之间。笔者基于T因子与不同坡度条件的关系对T因子赋值^[38]。

3 结果与分析 3.1 CSLE模型各因子计算结果 3.1.1 土壤可蚀性因子K

1) 描述性统计特征。基于K值估算结果，运用数理统计分析得到集水区K值的描述性统计特征(表 1)。

由表 1可知，K_EPIC、K_NOMO、K_M-NOMO、K_Torri和K_Shirazi的变化范围分别为0.032~0.060、0.046~0.092、0.047~0.088、0.009~0.066和0.018~0.044，最大值分别是相应最小值的1.875、2.000、1.872、7.333和2.444倍，可见K值的变动幅度较大。K值均值与中值之差十分接近于0，表明K值在集水区内分布较为均匀，未受到K值特异值的影响。K_M-NOMO的C_v值为0.067＜10%，表明K_M-NOMO在研究区内存在较弱强度的空间变异性；除此之外，K_EPIC、K_NOMO、K_Torri和K_Shirazi的C_v值分别为0.109、0.110、0.113和0.182，介于10%~100%之间，表明四者在集水区内存在中等强度的空间变异性。由于可蚀性受地形、气候、植被和土地利用等多种因素的综合影响，其空间变异性的存在是采样点土壤特性、植被生长状况、地形、气候和土地利用等多种因素综合作用的结果。

2) 正态分布检验。利用SPSS 20.0得到K值频数统计图(图 3)。可知，基于5种方法所得K值均呈倒钟形频数分布，与正态分布基本吻合^[40]。而非参数K-S检验的结果表明K-SP值均＞0.05，进一步证实K值在研究区内呈正态分布。

3) 半方差函数分析。半方差函数是地统计学解释土壤可蚀性空间变异结构的理论基础，包括块金值、变程和基台值3个重要参数^[40]。变程反映K值的空间变异特性，在变程范围内K值成空间相关；块金值反映的是最小取样间距内下K值的变异性及测量误差；而基台值则反映K值受土壤性质、地形、气候等结构性因子影响的程度。笔者基于ArcGIS 10.1地统计分析模块选择最佳半方差函数理论模型(表 2)。

由表 2可知，基于不同方法得到的K值，其最适半方差函数理论模型不同。K_EPIC的最适理论模型为指数模型，C₀/(C₀+C)值为11%＜25%，表明K_EPIC在变程内的空间相关性较为强烈；K_NOMO、K_M-NOMO、K_Torri和K_Shirazi的最适理论模型为高斯模型，C₀/(C₀+C)值分别为52%、49%、54%和42%，介于25%~75%之间，表明在变程内具有中等的空间相关性。步长为61.09 m，变程为733.06 m，步长小于变程，说明K值在集水区尺度下具有较好的空间相关性，通过Kriging插值能够得到较为准确的插值结果^[40]。

基于遴选出的最适半方差函数理论模型，经Kriging插值得到K因子空间分布图(图 4)。

由图 4可知，K_EPIC、K_NOMO、K_Torri和K_Shirazi的空间分布大致相同，表现为集水区东南部和中部K值较高而西北部较低的分布态势。东南部和中部地形较为破碎，土壤颗粒稳定性差，土壤结构等级系数和渗透等级系数相对较高是产生这一态势的主要原因。除东南部和中部地区以外，K_M-NOMO高值区域还分布于集水区的中北部、西南部和东部的局部地区，主要原因在于这些地区地势较高，地形坡度起伏较大，破碎度较高，植被生长较差，土壤结构松散。

3.1.2 降雨侵蚀力因子R

基于集水区周边20个雨量站点的日降雨数据，计算降雨侵蚀力，经Kriging插值得到R因子空间分布图(图 5)。

2006—2014年R值均值分别为1 365.06、1 416.39、820.89、2 560.23、1 251.06、1 072.43、1 811.11、4 416.59和1 765.42 MJ·mm/(hm²·h·a)。由图 5可知，除2008和2014年外，研究区R因子空间分异十分显著且不同年份的空间分异特征不同。如2009和2012年集水区R值大致呈现由东南向西北逐渐增大的态势；2011和2013年R值则表现为自西南向东北逐渐减小的分布特征；2008和2014年R值空间分异不明显，是由于相应年份降雨较为均匀所致；其余年份R值空间分异特征的出现，主要是由于集水区降水时空分布不均引起的。

3.1.3 坡度坡长因子LS

基于前文所述方法计算坡度和坡长因子，经坡度因子修正后，得到研究区坡度坡长因子LS(图 6)。

由图 6可知，LS因子最小值为0.02，最大值为78.12，平均值为11.76。集水区内部，LS值小于5的土地面积约占集水区总面积的20%；LS值小于平均值的土地面积约占总面积的60%；此外，大约95%的土地LS值低于27，这与谢红霞等^[38]在延河流域计算LS因子的结果基本一致。

3.1.4 植被覆盖与生物措施因子B

基于赋值法计算不同土地利用类型和不同植被盖度下的B值，利用ArcGIS 10.1得到集水区B因子空间分布图(图 7)。

由图 7可知，B值在集水区内呈现一定的空间分异特征。如2006—2011年，B值大致呈集水区东南部较低而西北部较高的分异特征；2012和2013年B值空间分布较为均匀，分异特征不明显；2014年集水区东北部植被盖度明显降低，B值呈减小的变化趋势，而西南部B值有所增加。此外，B值的年际变化较为明显，2006和2014年B值均值分别为0.11和0.15。集水区退耕还林(草)政策的实施使得坡耕地向林地和草地转化，是导致B值年际变化的主要原因。研究时段内，集水区耕地面积由2006年的139.14 km²减少至2014年的115.18 km²，草地面积由975.80 km²增加至1 041.53 km²。

3.1.5 工程措施因子E

基于研究区统计年鉴查找梯田、淤地坝数据，计算研究区2006—2014年E值(表 3)。

由表 3可知，集水区E值的年际变化不大，介于0.84~0.88之间。其中，2011年E值最小，为0.84，2013和2014年的E值最大，为0.88。

3.1.6 耕作措施因子T

基于研究区DEM提取坡度，根据坡度赋值法计算耕作措施因子，得到T因子空间分布图(图 8)。

由图 8可知，集水区T因子最小值为0.10，最大值为0.80，均值为0.64。且T值在集水区内空间分异不明显，除主要道路区域T值较低外，其余区域T值均较高。

3.2 基于不同K值的土壤侵蚀量计算

将基于5种方法得到的K因子及CSLE模型其他各因子的栅格图层加载到ArcGIS 10.1中，计算各因子的乘积得到集水区土壤侵蚀量。按照土壤侵蚀强度分级标准SL190—2007^[41]，将研究区土壤侵蚀划分为微度(≤500)、轻度(500~2 500)、中度(2 500~5 000)、强度(5 000~8 000)、极强度(8 000~15 000)和剧烈(＞15 000)6个等级，重分类后得到集水区土壤侵蚀空间分布图(图 9)。

由图 9可知，基于不同K值得出的土壤侵蚀量在同一年内，大致呈现相同的空间分布态势。例如：2006年大致呈集水区北部较低而南部和西北部较高的分布态势；2008年土壤侵蚀空间分异不明显，总体侵蚀程度较轻；2011年集水区东部侵蚀较轻而北部侵蚀有所加重。另外，从图中可以明显看出，2009、2013和2014年的土壤侵蚀量明显高于其他年份，且2013年侵蚀的空间分异程度相对较弱。其原因主要因为2009和2013年高强度降水事件的影响。此外，2013年极端降水不仅导致当年土壤侵蚀明显加重，而且造成了植被、梯田等水土保持措施的大面积损坏，也因此减弱了2014年的水土保持能力。

3.3 集水区K值最优估算方法遴选

利用ArcGIS 10.1计算集水区2006—2014年历年平均土壤侵蚀量，对比相应年份输沙量监测值，采用平均绝对误差(MAE)、平均相对误差(MRE)、均方根误差(RMSE)和精度因子(A_f)多指标相结合的方法，评价5种K值估算方法的区域适用性。MAE、MRE和RMSE越趋近于0，A_f越接近于1，则K值预测的不确定性越小，模型适用性越高^[23-42]。各评价指标的具体计算方法见文献[24]，评价结果如表 4所示。

由表5可知，基于5种K值估算方法的集水区土壤侵蚀模数模拟值差异显著。研究时段内，年均土壤侵蚀模数依次为65.59、106.00、108.47、76.69和47.68 t/hm²。相较于安塞水文站实际年均输沙模数17.73 t/hm²，由MAE、MRE和RMSE越趋近于0，A_f越接近于1，模型区域适用性越高的原则可知，Shirazi模型的适用性最高，MAE、MRE和RMSE值依次为30.93、3.25和43.66，A_f值为4.41，EPIC法次之，A_f值为5.80，Torri模型的适用性处于中间水平，修正诺谟方程和诺谟方程的适用性最低，与实际情况相差最大，A_f值分别为7.99和7.88。由此可知，相对于其他方法而言，Shirazi模型和EPIC公式更加适用于黄土丘陵沟壑区集水区尺度的K值模拟。目前，在黄土高原土壤侵蚀评价中，研究者大多采用EPIC方法计算K值^{[12, 21, 38, 43]}，因此，在未来该区域K值模拟乃至土壤侵蚀评价研究中，应考虑采用Shirazi模型，以期提高K值预测精度。另外，史东梅等^[23]在紫色丘陵区的研究中发现，EPIC法和诺谟法更加适用于紫色丘陵区的K值模拟；张文太等^[24]在我国亚热带地区的研究发现，Torri模型K值预测的不确定性最小，预测效果最佳。这也说明，由于不同地区土壤性质差异的存在，致使不同K值估算方法在不同地区的适用性不同，因此在某一特定区域的K值估算和侵蚀预报过程中，优先遴选适用于研究区的K值计算方法是十分必要的。

4 结论

1) 基于5种方法所得K值均值的大小关系为K_M-NOMO＞K_NOMO＞K_Torri＞K_EPIC＞K_Shirazi。从空间分布情况来看，K_EPIC、K_NOMO、K_Torri和K_Shirazi的分布状况大致相同，表现为集水区东南部和中部K值较高而西北部较低的分布态势。东南部和中部地形较为破碎，土壤颗粒稳定性差，土壤结构等级系数和土壤渗透等级系数均相对较高是产生这一态势的主要原因。除东南部和中部地区以外，K_M-NOMO高值区域还分布于集水区的中北部、西南部和东部的部分地区，主要原因在于这些地区地势较高，地形坡度起伏较大，破碎度较高，植被生长较差，土壤结构松散。

2) 基于5种K值的集水区土壤侵蚀模数模拟值差异显著。基于EPIC、诺谟方程、修正诺谟方程、Torri模型和Shirazi模型的2006—2014年平均土壤侵蚀模数依次为65.59、106.00、108.47、76.69和47.68 t/hm²。对比安塞水文站相应年份实际年均输沙模数17.73 t/hm²可知，Shirazi模型的区域适用性最高，EPIC法次之，Torri模型的适用性处于5种模型的中间水平，而诺谟方程和修正诺谟方程的适用性较低，与实际情况相差最大；因此，本研究认为相对于其他4种K值估算方法而言，Shirazi模型的估算结果更贴近实际情况，更加适用于集水区的K值模拟，在未来黄土丘陵沟壑区土壤可蚀性K值估算乃至土壤侵蚀评价工作中，应优先选用基于几何平均粒径的Shirazi模型。