一种等高线生成DEM的内插方法

摘要

一种等高线生成DEM的内插方法，它包括以下步骤：步骤一，获取并加载满足特定条件的矢量等高线数据，所述特定条件为：每一条等高线在空间上连续分布、不会中断，且奇数等高线和偶数等高线间隔出现；步骤二，进行形态学腐蚀重建和形态学膨胀重建分别生成下高原图像和上高原图像；步骤三，进行障碍欧氏距离变换获取测地路径和测地距离，并进行沿流水线的线性内插生成DEM。本发明生成的DEM具有高精度、高保真的特点，通过线性内插方法生成DEM，并且在只使用等高线数据的情况下获取的DEM精度等于或高于用地图代数方法获取的DEM精度，不但无需人工干预，而且不受到迭代次数的限制，真正实现了自适应算法，各行各业都可以应用。

说明书

技术领域

本发明涉及一种数字高程模型(DEM)的内插方法，具体地说是 一种等高线生成DEM的内插方法，属于地图学与地理信息科学技术 领域。

背景技术

数字高程模型(Digital Elevation Model，DEM)是对地球表 面地形地貌的一种离散的数字表达。自20世纪50年代后期被提出 以来，DEM受到极大的关注，并在测绘、土木工程、地质、矿山工程、 景观建筑、道路设计、防洪、农业、规划、军事工程、飞行器与战 场仿真等领域得到了广泛的应用。与此同时，对DEM的研究始终方 兴未艾，其中生成DEM的内插方法是研究的重点之一。一般而言， 不同数据源需要不同的内插方法来生成DEM。目前，生成DEM的数据 主要来源于地形图、遥感数据(既包括航天航空影像，又包括合成 孔径雷达干涉测量数据和激光雷达数据)、地面测量、既有DEM等； 其中，从地形图上获取DEM是目前应用最为广泛的一种方法。我国 的测绘部门就分别利用1：1万、1：5万和1：25万比例尺的数字线 划图(地形图的一种载体)生成了多种分辨率的DEM。

通常，由地形图获取DEM时，基于等高线的分布特征，有三种 方式生成DEM：等高线离散化、等高线内插和等高线构建狄洛尼不规 则三角网(Triangulated Irregular Network，TIN)。等高线离散 化方法实质是将等高线看作不规则分布的数据，并没有考虑等高线 本身的特性，这导致生成的DEM可能会出现一些异常；Lebel和Olson 提出的基于等高线数据生成DEM的最陡坡度(流水线)内插算法， 内插原理比较简单，但由于数字化的等高线远远没有纸质地形图等 高线直观，因此，该方法实现起来还存在许多问题；由于直接由等 高线构建的TIN存在“平坦三角形”(即，水平三角形)问题，因此， 目前工程生产中普遍采用基于等高线和附加的“特征数据”(如地形 结构线和特征数据点诸如山顶点、凹陷点、鞍部点等)构建TIN的 方法。

近几年提出了很多新的等高线生成DEM的内插方法，其中胡鹏、 胡海等人的研究成果比较具有代表性。他们的立论依据：上述“特 征数据”本质上是等高线的对偶形式，并不是必须的；而且在工程 生产中，很难控制特征数据的密度以平衡DEM的精度和工作量。基 于此，胡鹏、胡海等先后论述了利用地图代数直接由等高线内插生 成DEM的方法，即：MADEM。地图代数，本质上讲，是建立在距离变 换运算基础上的一种图像操作；它用来内插生成DEM时，不仅不需 要额外的辅助特征数据，而且生成的DEM具有较高的精度，满足“高 程序同构”的DEM精度评价标准。目前，基于地图代数开发的DEM 生成软件已经开始部分地应用于我国1：5万DEM的生产过程中，其 效率和精度逐渐获得科研和生产部门的认可。

