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

第４１卷第７期２０１６年７月武汉大学学报·信息科学版Ｇｅｏｍａｔｉｃｓ　ａｎｄ　Ｉｎｆｏｒｍａｔｉｏｎ　Ｓｃｉｅｎｃｅ　ｏｆ　Ｗｕｈａｎ　ＵｎｉｖｅｒｓｉｔｙＶｏｌ．４１Ｎｏ．７Ｊｕｌｙ　２０１６收稿日期：２０１５－０１－２６项目资助：国家自然科学基金（４１３７１４０５）；国家测绘地理信息局基础测绘项目（Ａ１５０６）；中央级公益性科研院所基本科研业务费专项资金（７７７１４１３）。

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

ｌｉｎｘｉａｎｇｇｕｏ＠ｇｍａｉｌ．ｃｏｍＤＯＩ：１０．１３２０３／ｊ．ｗｈｕｇｉｓ２０１４００９７文章编号：１６７１－８８６０（２０１６）０７－０８９６－０７集成形态学重建和测地距离变换的ＤＥＭ内插方法林祥国１１　中国测绘科学研究院摄影测量与遥感研究所，北京，１００８３０摘　要：等高线是获取数字高程模型（ＤＥＭ）常用的数据源之一，但内插方法对ＤＥＭ生成精度有显著的影响。

基于形态学重建和测地距离变换运算，提出一种等高线数据生成ＤＥＭ的内插方法。

形态学重建用于获取与空间一点对应的最邻近的上等高线和下等高线的高程值，测地距离变换用于获取该点到上下两条等高线的测地距离；使用沿流水线的线性内插获取该点的高程值。

实验表明，在只使用等高线数据生成ＤＥＭ的情况下，本文提出的内插方法获取的ＤＥＭ精度更高。

关键词：形态学重建；测地距离变换；测地距离；ＤＥＭ；内插中图法分类号：Ｐ２０８；Ｐ２３２ 文献标志码：Ａ 数字高程模型（ｄｉｇｉｔａｌ　ｅｌｅｖａｔｉｏｎ　ｍｏｄｅｌ，ＤＥＭ）是对地球表面地形的一种离散的数字表达［１］。

自２０世纪５０年代后期被提出以来，ＤＥＭ受到极大的关注，并在测绘、土木工程、地质、矿山工程、景观建筑、道路设计、防洪、农业、规划、军事工程、飞行器与战场仿真等领域得到了广泛的应用。

一般而言，不同数据源需要不同的内插方法来生成ＤＥＭ。

目前，生成ＤＥＭ的数据主要来源于地形图、遥感数据（既包括航天航空影像数据，又包括合成孔径雷达干涉测量数据和激光雷达数据）、地面测量、既有ＤＥＭ等［２］；从地形图上获取ＤＥＭ是目前应用最为广泛的一种方法。

我国测绘部门就分别利用１∶１万、１∶５万和１∶２５万比例尺的数字线划图生成了多种分辨率的ＤＥＭ。

通常，由地形图获取ＤＥＭ时，基于等高线的分布特征，有三种方式生成ＤＥＭ［１］：等高线离散化、等高线内插和等高线构建Ｄｅｌａｕｎａｙ不规则三角网（ｔｒｉａｎｇｕｌａｔｅｄ　ｉｒｒｅｇｕｌａｒ　ｎｅｔｗｏｒｋ，ＴＩＮ）。

等高线离散化方法实质是将等高线看作不规则分布的数据，并没有考虑等高线本身的地形特性［１］，这导致生成的ＤＥＭ可能会出现一些异常；基于等高线数据生成ＤＥＭ的最陡坡度（流水线）内插算法的内插原理比较简单，但由于数字化的等高线远远没有纸质地形图等高线直观，因此，该方法实现起来还存在许多问题［２］。

由于直接由等高线构建的ＴＩＮ存在“平坦三角形”（即水平三角形）问题［３］，因此，目前工程生产中普遍采用基于等高线和附加的“特征数据”（如地形结构线和特征数据点诸如山顶点、凹陷点、鞍部点等）构建ＴＩＮ的方法。

近几年提出了很多新的内插方法，胡鹏［４］、胡海［５］等人的研究成果比较具有代表性。

