留言板

尊敬的读者、作者、审稿人, 关于本刊的投稿、审稿、编辑和出版的任何问题, 您可以本页添加留言。我们将尽快给您答复。谢谢您的支持!

姓名
邮箱
手机号码
标题
留言内容
验证码

利用特征显著度提取地形特征线的方法

邹昆 沃焱 徐翔

邹昆, 沃焱, 徐翔. 利用特征显著度提取地形特征线的方法[J]. 武汉大学学报 ● 信息科学版, 2018, 43(3): 342-348. doi: 10.13203/j.whugis20150373
引用本文: 邹昆, 沃焱, 徐翔. 利用特征显著度提取地形特征线的方法[J]. 武汉大学学报 ● 信息科学版, 2018, 43(3): 342-348. doi: 10.13203/j.whugis20150373
ZOU Kun, WO Yan, XU Xiang. A Feature Significance-Based Method to Extract Terrain Feature Lines[J]. Geomatics and Information Science of Wuhan University, 2018, 43(3): 342-348. doi: 10.13203/j.whugis20150373
Citation: ZOU Kun, WO Yan, XU Xiang. A Feature Significance-Based Method to Extract Terrain Feature Lines[J]. Geomatics and Information Science of Wuhan University, 2018, 43(3): 342-348. doi: 10.13203/j.whugis20150373

利用特征显著度提取地形特征线的方法

doi: 10.13203/j.whugis20150373
基金项目: 

国家自然科学基金 61502088

国家自然科学基金 61300095

广东省自然科学基金 S2013010012307

广东省自然科学基金 S2012040011123

广东省高等学校优秀青-教师培养计划 Yq2013206

详细信息
    作者简介:

    邹昆, 博士, 副教授, 现从事地形特征提取、地形合成、纹理合成等研究。cszoukun@foxmail.com

  • 中图分类号: P208;P283

A Feature Significance-Based Method to Extract Terrain Feature Lines

Funds: 

The National Natural Science Foundation of China 61502088

The National Natural Science Foundation of China 61300095

the Natural Science Foundation of Guangdong Province S2013010012307

the Natural Science Foundation of Guangdong Province S2012040011123

the Cultivation Plan of the Outstanding Young Teachers in Guangdong Province Yq2013206

More Information
    Author Bio:

    ZOU Kun, PhD, associate professor, specializes in terrain feature extraction, terrain synthesis and texture synthesis. E-mail:cszoukun@foxmail.com

图(5)
计量
  • 文章访问数:  3232
  • HTML全文浏览量:  105
  • PDF下载量:  445
  • 被引次数: 0
出版历程
  • 收稿日期:  2016-01-15
  • 刊出日期:  2018-03-05

利用特征显著度提取地形特征线的方法

doi: 10.13203/j.whugis20150373
    基金项目:

    国家自然科学基金 61502088

    国家自然科学基金 61300095

    广东省自然科学基金 S2013010012307

    广东省自然科学基金 S2012040011123

    广东省高等学校优秀青-教师培养计划 Yq2013206

    作者简介:

    邹昆, 博士, 副教授, 现从事地形特征提取、地形合成、纹理合成等研究。cszoukun@foxmail.com

  • 中图分类号: P208;P283

摘要: 针对现有的地形特征线提取方法不能根据特征强弱对提取结果进行方便的筛选控制的问题,提出了一种基于特征显著度的山脊和山谷线提取方法。首先,对数字高程模型数据进行高斯滤波和下采样;然后,通过全局断面扫描确定特征点,并根据高度、高度落差和坡度计算特征显著度,在进行特征延伸后构造特征图,利用最小生成树算法得到特征森林;随后进行枝干分解并计算各枝干的特征显著度,根据特征显著度进行枝干裁剪;最后,通过尾钩剔除、位置修正和平滑等后处理得到最终特征线。实验结果表明,该方法能提供给用户方便有效的手段以根据特征显著度对特征进行筛选控制,提取的特征线与实际地形相符,有良好的抗噪能力且能较好地处理平坦地形特征。