但基于地图代数的内插方法也存在亟待改进之处。具体而言， 由于该方法是通过迭代求取半距等高线(即到两相邻等高线距离相 等的线)C_l/2,C_l/4,C_l/8,C_l/16,C_l/32…(C_l为地形图上等高线的 基本等高距)来生成DEM的，即迭代地求取两相邻等高线的Voronoi 图的边界、并将两等高线的平均值赋予该边界；至再分已无必要时， 以1/2n+1Voronoi图为界(其中，n为最大迭代次数)，分层赋相应 高程。本质上这也是一种线性内插方法。但是，该方法需要预先输 入最大迭代次数，超过了最大迭代次数时，未插值到的点均赋予最 后一个半等高距的增量。可见，最大迭代次数是该方法内插生成DEM 中的一个关键的指标。但由于该指标采用了人工干预的方式，没有 做到自适应，这在一定程度上限制了该方法的工程性应用。

发明内容

针对现有技术的不足，本发明提供了一种等高线生成DEM的内 插方法，它无须任何人工干预，通过线性内插方法生成DEM，并且 获取的DEM精度等于或高于用地图代数方法获取的DEM精度。

本发明解决其技术问题采取的技术方案是：一种等高线生成DEM 的内插方法，其特征是，包括以下步骤：

步骤一，获取并加载满足特定条件的矢量等高线数据，所述特 定条件为：每一条等高线在空间上连续分布、不会中断，且奇数等 高线和偶数等高线间隔出现；

步骤二，进行形态学腐蚀重建和形态学膨胀重建分别生成下高 原图像和上高原图像；

步骤三，进行障碍欧氏距离变换获取测地路径和测地距离，并 进行沿流水线的线性内插生成DEM。

优选地，在步骤一中，通过数字化纸质地形图中的等高线、或 者从数字线画图中剥离出等高线专题数据以获取矢量形式的等高线 数据，并给每一条等高线赋予相应的高程值，要求每一条等高线在 空间上连续分布且不得无故中断；获取的等高线中任意两个相邻等 高线中一个等高线为奇数等高线、另一个等高线为偶数等高线，且 奇数等高线和偶数等高线间隔出现。

优选地，在步骤二中，对矢量等高线进行栅格化操作生成一个 等高线影像，对等高线影像进行形态学腐蚀重建和形态学膨胀重建 分别生成一个下高原图像和上高原图像。

优选地，在步骤三中，分别由以偶数等高线为障碍物的奇数等 高线的距离变换和以奇数等高线为障碍物的偶数等高线的距离变换 两个障碍欧氏距离变换得到等高线影像上任意一个像素到其上等高 线和下等高线的测地路径的长度，进行沿流水线的线性内插生成 DEM。

进一步地，所述步骤一具体包括以下步骤：(1)获取矢量等高 线数据；(2)处理断裂的等高线，并给等高线赋予其相应的高程值； (3)检查相邻等高线高程值的奇偶特征。

进一步地，所述步骤二具体包括以下步骤：(1)矢量等高线栅 格化；(2)生成辅助影像；(3)进行形态学腐蚀重建和形态学膨胀 重建。

进一步地，所述步骤三具体包括以下步骤：(1)奇、偶等高线 的障碍欧氏距离变换；(2)进行沿流水线的线性内插。

本发明的有益效果如下：

(1)根据数学形态学的测地距离的概念，基于形态学重建和障 碍距离变换两种图像运算，本发明提出了一种基于形态学重建和障 碍欧氏距离变换集成的DEM内插方法。内插方法是影响数字高程模 型精度的一个重要因素，从等高线获取DEM的过程中，内插方法的 选择是一个十分重要的技术环节，它对生成DEM的精度有着重要的 影响。本发明中，形态学重建用于获取任一像素点对应的较高等高 线和较低等高线的高程值，障碍欧氏距离变换用于获取该点到上下 两条等高线的测地距离，使用线性内插获取该点的高程值，整个流 程相当的严密、且具有科学性。

(2)与胡鹏、胡海提出的“MADEM”内插方法相比，在只使用 等高线数据的情况下，本发明生成的DEM精度等于或者高于“MADEM” 获取的结果。但是，胡鹏、胡海提出的“MADEM”内插方法需要人工 干预，受到迭代次数的限制，限制了该方法的工程性应用；而本发 明采用DEM内插方法不但无需人工干预，而且不受到迭代次数的限 制，真正实现了自适应算法，各行各业都可以应用。

附图说明

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

