集成形态学重建和测地距离变换的DEM内插方法

第４１卷第７期２０１６年７月武汉大学学报·信息科学版

Ｇｅｏｍａｔｉｃｓ　ａｎｄ　Ｉｎｆｏｒｍａｔｉｏｎ　Ｓｃｉｅｎｃｅ　ｏｆ　Ｗｕｈａｎ　Ｕｎｉｖｅｒｓｉｔｙ

Ｖｏｌ．４１Ｎｏ．７

Ｊｕｌｙ　

２０１６收稿日期：２０１５－０１－

２６项目资助：国家自然科学基金（４１３７１４０５）；国家测绘地理信息局基础测绘项目（Ａ

１５０６）；中央级公益性科研院所基本科研业务费专项资金（７７７１４１３

）。第一作者：林祥国，副研究员，主要从事遥感数据信息提取的理论与方法研究。ｌｉｎｘｉａｎｇｇ

ｕｏ＠ｇｍａｉｌ．ｃｏｍＤＯＩ：１０．１３２０３／ｊ．ｗｈｕｇｉｓ２０１４００９７文章编号：１６７１－８８６０（２０１６）０７－０８９６－

０７集成形态学重建和测地距离变换的ＤＥＭ内插方法

林祥国１

１　中国测绘科学研究院摄影测量与遥感研究所，北京，１

００８３０摘　要：等高线是获取数字高程模型（ＤＥＭ）常用的数据源之一，但内插方法对ＤＥＭ生成精度有显著的影响。基于形态学重建和测地距离变换运算，提出一种等高线数据生成ＤＥＭ的内插方法。形态学重建用于获取与空间一点对应的最邻近的上等高线和下等高线的高程值，测地距离变换用于获取该点到上下两条等高线的测地距离；使用沿流水线的线性内插获取该点的高程值。实验表明，在只使用等高线数据生成ＤＥＭ的情况下，本文提出的内插方法获取的ＤＥＭ精度更高。关键词：形态学重建；测地距离变换；测地距离；ＤＥＭ；内插中图法分类号：Ｐ２０８；Ｐ２３２ 文献标志码：Ａ

数字高程模型（ｄｉｇ

ｉｔａｌ　ｅｌｅｖａｔｉｏｎ　ｍｏｄｅｌ，ＤＥＭ）

是对地球表面地形的一种离散的数字表达［１］

。自２０世纪５０年代后期被提出以来，Ｄ

ＥＭ受到极大的关注，并在测绘、土木工程、地质、矿山工程、景观建筑、道路设计、防洪、农业、规划、军事工程、飞行器与战场仿真等领域得到了广泛的应用。一般而言，不同数据源需要不同的内插方法来生成ＤＥＭ。目前，生成ＤＥＭ的数据主要来源于地形图、遥感数据（既包括航天航空影像数据，又包括合成孔径雷达干涉测量数据和激光雷达数

据）、地面测量、既有ＤＥＭ等［２］

；

从地形图上获取Ｄ

ＥＭ是目前应用最为广泛的一种方法。我国测绘部门就分别利用１∶１万、１∶５万和１∶２５万比例尺的数字线划图生成了多种分辨率的ＤＥＭ。

通常，由地形图获取ＤＥＭ时，

基于等高线的分布特征，有三种方式生成ＤＥＭ［１］

：

等高线离散化、等高线内插和等高线构建Ｄｅｌａｕｎａｙ不规则三角网（ｔｒｉａｎｇｕｌａｔｅｄ　ｉｒｒｅｇｕｌａｒ　ｎｅｔｗｏｒｋ，ＴＩＮ）。等高线离散化方法实质是将等高线看作不规则分布

的数据，并没有考虑等高线本身的地形特性［

１］

，这导致生成的ＤＥＭ可能会出现一些异常；基于等高线数据生成ＤＥＭ的最陡坡度（流水线）内插算法的内插原理比较简单，但由于数字化的等高线远远没有纸质地形图等高线直观，因此，该方法实

现起来还存在许多问题［

２］

。由于直接由等高线构建的ＴＩＮ存在“

平坦三角形”（即水平三角形）问题［

３］

，因此，目前工程生产中普遍采用基于等高线和附加的“特征数据”（如地形结构线和特征数据点诸如山顶点、凹陷点、鞍部点等）构建ＴＩＮ的方法。

近几年提出了很多新的内插方法，胡鹏［

４］

、胡海［５］

等人的研究成果比较具有代表性。“特征数

据”本质上是等高线的对偶形式，并不是必须的；而且在工程生产中，很难控制特征数据的密度以平衡ＤＥＭ的精度和工作量。因此，可利用地图代数直接由等高线内插生成ＤＥＭ，即ＭＡＤＥＭ。

地图代数是建立在距离变换［

６］

运算基础上的一种图像操作；它用来内插生成ＤＥＭ时，不仅不需要额外的辅助特征数据，而且生成的ＤＥＭ具有较

高的精度，满足“高程序同构”［７，８］

的ＤＥＭ精度评

价标准。

但是基于地图代数的内插方法也存在亟待改进之处。由于该方法是通过迭代求取半距等高线（即到两相邻等高线距离相等的线）Ｃｌ／２、Ｃｌ／４、Ｃｌ／８、Ｃｌ／１６、Ｃｌ／３２…（Ｃｌ为地形图上等高线的基本等高距）来生成ＤＥＭ的，即迭代地求取两相邻等高线的Ｖｏｒｏｎｏｉ图的边界、

