赵明伟，汤国安，李发源

(南京师范大学地理科学学院，江苏南京210046)

1 问题的提出

沟壑密度是反映黄土高原地区地面受侵蚀程度重要的区域水土流失因子，在该地区的土壤侵蚀与侵蚀产沙研究中得到了广泛的应用［1－3］。目前，不同学者提出了不同的提取沟壑网络的方法［4－7］，其中基于DEM的坡面流模拟方法因为数据来源广泛、方法简单、算法稳定性高而得到普遍应用。这一方法根据水文学坡面流模型来判别水流路径，首先对地面洼地进行填充，构建无洼地的连续地面［8］，并根据每一个栅格单元与相邻单元之间的最陡坡度，识别求取沿水流路径的汇流累积量值，再选择适当的汇流累积量阈值来确定河网。该方法可以直接生成连续的河网。但是，在计算格网单元流量累积量值时，本不会产生地面径流的平坦地区(例如黄土塬面等)也将产生很大的汇流累积量值，从而造成提取沟谷网络时会有大量的伪沟谷出现［9－10］。以上问题的出现，给黄土高原地区沟壑密度的正确提取以及水文分析的有效进行，都带来了很大的困难。

图1为实验样区水文分析中存在的伪沟谷示意图，(a)为实验样区的地形等高线图，(b)为该地区采用坡面流方法提取的沟谷网络，对比可以发现，原始提取的沟谷网络中存在大量的伪沟谷。

图1 伪沟谷示意

影响提取沟谷网络结果的关键因素是汇流累积量阈值的选取［11］，一般说来，当该阈值变大时，所提取的沟谷网络结果数量会变少。李俊、陈涛等以陕北黄土高原为研究样区，研究了沟壑网络提取中汇流累积量阈值的选取依据［12－13］。但是，对于黄土塬、黄土残塬等地区，由于大量平地的存在，试图仅仅通过调整汇流累积量阈值来消除伪沟谷是不可行的。图2为DEM数据的山体阴影晕渲图，可以看出A点处于沟谷中，而B点处于塬面平坦区;对应于水流路径汇流累积量图3，B点的汇流累积量值却大大高于A点。由此可见，在黄土塬、黄土残塬等存在大面积平坦区域的地区，不可能存在合适的阈值，使得提取的沟谷网络既包含较小的沟谷，又不存在伪沟谷。

如上所述，在平原台地地区(如黄土塬面)，由于大量平坦区域的存在，使得当采用坡面流模拟方法进行沟谷提取时会产生大量伪沟谷，且无法通过调整汇流累积量阈值来消除。因此，急需寻求一种简单而有效的解决方法。

图2 山体阴影晕渲图

图3 汇流累积量图

2 实验方法

2.1 实验样区与数据

实验样区位于陕西省淳化县城西北部泾河中游地区，黄土塬及残塬为主要地貌类型，塬面地形平缓，但被诸多大沟谷深切割裂，平均海拔1045.04 m，相对高差390 m，塬面平均坡度为8.1°，面积100 km2(10 km×10 km)。此外，实验中为了对比文中提出的实验方法消除伪沟谷的效果，还选取了地貌类型属于黄土台塬、黄土残塬的样区进行附加实验。附加实验样区的基本情况如表1所示。

表1 附加实验样区概况

实验数据为国家测绘部门生产的5 m分辨率的DEM数据，数据精度满足实验要求。

2.2 实验思路与方法

如前所述，基于DEM提取沟谷网络是在得到水流路径汇流累积量矩阵的基础上，根据一定的阈值提取沟谷网络。水流路径汇流累积量矩阵的产生示意图见图4，可以看出，某点汇流累积量值的确定既受水流方向的控制，同时也与分析区域面积紧密相关。因此，在大面积的平原台地地区，虽然地势变化缓慢，但是仍然会产生很大的汇流累积量值，因而采用常规方法提取水系时会产生很多伪沟谷。

图4 汇流累积量矩阵计算示意

根据以上分析，可以在计算汇流累积量矩阵时对每个栅格点增加权重值，使得平原台地地区汇流累积量减小，从而减小伪沟谷出现的可能性。实验发现，在黄土高原地区，将地形表面划分为正地形与负地形之后，沟谷网络属于负地形区域。因此，本实验根据地表面的正负地形属性来设置权重值，正地形区域的汇流累积量权重值为0，负地形区域的汇流累积量权重值为1;将权重矩阵与原始汇流累积量矩阵做点乘运算得到新的汇流累积量矩阵，再设置阈值提取沟谷网络。综合以上分析，在平原台地地区基于DEM数据提取正确的沟谷网络的步骤为:

(1)基于DEM数据得到汇流累积量矩阵［14］。首先，对原始DEM数据进行洼地填充处理，然后，利用GIS水文分析，得到提取区域的水流方向矩阵及汇流累积量矩阵。

(2)基于DEM数据获取实验样区地表的正负地形属性［15］，见图5。将正地形区域栅格属性设为0，负地形区域栅格属性设为1，作为权重矩阵。

(3)将步骤(1)得到的汇流累积量矩阵与权重矩阵做点乘运算，得到新的汇流累积量矩阵。

(4)基于新的汇流累积量矩阵设置阈值，得到正确的沟谷网络，其中阈值的选取要保证达到所要提取的最小沟谷的汇流累积量。

