CN104183017A - 基于地面三维激光点云实现地质体产状自动提取的方法 - Google Patents

Abstract

本发明公开了一种基于地面三维激光点云实现地质体产状自动提取的方法，利用三维点云数据空间位置和颜色信息，结合图像分割和拐点检测技术实现地质体出露面和未出露面产状自动提取。包括如下步骤：点云数据采集和预处理；获取点云数据的txt文本格式；采用最近邻法将点云数据重采样到二维图像；分割二维图像；评价图像分割精度；分别拟合地质体出露面和未出露面，并计算其产状；其中地质体未出露面在二维图像上表现为线状要素，且线条走向的大角度变化代表未出露面产状的变化，因此，计算未出露面产状较复杂，涉及拐点检测、赤平投影等技术。本发明在地质体产状提取方面，具有十分重要的实际应用价值，大量减少了危险且耗时耗力的实地测量工作。

Description

基于地面三维激光点云实现地质体产状自动提取的方法

技术领域

本发明涉及三维点云与影像分割技术和构造地质科学领域，特别是一种高精度、高密度、包含空间和颜色信息的三维激光点云实现地质体产状自动提取的方法。

背景技术

岩体结构面的优势产状和几何信息是评价岩体结构面发育的重要指标，特别是结构面中规模小、数据量大且随机分布的节理和裂隙，其优势产状是评价岩体稳定性的重要参数。野外实地测量不但耗时耗力，而且在一些危险地区，人员无法到达。三维激光扫描是一种远距离非接触的全自动立体扫描技术，又被称为“实景复制技术”。获取的点云数据精度高且数据量大，涵盖目标的三维空间位置和真彩色信息，目前已在数字地质露头、地质工程等领域得到应用。Gigliand Casagli(2011)和刘昌军等(2011)都提出了基于三维点云数据提取地质体产状的方法，他们的基本思想是利用结构面上多个点的空间坐标来拟合平面，并根据平面方程计算地质体产状。

