基于数字高程模型的淤地坝地区河网提取方法

张宏鸣 樊世豪 陈茹雪 鱼文虎 张 琦 袁琳琳

(西北农林科技大学信息工程学院,陕西杨凌 712100)

0 引言

从数字高程模型(Digital elevation model,DEM)中提取河网是数字流域[1]、分布式水文模拟[2-3]、土壤侵蚀模型建立[4]等的关键技术环节,可为进一步划分流域和计算汇流面积等其他水文参数提供基础[5]。但在黄土高原土流失严重的地区,淤地坝作为重要而独特的滞洪、拦淤、蓄水工程被广泛修筑于开阔的沟谷区域,在该地区利用DEM进行流域水系分析时,直接提取河网会受到干扰[6]。原始DEM中存在许多中心栅格高程比周围栅格低的洼地[7]。在这些洼地内,水流不能确定正确的流向,导致水流路径不连续甚至出现逆流现象[7],进而影响河网的正常提取。因此,国内外研究人员提出了许多DEM预处理方法以改善河网的流动路径,从而实现河网提取[8-11]。

DEM数据的预处理是通过修正洼地,使得每个栅格都有明确的水流方向。目前,洼地的处理方法主要分为填充法和裂开法[7-8]。填充法通过抬升洼地中的栅格来修改地形,直到匹配最低出口栅格的高程[8,11-13]。而裂开法则是通过降低洼地边缘处阻挡物的栅格高程,来模拟水流下切裂开洼地[8,14]。虽然这些方法能保证河网的连续性,但在淤地坝地区的河网提取结果存在河网偏移问题。由于淤地坝两侧道路高程较低,河网将其作为洼地出水口穿过淤地坝,而到达出水口位置之前,河网流向与淤地坝的横截面方向基本保持平行。因此针对淤地坝流域,需开发一种新的DEM预处理方法提取河网。

CARLSON等[15]提出一种桥梁识别并修正其栅格高程的河网提取方法,与填充法相比,这种先识别后处理的方法能建立更真实的河网。然而,位于山地的梯田具有相似的地形特征,对淤地坝检测及河网提取造成负面影响[16]。ZHU等[17]提出利用区域生长和形态学处理提取黄土高原沟道,虽然该方法并没有考虑黄土高原地区的淤地坝,但能将淤地坝的检测范围缩小至河道区域,从而避免梯田的影响。此外,SOFIA等[18]提出使用直线段检测(Line segment detector,LSD)算法[19]提取具有线性特征的梯田,为检测具有类似线性特征的淤地坝提供了一个新的思路。

综上所述,本文提出一种基于线性特征检测淤地坝并修正其所在高程的方法。首先综合淤地坝和河道的地形特征进行区域生长和形态学操作提取河道中心线,并结合改进LSD算法检测结果保留淤地坝轮廓直线,最后构建十字模型检测淤地坝,修正所在高程并填充剩余的洼地以提取河网。

1 材料和方法

1.1 研究区域

如图1所示,本文选择的两块实验区均位于黄土高原地区,且具有不同的地形特征。实验区1为王茂沟流域,位于陕西省绥德县韭园沟中游左岸,是我国最早的治理实验小流域之一[20]。该流域地貌复杂,以梁、峁为主,沟壑纵横,坡陡沟深,属典型的黄土丘陵沟壑地貌。多年平均降水量513.1 mm,汛期7—9月降水量占年降水量的73.1%。经过近70年治理,坝系布局较为完整,存在淤地坝淤满的情况,坡耕地大幅减少,梯田面积和植被覆盖度提高。实验区2为周屯沟流域的一条支流,位于陕西省延安市安塞区。流域内地貌复杂多样,境内川道狭长、梁峁遍布,土壤以黄绵土为主,物理风化作用强,土壤疏松,多粉沙,易受水力侵蚀。多年平均降水量500 mm左右,汛期7—9月降水量占年降水量的74%。该流域淤地坝和梯田相对较少,浅沟区域治理以种植草灌为主要措施。图1中两流域对于淤地坝区域的河网提取和分析,具有典型的代表意义。

图1 研究样区

1.2 实验方法

技术路线如图2所示,淤地坝地区的河网提取主要为:将无人机航测分辨率0.15 m DEM采用最邻近插值法生成高分辨率DEM和低分辨率DEM。将低分辨率DEM提取河道中心线,高分辨率DEM导出山体阴影图并提取轮廓直线。依据河道中心线进行角度滤波保留淤地坝轮廓直线,构建十字模型定位淤地坝。修正淤地坝栅格高程并填充剩余洼地来提取河网。