图2(a)至图2(g)为等高线的高程分布在630-920m之间， 长和宽分别为1400m×1100m时基于形态学重建和障碍欧氏距离变换 集成的等高线内插的示意图；

图2(a)为等高线影像的示意图；

图2(b)为下高原图像的示意图；

图2(c)为上高原图像的示意图；

图2(d)为奇等高线到偶等高线的测地距离的示意图；

图2(e)为偶等高线到奇等高线的测地距离的示意图；

图2(f)为内插生成的DEM的示意图；

图2(g)为DEM生成的等高线与原等高线套合图的示意图；

图3为等高线影像上某像素p的上下等高线的示意图；

图4为通过障碍欧氏距离变换获取测地距离的示意图；

图4(a)为未进行障碍欧氏距离变换的原图(白色为实体，黑 色为障碍物)；

图4(b)为进行障碍欧氏距离变换后的结果示意图；

图5(a)和(b)为基于胡鹏的“MADEM”内插方法生成的DEM 的示意图，图5(a)迭代次数为10次时的DEM的示意图，图5(b) 迭代次数为30次时的DEM的示意图；

图6本发明内插方法和胡鹏的“MADEM”内插方法获取的DEM 对比的差值示意图；

图7(a)至图7(c)为等高线的高程分布在240-940m之间， 长和宽分别为4400m×3100m时基于形态学重建和障碍欧氏距离变换 集成的等高线内插的示意图；

图7(a)奇等高线到偶等高线的测地距离(特征值在0至15.11m 之间)的示意图；

图7(b)内插生成的DEM(亮度代表高程的大小)的示意图；

图7(c)内插DEM生成的等高线(黑色)与原等高线套合图及 其局部放大图。

具体实施方式

为能清楚说明本方案的技术特点，下面通过具体实施方式，并 结合其附图，对本发明进行详细阐述。下文的公开提供了许多不同 的实施例或例子用来实现本发明的不同结构。为了简化本发明的公 开，下文中对特定例子的部件和设置进行描述。此外，本发明可以 在不同例子中重复参考数字和/或字母。这种重复是为了简化和清楚 的目的，其本身不指示所讨论各种实施例和/或设置之间的关系。应 当注意，在附图中所图示的部件不一定按比例绘制。本发明省略了 对公知组件和处理技术及工艺的描述以避免不必要地限制本发明。

针对基于等高线内插生成高精度、高保真DEM的需要，本发明 提供了一种基于形态学重建和距离变换的等高线生成DEM的内插方 法，尤其提出了一种基于形态学重建的上高原、下高原图像生成方 法、一种基于障碍欧氏距离变换的测地路径和测地距离生成方法， 以利用等高线内插生成高精度、高保真的DEM。

如图1所示，本发明的一种基于形态学重建和距离变换的等高 线生成DEM的内插方法，其特征是，任一等高线空间上连续分布、 不会中断、且奇数等高线和偶数等高线间隔出现，所述方法包括以 下步骤：

步骤一，获取并加载满足特定条件的矢量等高线数据；

步骤二，进行形态学腐蚀重建和形态学膨胀重建分别生成下高 原图像、上高原图像等两幅高原图像；

步骤三，进行障碍欧氏距离变换获取测地路径和测地距离，并 进行沿流水线的线性内插生成DEM。

进一步地，所述步骤一具体包括以下步骤：

(1)获取矢量等高线数据

原始的地形图如果是纸质地形图，可以通过数字化获取矢量形 式的数据；如果原始数据本身是以数字线划图的形式提供的，则可 以提取其等高线。

(2)处理断裂的等高线，并给等高线赋予其相应的高程值

对矢量等高线需要满足空间连续性的要求。即，在一个图幅范 围内，一条等高线不得无故中断。否则，对于断裂的等高线，需要 通过人工的连接操作连接断裂的等高线，以获取完整的等高线。给 每条等高线赋予其相应的高程值。

(3)检查相邻等高线高程值的奇偶特征

按照高程值的大小，对所有等高线进行排序，并从最小值开始 进行奇数或者偶数等高线的类型划分。原则上，两个相邻的高程值 中必有一个奇数、一个偶数，且奇数等高线和偶数等高线间隔出现。

