下载:
-
数字高程模型(digital elevation model，DEM)是对地球表面地形的一种离散的数字表达[1]。自20世纪50年代后期被提出以来，DEM受到极大的关注，并在测绘、土木工程、地质、矿山工程、景观建筑、道路设计、防洪、农业、规划、军事工程、飞行器与战场仿真等领域得到了广泛的应用。一般而言，不同数据源需要不同的内插方法来生成DEM。目前，生成DEM的数据主要来源于地形图、遥感数据(既包括航天航空影像数据，又包括合成孔径雷达干涉测量数据和激光雷达数据)、地面测量、既有DEM等[2]；从地形图上获取DEM是目前应用最为广泛的一种方法。我国测绘部门就分别利用1∶1万、1∶5万和1∶25万比例尺的数字线划图生成了多种分辨率的DEM。
通常，由地形图获取DEM时，基于等高线的分布特征，有三种方式生成DEM[1]：等高线离散化、等高线内插和等高线构建Delaunay不规则三角网(triangulated irregular network，TIN)。等高线离散化方法实质是将等高线看作不规则分布的数据，并没有考虑等高线本身的地形特性[1]，这导致生成的DEM可能会出现一些异常；基于等高线数据生成DEM的最陡坡度(流水线)内插算法的内插原理比较简单，但由于数字化的等高线远远没有纸质地形图等高线直观，因此，该方法实现起来还存在许多问题[2]。由于直接由等高线构建的TIN存在“平坦三角形”(即水平三角形)问题[3]，因此，目前工程生产中普遍采用基于等高线和附加的“特征数据”(如地形结构线和特征数据点诸如山顶点、凹陷点、鞍部点等)构建TIN的方法。
近几年提出了很多新的内插方法，胡鹏[4]、胡海[5]等人的研究成果比较具有代表性。“特征数据”本质上是等高线的对偶形式，并不是必须的；而且在工程生产中，很难控制特征数据的密度以平衡DEM的精度和工作量。因此，可利用地图代数直接由等高线内插生成DEM，即MADEM。地图代数是建立在距离变换[6]运算基础上的一种图像操作；它用来内插生成DEM时，不仅不需要额外的辅助特征数据，而且生成的DEM具有较高的精度，满足“高程序同构”[7, 8]的DEM精度评价标准。
但是基于地图代数的内插方法也存在亟待改进之处。由于该方法是通过迭代求取半距等高线(即到两相邻等高线距离相等的线)Cl/2、Cl/4、Cl/8、Cl/16、Cl/32…(Cl为地形图上等高线的基本等高距)来生成DEM的，即迭代地求取两相邻等高线的Voronoi图的边界、并将两等高线的平均值赋予该边界; 至再分已无必要时，以1/2n+1Voronoi图为界(n为最大迭代次数)，分层赋相应高程[9]，本质上这也是一种线性内插方法。但是，该方法需要预先输入最大迭代次数，超过了最大迭代次数时，未插值到的点均赋予最后一个半等高距的增量。可见，最大迭代次数是该方法内插生成DEM中的一个关键的指标。但由于该指标采用了人工干预的方式，没有做到自适应，这在一定程度上限制了该方法的工程性应用。
本文根据数学形态学的测地距离的概念，基于形态学重建[10, 11]和障碍距离变换[6, 12]两种图像运算，提出一种等高线内插方法。它无须任何人工干预，通过与胡海等[5]类似的线性内插方法生成DEM比较，获取的DEM精度等于或高于用地图代数方法获取的DEM精度。
-
测地变换的基本运算包括测地膨胀和测地腐蚀两个算子[13, 14]。测地膨胀涉及标记图像和掩膜图像，两幅图像的大小和定义域相同，但掩膜图像每个像素的值必须大于或者等于标记图像对应的同名像素的值。测地膨胀的实现过程是：利用基本的各向同性结构元素对标记图像作膨胀运算，运算过程中要求获取的结果图像必须保持在掩膜图像之下，即掩膜图像起着限制标记图像膨胀蔓延的作用[13]。 同理，测地腐蚀要求掩膜图像必须小于或者等于标记图像，其实现过程是：利用基本的各向同性结构元素对标记图像作腐蚀运算，运算过程中要求获取的结果图像必须保持在掩膜图像之上，即掩膜图像起着限制标记图像收缩的作用。
有界图像的测地膨胀或测地腐蚀变换，经过一定次数的循环总会收敛，即直至标记图像扩张或收缩完全被掩膜图像阻止；此时再循环一次，标记图像的任一像素的值不再发生改变，从标记图像中对掩膜图像进行形态学重建正是基于这种原理。从标记图像f中对掩膜图像g(f≤g)进行膨胀重建表示为Rδg(f)，其定义为f相对于g的测地膨胀直至稳定[10]：
(1) 式中，i为δg(i)(f)=δg(i+1) (f)时循环的次数。同理，从标记图像f中对掩膜图像g(f≥g)进行腐蚀重建表示为Rεg(f)，其定义为f相对于g的测地腐蚀直至稳定：
(2) 式中，i为εg(i)(f)=εg(i+) (f)时循环的次数。
本文利用形态学重建来生成两幅高原图像，以确定空间任意一点所毗邻的两条等高线(较高的等高线(上等高线)和较低的等高线(下等高线))的高程值。
-
相比欧式距离而言，数学形态学的测地距离[13, 14]是指用于连接图像平面某个子集内的两个像素的最短路径。其中，这个子集区域被称为测地掩膜。假设图像I中有一子集A(A为连通的集合)，则A中两个像素p和q间的测地距离为连接A中p和q间路径P=(p1，p2，…，pn)长度N的最小值：
(3) 式中，集合A为测地掩膜。具有最小长度的路径成为测地路径。
进一步，集合A中的像素p和A中的另一个子集Y间的测地距离dA(p,Y)为p和Y中任意像素q间的最小测地距离：
(4) 式中，集合Y为标记集合。
两像素间、像素与集合间的测地路径如图 1所示。测地掩膜中两个独立像素间的测地距离和测地路径主要取决于掩膜的形状：如果掩膜是凸形的，则测地距离等价于欧氏距离，且测地路径为直线段；否则，测地距离受像素的位置和测地掩膜形状的影响。本文中的测地距离的计算是使用障碍空间的欧氏距离变换[12]来实现的。
图 1 测地掩膜A以及p和Y间的测地路径
Figure 1. Geodesic Paths Between p and q in Geodesic Mask A,and Between p and Y in the Same A,and Y is a Subset of A
-
距离变换是计算并标识空间点集各点到参照体的距离的变换或过程[12, 15]，分为欧氏距离变换和出租车距离变换。前者标识的是欧氏距离，适用于自然形态的图形；后者标识的是曼哈顿距离，适用于规则形态的图形，本文中的距离变换特指欧氏距离变换。
欧氏障碍空间G(Ω)是指欧氏空间中存在障碍物集合(图 2)，即：Ω=w1,w2,…,wk，其中，wk等是障碍物图形，并且互不相交。设P1、P2为平面上的两个点。从其中的一个点出发，至多仅与Ω=w1,w2,…,wk的边界邻接，而到达另一点所得到的最短距离，称为欧氏障碍空间下P1、P2间的穿越障碍物集合Ω=w1,w2,…,wk的距离，简称P1、P2之间的障碍距离，记做Db(P1,P2)。
图 2 通过障碍欧氏距离变换获取测地距离
Figure 2. Geodesic Distance is Calculated by Euclidian Distance Transformation with Obstacles
可见，当将测地掩膜之外的图像空间记为障碍物时，对待处理的某一实体进行障碍欧氏距离变换，就可以获取该实体到掩膜内任意点的测地距离，如图 2所示。图 2中，黑色代表障碍物，白色代表当前的实体，灰色代表测地掩膜，其他色彩代表距离实体的测地距离大小。障碍欧氏距离变换的技术流程参见文献[12]。
本文利用障碍欧氏距离变换来确定空间任意一点到所毗邻的上等高线、下等高线的测地距离。
-
基于上述原理，本文设计了由等高线生成DEM的内插方法。
1) 获取矢量等高线
原始的地形图如果是纸质地形图，可以通过数字化获取矢量形式的数据；如果原始数据本身是以数字线划图的形式提供的，则可以提取其等高线。对矢量等高线有下述要求：在一个图幅范围内，一条等高线不得无故中断，否则通过连接操作获取完整的等高线数据；给每条等高线赋予其相应的高程值。对所有等高线的高程值进行排序，从最小值开始进行奇数或者偶数等高线的类型划分；原则上，两个相邻的高程值中必有一个奇数、一个偶数，奇数和偶数间隔出现。
2) 生成两幅高程图像
首先对矢量等高线进行栅格化操作，在栅格影像上，每个等高线上的像素用该等高线相应的高程值进行填充，对于不在任何等高线上的像素填充为一个无效值，如图 3所示，记为等高线图像CL。图 3中，白色区域间的像素点p被两条等高线所包围，高程值较大的等高线表示为Cu(p)，高程值较小的等高线表示为Cl(p)；同时，p为山顶点或凹点时可以只被Cu(p)或Cl(p)所包围，此时Cu(p)=Cl(p)。
图 3 某像素p的上下等高线
Figure 3. Upper Contour Line and the Lower Contour Line Associated with Pixel p
对于每个像素点来说，如果从等高线图像CL能够确定下面两幅高程图像Pl(下高程图像)和Pu(上高程图像)，则两条等高线的高程是可知的：
(5) 式中，CL[]为对应等高线的高程值。
利用形态学重建的原理来生成这两幅高程图像的方法如下：
(1) 生成一个辅助影像M，且M与栅格化的等高线影像大小和值域相同，生成的规则为：
(6) 式中，hmax为大于该等高线影像最大值的一个数。
(2) 以M为标记图像，以CL为掩膜图像，进行形态学腐蚀重建获取下高程图像Pl；同理，以CL为标记图像，以M为掩膜图像，进行形态学膨胀重建获取上高程图像Pu。使用上述方法获取的某等高线数据的高程图像如图 4(b)、图 4(c)所示。
3) 障碍欧氏距离变换和沿流水线的线性内插
沿着从上等高线Cu(p)到下等高线Cl(p)且过p的测地路径，利用线性插值计算每一像素p的插值。如图 4(d)、图 4(e)所示，图中位于上下等高线间的像素p对应于连接两等高线的测地路径可以由下述方法确定：从上等高线Cu(p)到下等高线Cl(p)且过p的最短路径由从p到Cl(p)以及p到Cu(p)的测地路径组成，该最短路径也即“流水线”、或称为“最速下降线”。地学研究中通常假定上等高线Cu(p)和下等高线Cl(p)间的高度差是沿着测地路径均匀分布的，p的插值H(p)可以由Cu(p)的高程值Pu(p)和Cl(p)的高程值Pl(p)加权得到：
(7) 式中，dM[(p,Cu(p)]为从p到Cu(p)的测地距离；dM[(p,Cl(p)]为从p到Cl(p)的测地距离。此时，连接每一个像素到上、下等高线的测地路径的长度由两个障碍欧氏距离变换得到：即以偶数等高线为障碍物的奇数等高线的距离变换，以奇数等高线为障碍物的偶数等高线的距离变换，如图 4(d)、4(e)所示。
图 4 基于形态学重建和障碍欧氏距离变换集成的等高线内插
Figure 4. Contours’ Interpolation by Integration of Morphological Reconstruction and Euclidean Distance Transformation with Obstacles
需要特别指出，对于等高线数据内插，本文的方法与胡鹏等提出的MADEM方法的原理是完全一致的：沿最速下降线的线性内插，而本文的方法仅仅是从另一途径实现相同的目的，且只有部分点(指“至再分已无必要时，以1/2n+1Voronoi图为界，分层赋相应高程[9]”的那些特殊点)的高程精度可能高于MADEM方法获取的高程的精度。其误差分析完全采用MADEM的误差分析模型[3]。
-
实验数据采用山东省泰安市境内泰山山脉2个局部区域1∶1万的等高线数据，等高距为10 m。为了验证两者的区别，对于其中的一份等高线数据分别使用胡鹏的内插方法和本文的内插方法生成分辨率为1 m的DEM，另外一份数据只采用了本文提出的内插方法。
第一场景等高线的高程分布在630～920 m之间，长和宽分别为1 400 m×1 100 m。使用本文方法获取的DEM如图 4(f)所示，其分辨率为1 m。另外，由此DEM生成了等高距为10m的等高线，见图 4(g)，该等高线能精确地逼近原始等高线，证明了本方法的正确性和生成DEM的高精度。使用胡鹏等的内插方法，分别经过10次和30次迭代，生成了两张DEM，如图 5(a)、5(b)所示。图 5说明，使用胡鹏等的内插方法获取的DEM，迭代的次数越多，越逼近本文的内插方法获取的结果。
图 5 基于胡鹏的内插方法生成的DEM(图例同图 4(f))
Figure 5. DEM Generation by Hu Peng’s Interpolation Method
为了比较两种内插方法获取的DEM的质量，对获取的DEM进行了减操作运算，如图 6所示。两种内插方法获取的DEM在大部分像素上，高程值相等，但在部分像素上存在一定的高程差异，最大的差值约0.3 m，较大的差异点主要分布在远离等高线的位置。可见，胡鹏等的内插方法获取DEM的精度在部分点受迭代次数的限制；而本文提出的内插方法不受此约束条件的影响，而且生成的DEM的质量高于或者等于胡鹏等提出的内插方法获取的DEM。两种方法获取的DEM质量趋向一致，原因在于尽管两种内插方法的实现过程不同，但两种内插方法本质上都是沿流水线对相邻两等高线的线性内插。文献[7]中论述了对等高线进行线性内插具有维护“高程序同构”的优点，既然原理相同，所以本文提出的内插方法同样具有此优势。
图 6 两种内插方法获取的DEM的差值图
Figure 6. Differences of the DEMs Generated by the Two Comparable Interpolation Methods
图 7 基于形态学重建和障碍欧氏距离变换集成的等高线内插
Figure 7. Contours’ Interpolation by Integration of Morphological Reconstruction and Euclidean Distance Transformation with Obstacles
第二场景等高线的高程分布在240～940 m之间，长和宽分别为4 400 m×3 100 m。由此数据生成的DEM见图 7(b)，其分辨率为0.5 m；同时，图 7(a)显示了奇等高线的障碍距离变换的结果。由此DEM生成的高程间隔为10 m的等高线与原等高线套合情况见图 7(c)。通过观察我们发现，两者的吻合度相当好，这与第一个场景的情况是类似的。上述两个场景的实验说明，本文提出的DEM内插方法具有很高的保真度。
-
内插方法是影响数字高程模型精度的一个重要因素[16]。从等高线获取DEM的过程中，内插方法的选择是一个十分重要的技术环节，它对生成DEM的精度有重要的影响。本文提出了一种基于形态学重建和障碍欧氏距离变换集成的DEM内插方法。形态学重建用于获取任一像素对应的上等高线和下等高线的高程值，障碍欧氏距离变换用于获取该点到上下两条等高线的测地距离，使用线性内插获取该点的高程值。与胡鹏等提出的内插方法相比，在只使用等高线数据的情况下，本文方法生成的DEM精度等于或者高于其获取的结果。但本文内插方法要求很高的空间分辨率，且没有考虑三角点和控制点在DEM生成中的作用。而对于后一问题，胡鹏等提出的基于地图代数的MADEM方法已严密而实际地解决了，三角点和控制点的高程是一些特殊的等高线，因而必须使用“加权的距离变换”。因此，后续的研究将放在如何集成三角点和控制点以进一步保证内插的DEM的高精度和高程序同构上。
An Interpolation Method for DEM Generation by Integration of Morphological Reconstruction and Distance Transformation with Obstacles
-
摘要: 等高线是获取数字高程模型(DEM)常用的数据源之一,但内插方法对DEM生成精度有显著的影响。基于形态学重建和测地距离变换运算,提出一种等高线数据生成DEM的内插方法。形态学重建用于获取与空间一点对应的最邻近的上等高线和下等高线的高程值,测地距离变换用于获取该点到上下两条等高线的测地距离;使用沿流水线的线性内插获取该点的高程值。实验表明,在只使用等高线数据生成DEM的情况下,本文提出的内插方法获取的DEM精度更高。Abstract: Contour line is one of the commonly used data sources for DEM generation, and the interpolation method has a significant impact on the accuracy of the contours-derived DEM. However, most of the existing interpolation methods generate low-accuracy DEM due to unreasonable methodologies. In this paper, an interpolation method is presented to build DEM from contour lines by integration of morphological reconstruction and distance transformation with obstacles. Particularly, morphological reconstruction is used to get the elevation values of the higher contour lines and the lower contour lines of any a spatial point between two contour lines, and distance transformation with obstacles is used to get the geodesic distances of the spatial point to the higher contour lines and the lower contour lines respectively. At last, linear interpolation along steepest slope is used to get the elevation values of the pixels to be interpolated. The experiments demonstrate that the accuracy of the DEM generated by our method is higher than the one generated by Hu Peng's MADEM method if only contours are employed to interpolate a DEM.
-
图 1 测地掩膜A以及p和Y间的测地路径
Figure 1. Geodesic Paths Between p and q in Geodesic Mask A,and Between p and Y in the Same A,and Y is a Subset of A
图 2 通过障碍欧氏距离变换获取测地距离
Figure 2. Geodesic Distance is Calculated by Euclidian Distance Transformation with Obstacles
图 3 某像素p的上下等高线
Figure 3. Upper Contour Line and the Lower Contour Line Associated with Pixel p
图 4 基于形态学重建和障碍欧氏距离变换集成的等高线内插
Figure 4. Contours’ Interpolation by Integration of Morphological Reconstruction and Euclidean Distance Transformation with Obstacles
图 5 基于胡鹏的内插方法生成的DEM(图例同图 4(f))
Figure 5. DEM Generation by Hu Peng’s Interpolation Method
图 6 两种内插方法获取的DEM的差值图
Figure 6. Differences of the DEMs Generated by the Two Comparable Interpolation Methods
-
[1] 李志林, 朱庆. 数字高程模型[M]. 武汉: 武汉大学出版社, 2003) Li Zhilin, Zhu Qing. Digital Elevation Models [M]. Wuhan: Wuhan University Press, 2003 [2] 汤国安, 刘学军, 闾国年. 数字高程模型及地学分析的原理与方法[M]. 北京: 科学出版社, 2005) Tang Guoan, Liu Xuejun, Lv Guonian. Principles and Methods of Digital Elevation Models and Geoanalysis [M]. Beijing: Science Press, 2005 [3] 胡鹏, 杨传勇, 吴艳兰, 等. 新数字高程模型理论、方法、标准和应用[M]. 北京: 测绘出版社, 2007) Hu Peng, Yang Chuanyong, Wu Yanlan, et al. New Theories, Methods, Standards and Applications of Digital Elevation Models [M]. Beijing: Surveying and Mapping Press, 2007 [4] 胡鹏, 黄雪莲, 吴艳兰, 等. DEM若干理论问题思考[J]. 哈尔滨工业大学学报, 2006, 38(12):2143-2147) Hu Peng, Huang Xuelian, Wu Yanlan, et al. Some Problems of DEM Error Theory [J]. Journal of Harbin Institute of Technology, 2006, 38(12): 2143-2147 [5] 胡海, 杨传勇, 胡鹏. DEM最优线性生成技术——MADEM[J]. 华中科技大学学报(自然科学版), 2007, 35(6):118-121) Hu Hai, Yang Chuanyong, Hu Peng. Optimal Linear Method for DEM Generation: MADEM [J]. Journal Huazhong University of Science and Technology (Natural Science Edition), 2007, 35(6): 118-121 [6] Zhang J, Lin X, Liu Z, Shen J. Semi-automatic Road Tracking by Template Matching and Distance Transformation in Urban Areas [J]. International Journal of Remote Sensing, 2011, 32(23): 8331-8347 [7] 胡鹏, 白轶多, 胡海. 数字高程模型生成中的高程序同构 [J]. 武汉大学学报·信息科学版, 2009, 34(3):352-357) Hu Peng, Bai Yiduo, Hu Hai. Elevation Order Isomorphism Characteristics of DEM [J]. Geomatics and Information Science of Wuhan University, 2009, 34(3):352-357 [8] 胡海, 吴艳兰, 胡鹏. 数字高程模型精度标准、质量理论和科学观念讨论[J]. 武汉大学学报·信息科学版, 2011, 36 (6):713-716) Hu Hai, Wu Yanlan, Hu Peng. Discussion of DEM Standards, Quality Theory and Conceptions [J]. Geomatics and Information Science of Wuhan University, 2011, 36 (6):713-716 [9] Fan Qingsong. Research on Technologies of DEM Generation by Cartographic Generalization [D]. Wuhan: Wuhan University, 2007: 43-50 [10] Lantuejoul C, Maisonneuve F. Geodesic Methods in Quantitative Image Analysis [J]. Pattern Recognition, 1984, 17(2): 117-187 [11] 沈晶, 刘纪平, 林祥国. 用形态学重建方法进行机载LiDAR数据滤波[J]. 武汉大学学报·信息科学版, 2011, 36(2):167-170) Shen Jing, Liu Jiping, Lin Xiangguo. Airborne LiDAR Data Filtering by Morphological Reconstruction [J]. Geomatics and Information Science of Wuhan University, 2011, 36(2):167-170 [12] 胡鹏, 游涟, 杨传勇, 等. 地图代数[M]. 武汉: 武汉大学出版社, 2007:127-133) Hu Peng, You Lian, Yang Chuanyong. Map Algebra [M]. Wuhan: Wuhan University Press, 2007: 127-133 [13] Lantuejoul C, Beucher S. On the Use of Geodesic Metric in Image Analysis [J]. Journal of Microscopy, 1981, 121: 39-49 [14] Soille P. Morphological Image Analysis: Principles and Applications [M]. Berlin/Heidelberg/New York: Springer-Verlag, 2008 [15] 周培德. 计算几何[M]. 北京: 清华大学出版社, 2002) Zhou Peide. Geometric Algebra [M]. Beijing: Tsinghua University Press, 2002 [16] 张锦明, 游雄, 万刚. 径向基函数算法中插值参数对DEM精度的影响 [J] . 武汉大学学报·信息科学版, 2013, 38 (5): 608-612) Zhang Jinming, You Xiong, Wan Gang. Effects of Interpolation Parameters in Multi-log Radial Basis Function on DEM Accuracy [J] . Geomatics and Information Science of Wuhan University, 2013, 38 (5): 608-612 -
点击查看大图
计量
- 文章访问数: 1076
- HTML全文浏览量: 35
- PDF下载量: 235
- 被引次数: 0