English Abstract

邹昆, 沃焱, 徐翔. 利用特征显著度提取地形特征线的方法[J]. 武汉大学学报 ● 信息科学版, 2018, 43(3): 342-348. doi: 10.13203/j.whugis20150373
引用本文: 邹昆, 沃焱, 徐翔. 利用特征显著度提取地形特征线的方法[J]. 武汉大学学报 ● 信息科学版, 2018, 43(3): 342-348. doi: 10.13203/j.whugis20150373
ZOU Kun, WO Yan, XU Xiang. A Feature Significance-Based Method to Extract Terrain Feature Lines[J]. Geomatics and Information Science of Wuhan University, 2018, 43(3): 342-348. doi: 10.13203/j.whugis20150373
Citation: ZOU Kun, WO Yan, XU Xiang. A Feature Significance-Based Method to Extract Terrain Feature Lines[J]. Geomatics and Information Science of Wuhan University, 2018, 43(3): 342-348. doi: 10.13203/j.whugis20150373
  • 山脊线和山谷线是非常重要的地形特征线,描述了最主要的地形骨架结构,其自动提取对于地形地貌研究有着重要意义,近年来也被用于基于样本的地形合成[1-2]

    从规则格网数字高程模型(digital elevation model,DEM)中提取山脊线和山谷线一直是一个研究热点,主要方法包括地表流水模拟分析法[3-5]、几何形态分析法[6-9]和基于图像处理技术的方法[10]等。地表流水模拟分析法提取结果整体性较好,但效率较低,且较难判定平坦区域的水流方向[4];几何形态分析法通常采用局部窗口对断面进行分析以确定特征点,简单高效,提取精度相对较高,但结果易受窗口大小影响[7],且对噪声敏感,特征线连续性较差。文献[9]通过降低特征点入选条件增强了特征连续性,但使得后续筛选工作量大幅增加;文献[10]利用图像处理中的灰度形态学算子提取山脊线和山谷线,取得了较好效果,但在平坦山脊和山谷处易产生断裂。上述方法均不能根据特征强弱对提取结果进行方便的筛选控制,此外基于规则格网DEM的特征线提取方法总体上抗噪性较差,且特征点连接成线时较为困难[11]。因此,本文提出了一种新的地形特征线提取方法,通过基于特征显著度的裁剪实现对提取结果特征显著程度的控制,通过全局断面扫描避免窗口参数选择,提高抗噪能力及对平坦地形特征的妥善处理能力,通过基于特征方向的特征延伸提高连通性。

    • 高斯滤波与下采样是对DEM数据的预处理操作。下采样一方面是为了保证利用本文方法最终得到的是单像素宽的特征线,另一方面可降低后续步骤的计算量。下采样时,在原DEM中沿水平/垂直方向每n列/行取1个样本点,n值可根据DEM分辨率和提取精度要求选取。为了提高下采样后DEM的保真度,同时过滤细微噪声,在下采样前先进行高斯滤波,滤波器窗口大小设为2n-1。

    • 断面极值法是一种简单高效的确定候选特征点的方法,通常将断面上的高程极大值点作为山脊线上的候选点,将极小值点作为山谷线上的候选点[8]。然而现有方法通常采用局部窗口进行断面分析,分析时所利用信息有局限性,效果易受窗口大小和噪声影响。例如图 1(a)给出一DEM断面示意图,人眼易判断出A处为平原,B处为小平台,CE两处为山脊,D处为山谷,DF处的小突起为伪特征,G处可能为山谷,但由于过于靠近边界且右方只有微小上升,无法确定。但若使用局部窗口进行分析,如果窗口过小,则会丢失C处中间的山脊点,而D处和F处的小突起易被误判为山脊点,AB处的平坦特征也可能被误判。造成此种差异的原因在于人眼是利用全局信息进行判断的,受此启发本文设计了基于全局断面扫描的特征点确定方法。

      图  1  全局断面扫描示意图

      Figure 1.  Illustration of Global Scan on a Terrain Profile

      为了在追求高效的同时尽量不遗漏特征,分别沿水平、垂直及两个斜45°方向对DEM数据进行逐行扫描。虽然对方向进行进一步细分还可降低特征的漏检率,但随着断面数的增多,计算量和误检率也会增加[8]。针对每一行高程数据处理步骤如下。

      1) 提取关键点。逐个访问每一数据点pi,如果是首尾点则直接加入关键点集合K,否则通过计算该点与相邻点的高程差值判断该点前后的地形走势,有3种情况:

      (1) 前后地形走势同为上升或下降,则不作任何处理。

      (2) 前后地形走势相反,则将此点加入K

      (3) pi位于平坦地势起点,则不断访问后续点直至高程差值不为0,从而获知平坦地势长度l及后续走势。如果l大于阈值W,或是前后走势相反,则将pi加入K。其中,W用于鉴别平坦线型特征(例如平顶山脊)与平坦区域特征(例如平原)。

      由上述步骤得到的关键点主要是地形断面上的折点,但剔除了与图 1(a)B处类似的平坦地势点,此外将断面首尾点加入K主要是为了方便后续分析。针对图 1(a)提取的关键点ki图 1(b)所示,绿色线给出了对应关键点的长度(非平坦地势起点的长度为0)。

      2) 剔除噪声和伪特征,提取候选特征点。方法是遍历K中各点,根据当前地形总体走势T(分为不确定走势TU、上升走势TA、下降走势TD三种)及相邻关键点的高度差进行分析判断。首先将首个关键点直接加入候选特征点集合C,初始化TTU,然后逐个访问后续关键点,对每一关键点ki进行如下操作(设前一候选特征点为c,前一关键点为ki-1,噪声阈值为V):

      (1) 如果kic高度相同,则视二者处于同一平坦地势,修改c的长度,如果长度大于W,置TTU,转访问下一关键点。

      (2) 如果ki的长度大于W,置TTU,将ki加入集合C,转访问下一关键点。

      (3) 判断当前T的值,有3种情况:

      T等于TU。如果kic的高度差Δh的绝对值大于V,将ki加入集合C,并根据Δh符号修改Th>0则置TTA,Δh < 0则置TTD)。

      T等于TA。如果ki的高度大于ki-1的高度,且也大于c的高度,说明c为伪特征点,将其替换为ki;如果ki的高度小于ki-1的高度,且cki的高度差大于V,则将ki加入集合C,并置TTD

      T等于TD。如果ki的高度小于ki-1的高度,且也小于c的高度,说明c为伪特征点,将其替换为ki;如果ki的高度大于ki-1的高度,且kic的高度差大于V,则将ki加入集合C,并置TTA

      对于末尾关键点,可直接加入集合C图 1(c)给出了提取的候选特征点ci图 1(b)中的k5k8k9被当作噪声剔除,而k6由于和前一候选特征点c4高度相同,合并到c4,延长了c4的长度。

      3) 提取特征点。对C进行一次遍历,将长度小于W且与两侧相邻候选特征点的高度差绝对值均大于V的候选特征点加入特征点集合F。不用判断两侧高度差符号相反是因为步骤2已保证这一点。图 1(d)给出了提取的特征点fi

      4) 计算特征显著度。特征显著度即山脊/山谷特征点的特征显著程度,其值越大,表明该点越可能在一条更显著的山脊/山谷线上。海拔越高、与两侧区域高度差越明显的地形特征越显著,结合人眼对地形特征线显著程度的这一感知特性,通过考虑高度、与两侧区域高度落差以及两侧的坡度来计算特征点的显著度。

      具体而言,对于山脊点的显著度计算,考虑以下因素。

      (1) 归一化高度。高度越高,山脊特征越显著。

      (2) 平均归一化高度落差。两侧高度落差越大,山脊特征越显著。高度差取的是与相邻候选特征点的高度差,而非与相邻特征点的高度差,因为高原上山脊的高度落差应计算山脊点与高原的海拔差值,例如图 1(d)f1的左侧落差Δhl1是其与图 1(c)c2的高度差值,如图 2所示。

      图  2  特征点属性及坡面角示意图

      Figure 2.  Illustration of Feature Point's Attributes and Slope Angle

      (3)平均归一化坡度。沿断面扫描方向的特征点两侧坡度,坡度越大,山脊特征越显著。为了简化计算,采用坡面角的正弦值作为归一化坡度,例如图 2f1两侧的坡面角分别为αl1αr1

      由此得到山脊点特征显著度S计算公式:

      $$ \begin{array}{l} {\rm{S}}({f_i}) = {k_h}\frac{{{h_i}}}{{{h_{{\rm{max}}}}}} + {k_d}\frac{{\left| {\Delta {h_{li}}} \right| + \left| {\Delta {h_{ri}}} \right|}}{{2\left| {\Delta {h_{{\rm{max}}}}} \right|}} + \\ \frac{{{k_s}}}{2}(\frac{{\left| {\Delta {h_{li}}} \right|}}{{\sqrt {\Delta {h_{li}}^2 + {d_{li}}^2} }} + \frac{{\left| {\Delta {h_{ri}}} \right|}}{{\sqrt {\Delta {h_{ri}}^2 + {d_{ri}}^2} }}) \end{array} $$ (1)

      式中, hi为特征点fi的高度,; Δhli和Δhrifi与相邻候选特征点的高度差; dlidrifi与相邻候选特征点的距离; hmaxhmin分别为最大和最小高度; Δhmax为最大高度落差,取值为hmax-hmin; khkdks为权重系数。特征点的相关属性可参见图 2

      山谷点的显著度计算考虑的因素类似,计算公式为:

      $$ \begin{array}{l} {\rm{S}}({f_i}) = {k_h}\frac{{{h_{{\rm{max}}}}-{h_i}}}{{{h_{{\rm{max}}}}}} + {k_d}\frac{{\left| {\Delta {h_{li}}} \right| + \left| {\Delta {h_{ri}}} \right|}}{{2\left| {\Delta {h_{{\rm{max}}}}} \right|}} + \\ \frac{{{k_s}}}{2}(\frac{{\left| {\Delta {h_{li}}} \right|}}{{\sqrt {\Delta {h_{li}}^2 + {d_{li}}^2} }} + \frac{{\left| {\Delta {h_{ri}}} \right|}}{{\sqrt {\Delta {h_{ri}}^2 + {d_{ri}}^2} }}) \end{array} $$ (2)

      与式(1)的区别仅在于归一化高度,此时高度越低,特征越显著。

      可根据Δhli或Δhri的符号判定fi是山脊点还是山谷点。对于长度l大于0的平坦地势特征点f,将其后的l个点均作为特征点,并将它们的特征显著度设为与f相同,但将处于中间位置的特征点显著度加上一个微小增量,从而引导最终生成的单像素特征线的位置处于平顶山脊或平底山谷的中央。

      由于需要进行4次全局扫描,对于每个DEM样本点,保存其特征显著度最大值,并记录对应的扫描方向。

    • 虽然沿4个方向进行了全局断面扫描,但仍可能存在局部特征断裂,尤其是在显著特征线与其分支特征线(例如主脉与支脉)的连接处附近,这部分区域地势起伏相对平缓,特征点容易被当作噪声过滤掉。

      对于每一个特征点尝试沿该点的特征方向进行延伸,特征方向与断面扫描方向垂直。如图 3所示(其中pi表示3×3邻域内的DEM数据点),对于特征点p0(假设为山脊点)尝试沿左上和右下两个方向延伸,以沿左上方向延伸为例:

      图  3  特征延伸示意图

      Figure 3.  Illustration of Feature Extension

      1)考察点p5p6p7,如果其中有一个点已为特征点,则无需沿此方向延伸,否则继续。

      2) 依次考查p6p5p7,在p0的8邻域内判断其高度是否不小于相邻点(即为局部最大值),如是则将pi设为特征点,并将其显著度设为0,特征方向设为p0pi连线方向。

      3) 将新增特征点pi设为待延伸点,并根据其特征方向尝试进行进一步延伸。

      对于山谷点的延伸,将步骤2)中的局部最大值判断改为局部最小值判断即可。

      本文采用的延伸方法与文献[10]类似,但是在特征点的筛选之前进行,以更好地保证连通性,而由于延伸得到的特征点的显著度为0,在后面基于显著度的裁剪中很容易将“过度延伸”的特征线(指主要由延伸特征点构成的特征线)裁剪掉。

    • 在得到特征点后,构造特征图:将特征点作为图的顶点,将相邻特征点连线构成边,边的权值赋为两个特征点的显著度之和,如出现两条边交叉的情况,则剔除权值低的边。本文构图方法与文献[9]类似,只是权值不是采用高度而是采用特征显著度。

      接着需要打破特征线环路,生成特征森林,文献[9]的多边形裂开算法效率很低,本文采用了文献[2]的方法:利用Kruskal最小生成树算法生成特征森林。

    • 特征森林中通常存在较多干扰特征线,文献[6, 9]通过剔除长度较短的特征线段来消除干扰,然而长度并不能客观地反映特征的强弱,显著特征线的两端也可能存在较短的显著特征段,因此本文采用基于特征显著度的特征线裁剪剔除干扰特征线,其中将特征线的显著度定义为其包含的所有边的权值之和。

      首先进行枝干分解,即从特征森林中提取特征线,图 4给出了枝干分解的示意图(其中pi表示叶子节点或分支节点,特征线越粗表示其显著度越高),分解结果是p0p1p2p3p4p7p6p54条特征线。分解方法是遍历所有叶子节点(度数为1的特征点),对每一叶子节点pi做如下处理。

      图  4  枝干分解示意图

      Figure 4.  Illustration of Branch Decomposition

      1) 如pi属于一个未访问过的特征树,通过两次深度优先遍历找到该树的最长路径(即显著度最大路径)[12],作为该特征树的主干,对应主脉或干流。

      2) 如pi不在主干上,则从其开始进行路径搜索,直至遇到分支节点pj(度数大于2),提取特征线pipj,并对pj进行考察:如果pj不在主干上,且与pj相连的所有叶子分支(分支一端为叶子节点)均已提取,则选取以pj为端点的显著度最大的叶子分支进行延伸,直至遇到下一分支节点pk,并对pk进行与pj相同的考察和处理。

      对于分解得到的所有特征线,根据其显著度降序排序,然后可根据用户指定的阈值进行特征线的裁剪。

    • 将特征线上特征点的坐标乘以下采样的n值即可得到原图中的特征位置。为了提高特征线的质量,对其进行以下后处理。

      1) 由于下采样丢失了部分信息,特征点可能并非断面上的极大值点,因此对其进行位置修正:沿与特征方向垂直的断面方向,以(n-1)/2为半径搜索高度最大值点代替当前特征点。

      2) 部分特征线的末尾存在不自然的短钩,可通过计算曲率和到尾端的距离进行判断并移除。

      3) 为了使特征线更加平滑和自然,采用文献[9]方法对特征线进行平滑处理,即用相邻特征点的坐标均值作为当前特征点的坐标。

    • 选用某山区30 m DEM数据验证本文方法的有效性,如图 5(a)所示,其含有丰富的山脊、山谷特征以及噪声,在左方主要山脊交汇处存在较宽的平顶山脊特征。实验中参数选择为:下采样n值取3,但在生成图 5(c)图 5(d)时为了使结果显示更加清晰,未进行下采样;阈值W取15,阈值V取地形最大高度差的5%;在计算特征显著度时,权重系数khkdks分别取3、2、1。

      图  5  实验结果

      Figure 5.  Experimental Results

      图 5给出了采用本文方法对实验区域进行特征线提取的各步骤结果,并与两种典型方法进行对比。由于山谷线提取与山脊线提取类似,实验结果展示以山脊线提取为主。

      图 5(b)图 5(c)分别在图 5(a)基础上叠加显示了文献[8]方法和本文方法提取的山脊特征点,分别用4种颜色显示了不同特征方向的特征点(红色对应南北方向,蓝色对应东西方向,绿色对应东南-西北方向,黄色对应东北-西南方向)。文献[8]方法基于梯度方向选取近似最佳断面方向(近似到本文所采用的4个扫描方向之一),然后将局部极大值点定为山脊特征点。而本文方法采用全局断面扫描,相对局部断面扫描更不易受噪声影响,且能更好地识别出平顶山脊,因此相对图 5(b)图 5(c)中的细碎特征要少很多,在显著山脊上的连通性也更好。此外也可看到,基于特征显著度得到的特征方向是符合地形特征实际的;在部分显著特征的两侧,存在不同特征方向的相邻特征,这是沿多个方向进行断面扫描的结果,其特征显著度相对较弱,在生成特征森林时会形成主要山脊旁的短枝(如图 5(g)所示),在后续的枝干裁剪中很容易被剔除。

      图 5(d)给出了本文方法在提取山脊特征点后进行特征延伸的结果(用黑色显示延伸得到的特征点),可以看到,特征延伸进一步增强了特征线的连通性。有部分特征线主要由延伸的特征点构成,由于这类特征点显著度为0,这些过度延伸的特征线在枝干裁剪中同样易被剔除。

      图 5(e)~5(g)给出了本文方法在提取山脊线时3种枝干裁剪的结果,分别根据特征显著度取前3、前100和全部特征线(即无裁剪)。由于特征线的显著度是由特征边的显著度累加而来,实际上考虑了高度、高度落差、坡度及长度4方面因素,从图 5中可以看到裁剪结果与人眼观测到的地形特征显著程度基本一致。然而,由于在进行枝干分解时追求主干显著度最大化,可能将主干两端的细小分支也加到主干上,从而使这些细小分支在枝干裁剪中无法裁剪掉,例如图 5(f)中左下角特征线的长钩状特征。可通过在进行枝干分解时考虑特征线段间的特征相似性来解决此问题,但会增加算法的复杂度。在实际应用中,所要保留的特征线的数量可由用户交互式选取或者根据需要预先设定,从而提供给用户根据特征强弱对提取结果进行筛选控制的手段。

      文献[9]方法是特征连通性保持得最好的方法之一。图 5(h)5(i)分别给出了文献[9]方法和本文方法提取的最终山脊线,其中文献[9]方法窗口大小设为11,本文方法提取的是显著度最大的前100条山脊线。可以看出,图 5(h)中多处存在近距离平行特征线,而实际上这些位置仅存在一条山脊,这说明文献[9]所采用的基于长度的裁剪等筛选操作并不足以剔除前期被纳入候选特征点的存在于真实特征点周围的大量伪特征点,而本文基于特征显著度的裁剪则很好地剔除了图 5(g)中的伪特征线。此外,本文方法所提取山脊线的连通性明显更好,尤其在显著山脊与其分支的连接处,虽然采用更大窗口可进一步提高文献[9]方法结果的连通性,但其平行伪特征线问题将更严重。

      图 5(j)给出了用本文方法提取的显著度最大的前100条山谷线,同样具有很好的连通性,且山谷线与实际地形比较吻合,其中有几处山谷线穿过了山脊,一方面是因为鞍点的存在,另一方面是预处理中下采样所带来的负面效果。

    • 传统的地形特征线提取方法不能根据特征强弱对提取结果的精细程度进行方便的控制,为此本文提出了一种新的规则格网DEM地形特征线提取方法。该方法利用高度、高度落差和坡度信息计算特征点的显著度,通过枝干分解及基于特征显著度的枝干裁剪实现对提取结果的控制。此外,通过全局断面扫描增强特征点识别的准确性与抗噪性,以及更好地处理平坦特征。实验结果证实了本文方法的有效性,在可对提取特征显著程度进行方便控制的同时,提取的特征线与人工识别结果有较强一致性。但枝干分解算法还可进一步改进,通过增强同一枝干上不同分段的相似性来提高分解结果的准确性。

参考文献 (12)

目录

    /

    返回文章
    返回