国家基本比例尺地形图图幅号查询系统的设计

国家基本比例尺地形图图幅号查询系统的 设计 领导形象设计圆作业设计ao工艺污水处理厂设计附属工程施工组织设计清扫机器人结构设计 摘 要 图幅号的快速查询对于企业和用户有着广泛的用途。本文介绍了国家基本比例尺地形图图幅号编号旧方法和新方法，通过比较新旧图幅号编号方法的差异，表明了两者在不同比例尺下的转换关系，并且介绍了高斯坐标与大地坐标之间的换算关系，进而可从目标点的高斯坐标推算出该地的图幅号。本文介绍的图幅号查询系统是用Visual C++编程语言设计实现的，重要功能模块包括高斯坐标与大地坐标的转换、经纬度查询新图幅号、大地坐标查询新图幅号、新旧图幅号转换、查询图幅信息等，系统设计界面友好，查询结果准确，可以使得用户方便快捷地查询关于图幅号的信息。 关键词：图幅号，查询系统，比例尺 Abstract Quick Query map sheet number for businesses and consumers have a wide range of uses. This article describes the basic national scale topographic map sheet number of old and new methods, by comparison between the old and new numbering map sheet number, indicating the relationship between the conversion at different scales, and introduced the Gaussian coordinates of the earth's conversion relationship between the land and thus can calculate the map sheet number from the Gaussian coordinates of the target point. Map sheet number inquiry system described in this article is to use Visual C + + Programming Language Design and Implementation of an important function modules include converting the Gaussian coordinates of the earth, the latitude and longitude query new map sheet number, geodetic coordinates query new map sheet number, old map sheet number conversion query map sheet information, system design and friendly interface, query results are accurate, so that the user can quickly and easily query information about the map sheet number. Keywords: map sheetnumber，query system，scale 目 录 摘 要 I Abstract I 第一章 绪论 1 1.1 地形图图幅号概述 1 1.1.1 地形图分幅与编号 1 1.1.2 图幅号的作用 1 1.2 图幅号查询系统设计意义 2 第二章 我国地形图图幅号编号方法 3 2.1 国家基本比例尺旧图幅号编号 3 2.2 国家基本比例尺新图幅号编号 3 2.2.1 新图幅编号规则 3 2.2.2 经纬度与图幅号转换关系 5 2.3 新旧图幅号转换关系 6 2.3.1 新图幅到旧图幅的转换 6 2.3.2 旧图幅到新图幅的转换 8 2.3.3 图幅理论面积的计算 11 第三章 系统设计与实现 12 3.1 系统设计原则 12 3.2 系统设计路线 12 3.3 系统界面 12 3.4 系统功能模块 15 3.4.1 系统功能模块图 15 3.4.2 功能模块代码 15 第四章 使用说明 23 4.1 编写目的 23 4.2 系统功能 23 4.3 运行说明 23 4.3.1 启动方法 23 4.3.2 查询 24 4.3.3 清除数据 26 4.3.4 保存 26 4.3.5 退出 27 4.4 异常说明 27 参考文献 28 致 谢 29 外文原文 30 中文翻译 41 1 第一章 绪论 1.1 地形图图幅号概述 1.1.1 地形图分幅与编号 为了保管和使用方便，我国对每一种基本比例尺地形图的图廓大小都做了规定，每一幅地形图给出了相应的号码标志，这就是地形图的分幅与编号。编号将划分的图幅，按比例尺大小和所在的位置，用文字符号和数字符号进行编号。编号是每个图幅的数码标记，它们应具备系统性、逻辑性和唯一性。常用的地图编号有行列式、自然序数式、行列－自然序数式和西南角图廓点坐标公里数编号等。地形图分幅有两种方法：一种是矩形分幅，另一种是经纬线分幅。 （1）矩形分幅 每幅地图的图廓都是一个矩形。适用于小区域地图分幅，易于拼接，主要用于工程建设大比例尺地形图。 （2）经纬线分幅（梯形分幅） 地图的图廓由经纬线构成。适用于大范围地图分幅，每个图幅都有明确的位置，但不易拼接，主要用于国家基本比例尺地形图。 1.1.2 图幅号的作用 地形图指的是地表起伏形态和地理位置、形状在水平面上的投影图。具体来讲，将地面上的地物和地貌按水平投影的方法，并按一定的比例尺缩绘到图纸上。随着计算机技术的应用和发展，许多企事业单位在生产管理中会用到大量的地形图，这些地形图以图号和图名来标识，人们通过查找索引表来找到所需要的图纸。总之，图幅号为了地形图的管理、保存而出现。 地图的分幅编号，在地图的生产、管理和使用方面都有重要意义。 首先是测制地图的需要，就测制某种比例尺地图而言，按每一分幅地图的范围和图号下达任务，不仅可以避免测制地图过程中遗漏或重复，节资增效，而且还能使所测地图的幅面控制在适当范围内，避免因幅面过大使绘图作业难以操作，影响绘图质量。 其次是印制地图的需要，若不分幅，地图幅面过大，一般印刷设备难以满足 要求 对教师党员的评价套管和固井爆破片与爆破装置仓库管理基本要求三甲医院都需要复审吗 ，势必要增加成本，而在复照时会给图面带来较大的边缘误差，影响地图的几何精度。 第三是管理和发行的需要，地图分幅编号后，便于分类分区有序地存贮；大小规格一致，易于包装。运输和存放；统一编号，有利于快速检索和发行。 第四是用图的需要，快速检索，有利于及时提供，提高工效；图面过大，折叠、携带和展阅都不便，只有将图面控制在一定大小范围内才便于在室内外的应用；分幅可以扩大地图的比例尺，便于更详细地表示各种地理要素，增加地图信息，以便更好地满足社会多方面的需求。 1.2 图幅号查询系统设计意义 国家基础地形图的使用单位大多已建立了数据库管理查询系统。查询地形图库依据的元素很多，如图名、档案号、图幅号、地区等，但最基本、最便捷的查询元素是坐标和经纬度，只要知道某点的坐标或经纬度，系统便可自动搜寻该点的图幅号，并进行生成。 国家基本比例尺地形图分幅编号的技术标准已实施多年，新绘制的地形图已采用新标准编制图幅号，而大量的旧地形图采用的是旧图幅号，新旧图幅号并存的局面将会继续存在较长一段时间。旧图幅号与新图幅号的换算较为烦琐，并且两者之间的转换在生产过程中具有重要意义，如土地二调采用新图幅编号方法，而一调时采用就图幅编号方法，为了调用后者的数据，新旧图幅号转换的实现就很有必要，本文设计的国家基本比例尺地形图图幅号查询系统实现了它们的转换与对照。 第二章 我国地形图图幅号编号方法 2.1 国家基本比例尺旧图幅号编号 （1）1:100万地形图分幅编号 1:100万地形图分幅和编号是采用国际标准分幅的经差6°、纬差4°为一幅图。从赤道起向北或向南至纬度88°止，按纬差每4°划作22个横列，依次用A、B、……V表示；从经度180°起向东按经差每6°划作一纵行，全球共划分为60纵行，依次用1、2……60表示。每幅图的编号由该图幅所在的“列号—行号”组成。例如，北京某地的经度为116°26’08’’、纬度为30°55’20’’，所在1:100万地形图的编号为J-50。 （2）1:50万、1:25万、1:20万、1:10万地形图分幅编号 1:50万、1:25万、1:20万、1:10万地形图均以1:100万图为基础分幅、编号。 将一幅1:100万地形图分为4幅,图幅经差为3°，纬差为2°，构成以A、B、C、D为代号的1:50万地形图，如：J-50-A； 将一幅1:100万地形图分为16幅，图幅经差为1°30’，纬差为1°，构成以[1]、[2]、…[16]为代号的1:25万地形图，如：J-50-[1]； 将一幅1:100万地形图分为36幅，图幅经差为60'，纬差为40'，构成以(1)、(2)、…、(36)为代号的1:20万地形图，如：J-50-(1)； 将一幅1:100万地形图分为144幅，图幅经差为30'，纬差为20'，构成以1、2、…、144为代号的1:10万地形图，如：J-50-92； (3)1:5万、1:2.5万、1:1万、1:5000地形图分幅编号 将一幅1:10万地形图分为4幅，构成以A、B、C、D为代号的1:5万地形图，图幅经差为15’，纬差为10’，如：J-50-92-A； 将一幅1:5万地形图分为4幅，构成以1、2、3、4为代号的1:2.5万地形图，图幅经差为7’30’’，纬差为5’，如：J-50-92-A-1； 将一幅1:10万地形图分为8行、8列，共64幅1:1万地形图，分别以(1)、(2)、…(64)表示，经差为3'45"、纬差为2'30"，如：J-50-5-(3)； 将一幅1:1万地形图分为2行、2列，共4幅1:5000地形图，分别以A、B、C、D表示，经差为1'52.5"、纬差为1'15"，如：J-50-5-(3)-A； 2.2 国家基本比例尺新图幅号编号 2.2.1 新图幅编号规则 1992年12月，我国颁布了《国家基本比例尺地形图分幅和编号GB/T139 89-92》新标准，1993年3月开始实施。新的分幅和编号方法如下： （一）地形图分幅 我国基本比例尺地形图均以1：100万地形图为基础，按规定的经差和纬差划分图幅。 在我国范围内（我国没有纬度60度以上的区域，60度以上区域经差纬差与低纬度不同）: 每幅1：100万地形图的经差是6度，纬差4度； 每幅1：50万地形图的经差是3度，纬差是2度，是将每幅1：100万地形图划分为两行两列，即一幅1：100万地形图包含四幅1：50万地形图； 每幅1：25万地形图的经差是1.5度，纬差是1度，是将每幅1：100万地形图划分为四行四列，即一幅1：100万地形图包含十六幅1：25万地形图； 每幅1：10万地形图的经差是30分，纬差是20分，是将每幅1：100万地形图划分为12行12列，即一幅1：100万地形图包含144幅1：10万地形图； 每幅1：5万地形图的经差是15分，纬差是10分，是将每幅1：100万地形图划分为24行24列，即一幅1：100万地形图包含576幅1：5万地形图； 每幅1：2.5万地形图的经差是7'30'，纬差是5'，是将每幅1：100万地形图划分为48行48列，即一幅1：100万地形图包含2304幅1：2.5万地形图； 每幅1：1万地形图的经差是3'45'，纬差是2'30''，是将每幅1：100万地形图划分为96行96列，即一幅1：100万地形图包含9216幅1：1万地形图； 每幅1：5000地形图的经差是1'52.5''，纬差是1'15''，是将每幅1：100万地形图划分为192行192列，即一幅1：100万地形图包含36864幅1：5000地形图。 （二）地形图编号 （1）1:100万地形图编号 1：1 00万地形图编号采用国际1：1 00万地形图编号标准。从赤道算起，每纬差4度为一行，至南北纬88度各分为22行，依次用大写拉丁字母（字符码）A、B、C、……V表示相应行号。从180度经线算起自西向东每经差6度为一列，全球分为60列，一次用阿拉伯数字1、2、3、……60表示其相应列号。由经线和纬线所围成的每一个小方格为一幅1：1 00万地形图，它们的编号由该图所在的行号和列号组合而成，例如北京所在的1：1 00万地形图图号为J 50。 （2）1:50万至1:5000地形图的编号 1：50万-1:5000地形图的编号均以1：100万地形图编号为基础，采用行列编号方法。即将1：100万地形图按所含各比例尺地形图的经差和纬差划分成若干行和列，横行从上到下、纵列从左到右按顺序分别用三位阿拉伯数字(数字码)表示，不足三位者前面补零，取行号在前、列号在后的排列形式标记；各比例尺地形图分别采用不同的字符作为其比例尺的代码(参见表格);1:50万-1：5000地形图的图号均由其所在1：1 00万地形图的图号、比例尺代码和各图幅的行列号，共十位码组成。如（图21） 图21 表21比例尺代码 比例尺 1:50万 1:25万 1:10万 1:5万 1:2.5万 1:1万 1:5000 代码 B C D E F G H 2.2.2 经纬度与图幅号转换关系 （一）已知某点经纬度，计算其编号 （1）按下式计算1:100万地形图图幅编号： a=[L/4°]+1 公式2-1 b=[B/6°]+31； 公式22 式中[ ]-------表示商取整； a-------1:100万地形图图幅所在纬度带字符码所对应的数字码； b-------1:100万地形图图幅所在经度带的数字码； B-------图幅内某点的经度； L-------图幅内某点的纬度； （2）按下式计算所求比例尺地形图在1:100万地形图图号后的行、列号： c=4°/△L-[(L/4°)/△L]； 公式23 d=[(B/6°)/△B]+1； 公式24 式中：[ ]-------表示商取整； ( )-------表示商取余； c-------所求比例尺地形图在1:100万地形图图号后的行号； d-------所求比例尺地形图在1:100万地形图图号后的列号； B-------图幅内某点的经度； L-------图幅内某点的纬度； △B-------所求比例尺地形图分幅的经差； △L-------所求比例尺地形图分幅的纬差； （二）已知图号计算该图幅西南图廓点的经纬度 按下式计算该图幅西南图廓点的经纬度： B=（b-31）*6°+（d-1）*△B； 公式25 L=（a-1）*4°+（4°/△L-c）*△L； 公式26 式中：B-------图幅内某点的经度； L-------图幅内某点的纬度； △B-------所求比例尺地形图分幅的经差； △L-------所求比例尺地形图分幅的纬差； a-------1:100万地形图图幅所在纬度带字符码所对应的数字码； b-------1:100万地形图图幅所在经度带的数字码； c-------所求比例尺地形图在1:100万地形图图号后的行号； d-------所求比例尺地形图在1:100万地形图图号后的列号； 2.3 新旧图幅号转换关系 2.3.1 新图幅到旧图幅的转换 （1）1:100万地形图从新图幅号换算旧图幅号 只要在新图幅号行列号之间加连接符‘-’即可，如新图幅号J50，其旧图幅号为J-50； （2）1:50万地形图从新图幅号换算旧图幅号 旧图幅号中的第3组代码——1:50万地形图在1:100万地形图中的位置代码“m”： m=（）*2+（）+1； 公式27 计算结果，m=1时取A，m=2时取B，m=3时取C，m=4时取D； x表示1:50万地形图新图幅号中的图幅行号，y表示1:50万地形图新图幅号中的图幅列号。 [ ]-------表示商取整； ( )-------表示商取余； （3）1:25万地形图从新图幅号换算旧图幅号 旧图幅号中的第3组代码——1:25万地形图在1:100万地形图中的位置代码“m”： m=（）*4+（）+1； 公式28 x表示1:25万地形图新图幅号中的图幅行号，y表示1:25万地形图新图幅号中的图幅列号。 [ ]-------表示商取整； ( )-------表示商取余； （4）1:10万地形图从新图幅号换算旧图幅号 旧图幅号中的第3组代码——1:10万地形图在1:100万地形图中的位置代码“m”： m=（）*12+（）+1； 公式29 x表示1:10万地形图新图幅号中的图幅行号，y表示1:10万地形图新图幅号中的图幅列号。 [ ]-------表示商取整； ( )-------表示商取余； （5）1:5万地形图从新图幅号换算旧图幅号 旧图幅号中的第3组代码——1:5万地形图所在的1:10万地形图在1:100万地形图中的位置代码“m”： m=[]*12+[]+1； 公式210 旧图幅号中的第4组代码——1:5万地形图在1:10万地形图的位置代码“n”： n=()*2+()+1； 公式211 计算结果，n=1时取A，n=2时取B，n=3时取C，n=4时取D； x表示1:5万地形图新图幅号中的图幅行号，y表示1:5万地形图新图幅号中的图幅列号。 [ ]-------表示商取整； ( )-------表示商取余； （6）1:2.5万地形图从新图幅号换算旧图幅号 旧图幅号中的第3组代码——1:2.5万地形图所在的1:10万地形图在1:100万地形图中的位置代码“m”： m=[]*12+[]+1； 公式212 旧图幅号中的第4组代码——1:2.5万地形图所在的1:5万地形图在1:10万地形图的位置代码“n”： n=()*2+()+1； 公式213 旧图幅号中的第5组代码——1:2.5万地形图在1:5万地形图的位置代码“p”： p=()*2+()+1； 公式214 计算结果，n=1时取A，n=2时取B，n=3时取C，n=4时取D； x表示1:2.5万地形图新图幅号中的图幅行号，y表示1:2.5万地形图新图幅号中的图幅列号。 [ ]-------表示商取整； ( )-------表示商取余； （7）1:1万地形图从新图幅号换算旧图幅号 旧图幅号中的第3组代码——1:1万地形图所在的1:10万地形图在1:100万地形图中的位置代码“m”： m=[]*12+[]+1； 公式215 旧图幅号中的第4组代码——1:1万地形图在1:10万地形图的位置代码“n”： n=()*8+()+1； 公式216 x表示1:1万地形图新图幅号中的图幅行号，y表示1:1万地形图新图幅号中的图幅列号。 [ ]-------表示商取整； ( )-------表示商取余； （8）1:5000地形图从新图幅号换算旧图幅号 旧图幅号中的第3组代码——1:5000地形图所在的1:10万地形图在1:100万地形图中的位置代码“m”： m=[]*12+[]+1； 公式217 旧图幅号中的第4组代码——1:5000地形图所在的1:1万地形图在1:10万地形图的位置代码“n”： n=()*8+()+1； 公式218 旧图幅号中的第5组代码——1:5000地形图在1:1万地形图的位置代码“p”： p=()*2+()+1； 公式219 计算结果，p=1时取A，p=2时取B，p=3时取C，p=4时取D； x表示1:5000地形图新图幅号中的图幅行号，y表示1:5000地形图新图幅号中的图幅列号。 [ ]-------表示商取整； ( )-------表示商取余。 2.3.2 旧图幅到新图幅的转换 （1）1:100万地形图从旧图幅号换算新图幅号 只要在新图幅号行列号之间取消连接符‘-’即可，如旧图幅号J-50，其新图幅号为J50； （2）1:50万地形图从旧图幅号换算新图幅号 新图幅号中的第4组代码——1:50万地形图的图幅行号“x”： x=[]+1； 公式220 新图幅号中的第5组代码——1:50万地形图的图幅列号“y”： y=()+1； 公式221 m表示1:50万地形图在1:100万地形图中的位置代码； m为A取1，m为B时取2，m=3时取C，m=4时取D； x,y不足3位前面补0； [ ]-------表示商取整； ( )-------表示商取余； （3）1:25万地形图从旧图幅号换算新图幅号 新图幅号中的第4组代码——1:25万地形图的图幅行号“x”： x=[]+1； 公式222 新图幅号中的第5组代码——1:25万地形图的图幅列号“y”： y=()+1； 公式223 m表示1:25万地形图在1:100万地形图中的位置代码； x,y不足3位前面补0； [ ]-------表示商取整； ( )-------表示商取余； （4）1:10万地形图从旧图幅号换算新图幅号 新图幅号中的第4组代码——1:10万地形图的图幅行号“x”： x=[]+1； 公式224 新图幅号中的第5组代码——1:10万地形图的图幅列号“y”： y=()+1； 公式225 m表示1:10万地形图在1:100万地形图中的位置代码； x,y不足3位前面补0； [ ]-------表示商取整； ( )-------表示商取余； （5）1:5万地形图从旧图幅号换算新图幅号 新图幅号中的第4组代码——1:5万地形图的图幅行号“x”： x=[]*2+[]+1； 公式226 新图幅号中的第5组代码——1:5万地形图的图幅列号“y”： y=()*2+()+1； 公式227 m表示1:5万地形图所在的1:10万地形图在1:100万地形图中的位置代码，n表示1:5万地形图在1:10万地形图的位置代码；计算中位置代码A代表1，B代表2，C代表3，D代表4； x,y不足3位前面补0； [ ]-------表示商取整； ( )-------表示商取余； （6）1:2.5万地形图从旧图幅号换算新图幅号 新图幅号中的第4组代码——1:2.5万地形图的图幅行号“x”： x=[]*4+[]*2+[]+1； 公式228 新图幅号中的第5组代码——1:2.5万地形图的图幅列号“y”： y=()*4+()*2+()+1； 公式229 m表示1:2.5万地形图所在的1:10万地形图在1:100万地形图中的位置代码，n表示1:2.5万地形图所在的1:5万地形图在1:10万地形图的位置代码；p表示1:2.5万地形图在1:5万地形图的位置代码，计算中位置代码A代表1，B代表2，C代表3，D代表4； x,y不足3位前面补0； [ ]-------表示商取整； ( )-------表示商取余； （7）1:1万地形图从旧图幅号换算新图幅号 新图幅号中的第4组代码——1:1万地形图的图幅行号“x”： x=[]*8+[]+1； 公式230 新图幅号中的第5组代码——1:1万地形图的图幅列号“y”： y=()*8+()+1； 公式231 m表示1:1万地形图所在的1:10万地形图在1:100万地形图中的位置代码，n表示1:1万地形图在1:10万地形图的位置代码； [ ]-------表示商取整； ( )-------表示商取余； （8）1:5000地形图从旧图幅号换算新图幅号 新图幅号中的第4组代码——1:5000地形图的图幅行号“x”： x=[]*16+[]*2+[]+1； 公式232 新图幅号中的第5组代码——1:5000地形图的图幅列号“y”： y=()*16+()*2+()+1； 公式233 m表示1:5000地形图所在的1:10万地形图在1:100万地形图中的位置代码，n表示1:5000地形图所在的1:1万地形图在1:10万地形图的位置代码；p表示1:5000地形图在1:1万地形图的位置代码，计算中位置代码A代表1，B代表2，C代表3，D代表4； x,y不足3位前面补0； [ ]-------表示商取整； ( )-------表示商取余。 2.3.3 图幅理论面积的计算 公式234 其中：A,B,C,D,E 为常数，按下式计算: ﹦(﹣)/ A﹦1+（3/6）+（30/80）+（35/112）+（630/2304） B﹦（1/6）+（15/80）+（21/112）+（420/2304） C﹦（3/80）+（7/112）+（180/2304） D﹦（1/112）+（45/2304） E﹦（5/2304） 式中：a—椭球长半轴(单位：米)，b—椭球短半轴(单位：米)； ΔL—图幅经差(单位：弧度)， (B2﹣B1)—图幅纬差(单位：弧度)，Bm﹦（B1﹢B2）/2； S—图幅理论面积（单位：平方米）。 第三章 系统设计与实现 3.1 系统设计原则 本系统的设计原则为简洁和准确。简洁要求系统界面友好，用户输入数据方便，系统代码简洁明了，效率高；准确指数据精度高，数据类型采用双精度，并且数据的格式符合国家规范中的要求，查询结果准确无误。 3.2 系统设计路线 图31 3.3 系统界面 本次图幅号查询系统是用VC++语言设计实现的，在主对话框中添加标签控件，设立四个标签：查询图幅号、查询图幅范围、新旧图幅号转换、关于。每个标签的界面如下： 图32查询图幅号界面 图33查询图幅范围界面 图34新旧图幅号转换界面 图35关于界面 3.4 系统功能模块 3.4.1 系统功能模块图 图36 3.4.2 功能模块代码 高斯坐标反算代码： GAOSI::GAOSI(double m1,double m2) //计算椭球参数 { a=m1; f=m2; b=(1-f)*a; c=a*a/b; ep=sqrt(a*a-b*b)/b; e=sqrt(a*a-b*b)/a; b0=1-3.0/4*ep*ep+45.0/64*mi(ep,4)-175.0/256*mi(ep,6)+11025.0/16384*mi(ep,8); b2=b0-1; b4=15.0/32*mi(ep,4)-175.0/384*mi(ep,6)+3675.0/8192*mi(ep,8); b6=-35.0/96*mi(ep,6)+735.0/2048*mi(ep,8); b8=315.0/1024*mi(ep,8); } double GAOSI::jtoh (double w,double w1,double w2) //角度转弧度 { double q0,q; q0=w+w1/60.0+w2/3600.0; q=q0*pi/180; return q; } double GAOSI::htoj (double m0) //弧度转角度 { double u; u=m0*180/pi*3600; return u; } double GAOSI::mi (double m1,int m2) //幂函数 { double q1=1; int s=1; while (s<=m2) { q1=m1*q1; s++; } return q1; } void GAOSI::an (double m) //计算a1---a6 { double n,v,N,sinh,cosh,t; sinh=sin(m); cosh=cos(m); t=tan(m); n=ep*cosh; v=sqrt(1+n*n); N=c/v; a1=N*cosh; a2=1.0/2*N*sinh*cosh; a3=1.0/6*N*mi(cosh,3)*(1-t*t+n*n); a4=1.0/24*N*sinh*mi(cosh,3)*(5-t*t+9*n*n+4*mi(n,4)); a5=1.0/120*N*mi(cosh,5)*(5-18*t*t+mi(t,4)+14*n*n-58*n*n*t*t); a6=1.0/720*N*sinh*mi(cosh,5)*(61-58*t*t+mi(t,4)); } double GAOSI::fun (double h) //计算Fx(B) { double h1; h1=c*sin(h)*(b2*cos(h)+b4*mi(cos(h),3)+b6*mi(cos(h),5)+b8*mi(cos(h),7)); return h1; } double GAOSI::funx (double h2,double h3) //计算Fx(B,l) { double h4,h5,h6,q2,q4,q6,h7; h4=ep*cos(h2); h5=sqrt(1+h4*h4); h6=c/h5; q2=1.0/2*h6*sin(h2)*cos(h2); q4=1.0/24*h6*sin(h2)*mi(cos(h2),3)*(5-tan(h2)*tan(h2)+9*h4*h4+4*mi(h4,4)); q6=1.0/720*h6*sin(h2)*mi(cos(h2),5)*(61-58*tan(h2)*tan(h2)+mi(tan(h2),4)); h7=q2*h3*h3+q4*mi(h3,4)+q6*mi(h3,6); return h7; } double GAOSI::funy (double h2,double h3) //计算Fy(B,l) { double h4,h5,h6,q3,q5,h7; h4=ep*cos(h2); h5=sqrt(1+h4*h4); h6=c/h5; q3=1.0/6*h6*mi(cos(h2),3)*(1-tan(h2)*tan(h2)+mi(h4,2)); q5=1.0/120*h6*mi(cos(h2),5)*(5-18*mi(tan(h2),2)+mi(tan(h2),4)+14*mi(h4,2)-58*mi(h4,2)*mi(tan(h2),2)); h7=q3*mi(h3,3)+q5*mi(h3,5); return h7; } void GAOSI::FSUAN (double m1,double m2,double m3) //反算迭代函数 { double B,aa,l,Bp,a1p,lp; B=m1/(c*b0); aa=a*cos(B)/sqrt(1-e*e*sin(B)*sin(B)); l=m2/aa; Bp=(m1-fun(B)-funx(B,l))/(c*b0); a1p=a*cos(Bp)/sqrt(1-e*e*sin(Bp)*sin(Bp)); lp=(m2-funy(Bp,l))/a1p; while ((Bp-B)>jtoh(0,0,0.0001)||(lp-l)>jtoh(0,0,0.0001)) { B=Bp; aa=a1p; l=lp; Bp=(m1-fun(B)-funx(B,l))/(c*b0); a1p=a*cos(Bp)/sqrt(1-e*e*sin(Bp)*sin(Bp)); lp=(m2-funy(Bp,l))/a1p; } fswd=htoj(Bp); fsjd=htoj(lp+m3); } 新旧图幅号转换代码： void prop3::OnButton1() { UpdateData(true); if(m_radio==0) //新图幅转旧图幅 { m_tf2=m_tf1.Left(1)+"-"+m_tf1.Mid(1,2); CString str; char *ch; ch=(LPSTR)(LPCTSTR)m_tf1; int m,x,y,n,p; x=atoi(m_tf1.Mid(4,3)); y=atoi(m_tf1.Right(3)); switch(ch[3]) { case 'B': m_blc="1:50万"; m=ykh((x-1),2)*2+ykh((y-1),2)+1; str.Format("%c",m+64); m_tf2=m_tf2+"-"+str; break; case 'C': m_blc="1:25万"; m=ykh((x-1),4)*4+ykh((y-1),4)+1; str.Format("%d",m); m_tf2=m_tf2+"-"+"["+str+"]"; break; case 'D': m_blc="1:10万"; m=ykh((x-1),12)*12+ykh((y-1),12)+1; str.Format("%d",m); m_tf2=m_tf2+"-"+str; break; case 'E': m_blc="1:5万"; m=fkh((x-1),2)*12+fkh((y-1),2)+1; str.Format("%d",m); m_tf2=m_tf2+"-"+str; n=ykh((x-1),2)*2+ykh((y-1),2)+1; str.Format("%c",n+64); m_tf2=m_tf2+"-"+str; break; case 'F': m_blc="1:2.5万"; m=fkh((x-1),4)*12+fkh((y-1),4)+1; str.Format("%d",m); m_tf2=m_tf2+"-"+str; n=ykh((fkh((x+1),2)-1),2)*2+ykh((fkh((y+1),2)-1),2)+1; str.Format("%c",n+64); m_tf2=m_tf2+"-"+str; p=ykh((x-1),2)*2+ykh((y-1),2)+1; str.Format("%d",p); m_tf2=m_tf2+"-"+str; break; case 'G': m_blc="1:1万"; m=fkh((x-1),8)*12+fkh((y-1),8)+1; str.Format("%d",m); m_tf2=m_tf2+"-"+str; n=ykh((x-1),8)*8+ykh((y-1),8)+1; str.Format("%d",n); m_tf2=m_tf2+"-"+"("+str+")"; break; case 'H': m_blc="1:5000"; m=fkh((x-1),16)*12+fkh((y-1),16)+1; str.Format("%d",m); m_tf2=m_tf2+"-"+str; n=ykh((fkh((x+1),2)-1),8)*8+ykh((fkh((y+1),2)-1),8)+1; str.Format("%d",n); m_tf2=m_tf2+"-"+"("+str+")"; p=ykh((x-1),2)*2+ykh((y-1),2)+1; str.Format("%c",p+64); m_tf2=m_tf2+"-"+str; break; case 0: m_blc="1:100万"; } } if(m_radio==1) //旧图幅转新图幅 { m_tf2=m_tf1.Left(1)+m_tf1.Mid(2,2); CString str; char *ch; ch=(LPSTR)(LPCTSTR)m_tf1; int m,x,y,n,p; if(m_tf1.GetLength()==4) { m_blc="1:100万"; } else if(ch[5]=='A'||ch[5]=='B'||ch[5]=='C'||ch[5]=='D') { m_blc="1:50万"; m=ch[5]-64; x=fkh((m-1),2)+1; y=ykh((m-1),2)+1; str.Format("B%03d%03d",x,y); } else if (ch[5]=='[') { m_blc="1:25万"; int first,last; first= m_tf1.Find('['); last= m_tf1.Find(']'); m = atoi(m_tf1.Mid(first+1,last-first-1)); x=fkh((m-1),4)+1; y=ykh((m-1),4)+1; str.Format("C%03d%03d",x,y); } else if (ch[5]=='(') { m_blc="1:20万"; m_tf2="国家基本比例尺新图幅号不包括1:20万"; } else if (chrn(ch,'-')==2) { m_blc="1:10万"; m=atoi(m_tf1.Mid(5)); x=fkh((m-1),12)+1; y=ykh((m-1),12)+1; str.Format("D%03d%03d",x,y); } else if(chrn(ch,'(')==0&&chrn(ch,'-')==3) { m_blc="1:5万"; str=m_tf1.Mid(5); int w0; w0=str.Find('-'); m=atoi(str.Left(w0)); n=ch[strlen(ch)-1]-64; x=fkh((m-1),12)*2+fkh((n-1),2)+1; y=ykh((m-1),12)*2+ykh((n-1),2)+1; str.Format("E%03d%03d",x,y); } else if(chrn(ch,'-')==4&&chrn(ch,'(')==0) { m_blc="1:2.5万"; str=m_tf1.Mid(5); int w0; w0=str.Find('-'); m=atoi(str.Left(w0)); p=atoi(m_tf1.Right(1)); n=ch[strlen(ch)-3]-64; x=fkh((m-1),12)*4+fkh((n-1),2)*2+fkh((p-1),2)+1; y=ykh((m-1),12)*4+ykh((n-1),2)*2+ykh((p-1),2)+1; str.Format("F%03d%03d",x,y); } else if(chrn(ch,'-')==3) { m_blc="1:1万"; str=m_tf1.Mid(5); int w0; w0=str.Find('-'); m=atoi(str.Left(w0)); int first,last; first= m_tf1.Find('('); last= m_tf1.Find(')'); n = atoi(m_tf1.Mid(first+1,last-first-1)); x=fkh((m-1),12)*8+fkh((n-1),8)+1; y=ykh((m-1),12)*8+ykh((n-1),8)+1; str.Format("G%03d%03d",x,y); } else { m_blc="1:5000"; str=m_tf1.Mid(5); int w0; w0=str.Find('-'); m=atoi(str.Left(w0)); int first,last; first= m_tf1.Find('('); last= m_tf1.Find(')'); n = atoi(m_tf1.Mid(first+1,last-first-1)); p=ch[strlen(ch)-1]-64; x=fkh((m-1),12)*16+fkh((n-1),8)*2+fkh((p-1),2)+1; y=ykh((m-1),12)*16+ykh((n-1),8)*2+ykh((p-1),2)+1; str.Format("H%03d%03d",x,y); } m_tf2+=str; } UpdateData(false); } 第四章 使用说明 4.1 编写目的 地图的分幅编号，在地图的生产、管理和使用方面都有重要意义。而地形图图幅号快速的查询，以及新旧图幅号之间的转换，对于生产单位的企业、用户来说就显得尤为重要。本手册可以让用户快速的了解和实现图幅号查询系统的各种功能，能让用户了解在处理过程中遇到问题的解决方法。《使用说明》适用与从事地形图管理的企业或者对地形图图幅号查询有需求的所有用户。 4.2 系统功能 本软件参考GB/T 13989-92国家基本比例尺地形图分幅和编号及旧地形图分幅和编号规则编制。主要功能如下： （1）根据经纬度或高斯坐标生成图幅号，包括1:100万、1:50万、1:25万、1:10万、1:5万、1:2.5万、1:1万、1:5000比例尺的图幅号； （2）由图幅号查询图幅经纬度范围，图幅理论面积等信息； （3）新旧图幅号之间的相互转换； （4）对每一模块的查询结果导出txt文件。 4.3 运行说明 4.3.1 启动方法 双击打开图幅号查询系统文件夹，3个DLL文件为动态链接库，双击文件名为“tfhcx”的exe文件即可打开图幅号查询系统。如（图41） 图41 4.3.2 查询 （一）由经纬度查询新图幅号 启动程序后，在主界面上，单击“由经纬度查询图幅号”复选框，选中后，输入目标地的经度和纬度，格式：90度7分34秒写成90.0734，然后选择右边的国家基本比例尺，点击生成图幅号，右下方的图幅号编辑框会出现生成的新图幅号。例如，山西省太原市市中心位于北纬37度52分，东经112度33分，比例尺选择1:1万，生成的图幅号为J49G052073。如（图42） 图42 （二）由高斯坐标查询图幅号 在主界面上，单击“由高斯坐标查询图幅号”复选框，选中后，选择参考坐标系（如北京54坐标系，或西安80坐标系等），输入目标地的高斯投影x,y坐标（y坐标加过500Km），并输入当地中央子午线经度，选择比例尺，点击生成图幅号。例如，北京市天安门中心在BJ-54坐标系中x=4419494.433933，y=447962.33797，中央子午线经度为117度，比例尺选择1:10万，生成的图幅号为J50D001005。如（图43） 图43 （三）查询图幅范围 点击查询图幅范围标签，出现查询图幅范围界面，在左边输入所要查询的新图幅号，图幅面积处选择一个参考系，点击查询图幅范围按钮，右边会出现图幅范围，图幅面积等信息。如输入上述图幅号J50D001005，结果如（图44）。 图44 （四）新旧图幅号转换 点击新旧图幅号转换标签，出现图幅号转换界面，在单选按钮处选择转换方式，输入原图幅号，点击转换图幅号按钮，右边的目标图幅号编辑框会出现转换后的图幅号。如选择转换方式为——>旧图号，输入上述图幅号J50D001005，转换的旧图幅号为J-50-5。如（图45） 图45 4.3.3 清除数据 在每个功能页面上查询完一次结果或者输入错误想重新查询，则点击每个页面上的清除数据按钮即可。 4.3.4 保存 本图幅号查询系统保存数据为txt文件，方法为点击每个页面上的保存按钮，弹出的对话框选择保存位置，输入文件名称，则会在相应位置生成一个txt文件。4.3.2（四）中图幅号转换的保存结果如（图46） 图46 4.3.5 退出 点击每个页面上的退出按钮即可退出该系统，或者点击界面右上角的“X”也可关闭该系统。 4.4 异常说明 由于某些点处在地形图分幅的边缘，可能在两幅图的公共边上，也可能为四幅图的公共点，此时查询出的图幅号就不止一个，会出现多个，如经度114度，纬度34度，比例尺为1:5万的图幅号为I49E012024，I50E012001，I49E013024，I50E013001，目标点在四幅图的公共点上。如（图47） 图47 参考文献 [1] 王永成,杨保山,侯来福. 浅谈地形图分幅与编号的计算方法[J]. 西部探矿工程. 2010(06) [2] 刘宏林. 国家基本比例尺地形图新旧图幅编号变换公式及其应用[J]. 测绘通报. 1998(08) [3] 项仲贞,许承权,方子岩. 中大比例尺地形图标准分幅与编号查询系统的实现[J]. 铁道勘察. 2008(01) [4] 张保纲北京市大比例尺地形图图幅号查询信息系统的研制[J] 解放军测绘学院学报1997-03-20 致 谢 光阴似箭，转眼间大学生活已接近尾声。在大学的最后阶段，毕业 论文 政研论文下载论文大学下载论文大学下载关于长拳的论文浙大论文封面下载 是我的重中之重，在论文完成的过程中，除了我自己的努力之外，也凝聚了很多人的心血。所以我要对帮助我完成论文的所有人表示感谢。 首先，我要对我的指导老师，表示诚挚的感谢。感谢我的老师，他以朋友的身份告诉我怎样完成这篇论文，怎样搜集资料。他在忙碌的工作中挤出时间来审查、修改我的论文，严格把关，循循善诱，在此我表示衷心感谢。 其次，我的同学也给予了热情的帮助，感谢他们为我提出有益的意见和建议，感谢他们对我无私的帮助。 最后，我要感谢我的父母，他们一直是我努力和学习的动力，感谢他们为我所付出的一切。 在写论文的过程中，尽管我付出了很大的心血，但由于知识和经验的不足，不当之处恳请各位老师给予指正。 外文原文 （出自：Wei Z, Yang B S and Li Q Q. 2012. Automated extraction of building footprints from mobile LIDAR point clouds. Journal of Remote Sensing, 16(2): 286–296） Automated extraction of building footprints from mobile LIDAR point clouds WEI Zheng, YANG Bisheng, LI Qingquan Abstract This paper presents a novel method for automated extraction of building footprints from mobile LIDAR point clouds.We rst generate the georeferenced feature image of mobile LIDAR point clouds using an interpolation method, and adopt imagesegmentation and contour extraction and tracing to extract building boundaries in the geo-referenced feature image as the coarselevel of building footprints in Two-dimensional imagery space. Then, the coarse level of building footprints is further rened byapplying planar segmentation on the extracted point clouds in the building boundaries. Finally, the triangulated irregular network (TIN) is used to achieve the ne level of building footprints. Dataset of residential areas captured by Optech’s LYNX mobile mapping system was tested to check the validities of the proposed method. Experimental results show that the proposed method provides a promising and valid solution for automatically extracting building footprints from mobile LIDAR point clouds. Key words: mobile LIDAR, footprints extraction, image segmentation, point cloud segmentation, triangulated irregular network 1 INTRODUCTION Mobile Laser Scanning (MLS) provides an efcient solution forcapturing surface data of roads, buildings and trees fast, and dy-namic. It has been widely used in many elds, such as digital city,basic surveying and mapping, urban planning, transportation, andenvironment protection. Some MLS systems have been developed,such as VLMS (Zhao & Shibasaki, 2005) of Center for Spatial In-formation Science, The University of Tokyo; StreetMapper (Street-Mapper, 2010) of Three-dimensional Laser Mapping and IGI; Lynx(LYNX, 2010) of Optech Inc. in Canada. The MLS system has many advantages in Three-dimensional reconstruction for large-scale urban scene. As the basic element of urban reconstruction, buildings can be detected and reconstructed from Aerial Laser Scanning (ALS) data, satellite images and aerial images, after their Two-dimensional boundaries or building foot prints have conrmed. For aerial LiDAR point clouds, footprints of buildings can be determinated according to roof surfaces detected from LIDAR point clouds; For image data, the building corners can be detected directly from images to form Two-dimensional boundaries of buildings, or using Two-dimensional cadastral map to recon struct buildings. Many researchs have been made to address this problem. Haithcoat, et al. (2001) used different thresholds of features (size, height, shape) to distinguish buildings and cars, and trees were eliminated by means of differential geometry. And, the Two-dimensional footprints of buildings were extracted from classied building objects in LIDAR point clouds. Wang, et al. (2006)firstly detected the boundary points of buildings using airborne LIDAR data and then simply linked these points to obtain a coarserepresentation of building footprints. To rene the extracted building footprints, they used Bayesian maximum estimation to embed the distances of the boundary points into the coarse polygon using prior knowledge of building edge directions. Ahmadi, et al. (2010)adopted modified active contour model to extract building footprints from high resolution aerial images. Unlike the traditional active contour model, the modied method can extract building footprints more precisely and control the affection from roads and trees around the buildings. Ortner, et al. (2007), Lafarge, et al. (2008)and Tournaire, et al. (2010) all used the marked point process to extract building footprints from aerial LIDAR point clouds. In this method, (Reversible-Jump) Monte Carlo Markov Chain (RJMCMC) and simulated annealing were iteratively used to extract building footprints. Due to the Digital Elevation Model (DEM) sampling space and image resolution, extracted footprints using above methods are deviated from the ground truth. However, it is fairly difficult to extract building footprints directly from MLS point clouds because of noise in the data, the huge volume of the data, and the absence of explicit or ancillary data (e.g.,image intensities). Manandhar and Shibasaki (2001, 2002) used the information of individual scan lines, geometric structures, and densities of points to classify mobile LIDAR points into roads and buildings. Abuhadrous, et al. (2004) classied points into buildings,roads, and trees by histograms of Z and Y directions in one prole.These approaches need the prior knowledge such as profiles andscan lines for classifying mobile LIDAR point clouds. However,for large scale MLS point clouds, these methods are hard to work and less efficient, due to topographic relief and multiple objects. We propose a coarse-to-ne approach which automatically and effectively extracts building footprints from mobile LIDAR point clouds in urban environment. The proposed method firstly generates the georeferenced feature image of mobile LIDAR point clouds,and then extracts the coarse boundaries of buildings by applying the method of image segmentation and contour tracing on the georeferenced feature image generated. Once the boundaries of buildings are extracted, the Three-dimensional points associated with buildings could be easily isolated and extracted from mobile LIDAR point clouds. Then building facades are extracted to determinate the fine footprints of buildings. Fig. 1 elaborates the flow chart of the proposed method. MLS Point Cloud Building footprint extaction Georeferenced feature image generation Facade TIN model Threshold segmentation and object classification Ext raction of building point clou ds Contour extraction and edge tracing Fig. 1 The owchart of the proposed method 2 BUILDING OBJECT EXTRACTION Because dealing directly with large-volume point clouds is time-consuming, mobile LIDAR data points are projected into an X-Y plane in Two-dimensional space to generate a georeferenced feature image of the point clouds. More specically, the georeferenced feature image generated should preserve the spatial distribution and density of the scanned points and it must be able to reect the relative intensities of buildings. Therefore, existing operators for extracting key line or point features from imagery can be easily adapted to extract building footprints from the generated feature image. The efficiency of building footprint extraction is thereby improved. To fulfill these requirements, the grey values of each pixel and the resolution of the georeferenced feature image should be determined. 2.1 Georeferenced feature image generation The proposed method first projects the scanned points onto the X-Y plane. The minimum bounding box can be easily determined from the coordinates of the scanned points. Let (Xmin, Ymin) and (Xmax, Ymax) be respectively the coordinates of the bottom left and top right corners of the georeferenced feature image, and let the Ground Sample Distance (GSD) be the spatial resolution. The width and height of the georeferenced feature image can be calculated as below W=(Xmax-Xmin)/GSD; H=(Ymax-Ymin)/GSD; (1) Once the width and height of the georeferenced image have been determined, the mobile LIDAR point clouds can be located in each cell (see Fig. 2). To generate a georeferenced feature image, the grey value of each cell must be calculated. The spatial distribution pattern of the mobile LIDAR point clouds implies the local geometric features of street-scene objects. Therefore, the grey values of the georeferenced feature image should reect the spatial distribution pattern of the mobile LIDAR point clouds. Fig. 2 Sketch map of cells and data points Once the weight of each point has been determined, the grey value of each cell can be calculated as follows using the Inverse Distance Weighted (IDW) interpolation method: Wijk=α*WXY ijk+β* WH ijk (2) WXY ijk=*GSD/Dk ij α+β=1.0 (3) Where nij is the number of scanned points in cell (i, j), H is the height difference from the point in cell (i, j) to the lowest point in the cell, hmin(ij) and hmax(ij) are the lowest and highest elevations respectively in cell (i, j), and Zmax and Zmin are the lowest and highest elevations respectively in the whole scanned area. δ is a positive innitesimal value to make the denominator of the equation nonzero.The point weights calculated by the above equations emphasize high objects and objects with a large standard deviation of height values. (x0 , y 0 ) is the coordinate of the center of cell (i, j), and α and β are weighting coefcients that require tuning. Eq. (2)—Eq. (3) demonstrate that the grey value of each cell in the georeferenced feature image is determined both by the planar distance between each point and the center of the cell and by the height difference between each point and the lowest-height point in the cell. Thus, a georeferenced feature image can be generated for specied values of α and β. 2.2 Coarse footprints and building extraction In light of the generation of the georeferenced feature image,the areas with points of higher elevations have larger gray values in the georeferenced feature image. Hence, the areas of buildings and trees have larger gray values compared with their neighbor areas(e.g., ground, roads). A threshold of gray value can be specied to extract the areas of larger gray values according to the histogram analysis on the georeferenced feature image. To eliminate the noise in the segmented georeferenced feature image and extract street-scene objects (e.g., buildings, trees), we dene the following operators. (1) Threshold segmentation In this study, we selected the method of discrete discriminant analysis (Goldstein & Dillon, 1978) to automatically specify the threshold because of its robustness. First, the method of discrete discriminant analysis nds the maximum and minimum gray value of the feature image. Then, it calculates the between-cluster variance under an assumed threshold. Finally, it nds the threshold with the highest value of between-cluster variances as the optimal threshold, which is used to classify the pixels into different categories. (2) Boundary extraction As we deal with segmented binary feature image (the background is black), we implement contour extraction in a simply way. A pixel is changed into white if it is a black pixel and all the 8-neighbor pixels are black as well. In this way, the black pixels whose 8-neighbor pixels are not all black will not change and be labeled as boundary pixels. After the contours being extracted, a contour tracing algorithm proposed by Pavlidis (1982) is applied to identify each contour as a sequence of edge points. (3) Size constraint Considering the sizes of trees and buildings, the extracted boundary will be eliminated if its size is less than a threshold. This threshold should be tuned with different GSD. Since all the contours are closure, we count the length of each contour as its size in our experiment. (4) Shape constraint Generally, tree shows an ellipse or circle shape and building shows a rectangle shape. Compactness of polygon shape can be applied to distinguish buildings and trees. We use Eq. (4) in Touya (2007) to calculate the compactness of each contour shape. Compactness=4π (4) where S and P are the area and perimeter of a contour, respectively. In this equation, the compactness would be 1 for a circle contour. Hence, the compactness of trees would be near 1, while the one of buildings not. Once the generated georeferenced feature image is segmented and the above three operators are invoked, the boundaries of streetscene objects (e.g., buildings) can be extracted. Hence, the spatial extension of the extracted street-scene objects can be easily calculated according to the extracted boundaries associated. Therefore,the point clouds corresponding to each extracted object can be extracted from the mobile lidar point clouds according to the spatial extension calculated. 3 BUILDING FOOTPRINT EXTRACTION Due to the DEM sampling space and image resolution, extracted footprints using above methods are deviated from the ground truth.Therefore, point clouds of extracted building should be used to determinate the ne locations of building footprints. The operations are as follows. (1) Using RANSAC (von Hansen, et al., 2006) to detect planar patches from building point clouds, these patches could be facades,and roofs of a building. (2) Build up a TIN model for each extracted facade of the building, and then extract the boundary of the TIN model. (3) Calculate the intersection line of different facade planes, and determinate the exact building corners of facades according to the intersection lines and the TIN boundaries. (4) Using prior knowledge (rectangle boundary) of building skeleton model, calculate the boundary behind the building front to refine the exact building footprints. 4 EXPERIMENTS AND RESULTS One dataset captured by Optech’s LYNX Mobile Mapper, was selected for assessing the performance of the proposed method. The test dataset consists of about 9 million points covering a 400 m by 350 m street scene. The density of point clouds is about 1—5 cm.The mobile lidar point clouds contain a lot of street-scene objects(e.g., buildings, trees and power lines). 4.1 Georeferenced feature image In the proposed approach, a georeferenced feature image is generated from mobile lidar point clouds before the street-scene objects extraction operator is invoked. For the georeferenced feature image generation, two parameters should be determined, namely,the spatial resolution, also sometimes designated as the GSD, and the distance weight α (height difference weight β=1–α). As the density of the mobile lidar point clouds is about 1 cm to 5 cm, we generated the georeferenced feature image with cell size of 0.25 m. As shown in Fig. 3, the areas of buildings, trees, and power lines have higher gray values compared with other areas (e.g., roads and ground). Moreover, the areas of buildings, trees, and power lines have prominent boundaries that are well preserved in the feature images. Simultaneously, the differences between different streetscene objects are more notable. It implies that the generated georeferenced feature image is able to well depict geometric features of street-scene objects and is suitable to extract them. Fig. 3 Georeferenced feature image (GSD = 0.25, α = 0.2, β = 0.8) 4.2 Coarse footprints extraction The generated feature image was firstly segmented by the method of discrete discriminant analysis to maintain the point clouds of street-scene objects. Then, a contour extraction and tracing operator was invoked to extract the coarse boundaries of streetscene objects. The size constraint was set to be 100 pixels, and the closed red boundaries were extracted coarse footprints of objects,as shown in Fig. 4. Fig. 4 Registration of coarse footprints and segmented feature image 4.3 Plane segmentation of building point clouds According to the extracted coarse footprints of buildings, the point clouds corresponding to each building can be extracted from the mobile LIDAR point clouds. In this experiment, 53 objects (50 buildings and 3 trees) were extracted, and each building was segmented using RANSAC. Fig. 5 shows the plane segmentation result of a selected building. Fig. 5 The plane segmentation result of selected building 4.4 Fine footprints extraction In this paper, 50 buildings were segmented and their facades were extracted to build up TIN models. The TIN boundary was considered as the coarse skeleton of each facade. Then intersection lines of different facades were calculated to determinate the exact footprints of buildings, using orthogonal relationship of footprints. Fig. 6 illustrates the extracted fine footprints and skeleton of the building in Fig. 5. Fig. 7 shows the extracted footprints of 47 building objects (three buildings were missed due to the lack of scanning data. The extracted footprints were transformed into KML le and output into Google Earth software. Fig. 6 The extracted ne footprints and skeleton of the building in Fig. 5 Fig. 7 The extracted footprints of 47 building objects 5 CONCLUSION We propose a promising solution for extracting point clouds of street-scene object from mobile LIDAR point clouds. The proposed method transforms the extraction of street-scene objects from Three-dimensional point clouds to the object extraction from Two-dimensional imagery. The image segmentation and contour extraction techniques were applied to extract the boundaries of street-scene objects in image space. Hence, the point clouds of the street-scene objects (e.g., buildings and trees) can be extracted in Three-dimensional according to the boundaries associated.Then building facades were extracted to determinate the fine footprints of buildings. Our method reduces computational cost in the street-scene object extraction from mobile LIDAR point clouds. The georeferenced feature image generated from mobile LIDAR point clouds is able to represent the spatial distribution of scanned points and preserve the local geometric features of street-scene objects. The results achieved appear encouraging and demonstrate that the proposed method provides an effective solution for extracting building facade footprints from mobile LIDAR data. Some issues remain to be addressed in future research. Reection intensities of point clouds and associated imagery will be incorporated for better point-cloud mapping and understanding, which will improve the extraction of facade footprints. 中文翻译 （译自：Wei Z, Yang B S and Li Q Q. 2012. Automated extraction of building footprints from mobile LIDAR point clouds. Journal of Remote Sensing, 16(2): 286–296） 车载激光扫描点云中建筑物边界的快速提取 魏征，杨必胜，李清泉 摘 要 以车载激光扫描点云数据为研究对象，提出一种由粗到细且快速获取点云中建筑物3维位置边界的方法。首先，通过分析格网内部点云的空间分布特征(平面距离、高程差异和点密集程度等)确定激光扫描点的权值，采用距离加权倒数IDW(Inverse Distance Weighted)内插方法生成车载激光扫描点云的特征图像。然后，采用阈值分割、轮廓提取与跟踪等手段提取特征图像中的建筑物目标的粗糙边界。最后，对粗糙边界内部的建筑物目标点云进行平面分割，提取建筑物的立面特征并构建立面不规则三角网TIN(Triangulated Irregular Network)，并在建筑物先验框架知识条件下自动提取建筑物的精确3维位置边界。 关键词：车载激光扫描，边界提取，图像分割，点云分割，不规则三角网 1 引言 车载激光扫描系统能够在高速移动状态下获取道路以及道路两侧建筑物、树木等地物的表面数据，已成为一种快速的空间数据获取手段，广泛运用于数字城市、基础测绘、城市规划、交通和环保等领域。国内外的研究机构与公司也相继推出了一些车载激光扫描系统，如日本东京大学空间信息科学中心研制的VLMS(Zhao和Shibasaki，2005)，three-dimensionalLaser Mapping和IGI公司合资开发的StreetMapper系统(StreetMapper，2010)，加拿大Optec公司的Lynx系统(LYNX，2010)，中国的山东科技大学、南京师范大学、首都师范大学、中科院深圳先进研究院也相继研制了车载激光扫描系统(卢秀山和黄磊，2007；江水 等，2007)。 车载激光扫描系统在大规模城市场景的三维重建中具有越来越明显的优势。建筑物作为城市三维建模的主要元素，无论是从机载激光扫描数据，还是从卫星影像以及航空影像数据来进行建筑物的三维建模，都需要首先确定建筑物的二维边界信息。对于机载点云数据来说,可以根据点云中提取的建筑物的顶面信息来确定建筑物的二维边界，即足印(footprint)。很多学者对此也作了深入的研究。 Haithcoat等人(2001) 采用大小、高程、形状等特征，通过不同的阈值区分建筑物和车辆目标，实现LIDAR点云中建筑物边界的提取。Wang等人(2006)提出了一种采用LIDAR点云数据区域增长自动提取建筑物边界的方法。Tournaire等人(2010) 采用基于目标的标记点过程方法，利用(可逆转)蒙特卡洛马尔科夫链(MCMC)和模拟退火相结合的方法对LIDAR点云DEM数据进行建筑物边界的提取。上述方法由于受到LIDAR点云DEM采样间隔以及影像分辨率的影响，其提取的建筑物边界对应于实际的建筑物地面边界还存在一定差距。 对于车载激光扫描数据来说，由于其存在大量竖直面上的点云信息，直接从点云数据中提取建筑物目标及其边界将面临一定困难。Manandhar和Shibasaki(2001，2002) 根据每个断面扫描点的点位空间分布特征，对扫描点进行分类，可以将建筑物、道路和树木等初步分离。Abuhadrous等人(2004) 通过断面直方图分析以及Z、Y方向投影来实现目标分类。但此类方法需要扫描数据的断面信息，因而难以处理散乱的点云。对于大范围内车载激光扫描数据而言，由于受到地形起伏和多目标的影响，上述方法难以奏效，而且效率低下。 针对上述问题，本文提出了一种车载激光扫描数据建筑物边界由粗到精的快速提取方法。该方法首先通过分析点云的空间分布特征生成点云的特征图像。其次，结合图像处理的方法对点云特征图像进行分割和特征提取，从而实现建筑物粗糙边界范围的确定，并利用该粗糙边界提取建筑物点云数据，对其进行平面特征提取，利用提取的立面信息确定建筑物的精确边界。 2 车载点云的建筑物目标提取 由于车载激光扫描获取的点云数据量大，难以直接对点云的几何数据进行分类和特征提取。本文从车载激光扫描散乱点云数据出发，将三维点云数据的分类与目标提取问题转移到二维图像中进行，通过生成的点云特征图像，利用图像阈值分割、轮廓提取、边界跟踪等方法提取建筑物的粗糙边界，然后反算到点云数据中获取粗糙边界范围内的建筑物点云数据，实现建筑物目标的提取。其处理 流程 快递问题件怎么处理流程河南自建厂房流程下载关于规范招聘需求审批流程制作流程表下载邮件下载流程设计 如图1所示。 原始 散乱点 云 点 云特征 图像生成 阙值分割 、目标分类 建筑物边界提取 轮廓提取、边界跟踪 构建立面三角网 建筑物点云提取 图1 车载点云建筑物边界提取流程 2.1 点云特征图像的生成 为了生成整个扫描区域的点云特征图像，将整个点云区域投影到XOY平面，需要首先确定生成图像的大小以及每个像元格网对应三维空间中的采样间隔，然后利用每个格网内的点云数据确定格网像元的灰度值，最后得到整个扫描区域某一地面采样间隔下的点云特征图像。 (1)确定点云特征图像的宽W、高H以及地面采样间隔GSD 通过点云数据，可以得知扫描区域的最大最小X、Y、Z坐标分别为：Xmin、Ymin、Zmin、Xmax、Ymax、Zmax，用户自定义地面采样间隔GSD如图2，则有： W=(Xmax-Xmin)/GSD; H=(Ymax-Ymin)/GSD; (1) 图2 格网示意图 (2)确定每个格网(i, j)的特征值Fij 假设落在第(i, j)个格网中的激光扫描点个数为nij，格网(i, j)的中心点为p0(x0, y0, 0)，利用nij个扫描点的三维坐标加权计算格网(i, j)的特征值Fij。显然，每个格网的特征值受落在该格网内的扫描点的个数、空间分布形式(即平面距离、高程差异等)决定。扫描点距离格网中心点越远则其权值越小；扫描点的高程值越大则其权值越大。 为了确定落在单元格(i, j)内每个点(如第k个点，0＜k≤nij)的权值Wijk。本文将格网的特征值Fij的计算分为两个部分：第一部分由格网中所有点与格网中心点之间的XOY平面距离Dk ij决定；第二部分由格网中所有点与格网中最低点之间的高程差异Hk ij决定。因此对于单元格(i, j)内的扫描点k而言，其权值可以用式(2)和式(3)表示： Wijk=α*WXY ijk+β* WH ijk (2) WXY ijk=*GSD/Dk ij α+β=1.0 (3) 式中，WXY ijk、WH ijk分别为扫描点k与格网点中心的距离的权值以及扫描点k高程的权值。WXY ijk反映了离散点与格网中心点的平面距离对格网中心点特征值的贡献。WH ijk反映了格网内离散点的高程对格网中心点特征值的贡献。 (3)生成点云特征图像 将格网特征值Fij归化到0—255灰度空间即可得到格网对应的点云特征图像的象元值G ij，从而得到整个扫描区域的点云特征图像。 2.2 建筑物粗糙边界确定与建筑物目标提取 根据点云特征图像的生成方法可知，点云特征图像反映了散乱激光扫描点云的空间分布规律与地物的几何特征，因此，可以充分利用图像处理的一些方法，如阈值分割等，对点云特征图像进行分割和特征提取，从而实现点云数据的分类以及目标提取。 进行阈值分割后的图像中蕴涵了大量的几何地物，如建筑物、树木等。对分割后的二值图像进行轮廓提取以及边界跟踪，从而实现建筑物和树木轮廓区域的粗糙边界特征提取。 建筑物2维轮廓边界提取主要采用如下4个步骤: (1)阈值分割 考虑到生成的点云特征图像中建筑物、树木与道路目标之间像元灰度值的差异，本文采用离散判别分析(Goldstein和Dillon，1978)方法，通过计算最大类间方差，自动确定图像分割的阈值。阈值分割后可以将建筑物和树木目标与道路目标分离。 (2)轮廓提取与跟踪 本文中阈值分割后的图像为二值图像，作者将像元灰度值为0且其八邻域像元灰度值均为0的像元灰度值都修改为255，实现轮廓提取。然后采用Pavlidis(1982) 提出的轮廓跟踪算法对提取的建筑物和树木轮廓进行跟踪。 (3)轮廓尺寸约束 提取的轮廓中会存在一些细小噪声轮廓，本文采用轮廓边界点的个数即轮廓周长对轮廓进行约束，从而删除轮廓尺寸小于某一阈值的轮廓，保留建筑物和树木轮廓。 (4)轮廓形状约束 经过尺寸约束后，图像中保留了树木和建筑物的轮廓，为了提取建筑物的粗糙边界，需要将树木和建筑物轮廓分开。由于树木轮廓的椭圆特征以及建筑物轮廓的块状特征，本文采用Touya(2007)提出的轮廓紧凑度公式如式(4)所示： Compactness=4π (4) 式中，S和P分别为轮廓的面积和周长，圆的紧凑度(Compactness)为1。因此树木轮廓的紧凑度会接近于1，而建筑物轮廓的紧凑度将与圆有较大差距。经过形状约束后，就可以提取特征图像中所需要的建筑物轮廓作为建筑物的粗糙二维边界。 如果在点云特征图像中找到建筑物目标的粗糙边界，即可根据点云特征图像与几何数据间的映射关系，实现粗糙边界内部建筑物点云数据的提取。 3 建筑物点云的边界提取 利用点云特征图像提取的建筑物粗糙边界，由于受到图像采样间隔以及图像分割结果的影响，与建筑物的真实地理位置的边界还存在较大的差距，因此需 要 根 据 提 取 的 建 筑 物 目 标 三 维 点 云 数 据 ，从中提取建筑物较精确的位置边界信息。具体步骤如下。 (1)对点云数据进行RANSAC平面分割(vonHansen 等，2006)，提取建筑物的不同立面、顶面等面片特征； (2)对建筑物的不同立面构建面片三角网，提取三角网的边界； (3)计算不同立面之间的交线，并根据三角网边界，获取准确的建筑物立面角点信息，并作为建筑物的精确边界角点； (4)根据建筑物框架模型的先验知识(矩形边界)，利用已有立面的地面边界，推算建筑物被遮挡面的轮廓边界，从而精化建筑物精确边界。 4 实验结果 本实验采用Optech公司的激光扫描系统LYNX Mo-bile Mapper(该系统由2个激光扫描仪、2个CCD相机、Applanix POS LV 420系统以及操作软件平台构成，扫描仪视场角360°，每秒获取100000个扫描点数据)采集的居民区扫描数据，其扫描范围大致为400 m×350 m，扫描点的密度大致为0.01—0.05 m，共8139716个数据点。 4.1 点云特征图像生成结果 根据式(1)可知，点云特征图像的生成需要特征图像的分辨率，即格网的大小。由于扫描点云的密度在0.01—0.05 m之间，为了保证落入格网内部的点数大于10，所以格网采样间距采用0.25 m进行实验。 本文采用杨必胜等人(2010)提出的α、β值的取值原则，生成如图3所示的点云特征图像。从图3可以看出，没有点数据的格网像元灰度均为0(采用背景色蓝色表示)，包含道路数据的格网像元灰度较低，而包含建筑物和树木以及电线数据的格网像元灰度较高，且图像较好地保持了建筑物、树木等的边界特征以及建筑物部分顶面信息。 图3 点云特征图像(GSD = 0.25，α = 0.2，β = 0.8) 4.2 建筑物粗糙边界提取结果 本文对图3的点云特征图像进行阈值分割，并对分割后的图像采用轮廓提取、边界跟踪的手段得到了建筑物区域以及树木区域的图像粗糙边界，其中轮廓尺寸约束阈值为100个像元，如图4所示。图4中红色闭合边界为提取的目标粗糙边界。 图4 目标粗糙轮廓边界与分割图像的套合图 4.3 建筑物点云平面分割结果 根据提取的建筑物粗糙边界，通过点云特征图像和区域点云之间的映射关系，可以提取不同建筑物边界内部的三维点云数据。实验中对图像中的53个目标进行了提取，并对提取的单个建筑物点云进行了RANSAC平面分割，图5为选取的某一建筑物目标的平面分割效果图。 图5 目标建筑物的平面分割结果图 4.4 建筑物精确边界提取 本文对50个建筑物目标依此进行了平面分割，并提取了不同目标的立面信息。对不同立面构建三角网并提取三角网的边界作为立面的初始框架边界，同时计算不同立面之间的交线位置及方向，根据建筑物自身的正交几何关系确定了建筑物的精确边界。图6为对图5中建筑物提取的边界，图7为47个建筑物目标提取的边界(其中3个建筑物因扫描数据的严重不足难以确定其边界)。本文通过构建立面TIN，寻找立面的三角网边界点，同时计算不同立面之间的交线信息，可以克服遮挡导致的立面扫描数据缺失问题。 图6 建筑物的边界与立面框架效果图 图7 提取的47个建筑物目标的边界结果 5 结论与展望 本文以车载激光散乱点云数据为研究对象，采取由粗到细的策略，提出了一种车载激光扫描数据建筑物边界的快速提取方法。 该方法，避免了大量的几何计算，而且充分借鉴了图像处理的一些方法，很大程度上提高了数据处理的效率，从而能够实现建筑物边界的快速提取。同时，本文通过构建立面TIN，计算立面交线，有效地克服了建筑物立面由于遮挡导致的扫描信息不足的问题。对Optech公司的扫描数据进行实验的结果表明，提取的边界很好地符合了建筑物平面位置轮廓，为后期的建筑物立面建模提供了很好的数据支持。本文提出的建筑物边界提取方法对于居民区比较实用，对于闹市区(存在很多高大建筑物)的建筑物边界提取可能还有一些局限，这也是本文需进一步改进的地方。 由于本文实验数据中缺乏实验区域的2维线划图信息，并不能准确地计算提取出来的建筑物边界与实际建筑物边界之间的位置偏差，因此才将提取的建筑物边界信息放在Google Earth上与影像图进行叠加，体现提取的建筑物边界的直观吻合效果。这也是本文提取方法论证上的不足之处。