1 引言 1.1 研究背景与意义

激光雷达即激光探测与测距技术(Light Detection And Ranging, LiDAR)是一种逐渐被广泛应用的测绘技术^[1]. 激光测距系统可以实现传感器向测量目标发射激光束, 并记录发射时刻, 激光束到达目标后按照同一光路反射, 再由传感器接收, 并记录接收时刻, 那么激光器至测量目标的距离即可通过公式(距离=光速×时间差/2)计算得出. 将此激光器安装在机电设备上, 控制其按一定角度水平或竖直摆动, 此时可称其为激光扫描仪. 再将其与高精度动态载体姿态测量系统(INS)、高精度动态GPS差分定位系统(DGPS)相结合, 即可构成一个三维激光扫描系统, 经此系统可获取到ASCII格式的三维点云数据, 记录每个点三维坐标以此表达目标地物三维空间信息. 同时, 按照载荷平台的不同可将其分为地面激光雷达、机载激光雷达、车载激光雷达和星载激光雷达. 激光雷达较于传统遥感测量手段具有几项明显的优势: (1) 直接获取三维坐标信息; (2) 采用主动式测量, 不受光照影响; (3) 复杂地形对雷达测量影响小; (4) 根据激光的同一脉冲的多次反射对数据进行分层可去除植被影响, 只获得地表数据, 提高精度. 因此, 激光雷达测量技术被广泛用于各地学研究领域, 机载雷达主要应用于林业、城市规划水利等, 车载雷达主要是对地物的侧面进行激光扫描, 星载雷达则是用于植被、极地冰川等研究, 地面雷达可用于工程测量、建筑物建模等领域.

将激光雷达测量技术用于测绘行业在近几年刚开始起步, 而将其用于地形测绘领域也是一个新的探究方向^[2]. 传统的地形测量必须携带GPS仪器、全站仪等众多负重设备, 还要在许多环境恶劣的野外进行测量, 再经内业处理后才能生成地形图, 制作周期长, 工作强度和危险性都较大. 若能将激光雷达这种遥感测量方式应用于现在测量行业不仅可以简化外业流程, 缩短工作时间, 也可提高内业自动化程度. 基于此, 本文探讨将某处地形的LiDAR点云数据作为原始数据, 获取此处地形的等高线.

1.2 国内外研究进展

利用点云数据绘制等值线在国内外已有一定的研究, 一些学者采用软件平台进行了一些实验; 冯梅等^[3]分别用GlobaMapper, GeoTIN, Tindem, JX-4C, VirtuoZo和TerraScan六种相关软件提取某块山区等高线; 李鸿轶等^[4]使用TerraSolid提取等高线; 白志远^[5]则将ArcGIS与LP360相结合提取等高线. 这些实验均基于一些山区地形进行, 且并没有选择大区域数据进行实验, 使用软件进行数据的处理虽有一定的便利性, 但对于等值线生成过程中的变量控制性较差, 且缺乏普适性. 而也有一些学者研究算法, 通过各种方法生成等值线, 如吴杭彬等^[6]采用迭代凸包的方法提取激光扫描数据的分层多细节等值线, 但要经过多次迭代凸包, 算法效率低. 王宗跃等^[7]建立三角网跟踪等高点, 采用样条函数内插等高线; Boissonnat等^[8]提出了基于TIN的等高线跟踪方法. 这些算法均有其各自的特色和优点, 但都只是对固有数据加工生成等值线, 而未生成不同比例尺下的等高线. 但在实际应用中, 随着比例尺的减小, 地形图势必会综合一些信息, 但为了体现地形的形态, 应保留等值线特征处, 才具有清晰性和精确性. 因此, 探究如何控制生成等值线的变量以期绘制表达更实际的等值线, 以及可以根据不同比例尺的设定, 绘制相应的等值线.

2 研究基础与方法 2.1 研究基础