然而目前研究仅限于出露地质体产状的人机交互测量，无法实现地质体未出露面产状的提取。野外，地质体经历了剥蚀和沉积等过程，其出露面不能真实代表其。例如岩层面产状的提取，表层的岩层(出露)因为外力作用已经变得破碎不平整，因此要准确计算岩层的产状，我们需要测量两地层之间的接触面，而该接触面是未出露的；又如作为火山作用与深部地幔性质、演化及壳幔相互作用的纽带的岩墙，工作者需要测量它与围岩接触面来计算岩墙产状，其产状有助于计算岩浆超压，研究岩浆动力学及侵位机制。因此在实际测量工作中，获取地质体未出露面产状具有更重要的意义。此外，目前对三维点云数据的研究大多未涉及其颜色信息RGB。鉴于以上情况，充分利用点云数据的空间位置和颜色信息，在出露地质体产状提取的基础上，结合影像分割技术，ENVIEX(http://www.esrichina.com.cn/2010/1029/123.html)是目前普遍运用的分割软件之一，边界拐点检测技术，类似的工作包括利用8k(k＝2，3，4…)邻域链码和道格拉斯-普克(Douglas-Peuker)抽稀算法(http://www.baike.com/wiki/抽稀)检测边界拐点，进一步实现地质体出露面和未出露面产状的自动提取，具有重要的实用价值。

Riegl VZ100(www.riegl.com)是一种目前应用广泛的三维激光扫描仪，其有效扫描距离最大可达1400m，距工作站点100m内点间距为5mm，水平扫描范围为0°～360°，垂直扫描范围为-40°～60°。通过不同站位的多角度扫描，可以获取覆盖实验区的全部点云数据。其配套软件RISCAN PRO可以对三维点云数据进行赋色、配准、去噪等预处理，从而为后续工作提供较好的实验数据。

发明内容

本发明的目的在于针对现有技术的不足，利用点云数据的空间位置和颜色信息，结合三维点云数据处理、影像分割和边界拐点检测技术提供一种地质体出露面和未出露面产状自动提取的方法。

本发明的目的是通过以下技术方案来实现的：一种基于地面三维激光点云实现地质体产状自动提取的方法，包括如下步骤：

(1)通过三维激光扫描仪采集三维激光点云数据；

(2)对步骤1采集的三维激光点云数据进行预处理，包括赋色、配准、去噪操作；

(3)获取点云数据的txt文本格式，该文本的每一行包含一个点云数据的空间坐标XYZ和颜色信息RGB；根据其空间坐标，采用最近邻法将点云数据重采样到二维图像，即将三维点投影到二维图像上，该图像包含6个图层数据，即点云数据的X、Y、Z、R、G、B；

(4)基于点云数据的颜色信息，对步骤3中获取的二维图像进行分割；采用监督评价方法评价图像分割精度，不断调整分割参数直到分割精度达到阈值；

(5)地质体出露面和未出露面在二维图像上分别表现为面状要素和线状要素，即分割对象和分割边界线，当分割精度达到设定阈值时，对于出露地质体，利用组成该对象的三维点，采用最小二乘法拟合结构面并计算其产状；对于未出露地质体，只能扫描到其边界线上的点云数据，且边界线走向的大角度变化代表其产状的变化，因此其产状的计算较出露面复杂，具体包括以下步骤：

(5.1)利用8k邻域链码和道格拉斯-普克抽稀算法检测边界拐点，其中，k为大于等于2的整数；

(5.2)以拐点为分隔点，根据边界线上点云数据的空间坐标XYZ，利用最小二乘法分段拟合地质体未出露结构面并计算其产状；

(5.3)将步骤5.2中计算得到的地质体产状进行赤平投影，对产状相近的地质体进行归类，并赋予同一标志；合并具有相同标志的相邻地质体，即相邻分界线，对于孤立的地质体，如果其相邻地质体属于同一类，则将该孤立地质体与其相邻地质体归为一类，然后分别根据每类边界点拟合平面并计算产状；重复步骤5.3，直到相邻地质体产状均不相同。

进一步地，所述步骤4中所述监督评价方法通过比较分割结果与参考对象之间的相似度或差异度来评价分割精度，具体为：采用基于面积和基于位置的组合度量来评价图像分割精度。

本发明的有益效果是：通过结合利用点云数据的空间位置信息和颜色信息，从而实现“三维→二维→三维”的技术路线，将研究较成熟的影像分割技术用于三维点云数据分割，并通过精度评价来提高其分割效果，减少对后续操作的影响。本发明基于点云数据的颜色信息获取了地质体在二维图像上的面状和线状表现，并根据三维点的空间坐标拟合地质体结构面，从而获取其产状信息。利用该发明，可以实现出露和未出露地质体产状的自动提取，从而减少了大量野外实地测量工作，不但节约了时间和资金，还减少了工作的危险性。

附图说明

图1是本发明基于地面三维激光点云实现地质体产状信息自动提取方法的流程图；

图2是本发明不同站点对同一目标的扫描示意图；

图3是本发明基于面积度量评价分割精度获取相交面积的流程图；

图4中(a)是具体实施例中对点云数据进行重采样得到的二维图像，其中白色线条是手动分割岩墙的边界，作为图像分割精度评价中的参考对象；(b)是ENVI EX软件分割岩墙得到的边界二值图；(c)是8k邻域链码检测的拐点图(黑色点为拐点)，其中k为2；(d)是道格拉斯-普克算法检测的拐点图(黑色点为拐点)，其中抽稀因子为3.5。

具体实施方式

如图1所示，本发明基于地面三维激光点云实现地质体产状自动提取方法，包括以下步骤：

1.通过三维激光扫描仪采集三维激光点云数据；

三维激光扫描仪Riegl VZ100在扫描过程中由于受激光扫描仪视场角的限制，单一站点的扫描范围往往不能覆盖整个研究区域，因此需要进行多站点扫描，且不同站点之间必须保证存在重叠区，以便实现不同站点间点云数据的配准；具体可参考文献《Simon J,Buckley J A,et al,2008.Terrestrial laserscanning in geology:data acquisition,processing and accuracyconsiderations.Journal of Geological Society,London,165:625-638》；

2.对步骤1采集的三维激光点云数据进行预处理，包括赋色、配准、去噪操作；

在三维激光扫描仪Riegl VZ100的配套软件RISCAN PRO中完成三维点云数据的预处理；Riegl VZ1000在扫描的同时会拍照记录点云数据的颜色信息，在RISCAN PRO中通过Color from images操作，自动完成相应照片对点云数据的赋色，为其增加RGB信息；配准是将不同站点坐标系下的点云数据转换到一个共同的基准坐标系下，主要包括步骤：1)确定同名点对；2)解算坐标转换数据函数，通过RISCAN PRO软件中的Coarse registration操作完成不同站点点云数据的配准，该步骤需要用户在两个站点中选取至少4组同名点对进行配准，RISCAN PRO软件会计算配准误差，根据该误差确定是否需要重新选点；去噪即是去掉那些受测量仪器或其他环境因素影响产生的噪声点，该步骤可以通过定义RISCAN PRO软件中的Filter操作完成，软件用户说明书中具体介绍了定义Filter操作的各参数；

