桑云云 罗明良 杨 慧 白雷超 姜佳艳 刘维明

（1西华师范大学生命科学学院，四川南充 637009；2西华师范大学地理科学学院，四川南充 637009；3西华师范大学四川省干旱河谷土壤侵蚀检测与控制工程实验室，四川南充 637009；4中国科学院水利部、成都山地灾害与环境研究所，四川成都 610000）

我国土壤侵蚀分布范围较广，是土壤侵蚀较为严重的国家之一[1-2]。坡长是土壤侵蚀领域研究中的重要因子之一[3]，当其他条件相同时，土壤侵蚀的强度由坡面的长度来决定，坡面越长，汇聚的流量越大，侵蚀力和冲刷力就越强[4-5]。目前，关于坡度因子的计算已经相对成熟，而坡长因子的精确、快速提取是目前研究的重点之一[6-7]。

对于坡长的计算，学者提出了多种方法，例如：Hickey 等[8]在水文分析的基础上提出了非累积流量的直接算法；Desmet 等[9]提出了基于累积流量的单位汇水面积算法；杨昕等[10]基于Arcview空间分析提出坡长快速近似算法。随着计算机和GIS 研究方法的不断创新，在DEM 数据和GIS 空间分析方法的支持下，坡长提取的方法得以不断改进。其中，非累积流量的直接算法可以通过按照不同的流向进行赋值来体现坡长的准确特征[11]，但此方法需要大量的计算和复杂地形因子的提取；快速近似算法虽然计算过程简便，但由于地形类型复杂，局部高程点无法准确识别，导致计算结果误差较大[8]。基于网格DEM 的坡长计算可用来模拟地表的水流路径，但要求水流在DEM 模拟的地形面上自由流动，因此地表的自然凹陷对坡长的计算结果有很大影响，需要在确定水流方向之前对DEM凹陷处进行填充和平整。此外，坡长的计算还需要考虑很多因素，不仅包括地形特征，还需要考虑具体地块的土地利用类型[12-13]。坡长作为土壤侵蚀计算中不可忽视的因子，迫切需要一个高度精确的计算方法。

为此，本研究以黄土高原甘泉和延川小流域为研究区，基于5 m分辨率数据，提出在反地形条件下不断逼近真实坡长结果的算法。该研究不仅降低了原有坡长计算方法中由于填洼造成的误差，也有效减轻了计算过程中的工作量，为坡长计算提供了一定的思路和借鉴。

1 研究区概况

本研究选取黄土高原侵蚀剧烈、地形起伏大的甘泉和延川2个典型的丘陵沟壑地貌区作为研究对象[14-15]。研究区海拔为600～1 800 m，甘泉位于陕西省延安市中部，属于陕西省北部黄土高原黄土丘陵，以黄土梁状丘陵沟壑地貌为主，沟壑密度约5.6 km/100 hm2。延川则以梁峁状丘陵地貌为主，梁峁起伏，残塬、梁、峁相间分布，沟壑纵横、河谷深切，沟壑密度在所选样区中最大。研究区内以大陆性季风气候为主，大陆性气候特征明显，年均降水量约447 mm，多短时强降雨，土壤侵蚀强烈，如图1所示。

图1 研究区地形

2 研究方法

2.1 研究原理

坡长一般定义为从坡面径流的起点到径流被拦截点或流路中断点的水平距离[8，16]。如图2 中A′B′是径流线在水平面上的投影，即坡长。

图2 坡长示意

本研究在反地形条件下计算坡长时，为确定水流的起点，将原始地形分解为若干个集水盆地。当反转地形时，原始正向地貌的山地反转成为负向流域及集水区系统，相应的山顶点也演化为各集水区中的洼地点，原有正向地形上的洼地成为向上凸起的小丘。反地形水流的起点在原始地形中集水盆地的凹陷点上。该方法虽省略了填洼过程，但并不会对原始地形流路和相关数据的提取产生影响。这种计算过程不仅简单易懂，而且可以减小多种数据融合后计算过程中产生的误差。