“特征数据”本质上是等高线的对偶形式，并不是必须的；而且在工程生产中，很难控制特征数据的密度以平衡ＤＥＭ的精度和工作量。

因此，可利用地图代数直接由等高线内插生成ＤＥＭ，即ＭＡＤＥＭ。

地图代数是建立在距离变换［６］运算基础上的一种图像操作；它用来内插生成ＤＥＭ时，不仅不需要额外的辅助特征数据，而且生成的ＤＥＭ具有较高的精度，满足“高程序同构”［７，８］的ＤＥＭ精度评价标准。

但是基于地图代数的内插方法也存在亟待改进之处。

由于该方法是通过迭代求取半距等高线（即到两相邻等高线距离相等的线）Ｃｌ／２、Ｃｌ／４、Ｃｌ／８、Ｃｌ／１６、Ｃｌ／３２…（Ｃｌ为地形图上等高线的基本等高距）来生成ＤＥＭ的，即迭代地求取两相邻等高线的Ｖｏｒｏｎｏｉ图的边界、并将两等高线的平均值赋予该边界；至再分已无必要时，以１／２ｎ＋１Ｖｏｒｏｎｏｉ图为界（ｎ为最大迭代次数），分层赋相应高程［９］，本质上这也是一种线性内插方法。

但是，　第４１卷第７期林祥国：集成形态学重建和测地距离变换的ＤＥＭ内插方法该方法需要预先输入最大迭代次数，超过了最大迭代次数时，未插值到的点均赋予最后一个半等高距的增量。

可见，最大迭代次数是该方法内插生成ＤＥＭ中的一个关键的指标。

但由于该指标采用了人工干预的方式，没有做到自适应，这在一定程度上限制了该方法的工程性应用。

本文根据数学形态学的测地距离的概念，基于形态学重建［１０，１１］和障碍距离变换［６，１２］两种图像运算，提出一种等高线内插方法。

它无须任何人工干预，通过与胡海等［５］类似的线性内插方法生成ＤＥＭ比较，获取的ＤＥＭ精度等于或高于用地图代数方法获取的ＤＥＭ精度。

１　基本原理本文提出的内插方法涉及测地变换［１３］和测地度量［１４］两大经典形态学领域。

１．１　测地变换测地变换的基本运算包括测地膨胀和测地腐蚀两个算子［１３，１４］。

测地膨胀涉及标记图像和掩膜图像，两幅图像的大小和定义域相同，但掩膜图像每个像素的值必须大于或者等于标记图像对应的同名像素的值。

测地膨胀的实现过程是：利用基本的各向同性结构元素对标记图像作膨胀运算，运算过程中要求获取的结果图像必须保持在掩膜图像之下，即掩膜图像起着限制标记图像膨胀蔓延的作用［１３］。

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

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

从标记图像ｆ中对掩膜图像ｇ（ｆ≤ｇ）进行膨胀重建表示为Ｒδｇ（ｆ），其定义为ｆ相对于ｇ的测地膨胀直至稳定［１０］：Ｒδｇ（ｆ）＝δ（ｉ）ｇ（ｆ）（１）式中，ｉ为δ（ｉ）ｇ（ｆ）＝δ（ｉ＋１）ｇ（ｆ）时循环的次数。

同理，从标记图像ｆ中对掩膜图像ｇ（ｆ≥ｇ）进行腐蚀重建表示为Ｒεｇ（ｆ），其定义为ｆ相对于ｇ的测地腐蚀直至稳定：Ｒεｇ（ｆ）＝ε（ｉ）ｇ（ｆ）（２）式中，ｉ为ε（ｉ）ｇ（ｆ）＝ε（ｉ＋１）ｇ（ｆ）时循环的次数。

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

１．２　测地距离相比欧式距离而言，数学形态学的测地距离［１３，１４］是指用于连接图像平面某个子集内的两个像素的最短路径。

其中，这个子集区域被称为测地掩膜。

假设图像Ｉ中有一子集Ａ（Ａ为连通的集合），则Ａ中两个像素ｐ和ｑ间的测地距离为连接Ａ中ｐ和ｑ间路径Ｐ＝（ｐ１，ｐ２，…，ｐｎ）长度Ｎ的最小值：ｄＡ（ｐ，ｑ）＝ｍｉｎ｛Ｎ（Ｐ）｜ｐ１＝ｐ，ｐｎ＝ｑ，且Ｐ Ａ｝（３）式中，集合Ａ为测地掩膜。

