Effect of Slope Cutoff in Extraction of Distributed Soil Erosion Slope Length

一般意义上, 沿坡度线方向, 由径流起点向下任一点, 其轨迹在平面的投影为坡长。如果用于计算坡长的数据为数字高程模型(digital elevation model, DEM), 则坡长与所对应单元格的等高线垂直且相等; 如果DEM分辨率足够高, 则可认为格网内等高线长度值与该点坡长值相等, 但方向垂直。基于此, 坡长计算公式可表示为^[13]:

式中, λ_{xi, yj}为坐标点(x_i, y_j)处的坡长; 积分下限(x₀, y₀)和上限(x_i, y_j)分别为径流源点和终点(i、j为对应的坐标位置); λ_c为单元坡长。根据格网的特性, 式(1)可被简写为:

式中, C为栅格边长; k为单元坡长系数, 用于计算单元坡长(由坡向决定); m为单元格编号。

因此分布式侵蚀学坡长可由以下步骤计算完成:①DEM无值点修复及洼地处理; ②获取坡度和坡向(水流方向), 生成坡度和坡向数组; ③依据坡向数据, 计算每一个栅格长度, 即单元坡长; 采用D8算法^[13], 如果该栅格单元的坡向是东、南、西、北, 则k=1, C为栅格长度; 如果是东南、东北、西南、西北, 则k=1.414;④遍历步骤③产生的坡度和坡向数组, 如果当前栅格满足截断条件, 则该点为坡度截断点; ⑤提取坡长, 沿着水流方向, 累计单元坡长, 遇到截断点, 停止累计, 同时作为新坡长的起点, 开始新坡长的累计。本文在前期研究的基础上, 只考虑坡度变化的截断作用, 使用C#语言开发了本文研究的测试工具(简写为LS_TOOL)^[18]。

为了方便分析和理解并充分考虑坡度变化过程, 本文研究采用了两个试验样例。

1) 试验数据1:数学曲面。模仿黄土高原高程及坡度变化, 将5段不同变化规律的抛物曲线连接成投影距离为2 000 m长、500 m宽的坡面(图 1(a)), 其中包含根据两个坡度变化的位置, 分别是坡度大于或等于2.86°(约占94.5%)时的坡度变化以及坡度小于2.86°(约占5.5%)时的坡度变化位置。根据McCool等的建议^[19], 地形因子的计算应该在0.7 m间距的地形图中生成, 能够与实际地形相符, 因此该数学曲面采用的分辨率为0.5 m。该数据不存在汇水, 不受沟道的影响, 适合于独立验证坡度变化引起的截断对坡长的影响。通过Matlab生成分辨率为0.5 m的数学曲面, 计算公式为:

Figure 1.  Ideal Mathematical Surface and Real DEM Data for Experiment

式中, z为曲面高程; x∈[0, 2 000] m; z∈[845, 1 579] m。高程分布情况见图 1(b), 坡度分布情况见图 2(c), 最大坡度为78.37°, 平均坡度为18.74°。当x∈[0, 400] m时, 坡度变缓; 当x∈(400, 800] m时, 坡度增加; 当x∈(800, 1 000] m时, 坡度衰减; 当x∈(1 000, 1 400] m时, 坡度减小; 当x∈(1 400, 1 860] m时, 坡度逐渐变陡; 当x∈(1 860, 2 000] m时, 坡度变缓。

Figure 2.  Slope Length Calculation Results of Mathematical Surface Under Different R₁

2) 试验数据2:县南沟DEM, 用于检验存在沟道的情况下, 只考虑坡度截断对坡长的影响。该流域地处黄土高原中部(见图 1(d)), 面积约44.85 km², 属典型黄土丘陵沟壑区, 地形破碎复杂, 高程和坡度变幅较大, 比较适合用作地形指标提取算法的测试。DEM数据是基于1:1万数字地形图(包括等高线、高程点和河流), 在ANUDEM软件下生成(见图 1(e)), 分辨率为2.5 m。1:1万数字地形图由国家基础地理信息中心提供, 等高距为5 m, 符合CH/T 1015.1-2007测绘标准^[20]。坡度小于2.86°的区域占0.66%, 主要集中在沟谷处; 坡度在20°~40°之间的区域占58.97%。

依据文献[18], 以2.86°为坡度变化点, 分别进行测试。①坡度大于或等于2.86°, 将下坡坡度变化的百分比(R₁)按照10%的比例递增, 从0增至1.00(对于坡度小于2.86°的位置, 截断值设置为1, 即不计算其截断); ②坡度小于2.86°, 将下坡坡度变化的百分比(R₂)按照10%的比例递增, 从0增至1.00(对于坡度大于或等于2.86°的截断值设置为1, 即不计算其截断); ③R₁、R₂同时变化对平均坡长、最大坡长的影响。为了检验不同情况下的坡度截断影响情况, 分别将R₁和R₂作为输入因子, 以平均坡长和最大坡长作为输出结果, 采用Canoco 5软件中的主成分分析(principal component analysis, PCA)进行评估^[21]。PCA方法通过线性变换将原始数据变换为一组各维度线性无关的表示, 可用于表示各特征分量的相关关系, 在Canoco软件中, 显示结果无量纲。