进一步地，所述步骤二具体包括以下步骤：

(1)矢量等高线栅格化

对矢量等高线进行栅格化操作，生成一个等高线影像，如图2(a)所示。 在等高线影像上，每个等高线上的像素用该等高线相应的高程值进行填充， 对于不在任何等高线上的像素填充为一个无效值(如图3所示)，记为等高线 图像CL。图3中，白色区域间的像素点p被两条等高线所包围，高程值较大 的等高线表示为C_u(p)，高程值较小的等高线表示为C_l(p)；同时，p为山顶点 或凹点时可以只被C_u(p)或C_l(p)所包围，此时C_u(p)＝C_l(p)。

对于每个像素点来说，如果从等高线图像CL能够确定下面两幅高原图像 P_l(记为：下高原图像，如图2(c)所示)和P_u(记为：上高原图像，如图 2(b)所示)，则任一像素对应的两条等高线的高程是可知的，即：

P_l(p)＝CL[C_l(p)]   (1)

P_u(p)＝CL[C_u(p)]

其中CL[]为对应等高线的高程值。

(2)生成辅助影像

以上述等高线影像为参考，生成一个辅助影像M，且M与栅格化的等高线 影像大小和值域相同，生成的规则为：

其中h_max为大于该等高线影像最大值的一个数。

(3)进行形态学腐蚀重建和形态学膨胀重建

以M为标记图像，以CL为掩膜图像，进行形态学腐蚀重建获取下高原图 像P_l；同理，以CL为标记图像，以M为掩膜图像，进行形态学膨胀重建获取 上高原图像P_u。

形态学腐蚀重建和形态学膨胀重建两种操作来自数学形态学的测地变 换。数学形态学中，测地变换的基本运算包括测地膨胀和测地腐蚀两个算子。 测地膨胀涉及两幅图像：标记图像和掩膜图像，且这两幅图像的大小和定义 域相同，但掩膜图像每个像素的值必须大于或者等于标记图像对应的同名像 素的值。测地膨胀的实现过程是：利用基本的各向同性结构元素对标记图像 作膨胀运算，运算过程中要求获取的结果图像必须保持在掩膜图像之下，即 掩膜图像起着限制标记图像膨胀蔓延的作用。

同理，测地腐蚀要求掩膜图像必须小于或者等于标记图像，其实现过程 是：利用基本的各向同性结构元素对标记图像作腐蚀运算，运算过程中要求 获取的结果图像必须保持在掩膜图像之上，即掩膜图像起着限制标记图像收 缩的作用。

有界图像的测地膨胀或测地腐蚀变换，经过一定次数的循环总会收敛， 即直至标记图像扩张或收缩完全被掩膜图像阻止；此时再循环一次，标记图 像的任一像素的值不再发生改变，从标记图像中对掩膜图像进行形态学重建 正是基于这种原理。从标记图像f中对掩膜图像g(f≤g)进行膨胀重建表 示为其定义为f相对于g的测地膨胀直至稳定：

$R_g^{\delta }\left(f\right)={\delta }_g^{\left(i\right)}\left(f\right)---\left(3\right)$

其中i为时循环的次数。同理，从标记图像f中对掩膜图像g (f≥g)进行腐蚀重建表示为其定义为f相对于g的测地腐蚀直至稳 定：

$R_g^{ϵ}\left(f\right)={ϵ}_g^{\left(i\right)}\left(f\right)---\left(4\right)$

其中i为时循环的次数。

使用上述方法获取的某等高线数据的高原图像如图2(b)、(c)所示。 进一步地，所述步骤三具体包括以下步骤：

(1)奇、偶等高线的障碍欧氏距离变换

等高线图像上，连接每一个像素到其上等高线和下等高线的测地路径的 长度由两个障碍欧氏距离变换得到：(1)以偶数等高线为障碍物的奇数等高 线的距离变换，如图2(d)所示；(2)以奇数等高线为障碍物的偶数等高线 的距离变换，如图2(e)所示。

