CN104899865A - 基于高斯尺度空间的山脉线提取方法 - Google Patents

Abstract

本发明公开了一种基于高斯尺度空间的山脉线提取方法，其步骤为：1.以高斯核函数对规则格网数字高程模型数据做多尺度分解表达；2.从分解后的各尺度数据中提取山脊线；3.融合各尺度山脊线，标记出初选山脉特征点；4.设定高程阈值对初选山脉特征点作进一步筛选；5.细化连接处理，获得山脉线。本发明先以不同参数的滤波器组对规则格网DEM数据进行多尺度分解，从中分别提取各尺度下的山脊线，再对山脊线进行融合获得山脉线，这样既能有效去除山脊特征的冗余，又能保留山体的主要脉络。

Description

基于高斯尺度空间的山脉线提取方法

技术领域

本发明属于数字地形分析领域的宏观地形特征提取技术，具体涉及一种以高斯尺度空间理论对DEM数据进行分解，分别提取各尺度下的山脊线，进而以高斯尺度参数构造各尺度山脊线权重，融合生成宏观山脉线的方法。

背景技术

数字地形分析是在DEM数据上进行地形属性计算和特征提取的数字信息处理技术。山脉作为一种主要的地貌构架，由具有明显走向的若干条山岭组成，常成为生物地理分布和自然地理区域的界限。山脉线是山区地形的主要脉络，分布于高海拔区域，连贯突出且长度较长，控制了地形的整体特征和态势，在科学利用地形进行基础设施修筑、区域规划、灾害防治、军事指挥决策等方面具有重要指导意义。

目前地形特征线提取方法主要针对山脊线、山谷线、沟缘线、坡脚线等微观地形特征线，如何运用上述研究成果得到宏观地形特征线成为亟需解决的问题。贺文慧等人在文献“贺文慧,汤国安,杨昕等.面向DEM地貌综合的山脊线等级划分研究——以黄土丘陵沟壑区为例[J].地理与地理信息科学,2011,27(2):30-33.”中采用模拟流水法提取山脊线，再综合考虑分水线与沟谷线的耦合关系，对山脊线进行等级划分，得到连续性较长、控制地形整体特征的主脉。由于从树状沟谷线经对偶提取出的山脊线存在大量的孤立、断裂分枝，该方法提取出的伪山脊线较多，适用性有限。李建军等人在文献“李建军,李钊,李俊山等.基于DEM的地形特征提取算法[J].信号与信息处理,2008,38(2):26-29.”中先用断面极值法提取山脊线，再从地形遮挡的角度出发，提取出主山脊。该方法以局部窗口法分析并筛选山脊点，因而对噪声敏感，结果中存在较多断裂和分支。罗寅在文献“罗寅.宏观山脉线抽取及其辅助技术研究[D].西安:西安建筑科技大学,2012.”中将山脉线视为符合某些山脉特征的山脊线，通过对山脊线矢量化和分级，提取出平均高程高、坡度缓、起伏度小的山脊线段，形成山脉线。该方法提取结果稳定性欠佳，对山脊线的分级依据还有待研究。易玮在文献“易玮.基于DEM的地形特征提取算法研究及应用[D].西安:西安建筑科技大学,2012.”中依据鞍部点、山顶点以及山脉线间的关联关系，设计了山脉结构线的提取方案，该方法需要借助地形特征中的鞍部点和山顶点，提取的山脉线存在不连续现象，且与实际地形偏差较大。

综上所述，目前对微观地形特征线自动提取的研究成果颇丰，但关于宏观地形特征山脉线的提取方法很少，仅有的几种方法提取的山脉线偏差较大、连续性差，无法支持实际应用。

发明内容

针对现有山脉线提取方法的不足，本发明结合多尺度分析与特征融合的思想，提出一种山脉线提取方法。它通过高斯核函数对DEM数据进行尺度分解，从中提取山脊线，再采用加权平均法融合各尺度山脊线，形成山脉线。本发明既有效去除了原始山脊特征的冗余，又保留了山体的主要脉络，提取的山脉线清晰、连续，与实际地形较为相符。