2.3 实验结果与分析

实验结果见图6，图7是基于该地区的DOM数据手动勾绘的沟谷网络。实验结果与图1(b)的原始结果比较可以看出，本文提出的消除伪沟谷的方法可以有效去除平坦塬面上的伪沟谷，并且细小的沟谷以及平坦的大沟谷都得到了很好的保留;根据处理结果计算的沟壑密度为1.06 km/km2，而根据原始数据计算的沟壑密度为2.05 km/km2，二者的差异是相当明显的。由手动勾绘得到的沟谷网络，相应的沟壑密度为0.92 km/km2，其值略小于伪沟谷消除后的结果。从图中可以看出，手动勾绘的沟谷大部分是直线，而消除伪沟谷后的结果中沟谷有一定的弯曲，相对于手动勾绘的沟谷网络，伪沟谷消除以后的沟谷网络更符合实际情况。

为了定量描述实验方法对伪沟谷的消除效果，如果将消除伪沟谷后的沟谷网络近似看作实际区域的沟谷网络，则定义伪沟谷消除率为伪沟谷总长度与实际沟谷总长度的比值，即

式中:R是沟谷消除率;∑L真是真实沟谷的总长度，这里把消除伪沟谷后的沟谷网络作为真实沟谷网络的近似值;∑L伪是伪沟谷的总长度，在计算时可以用原始提取的沟谷总长度减去真实沟谷的总长度。

表2为实验样区消除伪沟谷之前的沟壑密度(DS1)、消除伪沟谷后的沟壑密度(DS2)以及伪沟谷消除率(R)的计算值。

表2 实验前后沟壑密度以及伪沟谷消除率计算值

从表中可以看出，在各种地貌类型的实验样区中，黄土塬的伪沟谷消除率最高，黄土台塬次之，黄土残塬的伪沟谷消除率最低。这是因为黄土塬区存在大面积的平坦区域，因此在采用坡面流模拟方法进行沟谷网络提取时，这些平坦区域的汇流累积量值很大，造成了大量的伪沟谷;黄土台塬的情况类似，只是塬面平坦区域面积有所减小，因此伪沟谷消除率较之黄土塬区也有所下降;而黄土残塬地区由于破碎程度较高，不存在大面积的平坦区域，所以这些地区在提取沟谷网络时产生的伪沟谷数量较少，伪沟谷消除率在三类实验样区中也是最小的。

3 结论与讨论

本文根据黄土塬、黄土残塬以及黄土台塬等地区的地形特点，以该区域的正负地形属性作为权重矩阵，与原始汇流累积量矩阵做点乘运算产生新的汇流累积量矩阵，以此汇流累积量矩阵设置阈值所提取的沟谷水系有效地去除了原始方法结果中的伪沟谷。该方法的伪沟谷消除效果对于大面积平坦地区最为有效，随着样区地形破碎程度的增加，伪沟谷消除率呈递减的趋势。

［1］白占国.从地貌空间结构特征预测土壤侵蚀的研究——以神府、东胜煤田区为例［J］.中国水土保持，1993(12):23－24.

［2］卢金发.黄河中游流域地貌形态对流域产沙量的影响［J］.地理研究，2002，21(2):171 －178.

［3］廖义善，蔡强国，秦奋，等.基于DEM的黄土丘陵沟壑区流域地貌特征及侵蚀产沙尺度研究［J］.水土保持学报，2008，22(1):1 －6.

［4］O’Callaghan J F，Mark D M.The extraction of drainage networks from digital elevation data［J］.Computer Vision，Graphics，and Image Processing，1984，28(3):323 －344.

［5］韦中亚，周贵云，罗万勤.一种基于数学形态学的沟壑密度提取算法［J］.地理学与国土研究，2001，17(2):24 －27.

［6］陈永良，刘大有，虞强源.从DEM中自动提取自然水系［J］.中国图象图形学报，2002，7(1):91 －96.

［7］李昌峰，冯学智，赵锐.流域水系自动提取的方法和应用［J］.湖泊科学，2003，15(3):205 －212.

［8］谢顺平，都金康，罗维佳，等.基于DEM的复杂地形流域特征提取［J］.地理研究，2006，25(1):96－102.

［9］祝士杰，汤国安，贾旖旎.基于DLG水系的DEM修正方法研究［J］.现代测绘，2009，32(6):3 －4.

［10］徐涛，胡光道.基于数字高程模型自动提取水系的若干问题［J］.地理与地理信息科学，2004，20(5):11 －14.

［11］吴良超.基于DEM的黄土高原沟壑特征及其空间分异规律研究［D］.西安:西北大学，2005.

［12］李俊，汤国安，张婷，等.利用DEM提取陕北黄土高原沟谷网络的汇流阈值研究［J］.水土保持通报，2007，27(2):75－78.

［13］陈涛，朱金兆，马浩，等.黄土高原沟壑区沟道提取时阈值范围的研究［J］.长江科学院院报，2008，25(3):28 －31.

［14］周贵云，刘瑜，邬伦.基于数字高程模型的水系提取算法［J］.地理学与国土研究，2000，16(4):77 －81.

［15］周毅.基于DEM的黄土正负地形特征研究［D］.南京:南京师范大学，2008.