图2 技术路线图

1.2.1DEM平滑处理

高分辨率DEM提供详细的地形信息,但其中存在的大量噪声不仅会增加计算复杂性,还会对河道提取产生负面影响。本文主要根据河道两侧高、中间低的地形特征提取河道,不过多关注河道内的地形细节。有研究表明,分辨率10 m DEM的坡度统计百分比与分辨率5 m DEM相比变化不大,足以满足地形和水文建模应用的需求[21-23]。分辨率更低的DEM会导致更多地形信息的丢失,使得区域生长结果无法有效反映河道的实际情况,因此本方法使用分辨率10 m DEM数据进行区域生长。分辨率10 m DEM仍存在少部分噪声,为了去除这些噪声并保持河道和山谷的边缘信息,本文使用中值滤波窗口(3×3)进行滤波。如图3所示,滤波能有效去除河道内部噪声,且能降低部分淤地坝所在的栅格高程,从而有利于后续区域生长提取连续的河道。

图3 滤波前后坡度对比

1.2.2河道中心线提取

对低分辨率DEM的洼地进行填充,并用最陡坡降计算流向和流量。经过中值滤波和填充处理后,河道的整体地形相对平坦。利用坡度和高程之间的关系进行区域生长得到河道,坡度S定义为

(1)

式中H——相邻栅格的高程差

L——相邻栅格之间的水平距离

为每个栅格设置生长状态,将未标记为河道且大于流量累积阈值rflow的流量栅格作为种子点,以避免流量栅格的重复生长,提高算法效率。区域生长的生长准则为

(2)

ha>hb

(3)

式中h——栅格高程E——坡度阈值

栅格a向相邻栅格b生长,如果满足上述两条件,则将其标记为河道栅格,直到确定河道边界位置。式(2)能确保淤地坝坝顶以及河道等相对平缓地区的生长,而由于淤地坝两侧坡度变化较大,利用式(3)可保证位于坝顶的生长点向淤地坝上下游进行区域生长。不同坡度阈值的生长结果如图4所示,当坡度阈值设置为3°时,区域生长无法覆盖大型淤地坝,导致细化后的河道中心线发生偏移,进而影响后续角度滤波结果的准确性。当将坡度阈值设置为4°时,可以保证河道的连续性和完整性。然而,当阈值设置为5°时,沟头区域的生长结果可能不准确,因为较大阈值可能使生长扩展到坡度较缓的沟头区域。此时,沟头区域的高程较高,而式(3)会将生长过程扩展到所有低于沟头栅格高程的区域,进而导致生长结果的不准确。因此,为确保河网连续性和准确性,取阈值E为4°。

图4 不同坡度阈值的生长结果

对分辨率10 m DEM的河道生长结果进行形态学闭运算以填补河道中未生长的小区域,并平滑河道边界。然而,在分辨率10 m下,淤地坝的特征无法有效表示,因此需用高分辨率DEM来检测淤地坝。为与高分辨率DEM嵌套来定位淤地坝,将河道生长结果上采样成高分辨率河道数据,最后利用细化算法生成单栅格大小的河道中心线。

1.2.3淤地坝特征提取

如图5所示,为提高淤地坝的提取精度,本文总结淤地坝的4个特征:①淤地坝修建在截面两侧具有陡峭山地的河道上。②淤地坝与河道的流向近乎垂直。③淤地坝为横向线性结构坝工建筑物。④淤地坝的坝顶、上游和下游地势相对平坦。

图5 淤地坝区域示意图

由于淤地坝整体呈线性,本文使用LSD算法对分辨率1、2、5、10 m的山体阴影图进行轮廓直线提取。如图6所示,分辨率10 m山体阴影图只能检测到河道两侧的轮廓直线,分辨率5 m虽能检测到大型淤地坝的轮廓直线,但对于中小型淤地坝的检测并不理想。在分辨率2 m下,能很好地检测各种规模淤地坝的轮廓直线。而在分辨率1 m下,检测到的淤地坝轮廓直线与分辨率2 m相比并无较大区别,但河道内的噪声轮廓直线数量却大大增加,使得提取淤地坝轮廓直线的方法更为复杂,检测效果也会受到影响。因此,本文选择分辨率2 m DEM进行淤地坝检测,以获得较好的检测效果。

图6 不同分辨率LSD检测结果示意图