距离变换是计算并标识空间点集中各点到参照体距离的变换或过程，分 为欧氏距离变换和出租车距离变换。前者标识的是欧氏距离，适用于自然形 态的图形；后者标识的是曼哈顿距离，适用于规则形态的图形，本发明中的 距离变换特指欧氏距离变换。

欧氏障碍空间G(Ω)是指欧氏空间中存在障碍物集合(如图4所示)，即： Ω＝{w₁,w₂,…,w_k}，其中，w_k等是障碍物图形，并且互不相交。设P1，P2为 平面上的两个点。我们把从其中的一个点出发，至多仅与Ω＝{w₁,w₂,…,w_k}的 边界邻接，而到达另一点所得到的最短距离，称为欧氏障碍空间下P1，P2 间的穿越障碍物集合Ω＝{w₁,w₂,…,w_k}的距离，简称为P1，P2之间的障碍距离， 记做D_b(P1,P2)。

可见，当将测地掩膜之外的图像空间记为障碍物时，对待处理的某一实 体进行障碍欧氏距离变换，就可以获取该实体到掩膜内任意点的测地距离， 如图4所示。图4(a)中，黑色代表障碍物，白色代表当前的实体，灰色代 表测地掩膜；图4(b)中，其他色彩代表距离实体的测地距离大小。

本发明中，欧氏距离变换的技术流程如下：

①把所有空间点的距离值赋为一足够大的正数M₁，另外把所有实体点的 距离值赋为0。

②顺序访问各个栅格，也即行号从0,1,2,……依次访问，列号也从0,1,2,……依 次访问，并按照下式改写各点的平方值：

de²(i,j)＝min(de²₁(i,j),de²₂(i,j),de²₃(i,j),de²₄(i,j),de²(i,j))   (5)

③逆序访问各个栅格，并改写各栅格点的平方值：

de²(i,j)＝min(de²₅(i,j),de²₆(i,j),de²₇(i,j),de²₈(i,j),de²(i,j))   (6)

④把各点距离平方值转为距离值：

c(i,j)＝[de²(i,j)]^1/2   (7)

本发明中，障碍欧氏距离变换的技术流程可以参考上述的欧氏距离变换 的技术流程，但需要做下述改动：

①在原算法“①把所有空间点的距离值赋为一足够大的正数M₁，另外把 所有实体点的距离值赋为0。”的基础上增加：把所有障碍物点的距离值赋为 另一足够大整数M₂，并且设M₂>>M₁。

②在原有算法②、③、④步骤中对所有点距离平方值de²(i,j)为M₂的点均 不做处理，也即不改写它们的距离平方值。

③原有算法③、④步骤中所有de²(i,j)邻域中距离原点的距离平方值为M₂的点均不参与比较。

(2)进行沿流水线的线性内插

沿着从上等高线C_u(p)到下等高线C_l(p)且过p的测地路径、利用线性插值 计算每一像素p的插值。如图2(d)、(e)所示，图中位于上下等高线间的像 素p对应于连接两等高线的测地路径可以由下述方法确定：从上等高线C_u(p) 到下等高线C_l(p)且过p的最短路径由从p到C_l(p)以及p到C_u(p)的测地路径 组成，该最短路径也即“流水线”、或称为“最速下降线”。

地学研究中通常假定上等高线C_u(p)和下等高线C_l(p)间的高度差是沿着 测地路径均匀分布的，p的插值H(p)可以由C_u(p)的高程值P_u(p)和C_l(p)的高 程值P_l(p)加权得到：