LiDAR系统在提供数据更新便利的同时, 带来的是点云数据的海量性. 数据量越大对内存要求越高, 因此便会影响实际生成等值线的效率. 同时, 若想针对不同比例尺地形情况下, 均有其对应的等值线形态, 也需要对源数据进行处理, 判断哪些点对地形特征描述更为重要, 若前期对此重要性进行排序, 则后期即可根据比例尺确定某一个节点, 以此为界, 则即可在特定比例尺下, 绘制特征表述完全的等值线. 若要在保证其精度的同时能够删去冗余点, 费立凡、何津^[9]等曾在2006年提出了三维Douglas-Peucker算法, 使用此法综合随机的三维离散点, 无论进行多大尺度的综合, 保留的点均是地貌特征点, 相对次要的点将被删除. 2013年, 何津、费立凡等^[10]提出了基于三维Douglas-Peucker算法综合三维离散点, 再回放等高线, 但其根据某个确定的阈值选取点集合, 阈值的确定有一定的主观性. 在本文中, 对三维Douglas-Peucker算法进行改进, 第一步不是选取所需要的点子, 而是对点子进行重要性排队, 对于地形的描述而言, 越重要的点越排在前列. 在排序之后根据不同尺度下等高线所需要的综合程度确定所选点子的数量. 在这种方法里确定取舍点的分界线可以称其为“剁尾刀”的定位, 在多尺度输出的情况下可以提高几十倍的效率, 并且需要输出的尺度越多, “剁”得越多, 效率越高.

使用三维Douglas-Peucker算法可以实现三维点数据的综合, 但要绘制成等值线要经过点数据的栅格化, 进而绘制. 于彩霞等^[11]提出将点云数据先栅格化后提取海岸线防止海岸线的一些曲折及破碎问题. 基于此, 本文选择将点云数据先插值生成DEM, 进而绘制等值线.

2.2 研究方法

(1) 在原有三维Douglas-Peucker算法的基础上进行改进, 加速算法运行的速度以提高效率. 将原始算法中对点子的选取改为对点子的排序, 在排序之后对其进行选取提高了控制性. 传统的分裂点选取点子的方法是根据地形图制图要求确定阈值大小, 进而确定点到面的距离, 大于阈值的则保留, 小于阈值的删除. 而现今, 如何确定每个点子的特征重要性, 如何客观地对其排序, 到底保留哪些点可以回放出特征完全的等值线, 即为本文的研究内容.

(2) 在能够实现点子重要性排队后, 保留原始数据的多少比率下的数据才能生成符合要求的等值线, 也就是探究已知量比例尺与未知量数据保留比率之间的关系. 而在点子数据生成等值线过程时, 栅格的分辨率也会影响最终等值线绘制的综合程度. 所以探究生成栅格分辨率与等值线效果的关系也成为另一个内容. 　

2.3 前期数据准备

使用地面激光雷达扫描技术获取地形数据, 首先进行野外数据的采集, 考察测区环境并确定扫描仪、测站数和标靶的位置. 选定合适的测站数既要保证最终的数据能表达测区完整地形, 也要保证数据冗余程度低. 之后进行数据的初始处理, 要逐步进行点云匹配、消除噪声、坐标变换、图像分割. 当点云数据来自不同测站点时, 需要将其匹配到同一个坐标系, 通用的方法是通过GPS获取扫描站点的地理坐标, 并将点云数据转换到同一个地理坐标系中. 配准完成后, 手工对图像进行重采样, 消除噪声如电线杆、植被等, 最后可以对图像进行分割, 选取出感兴趣的区域.

本文所选用的ASCII格式的点云数据, 获取到的是使用单测站所测得的矿山数据, 因此其表现的是矿山的单侧面. 可将此文件在ArcMap中转化为Shapefile文件, 手动添加X、Y、Z字段, 并使用自动计算器计算每个点的三维坐标并赋值. 可将此类数据在ArcScene进行三维显示.

3 算法设计与实现 3.1 算法的设计

(1) 创建基面

遍历Shapefile文件中所有的点, 获取每个点的三维坐标, 进行依次比较后贮存X, Y, Z三方向的最大值与最小值. 记X_min, X_max分别为X方向上的最小值和最大值, Y方向与Z方向以此类推. 则O(X_min, Y_min, Z_min)、A(X_max, Y_min, Z_min)、B(X_min, Y_max, Z_min)分别构成OA和OB两个矢量线, 以此构成首基面OAB; C(X_max, Y_min, Z_max)、D(X_min, Y_min, Z_max)、A(X_max, Y_min, Z_min)分别构成CD和CA两个矢量线, 以此构成首基面CAD; E(X_max, Y_max, Z_min)、F(X_max, Y_max, Z_max)、A(X_max, Y_min, Z_min)分别构成EA和EF两个矢量线, 以此构成首基面EFA; OAB, CDA, EFA为三个正交的首基面, 如图2所示.

(2) 单首基面(一合一)3D点子重要性排序

单首基面对3D点子重要性排序的方法依据点到面的距离, 而三维点到面的距离与二维点到线的距离相对应, 因此利用二维曲线结构图更为直观的优点说明三维数据, 如图3所示.

二维下对点子进行重要性排队的步骤如下(流程见图4):

Setp 1. 无条件地将曲线的首末点(即点1和点2)选取(已被选取点标记为红色, 下同);

