三维岩体点云数据地质信息提取方法

王喆，刘鹏飞，王凯

(北京市勘察设计研究院有限公司，北京 100038)

1 引 言

三维激光扫描技术(3D Laser Scanning Technology)作为一项高新测绘技术兴起于20世纪90年代。相比于传统测量手段，该技术采用非接触的方式连续、自动采集目标对象的表面属性点信息，数据形式为三维点云[1]。目前，三维激光扫描技术已在多个领域，例如城市科学规划、设计、管理中普遍应用。

传统的地质研究中，通常使用实地勘测的方法获取岩体中结构面的信息，由于该方法依靠技术人员以罗盘、皮尺等工具进行实地测量，故受现场局限性较大。当受到场地限制(如高陡边坡)，通常难以获取全面的断面信息。由于此类基础断面数据不齐全，会对进一步的岩体裂隙分布规律判断产生干扰，最终影响岩体稳定性分析评价工作。并且此方法存在工作效率低、精度低、成本高、外业工作量大和安全隐患大等问题[2]。

利用三维激光扫描技术可实现无接触地直接采集目标岩体的点云数据，数据经过处理后得到采集目标的三维点云模型。由于该模型真实地逼近于现场情况，技术人员可在后期处理软件中对危岩体几何参数进行量测和计算，得到岩体结构面产状，从而完成一系列地质分析。故此方法在地质勘察工作中，尤其是大高差山区危岩地区，可基本代替传统人工勘测方法，对地质体信息提取工作具有极大的价值[3]。

2 PCL简介

PCL(Point Cloud Library)是一个大型跨平台(Windows、Linux、MacOSX、Android等)C++编程库，作为一个开源项目，其最早是由斯坦福大学的Radu博士等人维护和开发[4]。PCL是基于Boost、FLANN、Eigen、VTK、CUDA、OpenNI等第三方库集成的模块化C++模板库。基于PCL可建立高效的点云数据结构，依此结构处理点云大数据，可实现相关的通用算法，PCL均将其进行了封装并提供调用接口，其中包括点云的获取、分割、检索、配准、滤波、可视化、特征提取等。图1为PCL的架构图[5]。

图1 PCL架构图

2.1 1PCD文件

PCL作为针对点云数据开发的专业库，为了高效地处理点云大数据，此库采用了C++作为底层编程语言，并以此发布了PCL独有的点云文件数据格式：PCD格式。随着三维激光扫描技术在各行业领域的广泛应用，技术设备不断迭代更新，各大设备和软件厂商纷纷推出了各自专属的点云数据格式，其中较为主流的格式包括LAS、OBJ、PLY、X3D、STL等。但上述数据格式大都是根据本厂家特定的应用场景，在旧式传感器设备和算法基础上创建的。相比以上数据格式，PCD格式的主要优势如下[6]：

(1)采用mmap、munmap作为基础二进制数据类型，使其拥有高效的数据下载和存储能力。

(2)可存储多类C++基本数据类型，如int、float、double、char、short等，并针对无效的数据点设置特定的数据类型：NAN类型，以便于在PCL内部存储，从而提高PCL的点云处理能力。

(3)通过创建PCD文件格式，可以发挥C++强大的数据处理能力的优势，并且能最大程度上适应PCL与PCL的内部函数无缝衔接，从而让基于PCL开发的应用程序获得最佳的数据处理性能。

2.2 PCD文件操作

基于C++强大的底层支撑，PCL内部提供了丰富的点云文件操作模块，以模块提供的一系列接口开发程序，可以很轻松地实现针对点云的相关操作。以 I/O模块为例，该模块主要用来支持PCD文件的输入输出操作，共包含21个子类。例如PCD文件读写的接口就可以使用File-Reader类和File-Writer类分别定义。二者的继承关系如图2所示，其功能实现的关键代码为：

pcl：：PointCloud：：Ptrcloud(newpcl：：PointCloud)；

pcl：：io：：loadPCDFile(“test_pcd.pcd”，*cloud)；

pcl：：io：：savePCDFileASCII(“test_pcd.pcd”，cloud)；

图2 类File Reader和File-Writer的继承关系

3 平面分割算法

利用三维激光扫描获取的点云数据，具有连续、实时、数据量大，精度高等特点，经过矫正、降噪等预处理操作后，点云数据仍是一个庞大的完整岩体[7]。故欲进行地质体产状要素的计算与统计，需先将每一个结构面单独提出。本文采用RANSAC算法作为结构面提取的算法，RANSAC算法最早由Fischler和Bolles于1981年提出，该算法以一组包含异常数据的样本数据集为依据，计算出数据的模型参数，从而得到有效样本[8]。其基本思想是：每一个平面数据由该平面的“局内点”组成，其数据的分布可满足于特定的模型解释参数；刨除“局内点”，数据源中的其他点再分为“局外点”与“噪声点”，其中“局外点”是不能适应“局内点”模型的数据，剩余数据即为“噪声点”。该算法结果的合理性是有一定概率限制的，因此提高迭代次数才能将结果合理性的概率提高，即要保证在一定置信度下基本子集最小抽样数N与至少取得一个良好抽样子集的概率P满足如下关系[9]：