3.在RISCAN PRO软件中导出点云数据的txt文本格式，该文本的每一行包含一个点云数据的空间坐标XYZ和颜色信息RGB；根据其空间坐标，采用最近邻法将点云数据重采样到二维图像，即将三维点投影到二维图像上，该图像包含6个图层数据，即点云数据的X、Y、Z、R、G、B，具体的：将三维点云数据投影到xy平面上，并对该区域进行图像重采样，最近邻法可表达为：

$\mathrm{Dist}=\sqrt{{(X_j-x_i+Q_1)}^2+{(Y_j-y_i+Q_2)}^2}$

其中j和i是正整数，X_j和Y_j是点云数据的X和Y坐标，x_i和y_i是图像点坐标，Q₁和Q₂是常数，分别表示点云数据坐标中最小的X值和Y值，Dist是点云数据与重采样二维图像中某一确定像素点的距离，对于每一个图像点，即一个确定的i，我们需要遍历所有点云数据，即j＝1,2，…n(n是点云数据的个数)，计算每个点云数据与该像素点的距离，然后选取距离该点位置最近的点云数据(即Dist_min)的空间坐标XYZ和颜色信息RGB作为其波段属性值；

4.基于点云数据的颜色信息，利用ENVI EX软件中的Feature Extraction操作对步骤3中获取的二维图像进行分割；采用监督评价方法评价图像分割精度，该方法是通过比较分割结果与参考对象之间的相似度或差异度来评价分割精度的，在进行该操作前需要准备好参考对象和分割结果，其中参考对象可以通过目视解译手动分割获得；通过调整分割参数(ENVI EX软件中分割参数有Scale level和merge level)以达到分割精度阈值，该阈值由用户自己设定；利用基于面积和基于重心的组合度量来评价分割精度，组合度量避免了单一度量对分割不足或过度分割的偏好，主要度量公式是：

基于面积的度量：

${\mathrm{OverSegmentation}}_{\mathrm{ij}}=1-\frac{\mathrm{area}(X_i\∩Y_j)}{\mathrm{area}\left(X_i\right)}$

$\mathrm{UnderSegmentati}{\mathrm{on}}_{\mathrm{ij}}=1-\frac{\mathrm{area}(X_i\∩Y_j)}{\mathrm{area}\left(Y_i\right)}$

基于位置的度量：

D_ij＝dist(centroid(X_i),centroid(Y_j))

其中X是参考对象，Y是分割对象，OverSegmentation和UnderSegmentation分别表示过度分割和分割不足，D_ij表示参考对象中心和分割对象中心的距离；area(X_i)和area(Y_j)分别表示参考对象和分割对象的面积，area(X_i∩Y_j)是参考对象与分割对象的相交面积，centroid(X_i)和centroid(Y_j)分别表示参考对象和分割对象的中心；

OverSegmentation、UnderSegmentation和D的最优值均为0，当三者同时取最优值0时即为最优分割，即分割对象与参考对象完全重合；具体可参考文献《Clinton N,Holt A,Scarborough J,et al,2010.Accuracy assessment measures forobject-based image segmentation goodness.Photogrammetric Engineering andremote sensing,76(3),289-299》；

5.地质体出露面和未出露面在二维图像上分别表现为面状要素和线状要素，即分割对象和分割边界线，当分割精度达到设定阈值时，对于出露地质体，利用组成该对象的三维点，采用最小二乘法拟合结构面并计算其产状；

根据空间几何原理，结构面的方程可以写为：

ax+by+cz+d＝0