将连接试验曲线首末两点的直线作为首基线, 考察原曲线介于首末点之间的所有点, 可找到点线距为最大值的点(即点3), 故点3被选取; 如图4(a).

Setp 2. 介于已选点1和3之间、3与2之间, 分别考察原曲线介于局部首末点之间的所有点, 可找出最大值的点(即点4和点19标记为黑色), 将点4插入1和3之间, 将点19插入3和2之间, 如图4(b).

Setp 3. 在被插入的中间点点集中(此处指的是点4和点19), 由于点4的点线距大于点19的, 所以点4被选取(已被选取的点4被标成红色), 如图4(c).

Setp 4. 介于1与4之间、4与3之间, 分别考察原曲线介于局部首末点之间的所有点, 可找出点线距为最大值的点(即点5、15), 并插入各区间, 如图4(d).

Setp 5. 在被插入的中间点点集中(此处指的是点5、15、19等黑色标号所指), 由于点5的点线距为最大, 所以点5被选取(点5被标成红色), 如图4(e).

Setp 6. 如此循环往复, 其依照的口诀便是: ① 点子选上变红点; ② 红点之间插黑点; ③ 黑点之间要PK; ④ PK得胜变红点; ⑤ 全是红点就结束; 否则转Step 2不要停. 黑点逐个变成红点, 变成红点的先后顺序也就是其重要性顺序. 最终点子顺序如图4(f).

将上述排序方法转化为三维, 以OAB首基面为例, 则O点为原点分别与上述方法中的局部首末点相连构成一个矢量面, 此时计算点线距转化为计算点面距. 若以OAB为一个首基面, 设此面的一般方程式为 ${\rm{a}}X + bY + cZ + d = 0$ , 其中a、b、c、d均为已知量. 点到此基面的距离由公式(1)计算:

${{dis}} = \frac{{aX + bY + cZ + d}}{{\sqrt {{a^2} + {b^2} + {c^2}} }}$

(1)

其中, a、b、c、d为OAB基面的一般方程式的系数, X、Y、Z为文件中点的三维坐标.

(3) 三合一3D点子排序

按照上一步的方法, 获得3组分别以OAB、CDA、EFA为基面得到的顺序点队列I、II、III, 在这些队列中, 从左到右, 点的重要性依次降低, 对每个队列中的点依次编号. 则点子的序号越大, 说明此点的重要性越小. 以I为依据, 依次读取其中的点对应的序号, 并在II和III中找到相同的点, 获取其序号. 将同一点在此3个队列中的序号相加得到数值sum, 表示其重要性, 数值越小则表示其重要性越大. 按照值从小往大的顺序依次排列其对应的点, 将这些点按此顺序存入到队列中, 由此将体现3个方向特征重要性的点序列融合成一个点队列.

(4) 剁尾刀

以三合一过程获取的点队列为依据, 确定所要保留的点率, 由公式(2)得到最终保留点的数量.

${{n}} = {{r*count}}$

(2)

其中, count为源数据点的数量, r为保留比例, 取值范围为0~1, n为最终保留点的数量. 此时确定的数量即为要截取点子列表中的数据量, 从队列的左端开始往右获取剁尾刀的保留点界限, 直接得到左侧应保留的点集合.

3.2 算法实现

依托ESRI的ArcGIS平台, 使用C#语言, 在VisualStudio 2010环境下, 进行Add-In定制开发ArcMap扩展桌面应用程序功能. 将原始格式为ASCII点云数据经预处理后的得到shapefile格式的点云数据后创建X、Y、Z坐标字段. 将shapefile格式的点云数据作为源数据, 设定固定比例的特征保留率, 经上述算法处理后得到结果点云数据.

利用COM组件技术, 实现ArcObjects开放的特定接口, 实现3D Analyst Tools工具栏下的Topo To Raster和Contour功能, 绘制等值线, 再实现可视化接口使其直接在ArcMap中显示. 在实现点数据综合后运用二次开发技术实现等值线绘制的自动一体化过程.

4 算法设计与实现

测试数据区域为LiDAR系统单测站所采集的矿山点云数据, 其中共有29 677个点. 在点云数据经过本算法生成等值线的过程中, 点集合保留的比率和DEM生成的分辨率是算法的两个可控变量, 成为影响最终等值线的生成结果的主要因素. 点云数据在低精度DEM时生成等值线时, 海量数据的优势得不到体现, 绘制等值线时与传统方法绘制差异不显著, 在1 m至5 m精度范围内可以较好的反映地形信息^[12,13], 故选择在保留原始数据100%的基础上, 分别生成分辨率为1.0 m、1.5 m、2.5 m的DEM对比分析绘制等值线差异. 按梯度选择出对比分析保留10%、50%、100%的点数, 生成DEM分辨率为1.5 m, 绘制的等值线结果差异. 最终对比保留点率和分辨率这两个参数作为变量对回放等值线的影响.