具有最小长度的路径成为测地路径。

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

两像素间、像素与集合间的测地路径如图１所示。

测地掩膜中两个独立像素间的测地距离和测地路径主要取决于掩膜的形状：如果掩膜是凸形的，则测地距离等价于欧氏距离，且测地路径为直线段；否则，测地距离受像素的位置和测地掩膜形状的影响。

本文中的测地距离的计算是使用障碍空间的欧氏距离变换［１２］来实现的。

图１　测地掩膜Ａ以及ｐ和Ｙ间的测地路径Ｆｉｇ．１　Ｇｅｏｄｅｓｉｃ　Ｐａｔｈｓ　Ｂｅｔｗｅｅｎ　ｐａｎｄ　ｑｉｎ　ＧｅｏｄｅｓｉｃＭａｓｋ　Ａ，ａｎｄ　Ｂｅｔｗｅｅｎ　ｐａｎｄ　Ｙｉｎ　ｔｈｅ　Ｓａｍｅ　Ａ，ａｎｄ　Ｙｉｓ　ａ　Ｓｕｂｓｅｔ　ｏｆ　Ａ１．３　障碍欧氏距离变换距离变换是计算并标识空间点集各点到参照体的距离的变换或过程［１２，１５］，分为欧氏距离变换和出租车距离变换。

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

７９８武汉大学学报·信息科学版２０１６年７月欧氏障碍空间Ｇ（Ω）是指欧氏空间中存在障碍物集合（图２），即：Ω＝ｗ１，ｗ２，…，ｗ｛｝ｋ，其中，ｗｋ等是障碍物图形，并且互不相交。

设Ｐ１、Ｐ２为平面上的两个点。

从其中的一个点出发，至多仅与Ω＝ｗ１，ｗ２，…，ｗ｛｝ｋ的边界邻接，而到达另一点所得到的最短距离，称为欧氏障碍空间下Ｐ１、Ｐ２间的穿越障碍物集合Ω＝ｗ１，ｗ２，…，ｗ｛｝ｋ的距离，简称Ｐ１、Ｐ２之间的障碍距离，记做Ｄｂ（Ｐ１，Ｐ２）。

图２　通过障碍欧氏距离变换获取测地距离Ｆｉｇ．２　Ｇｅｏｄｅｓｉｃ　Ｄｉｓｔａｎｃｅ　ｉｓ　Ｃａｌｃｕｌａｔｅｄ　ｂｙ　ＥｕｃｌｉｄｉａｎＤｉｓｔａｎｃｅ　Ｔｒａｎｓｆｏｒｍａｔｉｏｎ　ｗｉｔｈ　Ｏｂｓｔａｃｌｅｓ可见，当将测地掩膜之外的图像空间记为障碍物时，对待处理的某一实体进行障碍欧氏距离变换，就可以获取该实体到掩膜内任意点的测地距离，如图２所示。

图２中，黑色代表障碍物，白色代表当前的实体，灰色代表测地掩膜，其他色彩代表距离实体的测地距离大小。

障碍欧氏距离变换的技术流程参见文献［１２］。

本文利用障碍欧氏距离变换来确定空间任意一点到所毗邻的上等高线、下等高线的测地距离。

２　内插方法技术流程基于上述原理，本文设计了由等高线生成ＤＥＭ的内插方法。

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

对矢量等高线有下述要求：在一个图幅范围内，一条等高线不得无故中断，否则通过连接操作获取完整的等高线数据；给每条等高线赋予其相应的高程值。

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

２）生成两幅高程图像首先对矢量等高线进行栅格化操作，在栅格影像上，每个等高线上的像素用该等高线相应的高程值进行填充，对于不在任何等高线上的像素填充为一个无效值，如图３所示，记为等高线图像ＣＬ。

图３中，白色区域间的像素点ｐ被两条等高线所包围，高程值较大的等高线表示为Ｃｕ（ｐ），高程值较小的等高线表示为Ｃｌ（ｐ）；同时，ｐ为山顶点或凹点时可以只被Ｃｕ（ｐ）或Ｃｌ（ｐ）所包围，此时Ｃｕ（ｐ）＝Ｃｌ（ｐ）。