1.2.4有效轮廓直线提取

利用淤地坝与河道近似垂直的位置关系,对轮廓直线进行角度滤波,以保留与河道中心线相交且近乎垂直的轮廓直线。如图7所示,淤地坝横跨河道而建,淤地坝的轮廓直线与河道的轮廓直线可能相邻。为了正确计算淤地坝轮廓直线的倾斜角,本文改进LSD算法分离倾斜角不同的直线。改进方法如图7a所示,由于检测到的轮廓直线按先后顺序输出,若差值在0°～5°范围内,则认为这两条直线属于同一直线,保留栅格F并接着检测下一栅格是否满足该条件;相反则将该栅格丢弃,以分离不同倾斜角的轮廓直线。

分离轮廓直线后,利用最小二乘法拟合轮廓直线l与河道中心线j的倾斜角,分别为θl和θj。如果这两倾斜角的差值满足

|θl-θj|∈[90°-α,90°+α]

(4)

则认为直线l属于淤地坝轮廓直线;反之则认为直线l为噪声轮廓直线并丢弃。由于淤地坝和河道中心线并非完全垂直,且基于栅格的倾斜角计算存在误差,本文将α设置为30°,尽可能避免淤地坝轮廓直线的丢失及噪声轮廓直线的影响。

1.2.5十字模型检测

角度滤波后仍存在少部分的噪声直线,需对每条轮廓直线构建十字模型以确定最后淤地坝的位置。由于多条轮廓直线表示同一个淤地坝,将其合并以避免十字模型的重复构建。合并的判断依据:若距离较近的两直线倾斜角差值小于5°,则合并两直线。如图8所示,确定合并之后的轮廓直线中点U,则有

(5)

k——垂直淤地坝横截面的直线斜率

根据淤地坝坝顶和上下游区域地势平坦的特征,设置平缓值SA以确定垂直淤地坝横截面直线的延伸长度,平缓值SA只根据相邻两栅格整型高程是否相等来变化。

建立十字模型检测淤地坝具体步骤为:

(1)将轮廓直线中点U作为垂直淤地坝横截面直线的起点,并存储在链表的第1个节点。利用式(5)计算垂直淤地坝横截面的直线斜率k,其中标准高程ch初始化为中心点U的高程,平缓值SA初始化为3,可根据分辨率进行调整。

(2)根据Bresenham算法沿斜率k确定下一点的坐标u1。判断u1与ch的整形高程值是否相等,若相等且SA=0,则停止向该方向的延伸并执行步骤(3),否则将u1存储至链表,并把ch更新为u1的高程;相反若两栅格高程相等且SA>0,则将SA减1,并重复步骤(2),计算下一个点u2。

(3)根据反向斜率-k,重复步骤(1)～(2)。

(4)根据垂直淤地坝横截面直线所对应的高程数据,确定高程最高的栅格坐标M(x,y)为淤地坝坝顶。

(6)根据垂直淤地坝横截面直线上的每个栅格,分别向横截面的对应方向延长lenl和lenr来构造淤地坝所对应的矩形。

(7)利用形态学闭运算来填补遗漏的栅格。

1.2.6淤地坝高程修正

确定淤地坝位置后,使用侵蚀操作符修正淤地坝高程。侵蚀操作符是一个9×9的窗口,以淤地坝所在栅格为窗口中心,窗口内的高程最小值替代为窗口中心的高程。根据经验,连续5个侵蚀操作符能完全处理淤地坝高程,这些参数可根据分辨率进行调整。除淤地坝外,一些因其他原因形成的洼地仍然存在,使用填充法处理剩余的洼地,最后采用D8算法提取河网。

1.3 结果对比

本研究将所述方法与LINDSAY[24]提出改进的填充法(填充法)和最小代价裂开法(裂开法)进行对比,其中裂开法的最大搜索距离设为100 m。填充法对填充洼地后的平坦地区进行微小抬升来避免平行流向问题,而裂开法提取的河网更自然,能够在合适的位置穿过洼地[24]。

为了评估淤地坝的检测效果,将本方法应用于两处典型的实验区,并与人工标记的淤地坝进行比较。采用精确率(P)、召回率(R)和F1值来验证淤地坝检测的有效性。此外,在相同河网阈值的条件下,将本文与填充法、裂开法提取的河网结果进行比较。

2 结果与讨论

2.1 淤地坝检测结果

在2个实验区上进行淤地坝检测,实验结果如表1所示。

表1 实验区1和实验区2的实验结果