5 结论与讨论 5.1 算法实现

由点子数据绘制等值线的过程中, 影响最终结果的主要因素为点子的数量和DEM的分辨率, 点云数据中的点子数量和DEM的分辨率与绘制等值线结果的精度成正比, 为对比这两个因素对等值线绘制的影响, 先固定一个变量再变化另一变量的不同情况进行对比分析得到其影响程度.

(1) 不同分辨率

导入源数据, 为保证点子数量不变, 故使用3D-DP算法保留100%的点集合分别转换成网格边长1 m、1.5 m、2.5 m的DEM数据, 进而绘制出等值线, 如图5、6、7所示.

通过实验可以看出, 利用重要性排序后的固定点子, 转换成边长变量为1.5 m的栅格DEM, 回放出的等高线比较符合工程测量与绘图的要求(见图6). 由前三幅图中可见, 根据综合后的点数据所绘制的等高线能够很好地表现所测范围的地形, 且大体趋势一致. 随着生成DEM的格网增大, 分辨率逐渐降低, 综合程度增大, 等值线趋于平滑, 且碎屑多边形逐渐减少, 提升主观视觉效果. 将这3种情况叠合在一起的结果如图8所示, 图幅的中间部位三者吻合度较高, 在图幅的上端和下端稍有差异, 但并不影响整体的地形表达效果.

(2) 不同比率保留点集合

输入源数据, 使用3D-DP算法, 根据设定的10%、50%、100%比率分别保留2968、14 839、29 677个点数据, 插值为分辨率为1.5 m的DEM, 进而绘制成等高距为5 m的等值线. 三类比率对应的等值线如图9、10、11所示, 将此三种数据叠加如图12所示.

通过实验可以看出, 利用保留率变量为100%的最重要的点子, 转换成固定边长为1.5 m的栅格DEM, 回放出的等高线比较符合工程测量与绘图的要求(见图10). 由图11等高线可以看出, 随着保留点子数量的减少, 即使是只保留至10%的点子, 依旧能够表现所测地区的基本地貌形态. 且由等高线对比分析时可以看出保留不同程度点子所绘制的等值线差异不大, 可说明通过3D-DP算法综合的点子能够保留下最能表达地形特征的点, 在一定限度的点数下, 可对表达同一地形的冗余点进行删除.

实验一和实验二分别降低分辨率和减少点子数量, 均加大了综合程度, 由实验可知减少点子数量相对于降低分辨率所绘制的等值线吻合度更高, 差异更小. 由此可以得出, 通过3D-DP算法保留特征点可以满足不同综合程度下的等值线绘制, 也因此验证3D-DP中点子的顺序排队和“剁尾刀”的方法易于通过给定点子保留百分比控制等高线的综合程度, 与之前国内外其他研究相比^[2–8], 本方法充分利用了点云数据的优势, 且体量更小, 易于调用, 具有更好的普适性, 在比例尺为1:500–1:20 000的范围内可以实现小范围内多尺度地图制图综合的自动化, 可输出不同比例尺下的等值线, 与已有平台兼容良好, 可嵌套在ArcGIS平台中, 结果可直接在ArcMap中显示, 对于多尺度的等高线输出更为方便、高效.

5.2 研究结论与展望

(1) 对于综合后的等高线, 与实测数据对比以后, 确定其回放的等高线可以满足对应比例尺下的地形图要求, 可以实现小范围内地图制图综合程度控制的自动化.

(2) 本文使用了3D D-P算法对数据进行综合, 并可以选择一定分辨率内的点子数据集合以减少运算时间, 提高等值线生成效率.但是, 若要面对大批量数据自动生成等值线时, 必须对数据进行分块处理, 除点位数据保留比率与DEM分辨率之外, 本算法还应对比不同类型的典型地形之间的等值线绘制效率与精度的差异进行比较, 来进一步提高算法的科学性和完整性, 以及引入并行运算技术, 来实现大范围高分辨率、多尺度等高线的地图制图自动化.

(3) 在仅有三维坐标的情况下生成等值线, 若要真实地描述地形, 在考虑地貌结构线的参与下, 以规则约束等值线的生成会使结果更贴合实际地貌.特别是处理山区中地貌与水系之间的协调套合, 使其真实体现地形.