为了实现上述任务，本发明采用以下技术方案：

一种基于高斯尺度空间的山脉线提取方法，包括以下步骤：

步骤一，以高斯核函数对规则格网DEM数据做多尺度分解表达；

步骤二，从分解后的各尺度数据中提取山脊线；

步骤三，融合各尺度山脊线，标记出初选山脉特征点；

步骤四，设定高程阈值对初选山脉特征点作进一步筛选；

步骤五，对筛选后的特征点进行细化、连接，获得山脉线。

进一步地，步骤一的具体过程包括：

步骤S10，确定分解层数n和第一层表达的尺度参数σ₀，并建立各尺度参数间的关系；为了体现尺度的连续性，令相邻两层DEM数据间的尺度参数呈k倍关系，即有：

σ_i＝kⁱσ₀

上式中，i＝0,1,...,n-1；σ_i为各尺度参数。

由上式可知，第一层尺度参数σ₀的选择是山脉线提取的关键。地形表面的粗糙程度与地貌类型密切相关，可选用地形的相对高度变化即地势起伏度来描述其变化剧烈程度。地势起伏度越大，山地越复杂，其细微变化的特征线就越多；为了得到宏观山脉线，此时应选用较大的初始尺度滤波参数σ₀。在对我国的主要山脉进行算法仿真实验后发现，σ₀的大致取值范围在之间。

步骤S11，依据各尺度参数σ_i(i＝0,1,…,n-1)对DEM数据进行n层分解。具体分解公式如下：