其中a，b，c，d为结构面空间常数，x，y，z是三维点空间坐标变量。结构面法线与X轴、Y轴、Z轴的夹角余弦(l,m,n)按以下公式计算：

$l=\frac{a}{\sqrt{a^2+b^2+c^2}}$

$m=\frac{b}{\sqrt{a^2+b^2+c^2}}$

$n=\frac{c}{\sqrt{a^2+b^2+c^2}}$

根据结构面的法向向量计算其倾向α和倾角β(倾向是倾斜平面上走向线的垂线在水平面上的投影所指的沿平面向下倾斜的方位，其中走向线是倾斜平面与水平面的交线，倾角是倾斜平面与水平面所成的角度)；

$\mathrm{\α}=\arctan \left(\frac{m}{l}\right)+Q$

$\mathrm{\β}=\arctan \left(\frac{n}{\sqrt{l^2+m^2}}\right)$

其中Q是常数，当l>0，m>0时，Q＝0°；当l>0，m<0时，Q＝360°；当l<0时，Q＝180°；具体可参考文献《Gigli G,Casagli N,2011.Semi-automaticextraction of rock mass structural data from high resolution LIDAR pointclouds.International Journal of Rock Mechanics and Mining Sciences,48(2),187-198》

对于未出露地质体，只能扫描到其边界线上的点云数据，因此利用其边界线上的三维点(只要有3个不在同一直线上的空间点便能拟合平面，三维激光仪的多角度扫描和野外露头的不平整可以保证边界线上的点不在同一直线上)来拟合结构面，且边界线走向的大角度变化代表其产状的变化，因此其产状的计算较出露地质体复杂，具体包括以下步骤：

(5.1)利用8k(k＝2，3，4…)邻域链码和道格拉斯-普克(Douglas-Peuker)抽稀算法检测边界拐点；8k邻域链码和道格拉斯-普克抽稀算法分别适用于栅格图像和矢量图像，针对不同的目标物，需要通过多次实验来确定最优k值或抽稀因子。

8k(k＝2，3，4…)邻域链码适用于具有闭合边界的二值栅格图，该算法首先利用8k邻域链码求得二值图像中物体边界个点的曲率，然后通过寻找曲率的局部峰值检测出拐点；具体可参考文献《尚振宏,刘明亚.二值图像中拐点的实时检测算法.中国图像图形学报,2005,10(3)》；

道格拉斯-普克(Douglas-Peuker)抽稀算法适用于矢量图，其基本思路为：

(a)对曲线的首末点虚连一条直线，求曲线上所有点与该直线的距离，并找出最大距离值dmax，用dmax与事先给定的阈值D相比：

(b)若dmax<D，则将这条曲线上的中间点全部舍去；若dmax≥D，保留dmax对应的坐标点，并以该点为界，把曲线分为两部分，对这两部分重复步骤1和2，直到所有dmax均<D，即完成对曲线的抽稀；

(5.2)以拐点为分隔点，根据边界线上点云数据的空间坐标XYZ，利用最小二乘法分段拟合地质体未出露结构面，然后按照与计算出露面产状相同的方法来计算地质体未出露面的产状；

(5.3)将步骤5.2中计算得到的地质体未出露面产状在软件DIPS中进行赤平投影，对产状相近的地质体进行归类，并赋予同一标志，合并具有相同标志的相邻地质体；对于孤立的地质体，如果其相邻地质体属于同一类，则将该孤立地质体与其相邻地质体归为一类，然后根据每一类地质体边界上的三维点重新拟合结构面并计算产状；并重复步骤(5.3)，直到相邻地质体产状均不同。该步骤通过合并产状相近的相邻地质体，可以增加拟合未出露地质体的三维点，从而避免因范围过小而拟合出错误的地质体未出露面。

下面以岩墙产状提取为例进行说明。

岩墙作为一种厚度比较稳定的板状侵入体，岩墙面几乎不出露，出露地表的仅是岩墙的厚度，利用出露面上多个点的空间坐标构建平面方程来获取结构面产状的方法不适用于岩墙。因此，利用本发明提出的基于地面三维激光点云实现地质体产状自动提取的方法来提取岩墙产状信息。