并将两等高线的平均值赋予该边界；至再分已无必要时，以１／２

ｎ＋１

Ｖｏｒｏｎｏｉ图为界（

ｎ为最大迭代次数），分层赋相应高程［

９］

，本质上这也是一种线性内插方法。但是，

　第４１卷第７期林祥国：集成形态学重建和测地距离变换的ＤＥＭ内插方法

该方法需要预先输入最大迭代次数，超过了最大迭代次数时，未插值到的点均赋予最后一个半等高距的增量。可见，最大迭代次数是该方法内插生成ＤＥＭ中的一个关键的指标。但由于该指标采用了人工干预的方式，没有做到自适应，这在一定程度上限制了该方法的工程性应用。

本文根据数学形态学的测地距离的概念，基

于形态学重建［

１０，１１］和障碍距离变换［６，１２］

两种图像运算，提出一种等高线内插方法。它无须任何

人工干预，通过与胡海等［

５］

类似的线性内插方法生成ＤＥＭ比较，获取的ＤＥＭ精度等于或高于用地图代数方法获取的ＤＥＭ精度。

１　基本原理

本文提出的内插方法涉及测地变换［１３］

和测地度量［

１４］

两大经典形态学领域。１．１　测地变换

测地变换的基本运算包括测地膨胀和测地腐

蚀两个算子［

１３，１４］

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

［１３］

。同理，测地腐蚀要求掩膜图像

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

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

进行膨胀重建表示为Ｒδ

ｇ（ｆ），其定义为ｆ相对于ｇ的测地膨胀直至稳定［

１

０］：Ｒδ

ｇ（

ｆ）＝δ（ｉ）ｇ（ｆ）（１）式中，ｉ为δ（ｉ）ｇ（ｆ）＝δ（ｉ＋１）

ｇ

（ｆ）时循环的次数。同理，从标记图像ｆ中对掩膜图像ｇ（ｆ≥ｇ）

进行腐蚀重建表示为Ｒε

ｇ（

ｆ），其定义为ｆ相对于ｇ的测地腐蚀直至稳定：

Ｒε

ｇ（

ｆ）＝ε（ｉ）

ｇ（ｆ）（２

）式中，ｉ为ε（ｉ）ｇ（ｆ）＝ε（ｉ＋１）

ｇ

（ｆ）时循环的次数。本文利用形态学重建来生成两幅高原图像，以确定空间任意一点所毗邻的两条等高线（较高的等高线（上等高线）和较低的等高线（下等高线）

）的高程值。１．２　测地距离

相比欧式距离而言，数学形态学的测地距

离［

１３，１４］

是指用于连接图像平面某个子集内的两个像素的最短路径。其中，这个子集区域被称为测地掩膜。假设图像Ｉ中有一子集Ａ（Ａ为连通的集合）

，则Ａ中两个像素ｐ和ｑ间的测地距离为连接Ａ中ｐ和ｑ间路径Ｐ＝（ｐ１，ｐ２，…，ｐｎ）长度Ｎ的最小值：

ｄＡ（

ｐ，ｑ）＝ｍｉｎ｛Ｎ（Ｐ）｜ｐ１＝ｐ，ｐｎ＝ｑ

，且Ｐ Ａ｝（３

）式中，集合Ａ为测地掩膜。具有最小长度的路径

成为测地路径。

进一步，集合Ａ中的像素ｐ和Ａ中的另一个子集Ｙ间的测地距离ｄＡ（ｐ，Ｙ）为ｐ和Ｙ中任意像素ｑ间的最小测地距离：

ｄＡ（

ｐ，Ｙ）＝∧ｑ∈Ｙ

ｄＡ（ｐ，ｑ）（４

）式中，集合Ｙ为标记集合。

两像素间、像素与集合间的测地路径如图１所示。测地掩膜中两个独立像素间的测地距离和测地路径主要取决于掩膜的形状：如果掩膜是凸形的，则测地距离等价于欧氏距离，且测地路径为直线段；否则，测地距离受像素的位置和测地掩膜形状的影响。本文中的测地距离的计算是使用障

碍空间的欧氏距离变换［

１２］

来实现的

。图１　测地掩膜Ａ以及ｐ和Ｙ间的测地路径Ｆｉｇ．１　Ｇｅｏｄｅｓｉｃ　Ｐａｔｈｓ　Ｂｅｔｗｅｅｎ　ｐａｎｄ　ｑｉｎ　ＧｅｏｄｅｓｉｃＭａｓｋ　Ａ，ａｎｄ　Ｂｅｔｗｅｅｎ　ｐａ

ｎｄ　Ｙｉｎ　ｔｈｅ　Ｓａｍｅ　Ａ，ａｎｄ　Ｙｉｓ　ａ　Ｓｕｂｓｅｔ　ｏｆ　

Ａ１．３　障碍欧氏距离变换

距离变换是计算并标识空间点集各点到参照

体的距离的变换或过程［

１２，１５］

，分为欧氏距离变换和出租车距离变换。前者标识的是欧氏距离，适用于自然形态的图形；后者标识的是曼哈顿距离，适用于规则形态的图形，本文中的距离变换特指欧氏距离变换。

７

９８