基于格网划分的海量DEM数据生成

第 16 卷 第 1 期 2004 年 1 月 计算机辅助 设计 领导形象设计圆作业设计ao工艺污水处理厂设计附属工程施工组织设计清扫机器人结构设计 与图形学学报 JOURNAL OF COMPU TER2AIDED DESIGN & COMPU TER GRAPHICS Vol116 , No11Jan1 , 2004 　 原稿收到日期 :2002210215 ;修改稿收到日期 :20032082181 本课 题 快递公司问题件快递公司问题件货款处理关于圆的周长面积重点题型关于解方程组的题及答案关于南海问题 得到国家“八六三”高技术研究发展 计划 项目进度计划表范例计划下载计划下载计划下载课程教学计划下载 (2001AA135180) 资助1胡金星 , 男 ,1974 年生 ,博士后研究人员 ,主要研究方向为 GIS、智能交通等1吴焕萍 ,男 ,1977 年生 ,博士研究生 ,主要研究方向为 GIS1 潘　懋 ,男 ,1954 年生 ,教授 ,博士生导师 ,主要研究方向为 GIS、信息地质等1马照亭 ,男 ,1976 年生 ,博士研究生 ,主要研究方向为 GIS1 　 基于格网划分的海量 DEM 数据生成 胡金星1 ,2) 　吴焕萍3) 　潘　懋3) 　马照亭3) 1) (上海交通大学电子信息与电气工程学院　上海　200030) 2) (上海通用卫星导航有限公司　上海　200040) 3) (北京大学地球与空间科学学院　北京　100871) 摘 要 　在自适应格网划分的分割 - 合并 Delaunay 三角剖分算法、格网线性内插方法的基础上 ,提出基于格网划分 的海量 DEM 数据生成算法1 该算法执行效率较高 ,对计算机硬件配置要求较低 ,并适合于并行处理1 关键词 　GIS ;数字高程模型 ;Delaunay 三角剖分 ;规则格网 ;内插 中图法分类号 　P208 ; TP391 Massive DEM Creation with Grid Partitioning Approach Hu Jinxing1 ,2) 　Wu Huanping3) 　Pan Mao3) 　Ma Zhaoting3) 1) ( School of Elect ronics and Elect ric Engineering , S hanghai Jiaotong U niversity , S hanghai 　200030) 2) ( S hanghai General S atellite Navigation Co L t d , S hanghai 　200040) 3) ( School of Earth and S pace Science , Peking U niversity , Beiji ng 　100871) Abstract 　　Applying divide2and2conquer Delaunay triangulation and GRID linear interpolation , a faster algorithm of massive DEM data creation is presented , which requires low computer configuration and fits for parallel processing1 Key words 　GIS ; DEM ; Delaunay triangulation ; grid ; interpolation 1 　前　　言 数字高程模型 (Digital Elevation Model ,DEM) 是构建地形表面空间位置与其高程信息的数字化表 示 ,是对地形表面在地形采样数据基础上的表面重 构1 从 DEM 概念的提出至今 ,DEM 的理论不断得 到发展和完善 ,特别是自 1994 年美国推行的地球空 间数据框架和 1998 年提出的数字地球概念以来 ,以 DEM、数字正射影像和 GIS 的集成组成了数字地 球、数字城市建设的核心 ,而 DEM 成为数字地球海 量数据库的重要组成部分、数字地球的数学基础和 三维虚拟现实的基本框架[122 ]1 在当今的信息社会 中 ,DEM 正在通过数字地球、数字城市的应用深入 到各行各业 ,在城市规划、军事、农业、水利、交通、航 空航天等领域的重要性日益明显1 DEM 基本建模方法主要有三角网建模和规则 格网 ( GRID)建模两种方式1 其中 ,三角网模型具有 精度高、速度快、效率高和顾及地性线 (如断裂线、构 造线)等特点 ,并且 Delaunay 三角网形态良好 ,在表 达地表形态方面表现较为出色 ; GRID 模型用带有 高程值的二维矩阵表示地形 ,其数据排列规则、结构 简单 ,用其计算 分析 定性数据统计分析pdf销售业绩分析模板建筑结构震害分析销售进度分析表京东商城竞争战略分析 每一点的高程值具有良好的实 时性 ,并能比较方便地存储、计算 (坡度、坡向、阴影、 轮廓) 、分析、处理和可视化 ,是目前应用最普遍的 DEM 1 从数字城市到数字中国再到数字地球 ,地形数 据呈几何级数增长 ,科技的高度发达使得海量空间 数据的廉价获取已成为可能 ,针对海量空间数据进 行分析和决策的需求也越来越多1 而现有的理论与 方法只能对整体的地形数据进行处理 ,不但普通计 算机无法满足海量数据的需求 ,即使是硬件配置较 高的计算机 ,内存仍然会受到限制 ,数据量达到一定 级别之后按照常规方法仍然无法处理1 如何在普通 计算机上实现海量地形数据的高效处理 ,是 DEM 的一个关键问题 ,而海量 DEM 数据的高效生成是 基础1 2 　基于格网划分的海量 DEM 生成算 法原理 　　DEM 原始数据可以通过解析航测遥感影像、实 测离散点数据或数字化等高线等方式得到 ,由于这 些数据是以不规则离散点的形式来表示一个规则平 面的空间高程 ,因而需要对它们进行空间内插 ,以获 得 GRID1 生成 GRID (即非格网数据形成格网数 据 ,也称从随机到栅格的内插) 的方法一般有两种 : (1)对原始数据直接进行内插 ; (2) 对原始数据先构 建 Delaunay 三角网 ,然后应用从随机到栅格的转化 方法 ,从 Delaunay 三角网内插成 GRID1 由于第二 种方法在构建 Delaunay 三角网时能充分考虑原始 数据自身的特性、可以灵活地适应任意复杂的地形 情况、能顾及地形特征、在精度和效率方面都较优 , 因此本文采用该方法生成 GRID1 面向海量数据构建 Delaunay 三角网的基本工 作原理是部分数据驻入内存1 首先对原始数据集进 行格网划分 ,对每块格网内的数据构建 Delaunay 子 三角网 ,并将结果存入空间数据库 ,仅将按照 Delau2 nay 法则受邻接格网内的点影响的边界三角形及其 临近三角形常驻内存 ;然后两两 Delaunay 子三角网 进行边界衔接处理 ,将改动和新生成的三角形存入 空间数据库1 由于只需将受影响的边界三角形及其 邻近三角形驻入内存 ,因此按顺序构建 Delaunay 三 角网时 ,驻入内存的数据只是若干格网的边界三角 形及其邻近三角形和当前格网的 Delaunay 三角网 数据 ,对计算机内存的要求较低1 按照格网划分构建的格网块 Delaunay 三角网 包含格网内部与边界三角形 ,对这些三角形按照一 定的内插方法进行内插就可以得到格网块的 GRID 把格网块的 GRID 存入到空间数据库或数据文件 中 ,然后依次对其余格网块进行处理 ,直至所有格网 块都处理完毕1 这样 ,每次处理的只是当前一个格 网块 ,对计算机内存的要求较低1 3 　基于格网划分的海量 DEM 生成 311 　数据预处理 数据预处理是构建 Delaunay 三角网、GRID 生 成之前的准备工作 ,是整个数据处理的一部分 ,一般 包括数据格式的转换、坐标系统的变换、数据的编辑 及数据分块处理等1 由于地形原始数据采样方式的 不同、各种数据的排列顺序或坐标系统各不相同 ,因 此为处理方便需要对数据分块1 目前 ,格网划分的 原则视当前计算机的软硬件水平而定 ,如果格网块 内的点太多 ,则再分割 ,直至所有格网块内的点数都 在给定的范围内 ;但格网内的点也不要太少 ,否则会 由于合并次数太多而影响整体效率1 把经过转换的 数据按照分块的方式进行存储 ,为后续的数据处理 做好准备1 312 　构建 Delaunay 三角网 如何快速、高效地构建 Delaunay 三角网 ,一直 是众多学者研究和关注的焦点1 迄今为止出现了不 少成熟的算法 ,主要有分割- 合并算法、逐点插入法 及三角网生长算法等[324 ]1 其中 ,三角网生长算法目 前较少采用 ;逐点插入法虽然实现较简单 ,占用内存 较小 ,但它的时间复杂度差 ,运行速度慢 ;分割- 合并 算法最为高效 ,但相对复杂1 为充分实现 Delaunay 三 角网的快速、高效的构建 ,本文在分割- 合并算法的 基础上 ,采用自适应块划分方法分割点集 ,然后依序 合并、构建 Delaunay 三角网 ,效率较高1 在对地形采样数据分块存储的基础上 ,按照 预处理分割的逆序逐块从空间数据库下载或从文件 直接读取格网数据 ;然后针对该格网块数据构建 Delaunay 三角网 ,使用基于自适应块划分方法分割 点集到每个子集只包含 2 或 3 个点为止 ,按照分割 顺序的逆顺序进行子块 Delaunay 子三角网的合并 ; 最后形成该格网块的 Delaunay 三角网1 由于相邻格网块内点的影响 ,在格网块相邻边 界周围会出现狭长的非 Delaunay 三角形 ,因此在衔 接格网块时需要对其进行处理1 首先找出两格网块 24 计算机辅助设计与图形学学报 2004 年 Delaunay 子三角网衔接的起始线与终止线 ,然后从 起始线开始到终止线结束 ,伴随着非 Delaunay 三角 形的删除和新 Delaunay 三角形的生成 ,在两个格网 块相邻处形成一个新三角形条带 ,并对它们进行 存储1 上述方法最终构建的是全局 Delaunay 三角网1 可以采用点坐标存储于三角形内的存储方式 ,但这 种方式由于点的重复存储冗余很大 ,三角形拓扑关 系很难维护1 如果按照常规的点坐标单独存放的方 式 ,对三角形只存放点的索引号 ,比较适合于全部数 据驻入内存进行处理 ,对硬件要求较高1 在实际应 用中 ,一般是按照分块的方式进行存储 ,一个格网块 内的点、三角形都单独存放 ,三角形只记录点的索引 号1 所以 ,本文在进行块之间的合并时 ,提供了另外 一种构建类全局 Delaunay 三角网的处理方法 :在边 界处为两格网块都增加点 ,维持全局形状不变 ,这样 使各格网块内的数据相互独立1 图 1 所示为两格网 块 Delaunay 三角网合并处理后形成的两种结果1 图 1 　两 Delaunay 三角网合并处理 313 　GRID 内插 在常规三角网内插方法中 ,一般是给定一点的 平面坐标 ,先检索出该点所在三角形 ,再进行内 插[5 ]1 虽然三角形的检索可以应用拓扑关系按方向 进行快速定位 ,但仍然会占用较多的时间1 本文采 用的方法是从格网块中依次取出一个三角形 ,对该 三角形内的格网点计算其高程值 ,计算完该格网块 中所有三角形内的所有格网点高程值 ,也就生成了 该格网块的 GRID1 GRID 内插方法主要有线性内插、多项式内插、 样条函数内插、傅立叶级数内插、最小二乘趋势内插 及平均移动内插等1 目前大多数 DEM 软件均采用 多项式内插和线性内插 ,本文应用线性内插方法生 成 GRID 1 如图 2 所示 ,假设 △abc 的三个顶点坐标分别 为 a ( x a , ya , z a) , b ( x b , yb , z b) , c ( xc , yc , z c) 三角 形内一点 e ( x , y) 的高程 z 的计算公式为 x y z 1 x a ya z a 1 x b yb z b 1 xc yc z c 1 = 0 或 x - x a y - ya z - z a x b - x a yb - ya z b - z a xc - x a yc - ya z c - z a = 01 令 A = ( yb - ya) ( zc - za) - ( yc - ya) ( zb - za) , B = ( zb - za) ( xc - xa) - ( zc - za) ( xb - xa) , C = ( xb - xa) ( yc - ya) - ( xc - xa) ( yb - ya) , 则 z = z a - A ( x - x a) + B ( y - ya) C (1) 图 2 　三角形内插示意图 　　对于一个三角形 ,参数 A , B , C 是相同的 ,只需 计算一次1 获取在该三角形内的格网点的高程值只 需将 A , B , C 代入式 (1) 即可1 Delaunay 三角网内 插 GRID 的实验统计数据如表 1 所示 ,实验中内插 时间不包含 I/ O 时间 ,计算机硬件配置为 PIII 933 CPU , 256M RAM1 从表 1 中可以看出 ,内插效率 较高1 表 1 　Delaunay 三角网内插 GRID 效率统计表 三角网信息 (个) GRID 信息 离散点 三角形 行 列 内插时间 (s) 75 622 151 213 500 500 01669 1 000 1 000 11025 2 000 2 000 11819 4 000 4 000 31937 152 622 305 201 500 500 11147 1 000 1 000 11164 2 000 2 000 21669 4 000 4 000 81286 基于格网划分的 GRID 生成是在基于格网划分 341 期 胡金星等 :基于格网划分的海量 DEM 数据生成 构建 Delaunay 三角网、并按格网块存储 Delaunay 三 角网的基础上进行的1 生成过程如下 : Step11 从空间数据库或从文件中依次获取一个格网块 的 Delaunay 三角网数据 ; Step21 根据格网划分规则 (格网的 x , y 方向间距) 、该 格网块在所有格网块中的位置 (或格网范围) ,计算该格网的 起始格网点和终止格网点 ; Step31 从该格网块三角网中依次取出一个 △abc ,求出 其三个顶点 a , b , c 坐标中 x 和 y 的最小/ 最大值 xmin = min ( xa , xb , xc) , ymin = min ( ya , yb , yc) , xmax = max ( xa , xb , xc) , ymax = max ( ya , yb , yc) 1 　　Step41 根据该格网块的格网最小/ 最大值和格网的 x , y 方向间距 ,计算该三角形格网的起始格网点和终止格网点 (如图 2 中的 p 和 q) ; Step51 从三角形格网起始点 p 开始 ,按照 y 方向间距 逐行、在行内按照 x 方向间距逐列进行扫描 ,如果当前格网 点在三角形内 ,则应用式 ( 1) 计算其高程 ,直至终止点 q 结 束 ,则该三角形内格网点的高程值计算完毕 ; Step61 转 Step3 重复执行 ,直至该格网块内所有三角形 都处理完毕 ,即该格网块 GRID 内插完毕1 把它存入空间数 据库 (存储在 SQL Server 或 Oracle 的 BLOB 类型字段中) 或 数据文件中 ,释放该格网块的 Delaunay 三角网和 GRID 占用 的内存资源 ; Step71 转 Step1 重复执行 ,直至所有格网块都处理完 毕 ,则 GRID 生成完毕1 这样 , GRID 存储为若干条记录的格网块 GRID ,其记录数与格网划分块数相同1 4 　结　　论 本文算法按照格网划分方式构建 Delaunay 三 角网、GRID 生成 ,该算法只把当前处理格网块的数 据驻入内存 ,对计算机内存要求较低 ,完全可以在普 通计算机器上实现海量 DEM 的数据生成1 另外 ,格 网划分方式可以充分利用并行处理 ,针对各格网块 并行构建 Delaunay 三角网 ,然后再集中进行边界衔 接处理 ;各格网块的 GRID 内插 ,则完全适合于并行 处理1 应用以上算法原理 ,基于 Visual C + + 610 编译 环境 ,本文高效地实现了海量数据的 Delaunay 三角 网构建、GRID 生成1 为算法实现的高效以及跨操作 系统 (如 Linux , Unix 等) ,程序代码用 标准 excel标准偏差excel标准偏差函数exl标准差函数国标检验抽样标准表免费下载红头文件格式标准下载 C + + 编 写1 实验表明 ,本文算法的执行效率较高 ,对计算机 硬件配置的要求较低1 参 　考 　文 　献 [ 1 ] Wang Hongwu , Dong Shihai1 A view2dependent dynamic mul2 tiresolution terrain model [J ]1 Journal of Computer2Aided De2 sign & Computer Graphics , 2000 , 12 (8) : 575～579 (in Chi2 nese) (王宏武 , 董士海1 一个与视点相关的动态多分辨率地形模 型[J ]1 计算机辅助设计与图形学学报 , 2000 , 12 (8) : 575～ 579) [ 2 ] Sun Min , Xue Yong , Ma Ainai1 3D visualization of large DEM data set based on grid division [J ]1 Journal of Computer2Aided Design & Computer Graphics , 2002 , 14 ( 6) : 566～570 (in Chinese) (孙　敏 , 薛 　勇 , 马霭乃1 基于格网划分的大数据集 DEM 三维可视化 [J ]1 计算机辅助设计与图形学学报 , 2002 , 14 (6) : 566～570) [ 3 ] Lee D T , Schachter B J1 Two algorithms for constructing a De2 launay triangulation [J ]1 International Journal of Computer and Information Sciences , 1980 , 9 (3) : 219～242 [ 4 ] Dwyer R A1 A faster divide2and2conquer algorithm for con2 structing Delaunay triangulations [J ]1 Algorithmica , 1987 , 2 (2) : 137～151 [ 5 ] Tan Bing , Jiang Dinghua , Xu Suqin1 The DEM interpolation based on Delaunay triangulations networks [J ]1 Journal of Insti2 tute of Surveying and Mapping , 2001 , 18 (Suppl) : 26～28 (in Chinese) (谭　兵 , 蒋定华 , 许素芹1 一种顾及地形特征的 DEM 内插 方法[J ]1 测绘学院学报 , 2001 , 18 (增刊) : 26～28) 44 计算机辅助设计与图形学学报 2004 年