图1是基于地面三维激光点云实现地质体产状自动提取方法的流程图。使用者按步骤1采集原始点云数据；按步骤2对点云数据进行预处理，包括赋色、配准、过滤操作；按步骤3导出点云数据的txt文本格式，其中该文本的每一行包含一个点云数据的空间坐标XYZ和颜色信息RGB，并根据最近邻法将其重采样为二维图像，实现三维到二维的转换；按步骤4分割研究对象岩墙，利用监督评价方法评价分割精度，图3是基于面积度量评价分割精度获取相交面积的流程图，通过不断调整分割参数使分割精度达到设定阈值(80％)；鉴于岩墙属于未出露地质体，在二维图像上表示为线状要素，因此按步骤5.1检测岩墙边界拐点；按步骤5.2分段计算岩墙产状，实现二维到三维的转换；按步骤5.3不断对分段岩墙进行归类，并计算产状。

参阅图4所示，图4a是根据最近邻法将点云数据重采样得到的二维图像，图4a中的白色线条是手动分割岩墙的边界，作为图像分割精度评价中的参考对象。图4b是利用ENVI EX软件对岩墙进行分割后获取的边界二值图，以图4a中的分割岩墙为参考对象，采用基于面积和位置的组合度量评价图4b岩墙的分割精度，总体精度达80％以上，但比较图4a和图4b可知：ENVI EX软件对岩墙的分割忽视了左侧岩墙分支。利用8k(k＝2，3，4…)邻域链码和道格拉斯-普克(Douglas-Peuker)抽稀算法分别检测岩墙边界拐点(黑色点)得图4c(k取2)和图4d(抽稀因子取3.5)，其中图4c和图4d分别包含35和49个拐点。对于该实验数据，D-P算法抽稀检测出的拐点更多，对岩墙的分段更细。最后通过分割线上的空间点分段拟合岩墙未出露面，并不断合并产状相近的相邻面，最终得到两个大平面，其产状(倾向、倾角)分别是N55.4°W，76.4°(图4b中近竖直方向的岩墙①)和N25.8°E，60.6°(图4b中右侧岩墙分支②)。

Claims (2)

1.一种基于地面三维激光点云实现地质体产状自动提取的方法，其特征在于，该方法包括如下步骤：

（1）通过三维激光扫描仪采集三维激光点云数据；

（2）对步骤1采集的三维激光点云数据进行预处理，包括赋色、配准、去噪操作；

（3）获取点云数据的txt文本格式，该文本的每一行包含一个点云数据的空间坐标XYZ和颜色信息RGB；根据其空间坐标，采用最近邻法将点云数据重采样到二维图像，即将三维点投影到二维图像上，该图像包含6个图层数据，即点云数据的X、Y、Z、R、G、B；

（4）基于点云数据的颜色信息，对步骤3中获取的二维图像进行分割；采用监督评价方法评价图像分割精度，不断调整分割参数直到分割精度达到阈值；

（5）地质体出露面和未出露面在二维图像上分别表现为面状要素和线状要素，即分割对象和分割边界线，当分割精度达到设定阈值时，对于出露地质体，利用组成该对象的三维点，采用最小二乘法拟合结构面并计算其产状；对于未出露地质体，只能扫描到其边界线上的点云数据，且边界线走向的大角度变化代表其产状的变化，因此其产状的计算较出露面复杂，具体包括以下步骤：

（5.1）利用8k邻域链码和道格拉斯-普克抽稀算法检测边界拐点，其中，k为大于等于2的整数；

（5.2）以拐点为分隔点，根据边界线上点云数据的空间坐标XYZ，利用最小二乘法分段拟合地质体未出露结构面并计算其产状；

（5.3）将步骤5.2中计算得到的地质体产状进行赤平投影，对产状相近的地质体进行归类，并赋予同一标志；合并具有相同标志的相邻地质体，即相邻分界线，对于孤立的地质体，如果其相邻地质体属于同一类，则将该孤立地质体与其相邻地质体归为一类，然后分别根据每类边界点拟合平面并计算产状；重复步骤5.3，直到相邻地质体产状均不相同。

2.根据权利要求1所述的一种基于地面三维激光点云实现地质体产状自动提取的方法，其特征在于，所述步骤4中所述监督评价方法通过比较分割结果与参考对象之间的相似度或差异度来评价分割精度，具体为：采用基于面积和基于位置的组合度量来评价图像分割精度。