1) 对于数学曲面, 在不同截断情况下, 坡长分布情况如图 2所示。由图 2可以看出:①当R₁=0(见图 2(a)), 下坡坡度变小时, 坡长截断; 当下坡坡度表现为增大时, 没有出现截断位置, 坡长累计增长, 结果符合坡度变化截断的定义, 当坡度减小到一定程度时, 坡面侵蚀过程主要表现为沉积, 坡长截断。②随着R₁增加, 坡度截断位置逐渐减少(见图 2(a)~2(f))。由于R₁变大, 当坡度下降比例较大时, 坡长截断, 坡度变化的位置随着R₁的增加逐渐减少, 截断点减少, 坡长增加。③当R₁=1.0时, 即下坡坡度变为零, 在该曲面中只有边界处的坡度变化达到100%, 因此坡长累计过程未出现截断(图 2(f))。

不同截断情况下, 坡长平均值与R₁的关系说明:①随着R₁增加, 截断位置减少, 坡长平均值增大, 由75.69 m(R₁=0)增至999.98 m(R₁=1.0)。②当R₁ < 0.5时, 坡长增大缓慢; 当R₁＞0.5时, 坡长增大明显。③最大和最小R₁值导致的坡长变化约为92.43%。R₁与坡长λ的关系可表达为如下的多元回归式:

2) 对于县南沟, 为方便说明问题, 截取部分区域(包括坡面、沟谷), 坡长分布情况见图 3, 参考坡度分布见图 3(a)。由图 3可以看出:①在坡度截断R₁=0时, 下坡坡度变小的位置均出现截断(图 3(b)中的红色区域), 坡度小于2.86°的位置多集中在沟谷处, 坡度大于或等于2.86°的位置多集中在坡面处, 因此沟谷处可见连续的坡长(蓝线), 坡面位置多截断。坡度增大的位置, 坡长继续累计, 因此部分区域坡长呈现增长趋势(图 3(b)中的黄色区域)。②随着R₁的增加, 坡度变化达到该值的位置逐渐变少, 因此截断点随之减少。从图 3(b)至图 3(d), 短坡长(红色)减少, 长坡长(黄、蓝)增多。③当R₁=1.0时, 下坡坡度变化不再引起坡长截断, 因此坡长连续地分布在该区域。

Figure 3.  Slope Length Calculation Results of Xiannangou Catchment Under Different R₁ and R₂

对于县南沟流域, 在不同截断情况下, 由坡长平均值与截断比例R₁的关系说明:①随着R₁增加, 截断位置减少, 因此平均坡长增加, 由7.83 m(R₁=0)增至43.15 m(R₁=1.0)。②当R₁＜0.4时, 平均坡长增加显著; 当R₁＞0.4时, 坡长增长缓慢; 在R₁=0.5前后, 坡长变化缓慢。③最大和最小R₁值导致的坡长变化约为81.86%。R₁与坡长λ的关系可表达为如下的多元回归式:

与数学曲面不同, 由于实际地形的坡度变化较数学曲面复杂, 因此研究区的坡长随不同R₁的变化而变化。

1) 对于数学曲面, 坡度小于2.86°且出现坡度衰减的情况, 只出现在x坐标为1 400~1 500 m的位置(图 1(c)中的红色区域)。不同截断条件下, 坡长分布情况如图 4所示。由图 4可以看出:①在R₂=0~0.9之间, 坡度变化导致坡长截断, 但由于该范围的坡度位置较短, 主要集中在x坐标为1 400~1 500 m之间(图 4(a)、图 4(b)), 坡长截断作用不如R₁明显。②不同截断情况下, 统计坡长平均值与截断比例R₂的关系。由于数学曲面的特殊性, 坡度截断值R₂在0~0.5之间, 坡长变化较明显, 平均坡长从531.49 m增加为552.53 m; R₂在0.5~1.0之间, 坡长基本无变化, 这主要是由于坡度变化从2.86°降至0.30°, 变化范围小且距离短, 因此R₂的变化对坡长截断效果不显著(图 4(c))。③最大和最小R₂值导致的坡长变化约为53.15%。R₂与坡长λ的关系可表达为如下的多元回归式:

Figure 4.  Slope Length Calculation Results of Mathematical Surface Under Different R₂

2) 对于县南沟, 截取与坡度大于或等于2.86°一致的位置进行对比, 结果发现:①坡度小于2.86°的位置主要集中于沟谷处(图 3(a)), 因此随着R₂的增加, 坡度截断的位置逐渐减少, 沟道处的坡长逐渐连续(图 3(d)、图 3(e)、图 3(f))。对于县南沟流域, R₂和R₁的截断效果互补, R₁的截断位置主要集中在坡面, R₂的截断位置主要集中在沟谷。②R₂导致的坡长总变化约为17.08%, 相对于R₁对平均坡长的影响, R₂的影响较低, 主要原因在于:小于2.86°的坡度点较少, 且集中于沟谷, 导致R₂对平均坡长的影响较小。坡长λ与R₂的关系表达为如下的多元回归式:

1) 对于数学曲面, 将不同截断值(R₁、R₂)所对应的平均坡长作曲面图(图 5(a))。其中, R₁为坡度小于2.86°变化率的截断值, R₂为坡度大于或等于2.86°变化率的截断值。随着R₁和R₂的增加, 坡度变化引起截断的比例减少, 平均坡长增长。当R₁=0, R₂=0时, 平均坡长最短; 当R₁>0.6, R₂>0.5时, 坡长增加迅速。但数学曲面坡面变化平稳, 因此并不是所有R₁、R₂都出现截断。

Figure 5.  Effect of Average and Max Slope Length Under Different R₁ and R₂

2) 对于县南沟, 将不同截断值(R₁、R₂)所对应的平均坡长作曲面图(图 5(b))。与数学曲面类似, 随着R₁和R₂的增加, 坡度变化引起截断的比例减少, 平均坡长增长。当R₁=0且R₂=0时, 坡度下降的位置均表现为截断(图 3(g)), 平均坡长最短; 当R₁=1.0, R₂=1.0时, 无坡度截断(图 3(d)), 坡长一直累计, 因此平均坡长最长; 当R₁、R₂介于0~1.0之间时, 对比图 3(b)至图 3(h), 由于实际地形变化起伏较规律, 坡度分布更加均匀, 因此平均坡长随截断值(R₁、R₂)的变化更加均匀。

3) 统计数学曲面和县南沟两个测试数据的坡长最大值与R₁、R₂的关系(图 5(c)), 最大值的变化过程与平均值类似, 当R₁>0.7, R₂>0.5时, 坡长增加迅速, 这与前人研究^{[3, 10]}将坡度截断分别设置为0.5和0.7的结果基本吻合。

4) 县南沟PCA分析结果如图 6所示(结果无量纲)。图 6(a)中, 坡长与R₁、R₂箭头同向, 表明正相关, 即随着R₁、R₂的增大, 坡长平均值亦增长; 同时在坐标轴上的投影长度R₁>R₂, 表明R₁对平均坡长的作用大于R₂对平均坡长的作用。图 6(b)中, 最大坡长与R₁、R₂箭头同向, 表明正相关, 即随着R₁、R₂的增大, 坡长最大值增长; 同时在坐标轴上的投影长度R₂>R₁, 表明R₂对最大坡长的作用大于R₁对最大坡长的作用。该结果可以理解为:①R₁和R₂值的增加导致坡度变化增大才能截断, 因此无论是平均坡长还是最大坡长都会增大; ②由于样地坡度分布的原因, R₂的作用主要集中在坡面, R₁的作用集中在沟道, R₂导致坡长累计过程截断, 而坡面汇集到沟道处, 坡长增加显著, 但沟道单元格数量远小于坡面单元格数量, 即使坡长增加显著, 但由于较多的单元格坡长值较小, 所以导致了R₂对平均坡长的影响大, 但坡长最大值随着R₁的变大, 增长较明显, 因此R₁对坡长最大值的影响较大。该结果可以从图 5(b)、5(c)中看到, 图 5(b)中, 平均坡长沿X轴(R₁)方向的变化不如沿Y轴(R₂)方向的变化明显; 图 5(c)中, 最大坡长沿X轴(R₁)方向的变化比沿Y轴(R₂)方向的变化明显。③本文未将沟道截断考虑进来, 如果考虑沟道截断, 那么R₁对沟道的截断作用将可以忽略, 因此黄土高原地区坡度小于2.86°的截断可以不予考虑。

Figure 6.  PCA of Average Slope Length of Xiannangou Catchment Under Different R₁ and R₂

本文对坡长的坡度截断进行测试与结果分析, 得到如下结论:①由于研究区的坡度分布特性, 截断参数R₁对坡长最大值的作用更大, 截断参数R₂对坡长平均值的作用更大。在黄土高原地区, R₁的作用主要集中在坡面上, R₂的作用主要集中在沟谷处, 且沟道截断使得R₂的作用效果不明显; ②截断参数R₁和R₂越大, 坡度截断点越密集, 坡长截断越明显, 平均坡长越短; ③坡长增长的最大值统计规律与截断参数的经验设定有一定关联关系, 初步建议黄土高原地区R₁设置为0.7, R₂设置为0.5。但分布式侵蚀学坡长的截断因素有两个, 本文研究了坡度截断, 给出了坡长提取的变化情况及原因, 沟道截断的影响结果需要进一步研究。同时, 在不同分辨率下, 截断参数该如何设置是侵蚀模型应用遇到的较棘手的问题之一, 亟待进一步研究。坡长提取过程中的不确定性较多, 尽管前人进行了算法、分辨率上的研究^[22], 但综合截断参数、分辨率变化的研究还需尽快展开。