根据坡长的定义，如果将坡上的水流路径（M）设为直角三角形的斜边，则坡长可定义为坡上流线的长度乘以cos（α），其中α为流向与水平方向的夹角[17-19]。为计算α，坡长逼近算法（ASL）首先计算从水流流动起点到流路上某一点的高差，定义为δh，然后用反正弦函数求解α，即α=arcsin（δh/M）。按照这种思路，坡长逼近算法ASL的设计步骤如下。

反地形和相关水文分析（图3），如流向、汇流累积量和分水岭等[20]。其中，默认的流向算法为D8流向算法，即单元梯度的最佳代表值是以其为中心的3×3窗口中周围8个方向的最大梯度，水流的流动方向与最大梯度相同。

图3 水文分析

分区统计后利用栅格计算器计算前一步分区统计结果与原始DEM的差值，利用重分类工具将上一步栅格计算的结果进行赋值显示后，复制为1 的区域即为区域最高点。接着使用水文分析中的水流长度工具基于所得到的流向数据计算出反地形下的水流长度（上坡长度：计算沿水流路径从每个像元到分水岭顶部的最长上坡距离），即为地形坡面长度。

利用网格计算器构建地图表达式α=arcsin（δh/M），并计算坡面与水平面之间的包含角α；使用三角函数cos构建函数表达式L=M×cos（α），计算结果L即坡长（图4）。

图4 坡长计算示意

2.2 坡长逼近算法验证

为验证ASL计算坡长的精度，以相对差值系数和线性函数为描述指标，对3种坡长计算结果进行分析。

2.2.1 LS_TOOL（LST）LS_TOOL 算法是基于GIS 技术和多方位地形系数的估算方法，用于获取区域尺度的地形坡长坡度参数等数据。以LS_TOOL 算法计算的坡长（L）作为实际坡长，验证本文提出的近似算法的可靠性[18]。

2.2.2 SAGA（SGL）SAGA软件中的坡长计算工具可以计算特殊情况下的坡长。首先用Zevenbergen &Thorne 方法计算每个单元的斜率，然后使用D8单流向算法和阈值标准，将坡长值向下累积。

2.2.3 相关性分析 利用线性回归方法拟合出算法间的一元线性函数，当2 种算法之间的坡长值完全一致时，则函数关系为X=Y，XY 散点分布图是一条斜率为1的直线；当2种算法之间的坡长值存在差异时，另一种函数关系为Y=a+bx，表明2种算法之间存在系统差异。接着计算相对差异系数，相对差异系数通过逐个栅格单元进行对比，可分区域分析算法之间的相似性。当2 种算法进行比较时，选择其中一个作为基准算法，另外一个则作为比较算法，它们之间的差值θ如公式 （1）所示，可以通过减去相应位置的结果而得到（是基准算法在网格单元i的坡长为整个DEM的平均坡长，是比较算法在相应位置的坡长。）

为能够更清楚地比较坡长算法之间的相关关系，将其转换为相对差异系数α，表示比较算法和基准算法之间的偏差，其公式如下：

当α=1 时，2 种计算方法的结果完全一致，没有差异；α值越小，则算法之间的差异越大；当α<0时，则2种算法基本没有可比性。

3 结果与分析

3.1 不同区域3种坡长算法

ASL、LST和SGL 3种计算方法得到的甘泉和延川样区的坡长结果如图5～6 所示。从图5～6 可以看出，同一样区3种算法的坡长变化大致相同，但坡长结果均有差异，主要表现在3 种坡长算法的最大值上，ASL算法的最大值分别为564和468，LST算法的最大值为493 和353，SGL 算法的最大值为399 和389，ASL 算法相对于其他2 种算法较大的坡长值主要集中在山脊处。

图6 延川样区坡长（aASL、bLST、cSGL）

经过抽稀提取2 个研究区20%的坡长结果，并将结果划分为3个部分绘制累积坡长（图7）。从图7可以看出，3种算法得到的坡长在0～300 m，占提取坡长总数的98%左右；而大于300 m的坡长较少，占2%左右。ASL受流向和网格分辨率的影响，在甘泉样区有超过500 m的离群值，但其累计数量不大，仅占坡长总数的0.001 7%；而延川样区超过400 m 的离群值占总数的0.002 4%，这些离群值均出现在山脊处。