$\left\{\begin{array}{c}{K}_i(x,y,{\mathrm{\σ}}_i)=\frac{12}{{2\mathrm{\π\σ}}^2}{e}^{-(x^2+y^2)/\left({2\mathrm{\σ}}_i^2\right)}\\ g_i(x,y,{\mathrm{\σ}}_i)=K_i(x,y,{\mathrm{\σ}}_i)\⊗f(x,y)\end{array}\right.$ i＝0,1,…,n-1

式中，K_i(x,y,σ_i)为第i+1层高斯核，g_i(x,y,σ_i)为分解后的第i+1层DEM数据，表示卷积运算。

进一步地，步骤二的具体过程包括：

步骤S20，若将规格网DEM的行列坐标(x,y)视为图像的像素点坐标，该坐标处的高程值f(x,y)视为图像的灰度值，则DEM数据可渲染为灰度图像；显然地势凸出的山脊部分对应于渲染灰度图中的带状高亮区域。

步骤S21，为了提取山脊线，利用“提取地形特征线的形态学新方法”，孔月萍等人在文献“孔月萍,方莉,江永林等.提取地形特征线的形态学新方法[J].武汉大学学报信息科学版,2012,37(8):996-999.”中利用数学形态学原理，设计了具有高通特性的“高帽”算子，运用它获得规格网DEM中具有正脉冲峰值的带状高亮区域，标记出山脊线。该方法在山脊特征提取、运算复杂度、抗噪鲁棒性等方面取得了较好的平衡，获得的山脊线清晰、连贯，伪信息少；以该方法分别作用于分解后的各层DEM数据，可提取出n个尺度下的山脊线标记图，记为F_i(x,y,σ_i)∈{0,1},i＝0,1,…,n-1。

进一步地，步骤三的具体过程包括：

步骤S30，由于σ_i越大，高斯滤波后的信号特征越宏观，从中提取的山脊线越能代表主脉，则该尺度山脊线对运用融合法形成山脉线的贡献应越大，因此以归一化的规范化尺度参数作为各尺度山脊线融合时的权重w_i：

$w_i=\frac{{\mathrm{\σ}}_i^2}{\underset{j=0}{\overset{n-1}{\mathrm{\Σ}}}{\mathrm{\σ}}_j^2}$

采用加权平均法融合各尺度山脊线：对各尺度山脊线标记图赋予不同的权重w_i，并按照下式进行加权求和，获得融合结果A(x,y)：

$A(x,y)=\underset{i=0}{\overset{n-1}{\mathrm{\Σ}}}{w}_i{F}_i(x,y,{\mathrm{\σ}}_i)$

步骤S31，设定融合阈值TA，将A(x,y)大于该阈值的点标记为初选山脉特征点：A(x,y)越大，说明该点在山脉线上的概率越大，选定融合阈值TA，将A(x,y)≥TA的点标记为初选山脉特征点。建议TA的取值范围介于最大和次大尺度权重之间，即w_n-1＜TA≤w_n。

进一步地，在步骤四中，由于山体的主要脉络通常具有较高高程，故选择适当高程阈值TD对初选山脉特征点作进一步筛选。

进一步地，步骤五中，采用文献“孔月萍,方莉,江永林等.提取地形特征线的形态学新方法[J].武汉大学学报信息科学版,2012,37(8):996-999.”中的后处理法对筛选后的特征点进行细化、连接，获得山脉线。

与现有技术相比，本发明具有以下优点：

第一，山脉由具有明显走向的若干条山岭组成，而山岭中贯穿山顶、山背的最高凸棱线为山脊线，可见山脊线中包含了山脉线的位置和走向信息，但冗余了较多的短小支脉。本发明充分考虑了山脉线和山脊线的差异性与相似性，设计了从山脊线中提取山脉线的合理方法。

第二，高斯滤波可消除DEM数据中的局部起伏，但对地形整体态势和宏观结构影响不大，本发明以不同参数的滤波器组对规则格网DEM数据进行多尺度分解，从中提取各尺度下的山脊线，再采用加权平均法进行融合，获得山脉线，这样既有效地去除了山脊特征的冗余，又保留了山体的主要脉络。

第三，本发明考虑了各尺度宏观特性与其尺度滤波参数间的联系，依尺度参数分配各尺度特征权重，具有充分的合理性与可行性。

第四，本发明方法简单易行，提取的山脉线清晰、连续，与实际地形较为相符。

附图说明

图1为本发明方法的流程图；

图2(a)为某山区规则格网DEM数据的三维透视图；

图2(b)为某山区规则格网DEM数据的灰度晕渲图；

图3(a)为某山区高斯尺度空间分解后第一层数据提取的山脊线(为了便于观察，图中叠加了等高距为50米的等高线)；

图3(b)为某山区高斯尺度空间分解后第二层数据(σ₁＝2)提取的山脊线；

图3(c)为某山区高斯尺度空间分解后第三层数据提取的山脊线；

图3(d)为某山区高斯尺度空间分解后第四层数据(σ₃＝4)提取的山脊线；

图3(e)为某山区高斯尺度空间分解后第五层数据提取的山脊线；

图3(f)为本发明提取的某山区山脉线；

图4(a)为采用罗寅的方法提取的某山区山脉线；

图4(b)为采用易玮的方法提取的某山区山脉线；

图5(a)为昆仑山脉部分SRTM(Shuttle Radar TopographyMission)-DEM数据的三维透视图；

图5(b)为昆仑山脉部分SRTM-DEM数据的灰度晕渲图；

图5(c)为本发明对昆仑山脉部分SRTM-DEM数据提取的山脉线；

图6(a)为秦岭山脉部分SRTM-DEM数据的三维透视图；

图6(b)为秦岭山脉部分SRTM-DEM数据的灰度晕渲图；

图6(c)为本发明对秦岭山脉部分SRTM-DEM数据提取的山脉线；

图7(a)为太行山脉部分SRTM-DEM数据的三维透视图；

图7(b)为太行山脉部分SRTM-DEM数据的灰度晕渲图；

图7(c)为本发明对太行山脉部分SRTM-DEM数据提取的山脉线；

图8(a)为天山山脉部分SRTM-DEM数据的三维透视图；

图8(b)为天山山脉部分SRTM-DEM数据的灰度晕渲图；

图8(c)为本发明对天山山脉部分SRTM-DEM数据提取的山脉线。

具体实施方式

参照图1，以1：25万某山区DEM数据(数据格网精度为50米，高程范围为1322～1980米，地势起伏度为498米，属山地地形中的小起伏中山类别，其三维透视图及灰度晕渲图分别如图2(a)、2(b)所示)为例，本发明的具体实施方式如下：

步骤一，以高斯核函数对规则格网DEM数据做多尺度分解表达

步骤S10，选定分解层数n＝5，相邻两层DEM数据间的尺度参数比例系数依据该地区的地势起伏度选定第一层表达的尺度参数依据σ_i＝kⁱσ₀得各层高斯分解的尺度参数σ₀～σ₄为 $(\sqrt{2},\mathrm{2,2}\sqrt{2},\mathrm{4,4}\sqrt{2});$

步骤S11，依据各尺度参数σ_i(i＝0,1,…,n-1)对DEM数据进行n层分解。具体分解公式如下：

$\left\{\begin{array}{c}{K}_i(x,y,{\mathrm{\σ}}_i)=\frac{1}{{2\mathrm{\π\σ}}^2}{e}^{-(x^2+y^2)/\left({2\mathrm{\σ}}_i^2\right)}\\ g_i(x,y,{\mathrm{\σ}}_i)=K_i(x,y,{\mathrm{\σ}}_i)\⊗f(x,y)\end{array}\right.$ i＝0,1,…,n-1

步骤二，从分解后的各尺度数据中提取山脊线

将规格网DEM的行列坐标(x,y)视为图像的像素点坐标，该坐标处的高程值f(x,y)视为图像的灰度值，将DEM数据渲染为灰度图像；采用文献“孔月萍,方莉,江永林等.提取地形特征线的形态学新方法[J].武汉大学学报信息科学版,2012,37(8):996-999.”中的方法提取各尺度下的山脊线，结果如图3(a)～图3(e)所示。

步骤三，融合各尺度山脊线，标记出初选山脉特征点

步骤S30，确定各尺度山脊线融合时的权重：对各层规范化尺度参数进行归一化处理，即得到各尺度山脊线融合时的权重值w₀～w₄为依据公式 $A(x,y)=\underset{i=0}{\overset{n-1}{\mathrm{\Σ}}}{w}_i{F}_i(x,y,{\mathrm{\σ}}_i)$ 计算各尺度山脊线的融合结果A(x,y)；

步骤S31，选定阈值TA＝32/62，将A(x,y)≥TA的点标记为初选山脉特征点；

步骤四，设定高程阈值对初选山脉特征点作进一步筛选

因山体的主要脉络通常具有较高的高程，可认为山脉线特征点均不低于平均高程，故采用该区域的平均高程作为TD，对标记出的初选山脉特征点作进一步筛选。

步骤五，对筛选后的特征点进行细化、连接，获得山脉线

采用文献“孔月萍,方莉,江永林等.提取地形特征线的形态学新方法[J].武汉大学学报信息科学版,2012,37(8):996-999.”中的后处理法对筛选后的特征点进行细化、连接，获得山脉线，如图3(f)所示。

对比图3(f)与图3(a)～图3(e)可知，本发明方法去除了分布在主脉两侧的短小支脉，得到了较长的连续主脉，符合山脉线特征。

为了更加清楚地说明本发明的效果，发明人进一步提供了以下实验：

一、对比实验

采用文献“罗寅.宏观山脉线抽取及其辅助技术研究[D].西安:西安建筑科技大学,2012.”和“易玮.基于DEM的地形特征提取算法研究及应用[D].西安:西安建筑科技大学,2012.”中的方法对上述山区DEM数据进行山脉线提取，提取结果分别如图4(a)、图4(b)所示。与图3(f)对比可知，本发明方法提取的山脉线更加准确、完整、连续，说明本发明方法优于现有方法。

二、应用验证实验

由于地质构造的影响，我国的山脉呈现以下几种结构：

①我国西部地壳由于受印度板块向北漂移的强烈挤压作用形成巨大弧形山系，从喜马拉雅山脉起，层层向北推进，像波浪一样，主要包括唐古拉山山脉、巴颜喀拉山脉、昆仑山脉、阿尔金山脉、祁连山脉等；②我国东部，由于太平洋板块的强大作用，形成一系列北东-南西或北北东-南南西走向的山脉，主要包括长白山、大兴安岭、太行山脉、巫山、雪峰山等；③东西走向的山脉，主要包括阴山、秦岭、南岭、天山等。

为此，本发明选择了弧形山系的昆仑山脉、北东-南西走向的太行山脉、东西走向的秦岭山脉和天山山脉这四组具有代表性的山地区域SRTM-DEM数据对本发明方法进行测试。其相关参数如表1所示，三维透视图如图5(a)、图6(a)、图7(a)、图8(a)所示，灰度晕渲图如图5(b)、图6(b)、图7(b)、图8(b)所示。依据这些地区的地势起伏度选择初始尺度参数σ₀分别为：4、4、3.4、4。采用本发明提取的山脉线如图5(c)、图6(c)、图7(c)、图8(c)所示，可见本发明提取出的山脉线清晰、连续，与实际地形较为相符，由此验证了本发明适用于我国主要山脉的山脉线提取。

表1

Claims (7)

1.一种基于高斯尺度空间的山脉线提取方法，其特征在于，包括以下步骤：

步骤一，以高斯核函数对规则格网DEM数据做多尺度分解表达；

步骤二，从分解后的各尺度数据中提取山脊线；

步骤三，融合各尺度山脊线，标记出初选山脉特征点；

步骤四，设定高程阈值对初选山脉特征点作进一步筛选；

步骤五，对筛选后的特征点进行细化、连接，获得山脉线。

2.如权利要求1所述的基于高斯尺度空间的山脉线提取方法，其特征在于，步骤一的具体过程包括：

步骤S10，确定分解层数_n和第一层表达的尺度参数σ₀，并建立各尺度参数间的关系，σ₀的取值范围在之间；

步骤S11，依据各尺度参数σ_i(i＝0,1,…,n-1)对DEM数据进行n层分解。

3.如权利要求2所述的基于高斯尺度空间的山脉线提取方法，其特征在于，步骤S11中，令相邻两层DEM数据间的尺度参数呈k倍关系，即有：

σ_i＝kⁱσ₀

4.如权利要求1所述的基于高斯尺度空间的山脉线提取方法，其特征在于，步骤二的具体过程包括：

步骤S20，将规格网DEM的行列坐标(x,y)视为图像的像素点坐标，该坐标处的高程值f(x,y)视为图像的灰度值，将DEM数据渲染为灰度图像；

步骤S21，利用“提取地形特征线的形态学新方法”，获取规格网DEM中具有正脉冲峰值的带状高亮区域，标记出山脊线，可提取出n各尺度下的山脊线标记图，记为F_i(x,y,σ_i)∈{0,1},i＝0,1,…,n-1。

5.如权利要求1所述的基于高斯尺度空间的山脉线提取方法，其特征在于，步骤三的具体过程包括：

步骤S30，对各尺度山脊线标记图赋予不同的权重w_i，并按照下式进行加权求和，获得融合结果A(x,y)：

$A(x,y)=\underset{i=0}{\overset{n-1}{\mathrm{\Σ}}}{w}_i{F}_i(x,y,{\mathrm{\σ}}_i)$

步骤S31，设定融合阈值TA，将A(x,y)大于该阈值的点标记为初选山脉特征点，其中w_n-1＜TA≤w_n。

6.如权利要求5所述的基于高斯尺度空间的山脉线提取方法，其特征在于，步骤S30中，权重w_i以如下归一化的规范化尺度参数表示：

$w_i=\frac{{\mathrm{\σ}}_i^2}{\underset{j=0}{\overset{n-1}{\mathrm{\Σ}}}{\mathrm{\σ}}_j^2}$

7.如权利要求1所述的基于高斯尺度空间的山脉线提取方法，其特征在于，步骤五中采用“提取地形特征线的形态学新方法”中的后处理法对筛选后的特征点进行处理。