该方法在2个实验区中检测淤地坝的F1值分别为86.67%和86.95%,结果表明该方法能较为准确地检测淤地坝。淤地坝的检测结果如图9所示,所述方法中改进LSD算法,确保在进行角度滤波时能正确计算淤地坝轮廓直线的倾斜角,且角度滤波操作能有效避免位于山地中梯田轮廓直线的干扰。此外,淤地坝存在淤满和未淤满的两种地形特征,如图8所示上游虚线部分表示淤地坝已蓄满,其上游高程近似等于坝顶高程。本文构建的十字模型能有效识别这两种情况的淤地坝,因为十字模型中设置的平缓值SA仅用于判断淤地坝上下游高程是否平缓,而不考虑其上下游和坝顶之间的高程关系,从而提高了淤地坝的整体检测效果。本方法利用线性特征检测黄土高原地区的淤地坝,但该方法也可拓展到具有相似线性特征的建筑物,如桥梁、道路等。

图9 淤地坝检测结果

由于淤地坝地区地形复杂,检测过程中可能存在淤地坝误检和漏检的情况。位于河道交汇处两侧的道路轮廓直线可能与河道中心线近似垂直,并且道路与淤满淤地坝的特征相似,导致误检情况的发生。此外,溃坝或淤地坝周围的植被可能会破坏淤地坝的地形特征,导致部分淤地坝漏检。当淤地坝上下游存在植被时,植被高度会使所在区域的高程整体高于真实地面高程,垂直误差高达10 m以上,这些高程更高的植被栅格会代替原本的坝顶,导致在十字模型建模时会因横截面高程差过大而将被舍弃。

2.2 河网提取结果对比

仅使用填充法提取的河网结果误差最大,尤其是位于淤地坝上游区域的河网会出现大角度的偏移。这是因为淤地坝两侧的道路栅格高程较低,可能被误判为洼地的出口栅格,导致河网沿着两侧的道路穿过淤地坝。在到达出口栅格之前,河网流向与淤地坝的横截面方向基本保持平行。裂开法提取的河网在地形平缓的实验区1中能较为准确地穿过淤地坝,但相较于填充法的河网提取结果更为曲折,此外,在复杂地形中仍存在河网偏移问题,这表明裂开法的效果并不理想。

如图10所示,本文通过先修正淤地坝所在高程,再使用填充法提取河网,能有效解决河网偏移问题。相较于仅使用填充法直接改变整个淤地坝上游高程,修正后的淤地坝上游不存在洼地,提取的河网结果能更好地反映淤地坝上游的真实地形。与裂开法的河网提取结果相比,本文方法的河网提取结果更为平滑,且能有效地穿过淤地坝。在平坦河道中,这3种方法都存在共有的问题,即当多条支流汇入主干河网时,支流可能会沿着河道穿过一段距离后才汇入主干河网,这主要由于平坦区域的河网流向分配不理想所导致的。

图10 河网检测结果对比

2.3 淤地坝检测误差分析

淤地坝存在误检和漏检的情况,需要使用填充法来填补剩余的洼地。如图11a所示,因溃坝而导致淤地坝漏检时对河网的影响最小,因为破损淤地坝上游并未形成洼地,河网会沿破损处穿过淤地坝。如图11b所示,当受淤地坝下游植被影响而导致漏检时,河网提取结果仍存在偏移问题,而由于淤地坝检测的精确率较高,整体的河网提取结果仍能得到有效的改善。图11c中误判为淤地坝区域的最小高程替代成整个区域的高程,改变了原始区域的地形,导致该区域的河网发生轻微偏移。将来可利用高程关系判断淤地坝的上下游,在其两侧确定一条合理的单栅格路径,并只修正路径所在的栅格高程,这样既能减少对DEM的修改,又能尽可能降低误判对河网提取结果的影响。

图11 淤地坝漏检、误检的影响

3 结束语

提出了一种通过对淤地坝进行检测并修正其高程的自动化方法来实现河网的有效提取。在2个典型的实验区上进行验证,结果表明,该方法能较为准确地检测不同规模的淤地坝,并能有效解决河道偏移问题。然而,该方法仍有一定的局限性,例如溃坝或淤地坝上下游植被的影响可能会导致漏检和误检,使提取的河网结果发生轻微偏移。与目前的河网提取方法相比,本文方法提取的河网结果更符合真实地形,为黄土高原淤地坝地区的水土保持规划以及流域数字地形分析提供了良好的补充与支持。