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

技术领域

本发明属于数字地形分析领域的宏观地形特征提取技术，具体涉 及一种以高斯尺度空间理论对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,{\sigma }_i)=\frac{12}{{2\mathrm{\pi \sigma }}^2}{e}^{-(x^2+y^2)/\left({2\sigma }_i^2\right)}\\ g_i(x,y,{\sigma }_i)=K_i(x,y,{\sigma }_i)\otimes 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{{\sigma }_i^2}{\underset{j=0}{\overset{n-1}{\Sigma }}{\sigma }_j^2}$

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

$A(x,y)=\underset{i=0}{\overset{n-1}{\Sigma }}{w}_i{F}_i(x,y,{\sigma }_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 Topography  Mission)-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,{\sigma }_i)=\frac{1}{{2\mathrm{\pi \sigma }}^2}{e}^{-(x^2+y^2)/\left({2\sigma }_i^2\right)}\\ g_i(x,y,{\sigma }_i)=K_i(x,y,{\sigma }_i)\otimes 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}{\Sigma }}{w}_i{F}_i(x,y,{\sigma }_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