$H\left(p\right)=\frac{P_u\left(p\right)*d_M\left[(p,C_l(p)\right]+P_l\left(p\right)*d_M\left[(p,C_u(p)\right]}{d_M\left[(p,C_l(p)\right]+d_M\left[(p,C_u(p)\right]}---\left(8\right)$

其中，d_M[(p,C_u(p)]为从p到C_u(p)的测地距离，d_M[(p,C_l(p)]为从p到C_l(p)的测 地距离。图2(f)展示了由图2(a)的等高线经过本发明的处理生成的DEM。

需要特别指出，对于等高线数据内插，本发明生成DEM的内插方法与胡 鹏等提出的MADEM方法的原理是完全一致的：沿最速下降线的线性内插，而 本发明生成DEM的内插方法仅仅是从另一途径实现相同的目的，且只有部分 点(指“至再分已无必要时，以1/2ⁿ⁺¹Voronoi图为界，分层赋相应高程”的 那些特殊点)的高程精度可能高于MADEM方法获取的高程的精度。

本发明创造性地提出了一种基于形态学重建和距离变换的等高 线生成DEM的内插方法，尤其提出了一种基于形态学重建的上高原、 下高原图像生成方法、一种基于障碍欧氏距离变换的测地路径和测 地距离生成方法，以利用等高线内插生成高精度、高保真的DEM。

内插方法是影响数字高程模型精度的一个重要因素,从等高线 获取DEM的过程中，内插方法的选择是一个十分重要的技术环节， 它对生成DEM的精度有着重要的影响。本发明中，形态学重建用于 获取任一像素点对应的较高等高线和较低等高线的高程值，障碍欧 氏距离变换用于获取该点到上下两条等高线的测地距离，使用线性 内插获取该点的高程值，整个流程相当的严密、且具有科学性。

为了验证基于形态学重建和距离变换的等高线生成DEM的内插 方法的有效性，使用两景实验数据进行了DEM的内插实验。实验数 据采用了山东省泰安市境内泰山山脉2个局部区域1:1万的等高线 数据，等高距为10米。同时，为了验证两者的区别，对于第一个实 验数据分别使用胡鹏的“MADEM”内插方法和本发明的内插方法生成 分辨率为1m的DEM，第二个实验数据只采用了本发明提出的内插方 法。

第一个实验数据等高线的高程分布在630-920m之间，长和宽分 别为1400m×1100m。使用本发明方法获取的DEM如图2(f)所示， 其分辨率为1m。另外，由此DEM生成了等高距为10m的等高线，见 图2(g)，该等高线能精确地逼近原始的等高线，这证明了本方法 的正确性和生成DEM的高精度。同时，使用了胡鹏的“MADEM”内插 方法，分别经过10次和30次迭代，生成了两张DEM，如图5(a)、 (b)所示。图5说明，使用胡鹏的内插方法获取的DEM，迭代的次 数越多，越逼近本发明生成DEM的内插方法获取的结果。

为了比较两种内插方法获取的DEM的质量，对获取的DEM进行 了减操作运算，如图6所示。通过图6可以看出，两种内插方法获 取的DEM在大部分像素上，高程值相等，但在部分像素上存在一定 的高程差异，最大的差值约0.3米，较大的差异点主要分布在远离 等高线的位置。可见，胡鹏的内插方法获取DEM的精度在部分点受 迭代次数的限制；而本发明生成DEM的内插方法不受此约束条件的 影响，而且生成的DEM的质量高于或者等于胡鹏提出的内插方法获 取的DEM。两种方法获取的DEM质量趋向一致的原因在于，尽管两 种内插方法的实现过程不同，但两种内插方法本质上都是沿流水线 对相邻两等高线的线性内插。

第二个实验数据等高线的高程分布在240-940m之间，长和宽分 别为4400m×3100m。由此数据生成的DEM见图7(b)，其分辨率为 0.5m；同时，图7(a)显示了奇等高线的障碍距离变换的结果。由 此DEM生成的高程间隔为10m的等高线与原等高线套合情况见图7 (c)。通过观察我们发现，两者的吻合度相当好，这与第一个场景 的情况是类似的。上述两个场景的实验说明，本发明生成DEM内插 方法具有很高的保真度。

可见，本发明专利的等高线内插方法取得了高精度、高保真的 DEM，与胡鹏提出的“MADEM”内插方法相比，在只使用等高线数据 的情况下，本发明专利方法生成的DEM精度等于或者高于胡鹏方法 获取的结果。但是，胡鹏、胡海提出的“MADEM”内插方法需要人工 干预，受到迭代次数的限制，限制了该方法的工程性应用；而本发 明采用DEM内插方法不但无需人工干预，而且不受到迭代次数的限 制，真正实现了自适应算法，各行各业都可以应用。

以上所述只是本发明的优选实施方式，对于本技术领域的普通 技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改 进和润饰，这些改进和润饰也被视为本发明的保护范围。