P=1-(1-εk)n

(1)

式中，k为计算模型参数需要的最小数据量；ε为局内点与数据集点数的比值；P一般取值为0.9～0.99。对上式两边取对数可得：

(2)

算法迭代的核心逻辑为：

(1)首先随机选择一个面为最佳拟合平面，比较当前迭代的点数与该最佳平面点数，如当前迭代点数更多，将当前点设置为最佳拟合平面；

(2)如当前点数量与最佳拟合平面点数量相同，但当前点满足更小阈值时，同样将当前点设置为最佳拟合平面。

总体上，RANSAC对噪声点有一定的抑制作用，并能较好地分割多面物体，是一种有效的稳健估计算法[9]。

具体算法流程如下(如图3所示)：

图3 Ransac算法流程

①根据式(2)计算最小循环次数N；

②计算拟合平面的标准差，如果其大于阈值δ0，重新确定平面的初值；否则继续下一步；

③计算面内点的数量，并与阈值点数比较，如大于阈值，则计算平面的标准差；否则返回上一步；

④重复步骤②、步骤③，直到根据判断准则得出最佳平面；

⑤重复步骤①～步骤④，直到模型点的个数小于给定的阈值Nmin为止，得出最终的分割结果。

4 地质产状算法

本文采用空间解析几何和向量代数法计算结构面产状，其原理为，过不在同一直线上三点确定的平面法向量的Z轴方位角为γ，其角度范围在0°～360°取值，计算时不考虑γ＇的方向性，规定其取值范围为0°～90°，由图4可以证明：岩层的倾角δ可通过法向量与Z轴的夹角γ＇获得(角度相等)，并且岩层倾向可以由法向量在水平面(xoy)上的投影向量的方向算出[10]。

图4 岩层产状计算示意

4.1 面法向量算法原理

在获取每个独立结构面点云数据后，本文采用最小二乘法对点云数据进行平面拟合[11]，本文假设误差只存在于Z轴方向上，从而计算该平面的平面法向量，设点云拟合的平面方程为：

Z=aX+bY+c

(3)

当取Z为观测值，n为观测点数目时，可推算出其观测值方程为：

Z+V=aX+bY+c

(4)

将上式改写为误差方程：

(5)

其中：

(6)

(7)

4.2 面产状算法原理

(8)

规定δ的取值范围为0°～90°，则：

(9)

(10)

由式(10)可简化得：

(11)

根据下表对X和Y值的正负进行修正即可得岩层倾向θ，详见表1[12]。

倾向基本角θ修正 表1

5 实例分析

本方法中涉及两类算法，一是基于Ransac的点云结构面分割算法，二是结构面产状算法，后者在计算过程中皆使用固定常量，利用计算机编程逐步计算即可。而前者Ransac算法中存在人为设定参数，直接影响到计算结果。故本章将依托龙庆峡景区边坡防护工程获取边坡点云数据，以此算法进行测试验证。龙庆峡距北京城区 85 km，主要为高山峡谷地貌，分布岩壁大都近直立，高约 200 m左右，节理发育明显，适宜进行成果对比。

RANSAC的点云分割效果主要由两阈值参数决定，分别为点到平面的距离阈值δ0和聚类面最小点云数量Nmin[13]，为取得最佳的分割结果，本文通过控制变量法试验求取δ0与Nmin的最佳参数值。此试验以龙庆峡边坡峭壁点云数据为样本，选择 15 m×20 m的典型岩石出露面两处，将距离阈值δ0设为 5 cm、7 cm、10 cm、12 cm、15 cm，将迭代条件中剩余点云数量Nmin设为全部点云的20%和固定数值100个两种。比较分析分割效果得到统计图如图5、图6所示。

图5 点云分割运行耗时统计

图6 点云分割节理面个数统计

软件耗时方面，当距离阈值δ0越大，运行时间越短。精度方面，在δ0值不超过 12 cm时，大部分点云可以被识别，当δ0值达到 10 cm和Nmin设置为100个点云，点云分割效果达到最佳，如图7所示。

图7 点云分割结果

当结构面分割完成，各结构面点云文件可依次输入产状计算算法，并最终输出其计算结果。为进一步评价此结果准确性，本文将此次结构面产状计算结果分组统计，并与人工现场地质调查结果进行比对，对比结果如表2所示：

计算结果比对 表2

根据统计结果及现场调查成果，区内发育的三组结构面中第一组为岩层的层面，其他两组为近垂直相交的陡倾节理，控制边坡岩体结构发育特征。统计结果与地质调查数据基本一致。

6 结 语

本文提出的基于PCL的地质信息提取方法，通过三维激光扫描仪所获取的点云数据，经过RANSAC算法、向量代数法，最终得到了岩体的相关地质参数，此技术的出现为后续地质调查与危岩体的治理和设计提供了可靠的资料。但本算法设置的阈值较多，而地质岩体的产状信息多与其特定的地理环境相关，故此算法需要在更多地区，更多岩性条件下进行测试，从而形成更加完善的阈值设定机制，并为下一步的三维地质的分析与建模奠定坚实的基础。