由等高线反演地形曲面的约束优化方法.pdf

1 由等高线反演地形曲面的约束优化方法由等高线反演地形曲面的约束优化方法由等高线反演地形曲面的约束优化方法由等高线反演地形曲面的约束优化方法 宋敦江1，岳天祥2，杜正平2(1 中国科学院科技政策与管理科学究所,北京,100190)(2 中国科学院地理科学与资源研究所,北京,100101)摘要摘要摘要摘要：基于HASM(High Accuracy Surface Modeling)，提出一种由等高线建立DEM的优化控制(Optimization Control)HASMOC方法。在区域离散化格网的基础上，结合HASM方程，通过建立等高线树，计算每个格网点高程值的上下界，建立不等式约束方程，求解该方程即得所求DEM。通过等高线间的约束优化控制，最小化HASM方程的函数值，既能保证地形曲面的整体光滑性，又能保证DEM结果对于原始等高线数据的忠实性。最后通过实际案例，将HASMOC方法与TIN方法的DEM模拟结果及回放的等高线进行比较分析.，结果显示，HASMOC方法能较好地克服TIN的很多缺点。关键词关键词关键词关键词：等高线;HASM 优化控制;约束不规则三角网(TIN);等高线回放 中图分类号中图分类号中图分类号中图分类号:文献标识码文献标识码文献标识码文献标识码:A new Method of DEM Generation from Contour line Using HASM and Optimization Control Song Dunjiang*a,Yue Tianxiangb,DuZhengpingb aInstitute of Policy and Management CAS,Zhongguancun East Road No.55,Beijing China 100190 bInstitute of Geographical Sciences and Natural Resources Research Chinese Academy of Sciences(CAS).11 A,Road DaTun,ChaoYang District,Beijing,China,100101 Abstract:DEM can be constructed from many kinds of data sources among which traditional contour lines from manual interpretation of landforms is still a nontrivial ones.This paper proposed a new method(HASMOC)to construct DEM from contour lines based on the Optimization Control.Presently,HASM(High Accuracy Surface Modeling)which is used to handle discrete points only.Firstly a regular grid points are set up,on which each elevation will be predicted using the new method,each grid point represents the square grid area surrounding itself.Secondly contour lines are used to construct contour trees,which are then used to get the bound constraints for each grid points.Finally,integrating the HASM equations,a numerical optimization problem is formed.By solving the optimization problem DEM can be secured,which is continuous and of high fidelity to original contours.A real contour lines example is given,and results from HASMOC is compared with that from TIN(Triangulated Irregular Network)method,the former is superior to the latter in DEM construction from real contour lines.Keywords:Contour Trees;HASMOC(Optimization Control);TIN(Triangulated Irregular Network);Derived Contours 基金项目基金项目基金项目基金项目:1.引言引言引言引言 尽管新一代地球空间信息技术如 InSAR、LiDAR 和数字摄影测量等技术可快速获取 DEM数据，但是传统的由人工解译地形得到的等高线仍然是一种十分重要的用于建立 DEM 的数据源1111。广泛应用于科研与人们日常生活中的中国 1：25 万、1：5 万及 1：1 万 DEM 数据就是由等高线和离散点通过约束 TIN(Triangulated Irregular Networks)方法构建得到 2 2 2 2-3 3 3 3 。世2 界上某些地区还只有等高线数据，等高线数据从摄影测量立体相对模型中自动生成出来，等高线是人们对地形(Landforms)的一种抽象理解和概括而获得的一种数据，蕴含着大量的地形特征信息 3 3 3 3 。等高线表示地形时有很多类似 TIN 的地方，在地形变化平缓的地方等高线稀疏，在地形变化急促的地方等高线密集。等高线蕴涵着大量地形特征的信息，所以在建立DEM 时应针对等高线选择特定的方法，而不是选择普通的建立 DEM 的算法。由于等高线的复杂性，或者确切地说实际地形的复杂性，使得近 20 多年来，国内外对于等高线构建 DEM方法的研究一直没间断过，由等高线建立DEM的方法主要有TIN方法 4 4 4 4-5 5 5 5 、薄板样条法(Thin Plate Spline,TPS)6 6 6 6-7 7 7 7 、等高线膨胀法(Contour Dilation)8 8 8 8 、最大中间等高线法(Max Intermediate Contours,MIC)9 9 9 9 、最陡坡度法11110 0 0 0 和地形骨架线法(Skeleton)11111 1 1 1 及综合方法(综合中间等高线、最陡方向及地形骨架线等方法的各自优点)11112 2 2 2 。但是以上这些方法各有优缺点，如 TIN 方法在谷地和山峰生成“平三角”，两等高线间利用线性插值，与实际地形的连续过渡的现象有些不符；最陡坡度法只考虑到两条等高线之间的线性关系，而实际上两条等高线之间的最陡方向线常常难以求得，使最陡坡度法成为一个方法，而不是一个具体的算法，限于篇幅，其他各种方法这里就不一一叙述。HASM(High Accuracy Surface Modeling)是我国学者21世纪初提出的基于微分几何曲面论的曲面建模新方法11113 3 3 3-1 1 1 14 4 4 4 。HASM 目前主要用于离散点模拟曲面的研究，由于需要对研究区域内每个待模拟网格点建立曲面方程，HASM 最后需要求解一个巨大的线性方程组，计算量十分巨大，提高运算速度是 HASM 进一步发展需要解决的重要问题。已有数值模拟结果表明，HASM 用于由离散点建立 DEM 得到的结果比地理信息系统中传统方法的精度高许多倍 1 1 1 13 3 3 3 ，但是除离散点外的其他类型数据的 HASM 应用研究却很少。等高线能较好的概括和表现地形地貌特征，等高线既蕴含着定量的信息，如等高线所经过位置点的高程信息；还蕴含着定性信息，如两条等高线间的网格点的高程上下界，等高线还蕴含着山峰、山谷和鞍部以及山脊线、沟谷线等定性信息 1 1 1 1 。等高线树是用于表示等高线空间拓扑关系的一种数据结构，等高线拓扑关系主要包括父子包含关系，以及兄弟并列关系，由闭合等高线构建等高线树的构建方法的研究有很多 11115 5 5 5 ，由未闭合等高线构建等高线树的研究却很少 1 1 1 16 6 6 6 。通过等高线树，可以获得区域内任意一个格网点高程的上下界（为包围该格网点的所有等高线高程的最小值与最大值，开区间）。尽管优化控制方法应用广泛，但是运用优化控制方法进行DEM 构建的研究几乎没有，本文通过等高线树来计算任意点的高程上下界，引入优化控制理论，提出构建 DEM 的 HASMOC 方法。2.等高线反演等高线反演等高线反演等高线反演 DEM 的的的的优化控制优化控制优化控制优化控制方法方法方法方法(HASMOC)2.12.12.12.1 等式约束的等式约束的等式约束的等式约束的 HASMHASMHASMHASM 算法算法算法算法 设(,)|*,*,01,01ijijx yxi h yj hiIjJ=+是计算区域进行均匀正交剖分的网格,其中h为网格分辨率，则 HASM 主方程 Dirichlet 边界问题的有限差分迭代形式可表达为,3()()+=+=+122)(2122)(2,1,1,222,1,1,122211,1,11,1,1,211,1,1,11121,11,1,1njinjinjinjinjinjinjinjinjinjinjinjinjinjinjinjinjinjinjinjinjinjinjinjiGENhffhffhfffGELhffhffhfff (1)njif,表示迭代过程中网格点),(jiyx上第n次迭代的模拟值，nji,(2,122,12,11=ji；)等变量的具体含义，可参见13-14。设 0,),(1,11,12,11,1,111,112,111,11,211,212,211,21,111,112,111,11=+nffffffffffffffffFTnJInJInInInJInJInInInJnJnnnJnJnnn 以上两个方程可以用矩阵形式表达为：=+nnnncBFbAF11 (2)其中，A和nb分别为方程组（1）中第 1 个方程的系数矩阵和右端项向量；B和nc分别为方程组（1）中第 2 个方程的系数矩阵和右端项向量。设：,=nnncbqBAP (3)则 HASM 方法可以表示为如下等式约束最小二乘问题：=+tFStsqPFnnn121.|min (4)其中)(JITRS和1TRt分别为采样矩阵和采样向量，T为采样点个数。建立等式约束最小二乘方法的目的是为了在保证采样点处模拟值接近采样值的条件下，曲面的整体模拟误差最小。由于有限差分方法的计算特点，对模拟区域的边界，迭代过程中其高程不变，均采用初始插值，即0,jinjiff=，(0n，0=i或0=j)。以上是离散点曲面建模的方法，对等高线反演 DEM，为了充分利用等高线信息，保证回放等高线的质量，可采用区域控制，构建不等式约束的方法。2.22.22.22.2 不等式约束的不等式约束的不等式约束的不等式约束的 HASMHASMHASMHASM 算法算法算法算法 等高线(结合图廓边界)，将研究区域分成各子区域。对于每个子区域内的格网点，根据包围它的等高线，可以确定格网点高程值的范围。如图图图图 1 所示，方格代表一个栅格单元，中间的点表示该栅格的中心点，数字 2089 的表示格网点的标志号，等高线的间距为 2.5m，等高线的高程值用红色数字标注，除了等高线数据外，还有 1 个采样点 A（绿色三角点）位于网格点 46 上。网格点 46，被一条高程为 132.5 的山顶等高线包围，同时其附近还存在一个采样点，而网格点 51 被等高值为 117.5 和 120 的等高线所包围。这样，对每个网格点，可以根据其所在的子区域，获得其高程范围。4 网格点 51 处于等高值为 117.5 和 120 这两条等高线之间，所以对网格点 51 的高程，应该满足不等式)51()51()51(ufl，其中5.117)51(=l和120)51(=u。对网格点 46，注意到其所在的子区域为山峰，其中有一个离散采样点 A，高程为 134，对这种情况，一般可以近似认为，这个采样点，是这山峰的高程最大值点。于是，对网格点 46，其高程应该满足不等式)46()46()46(ufl，其中5.132)46(=l，134)46(=u。若网格点 46 所在的这个山峰没有采样点，则其高程所满足的不等式中，上界值)46(u可以取为 135，即为围成这片区域等高线的等高值加上 1 倍的等高距。对于洼地区域，也做类似处理。对于每个网格点，其高程都可以建立不等式，从而得到不等式组，而且，在 HASMOC 迭代过程中，曲面应该一直满足这一不等式组，即：uFln+1 (5)TJIJIIIJIJIIIJJJJlllllllllllllllll),(,1,2,1,11,12,11,1,21,22,21,2,11,12,11,1=TJIJIIIJIJIIIJJJJuuuuuuuuuuuuuuuuu),(,1,2,1,11,12,11,1,21,22,21,2,11,12,11,1=这里jil,、jiu,分别表示网格点)10,10)(,(+JjIiyxji高程的两个控制值。当等高线正好穿过或者非常接近网格点的，可以把这个网格点当成采样点。实际上，可以设定一个参数，对某个网格点),(jiyx，当最近的等高线与它的距离小于时，则将这个网格点定义为采样点，如图中格网点编号为 21、22、27、43、62、56 和 57 等等格网点，高程值为对应的等高线的等高值。图 1 根据等高线确定离散格网点上下界约束示意图 Fig.1 Lower and upper bounds determination of grids from contour lines 这样，针对等高线反演 DEM 的 HASM 方法为如下等式和不等式约束优化：1211min.nnnnPFqSFtstlFu+=(6)5 2.32.32.32.3 任意点的值范围的确定方法任意点的值范围的确定方法任意点的值范围的确定方法任意点的值范围的确定方法 HASMOC 方法建立 DEM 的重要步骤之一是计算所有网格中心点高程值的范围。主要分为两步：第一步是建立等高线树（见图图图图 2），其中图图图图 2a)是等高线，共计有 13 条等高线，形成 3 个山峰，在图图图图 2b)对应的等高线树中则形成 3 个叶子结点；第二步是按某个阈值(一般为模拟栅格单元宽度的 1e-3 倍)对等高线进行离散化处理，如果格网中心点离等高线距离小于这个阈值，则该格网点值为等高线的高程值（该格网点也被称为采样数据）；如果格网中心点离等高线距离大于这个阈值，则通过等高线树，确定该点的高程值的上下界范围：1）两条父子关系的等高线之间的高程值上下界为包围该点的两条等高线的值；2）马鞍形处的点高程值上下界为包围该点的所有等高线高程值的极大值和极小值；3）山峰点的高程值上下界为包围该点的等高线（只有 1 条）加上 1 倍的等高距，类似地还可以获得山谷点的高程值上下界。将格网点的高程值约束代入方程组（6）中的约束条件中，求解方程组（6）即可获得忠实于原始等高线数据的 DEM。为方便求解，本文在利用 HASMOC 方法进行实际模拟DEM 时，将采样数据的等式约束方程转换为在采样值附近(上下界为采样值一个微小量)，这样所有格网点都转换为上下界约束。a)b)图 2 等高线及等高线树图，a)等高线，b)等高线树 Fig.2 Contour lines and the corresponding contour lines trees 3.计算案例计算案例计算案例计算案例 中国科学院千烟洲生态试验站位于江西省中部，隶属于吉安市泰和县，处于 11504 13E，2644 48N，是一个由约 80 个小山丘，9 条沟谷组成的小山村，总面积约为 2.04 km2，海拔高程 60-150m 左右，属典型的红壤丘陵区。地貌类型主要是低丘，丘顶浑圆，海拔多在 100m 以下，最高达 147m，相对高度 20-50m，丘坡坡度以 10-30居多。本实验数据是在江西省吉安市泰和县纸质地形图的基础上，针对千烟洲生态试验站附近的地形图扫描然后人工矢量化得到的数据，1954 年北京坐标系，1956 年黄海高程系，等高距 2.5 米，199?年（原图看不清），该区域内多植被，且大部分为水稻田。为验证 HASMOC 方法在稀疏等高线情况下的有效性，从 2.5 米等高距的等高线中抽取了 20 米等高距(共计 76 条)(见图 3)进行试验。6 图 3 千烟洲扫描矢量化原始等高线数据(20 米等高距)Fig.3 Origin Contour Data of QianYanZhou from Cartographic Scanning(20 m interval)a)b)图 4 DEM 结果比较图，a)TIN 方法，b)HASMOC 方法 Fig.4 Interpolated DEM of,a)TIN,b)HASMOC 利用图图图图 3 中的等高线数据，分别使用 TIN 和 HASMOC 方法建立 5 米分辨率 DEM(图图图图 4)，TIN 方法模拟得到的结果为等高线的凸包，HASMOC 模拟方法返回的是外接矩形区域(521 行546 列)，为了统一起见，只显示了与 TIN 相同范围的 DEM。从图 4 中可以看出，TIN 在地形中的山顶和谷地出现了很多“磨平”现象，遗弃了很多地形细节，而 HASMOC 方法则基本保留了地形细节，地形曲面连续光滑、自然。在多处“弯月”型山谷中，TIN 方法构建的DEM 严重失真，未保证建立的 DEM 忠实于原始等高线；在多个马鞍处，TIN 方法也出现磨平现象。通过回放等高线(图 5),可以再次看到，TIN 方法回放的等高线(35 条)与原始等高线的相差非常大，而 HASMOC 方法回放的等高线与原始等高线基本完全一致(都是 76 条)。类似地，我们还进行了从 2.5 米等高距的等高线中抽取了 10 米等高距实验，得到类似的结果，TIN 方法回放的等高线(110 条)与原始等高线(152 条)的相差非常大，而 HASMOC 方法则仍与原始等高线基本完全一致。7 a)b)图 5 回放等高线(蓝色+横条)与原始等高线(黑色)叠加分析图，a)TIN 方法，b)HASMOC 方法 Fig.5 Overlay of the derived contour lines with original contour lines from,a)TIN,b)HASMOC 708090100110120130140012345678x 104高高高高高高高高频频频频频频频频a)70809010011012013014000.511.522.53x 104高高高高高高高高频频频频频频频频 b)图 6 高程频率统计柱状图，a)TIN 方法，b)HASMOC 方法 Fig.6 Histogram of DEM from,a)TIN,b)HASMOC 根据17的研究，由等高线建立的 DEM 容易受到等高线数据本身的影响，如 DEM 中等高线高程值附近的值相对较多，这是很多 DEM 构建方法无法避免的，这与实际地形严重不符，于是 DEM 的高程频率分布柱状图也是建议构建 DEM 方法优劣的一个标准。从图 6 中可以看出，TIN 方法构建的 DEM 受等高线数据的影响很大，分别在等高线值为 80、100 和 120处的高程点分布较多，HASMOC 也似乎受到了一定的影响，但是影响很小。4.结论与展望结论与展望结论与展望结论与展望 HASM 是从数学领域发展起来的一种曲面建模新方法，目前主要用于离散点的曲面模拟。本文引入优化控制理论，通过等高线树来计算任意点的高程上下界，最小化 HASM 方程的函数值，提出一种由等高线构建 DEM 的 HASMOC 方法。HASMOC 方法模拟得到了较 TIN更加符合实际的 DEM，HASMOC 方法的回放等高线与原始等高线完全一致，TIN 的回放等高线与原始等高线差别很大；另外从模拟的 DEM 来看，TIN 产生大量的磨平现象，丢失了很多地形细节，而 HASMOC 方法却保留了大量的地形细节。8 理论上，通过对等高线间各子区域的约束优化控制，HASMOC 方法能保证回放等高线的最大偏离距离不超过一个栅格单元，这点在实际应用的得到了证明。由于由等高线提取地形特征线过程比较复杂繁琐，故本文只研究了等高线间的约束优化这种最简单的情形。下一步的研究方向包括从等高线中提取地形特征点与特征线，并将地形特征点与特征线作为HASMOC 方法新的(等式或不等式)约束条件，构建更高精度 DEM。由于需要对计算研究区域内每个待模拟网格点的曲面参数变量，HASMOC 方法计算量大，另外大规模约束优化求解也是计算数学领域一个难点问题，提高运算速度是 HASMOC进一步发展需要解决的重要问题。致谢致谢致谢致谢：感谢中国科学院千烟洲试验站的张文江老师提供实验数据。参考文献参考文献参考文献参考文献 1.C.Gallant,Michael F.Hutchinson John,Digital Elevation Models and Representation of Terrain Shape.Terrain Analysis:Principles and Applications M.J.C.G.John P.Wilson,John Wiley and Sons,2000 2.Zhilin Li,Qing Zhu,Christopher Gold.Digital terrain modeling:principles and methodology M.Boca Raton,Florida New York,N.Y.,Taylor&Francis.2005 3.El-Sheimy,Naser Valeo,Caterina.Digital Terrain Modelling:Acquisition,Manipulation and ApplicationsM,Artech House,Inc.:Norwood,c2005 4.Clarke,A L Gruen,A Loon,J C.Application of contour data for generating high fidelity grid digital elevation models.5th Int.Symp.Comp.Assisted Cartography&Int.Soc.Photogrammetry and Remote Sensing Commission IV,Crystal City VA(USA),1982 5.GARCIA,ANTONIO BELLO,et al.A contour line based triangulating algorithm.In Proceedings of the 5th International Symposium on Spatial Data Handling(Charleston SC USA,1992),P.Bresnahan,E.Corwin,and D.Cowen,Eds.,vol.2,IGU Commission of GIS,1992,411423.6.Hancock,P.A.and M.F.Hutchinson(2006).Spatial interpolation of large climate data sets using bivariate thin plate smoothing splinesJ.Environmental Modelling&Software,21(12):1684-1694.7.Hutchinson,Michael.Calculation of hydrologically sound digital elevation models.In Proceedings of the Third International Symposium on Spatial Data Handling,International Geographical Union,Columbus,Ohio,1988 8.Taud,Hind,Jean-Francois Parrot and Roman Alvarez.DEM generation by contour line dilation J.Computers&Geosciences,1999,25(7):775-783.9.Gousie,M.B.and Franklin,W.R.Constructing a DEM from Grid-based Data by Computing Intermediate Contours.In GIS 2003:Proceedings of the Eleventh ACM International Symposium on Advances in Geographic Information Systems,E.Hoel and P.Rigaux,Eds.,2003,71-77.10.Oky Dicky Ardiansyah,Prima and Ryuzo Yokoyama.DEM generation method from contour lines based on the steepest slope segment chain and a monotone interpolation function.J ISPRS Journal of Photogrammetry and Remote Sensing,2002,57(1-2):86-101.11.GOLD,DAVID THIBAULT AND CHRISTOPHER M.,Terrain Reconstruction from Contours by Skeleton ConstructionJ.GeoInformatica,2000,349-373.12.Hu Pen,Yang Chuanyong,Wu Yanlan,Hu Hai,New DEM theories,methods,standards and applicationM,Surveying and Mapping Press,2007(胡鹏等,新数字高程模型M,测绘出版社,2007)13.Yue,T.X.,Du,Z.P.,Lu,J.,High accuracy surface modeling and and error analysisJ,Progress in Natural Science,2004,14(2):83-89.(岳天祥,杜正平,刘纪远,高精度曲面建模与误差分析J,自然科学进展,2004,14（2）：83-89.)14.YUE T.X.,Du Z.P.,SONG D.J.,GONG J.,A new method of surface modeling and its application to 9 DEM construction,Geomorphology,2007,91(12),161-172 15.Jun Chen,Chaofei Qiao,Renliang Zhao,A Voronoi interior adjacency-based approach for generating a contour treeJ,Computers&Geosciences,30(2004),355367 16.WU Hehai,Construction of Contour Trees and its ApplicationJ，Science of Surveying and Mapping,1996(1),2-7,(毋河海,等高线树的自动建立及其应用J,测绘科技动态,1996(1),2-7)17.R.Carla,A.Carrara,and G.Bitelli.Comparison of techniques for generating digital terrain models from contour linesJ,International Journal of Geographic Information Science,1997,11(5):451-473 _ 收稿日期收稿日期收稿日期收稿日期：2010-5-4；修回日期：2010-#-#