图7 坡长累积频率分布（a甘泉坡长结果频率分布，b延川坡长结果频率分布）

3.2 相关性分析

为检验ASL算法的可靠性，利用ArcGIS 10.2中的抽稀工具提取了3 种算法坡长结果20%的数据，并分别计算ASL 与LST、ASL 与SGL、SGL 与LST 算法之间的相关关系。

3.2.1 线性回归分析 每2种坡长算法间的一元线性函数以及R2值见表1。

表1 3种算法的相关性分析结果

算法之间的线性函数的R2值均在0.900 以上，其中最大的R2值出现在延川样区的ASL 与LST 之间，为0.993；最小R2值出现在延川样区的SGL 与LST之间，为0.937，表明3种算法的坡长结果之间存在强相关性。

3.2.2 相对差异系数 利用抽稀得到的2个样区20%坡长数据计算其相对差异系数（a，b），结果见表2。

表2 3种算法的相对差异系数

从表2 可以看出，算法间的相对差异系数均在0.80 以上，说明每2 个样区的坡长结果高度相关。其中，甘泉区ASL 与LST 的相对差异系数最高，为0.96；延川样区ASL与SGL的相对差异系数最低，为0.80。说明整体上这2 个样区的坡长结果相似度很高，但不同样区算法间的相对差异系数不完全相同，表示算法之间存在一定的差异。

4 讨论和结论

4.1 讨论

影响坡长计算结果并产生差异的原因，首先是地形表达的有限分辨率限制了地形高度的表达精度[21-22]。受网格大小的限制，一些山峰或山脊表达为平坦地形，导致集水区的划定结果出现一定误差。由图8 中的等高线可以看出，地形变化比较明显，顶部有大面积的平坦区域。结合谷歌图片，可以看出坡长出现较大值时存在大量的梯田、黄土脊状山丘等平坦地貌类型，地形相对破碎，这对DEM反映真实地形带来了较大挑战，进而影响了坡长的计算。如样本区有较大的平坦区域，则会影响反地形下的分割精度，从而导致坡长计算结果误差较大[23]。本研究通过反地形计算坡长，去除了在原始地形中填充洼地带来的误差，同时利用函数方法计算坡长提高了坡长计算的准确度，降低了试验过程中可能带来的误差，同时也减轻了计算过程中的工作量。

图8 研究区坡度和等高线（a甘泉，b延川）

其次，D8流向算法的局限性可能会影响计算结果。流向是8 个可能的网格方向之一，它们之间的间隔为45°，流向是确定的且不连续的。因此，D8算法具有高度敏感性，即较小的高度误差也会导致水流方向的改变，进而导致坡长计算结果出现偏差。其他2种坡长算法同样采用D8流向算法，因此不同流向算法对坡长的具体影响，还需要在后续的研究中做进一步的分析。

4.2 结论

本研究基于5 m 分辨率高程数据、利用ASL、LST、SGL 3种坡长计算方法对黄土高原丘陵沟壑区2 个样区的坡长进行计算，并利用线性回归和相对差异系数对计算结果进行了相似性分析，得出以下结论。

相同精度的DEM数据得到2个样区内不同计算方法结果具有差异。本研究提出的逼近算法其坡长结果98%的数据主要集中在0～300 m，并且线性回归以及相对差异系数分析结果表明3种算法计算得到的坡长结果非常相似。3个坡长计算方法通过回归函数得到的R2值均在0.9 以上，同时2 个样区的3个算法间的相对差异系数也主要集中在0.9附近，表明3 种坡长计算结果拟合效果均较好。因此，为获取更高精度的坡长结果，可将3 种算法的坡长平均值作为最逼近真实的坡长用于分析地区的土壤侵蚀状况。

坡长因子作为估计土壤侵蚀状况的重要参数，寻求一个不断逼近真实坡长结果的算法对提高土壤侵蚀治理水平十分重要。本研究提出的坡长逼近算法为坡长计算提供了更多的思路。