融合激光点云和数字影像的岩体结构面产状测量方法

摘要

一种融合激光点云和数字影像的岩体结构面产状全自动数字化测量方法。本发明运用摄影测量技术，采用混合全局和局部阈值法提取出迹线的轮廓线，通过骨架提取、迹线标记提取出节理迹线，并结合其他图像处理方法，实现了岩石节理迹线产状的信息获取；最后根据迹线与结构面的距离判断迹线与结构面的位置关系，并将迹线和结构面进行归并、分组，以各组产状表征岩体结构面产状。本发明融合了激光点云和数字影像的测量方法，充分发挥了这两种测量方法快速、高效、非接触、无视地形的性质，同时又提高了测量的精确度，并具备自动化程度高、普适性强、表述全面等优点。

说明书

技术领域

本发明涉及一种融合激光点云和数字影像的岩体结构产状全自动数字化测量方 法，属于工程地质勘察技术领域。

背景技术

岩体结构面产状测量方法归纳起来主要有四类：

①、测线法、窗口统计法即通过皮尺和罗盘人工现场逐一测量结构面产状信息。 该方法在面对实际工程中岩体结构面分布广、数量多、随机性强的特点，低效、费 力、耗时使之难以满足工程需求，而且遇到人迹罕至陡崖险坡受现场条件的影响， 测量方法便难以实施，从而导致测量数据的全面性与代表性也受到限制。

②、钻孔定向取芯法(Rosengren,1968)、孔内照相法(Eoek&Pentz,1968)、 数字全景钻孔摄像法(王川婴,2001)即通过钻孔取芯后对岩芯中所含结构面产状的 量测或通过孔内摄像对孔壁内侧结构面纹理进行数字图像量测来获取钻孔深度范围 内相交的结构面产状信息。该方法可以获得深部岩体结构面产状信息，但是孔径范 围测量尺寸使得该方法对于长大结构面的测量精度较低，同时钻孔破岩与孔位垂直 度也会影响到测量精度。

③、摄影测量法(Ross-Brown&At-kinson,1972；Hagan,1978)即基于数字图像与 摄影测量的基本原理，应用计算机三维成像技术、影像匹配、图像插值、模式识别 等多学科理论与方法相融合，可以获取地物表面的相对坐标空间几何信息，通过计 算这些几何信息可以得到暴露地表岩体结构面产状信息。该方法采用非接触的测量 手段，可以不受现场地形条件的制约，但是该方法精度受摄影图像质量与双目摄像 间距控制，同时仍然需要一定的人工干预，单条岩体结构面的工作效率比精测线法 有所提高，而面对实际工程中众多结构面大量人工干预仍然是无法接受的。

④、地面三维激光扫描法(TerrestrialLaserScanning，简称TLS)，即通过对地 物目标的逐点激光扫描测距可快速获得描述地物表面相对几何坐标信息的高精度点 云数据(激光数据)，通过一定人工干预与数据截取可以计算出那些暴露在岩壁表 面上的岩体结构面产状。该方法取消了摄影测量法中相对繁琐的摄影测量算法，提 高了测量速度与精度，但是仍然需要人工干预在海量点云数据中手动截取描述暴露 岩体结构面产状的点云数据，而且该方法对于未暴露岩体结构面或以迹线形式显现 的岩体结构面产状是无法测量。

综上所述，面对实际工程中岩体结构面数量多、分布广和随机性强等特点，目 前技术瓶颈在于如何兼顾非接触式(取消场地制约)、全自动(取消人工干预)、 测量速度快(算法简洁)、普适性(适合多种赋存环境岩体结构面)。

发明内容

技术问题：本发明目的在于借助地面三维激光扫描技术的快速精准空间测量手 段和计算机视觉图像识别技术的目标信息自动化筛取，将含有精确几何信息激光点 云数据与目标信息自动化筛取后的数字图像数据进行数据融合与几何配准，提出一 种基于TLS三维激光扫描仪的全自动、非接触式岩体结构产状数字化测量方法，即 融合激光点云和数字影像的岩体结构面产状测量方法。

技术方案：主要针对野外岩体结构面的数量多、分布广和随机性强等特征，传 统岩体结构面产状测量方法工作效率低下，甚至遇到人迹罕至陡崖险坡根本无法实 施。本发明借助于地面激光激光的精准空间测量技术和数字影像的计算机视觉识别 技术，提出了一种全自动、非接触式的岩体结构面产状数字化测量方法。

融合激光点云和数字影像的岩体结构面产状全自动数字化测量方法，其测量内 容包括了岩体结构面产状的点云数据提取、岩体节理迹线产状的图像数据与点云数 据相融合的提取，其表述方式采用迹线与结构面匹配分组、共同表述的方法。

在点云数据的提取中，先用滤波法进行去噪，再对点云数据进行三角剖分从而 建立初步的岩体三维模型，之后充分发挥点云数据精确度高的特点，通过对三角面、 三角面法向量的计算构造出特征矩阵并采用基于高密度联通区域的DBSCAN聚类 分析方法对结构面进行初步分组，随后借由三角面的连通性判别进行结构面的二次 分组，最后计算出结构面的产状；岩石节理迹线需从图像数据中提取出来，先对数 字影像进行灰度化处理，再通过图像混合阈值分割法提取出岩石的轮廓线，经过像 素剥离和去除断肢提取出迹线的骨架信息，最后将迹线标记出来并进行产状计算； 结构面和迹线的产状测量完成后，将迹线与结构面合理分组，通过计算出迹线和结 构面之间的位置关系并进行判别归并，用归并后的最终结果表征岩体结构面的产状 信息。

由于通过点云数据提取迹线产状是较为困难的，因此本发明融合图像数据和点 云数据进行迹线产状的提取。岩石节理迹线需从图像数据中提取出来，本发明先对 数字影像进行灰度化处理，再通过图像混合阈值分割法提取出岩石的轮廓线，经过 像素剥离和去除断肢提取出迹线的骨架信息，最后将迹线标记出来并进行产状计算。

结构面和迹线的产状测量完成后，为更全面、准确地表述岩体结构面产状，本 发明将迹线与结构面合理分组，通过计算出迹线和结构面之间的位置关系并进行判 别归并，用归并后的最终结果表征岩体结构面的产状信息。

融合激光点云和数字影像的岩体结构面产状全自动数字化测量过程主要包括如 下步骤：

第一步：采集现场数据；

第二步：将点云数据和影像数据进行匹配、融合；

第三步：对点云数据进行三角剖分，计算出剖分后三角形法相量的参数，对法 向量进行聚类分析；

第四步：将聚类分析后的结果进行合并，对合并后的数据进行产状计算；

第五步：对影像数据进行矫正、灰度化，然后提取出节理的轮廓线和骨架并在 点云中标记出骨架；

第六步：计算迹线对应点云任意三点构成三角形的内角，对迹线进行产状判别 并进行产状计算；

第七步：近似拟合提取出的结构面并计算结构面与各迹线的距离，然后按照结 构面和迹线的位置关系将结构面和迹线进行分组；

第八步：测量工作完成。

有益效果：本发明主要针对野外岩体结构面的数量多、分布广和随机性强等特 征，传统岩体结构面产状测量方法工作效率低下，甚至遇到人迹罕至陡崖险坡根本 无法实施。本发明借助于地面激光激光的精准空间测量技术和数字影像的计算机视 觉识别技术，提出了一种全自动、非接触式的岩体结构面产状数字化测量方法。

在激光点云数据处理的过程中，本发明采用基于密度的法向量聚类算法以自动 化的形式提取出主要的结构面产状信息，并通过TIN中三角形的连通性判别将结构 面进行分组。在此过程中，没有现有的K值聚类、模糊聚类等需要人工假定聚类结 果组数的缺点，数据处理结果不受人工干预的影响，因此能自动化识别出各种形状 的岩体结构面。

在数字影像数据处理的过程中，本发明采用混合全局和局部阈值法提取出迹线 的轮廓线，并通过骨架提取、迹线标记提取出迹线的产状信息。经实践检验，混合 全局和局部阈值法在岩体节理迹线的轮廓线提取方面能够达到较高精度，骨架提取 所采用判定条件也完全适用于迹线处理，而随后的迹线标记更能同时发挥数字影像 数据信息丰富与激光点云数据精度高的优点。

综上所述，本发明在产状测量过程中无人工干预，采集数据后结合一系列数字 处理方法能获取高精度的岩体结构面产状信息，是一种全自动、非接触式的岩体结 构面产状数字化测量方法。

附图说明

图1为本测量方法的技术方案流程图。

具体实施方式

本发明的融合激光点云和数字影像的岩体结构面产状全自动数字化测量方法， 主要包括现场数据采集、数据融合与配准、点云数据特征提取、点云部分的产状信 息获取、影像特征提取、影像部分的产状信息获取、激光点云与数字影像融合，具 体步骤如下：

第一步：现场数据采集

1)仪器现场安装

根据测量区域位置与大小、岩体结构特征以及复杂工程环境对目标对象可能造 成光线遮挡与反射的影响，确定TLS三维激光扫描仪架设测点的站点数量、测站点 具体位置以及仪器架设的高度、扫描角度，并且按照上述要求组装和架设TLS三维 激光扫描仪，若需要对不同站点的扫描数据进行拼接，还需要确定合理的控制点布 设位置，以便对不同站点的点云数据进行配准；

2)现场扫描测量

根据拟测量区的岩体结构完整性，选择TLS三维激光扫描仪的合适扫描分辨率、 校正参数、扫描角度。根据仪器架设位置距测区的距离与光线，确定高分辨率数码 相机的曝光度、光圈、快门时间、焦距等参数，以期保证能够获得高质量数字图像 数据；

设置好参数后，用扫描仪对需要扫描的区域进行激光扫描，从而获取点云数据， 目标物体距离扫描仪越远，目标物体表面扫描点的间隔就越大，根据需要适当地调 整扫描精度；

在激光扫描的同时，利用置于扫描仪顶部的高分辨率数码相机进行岩体结构特 征的图像采集，从而获取图像数据；

第二步：数据配准与融合

1)相机内部参数获取

三维激光扫描仪的自带数码相机一般均为固定焦距相机，其焦距可以直接从镜 头标识中获得，数码相机的内部参数只与相机内部结构和镜头有关，可由相机及镜 头的出厂说明书与技术规格表中获得；

2)相机外部参数获取

利用三维扫描仪系统采集数据时，相机是随着扫描仪转动而转动的，因此在进 行相机外部参数确定时需要考虑相机相对于初始坐标系的旋转矩阵，扫描仪系统会 涉及到四个坐标系：相机坐标系CMCS、扫描仪坐标系SOCS、工程项目坐标系PRCS、 世界坐标系GLCS；三维扫描仪厂商提供了坐标旋转、平移相关矩阵参数，主要有 ①Mounting矩阵，相机坐标系与扫描仪坐标系之间的坐标转换参数矩阵；②COP矩 阵，相机拍摄瞬间相机坐标系和初始相机坐标系的旋转矩阵；③SOP矩阵，各扫描 站点坐标系相当于工程坐标系的旋转平移矩阵；

3)坐标系统的相互转换

①世界坐标系与相机坐标系

世界坐标P_w(X_w,Y_w,Z_w)，相机坐标P_u(X_u,Y_u,Z_u)，两者转换关系：

$\left(\begin{array}{c}{X}_u\\ Y_u\\ Z_u\\ 1\end{array}\right)=\left(\begin{array}{cc}R& t\\ 0^T& 1\end{array}\right)\left(\begin{array}{c}{X}_w\\ Y_w\\ Z_w\\ 1\end{array}\right)=\left(\begin{array}{cccc}{r}_{11}& r_{12}& r_{13}& t_1\\ r_{21}& r_{22}& r_{23}& t_2\\ r_{31}& r_{32}& r_{33}& t_3\\ 0& 0& 0& 1\end{array}\right)\left(\begin{array}{c}{X}_w\\ Y_w\\ Z_w\\ 1\end{array}\right)=M·{\mathrm{COP}}^{-1}·{\mathrm{SOP}}^{-1}\left(\begin{array}{c}{X}_w\\ Y_w\\ Z_w\\ 1\end{array}\right)---(2-1)$

其中，旋转矩阵t是世界坐标系原点在相机坐标系中的位置坐标，旋转矩阵R为正交 旋转矩阵，满足：

$\left(\begin{array}{c}{r}_{11}^2+r_{12}^2+r_{13}^2=1\\ r_{21}^2+r_{22}^2+r_{23}^2=1\\ r_{31}^2+r_{32}^2+r_{33}^2=1\end{array}\right)---(2-2)$

②相机坐标系与图像坐标系

相机坐标P_u(X_u,Y_u,Z_u)，图像坐标为(x,y)，两者转换关系：

$Z_u\left(\begin{array}{c}x\\ y\\ 1\end{array}\right)=\left(\begin{array}{cccc}f& 0& 0& 0\\ 0& f& 0& 0\\ 0& 0& 1& 0\end{array}\right)\left(\begin{array}{c}{X}_u\\ Y_u\\ Z_u\\ 1\end{array}\right)---(2-3)$

其中，f为相机焦距。

③世界坐标系与图像坐标系

世界坐标P_w(X_w，Y_w，Z_w)，图像坐标为(x，y)，两者转换关系：

$Z_u\left(\begin{array}{c}x\\ y\\ 1\end{array}\right)=\left(\begin{array}{cccc}f& 0& 0& 0\\ 0& f& 0& 0\\ 0& 0& 1& 0\end{array}\right)\left(\begin{array}{c}{X}_u\\ Y_u\\ Z_u\\ 1\end{array}\right)=\left(\begin{array}{cccc}f& 0& 0& 0\\ 0& f& 0& 0\\ 0& 0& 1& 0\end{array}\right)M·{\mathrm{COP}}^{-1}·{\mathrm{SOP}}^{-1}\left(\begin{array}{c}{X}_w\\ Y_w\\ Z_w\\ 1\end{array}\right)---(2-4)$

第三步：点云特征信息提取

1)点云数据去噪

本发明采用一般的中值滤波法对点云数据进行去噪处理，中值滤波法的原理就 是用点云目标点的邻域内诸点坐标的中间值来替代点云目标点的坐标值，

2)三角剖分

三角剖分是拓扑学中最基本的方法，它能将零散的点云剖分为无数曲边三角形， 在此将用Delaunay三角剖分法进行三角剖分，运用Matlab中的delaunay函数实现 三维点云的Delaunay剖分，delaunay函数输出的数据是完成划分后各个三角形的编 号和顶点坐标，delaunay函数在Matlab中的用法如下：

Ⅰ、输入所有点的x,y,z坐标；

Ⅱ、运用delaunay函数进行三角划分，得到各三角形的ID编号和顶点坐标：

Tri＝delaunay(x,y)；

Ⅲ、调用z坐标，运用trimesh或者trisurf指令画出三维曲面图。

3)三角形的法向向量计算

利用剖分后三角形的三个顶点的坐标，计算出三角形的法向量坐标和三角形的 形心坐标，并进一步计算出法向量与坐标Z轴、X-Z面的夹角φ_T、θ_T,形心距坐标原点 的距离r_T，

设Matlab中算得的某三角形的三个顶点为

A(x_T1,y_T1,z_T1),B(x_T2,y_T2,z_T2)，C(x_T3,y_T3,z_T3)，则

$\overrightarrow{AB}=(x_{T2}-x_{T1},y_{T2}-y_{T1},z_{T2}-z_{T1}),\overrightarrow{AC}=(x_{T3}-x_{T1},y_{T3}-y_{T1},z_{T3}-z_{T1}),$

三角形的形心坐标即为则形心O距坐标原点的距 离$r_T=\sqrt{{(\underset{k=1}{\overset{3}{\Sigma }}{x}_{Tk}/3)}^2+{(\underset{k=1}{\overset{3}{\Sigma }}{y}_{Tk}/3)}^2+{(\underset{k=1}{\overset{3}{\Sigma }}{z}_{Tk}/3)}^2}---(3-1)$

设ΔABC的法向量为通过构建如下方程组可得求法向量坐标：

$\left(\begin{array}{c}(x_{T2}-x_{T1})a_T+(y_{T2}-y_{T1})b_T+(z_{T2}-z_{T1})c_T=0\\ (x_{T2}-x_{T3})a_T+(y_{T2}-y_{T3})b_T+(z_{T2}-z_{T3})c_T=0\\ (x_{T3}-x_{T1})a_T+(y_{T3}-y_{T1})b_T+(z_{T3}-z_{T1})c_T=0\end{array}\right)---(3-2)$

随后利用三角函数关系求得法向量的φ_T、θ_T值：

${\theta }_T=\arccos \frac{c_T}{\sqrt{{a_T}^2+{b_T}^2+{c_T}^2}}$

第四步：产状信息获取即点云部分

1)法向量聚类分析

①数据中心化

数据中心化是聚类分析前的必要准备工作，本发明中数据中心化的方法，是用单个 三角形的φ_T、θ_T、r_T值减去所有三角形φ_T、θ_T、r_T值的平均值；

点云数据总共被剖分为n个三角形，那么能计算并得到对应n个法向量的坐标 (a_Ti,b_Ti,1)(i＝1,2,3...n)及φ_T、θ_T、r_T值φ_Ti、θ_Ti、r_Ti(i＝1,2,3....n)；将n个法向量的φ_T、θ_T、r_T值 用如下的n×3阶矩阵T表述，并在数据处理过程中保留ID：

$T=\left(x_{ij}\right)=\left(\begin{array}{ccc}{x}_{11}& x_{12}& x_{13}\\ x_{21}& x_{22}& x_{23}\\ \begin{array}{c}·\\ ·\\ ·\end{array}& \begin{array}{c}·\\ ·\\ ·\end{array}& \begin{array}{c}·\\ ·\\ ·\end{array}\\ x_{n1}& x_{n2}& x_{n3}\end{array}\right),$φ_Ti＝x_i1,θ_Ti＝x_i2，r_T＝(x_i3),(i＝1,2,…n),(j＝1,2,3)

记第j列的平均值为：

$\overline{x_j}=\frac{1}{n}\underset{i=1}{\overset{n}{\Sigma }}{x}_{ij}(i=1,2,...n),(j=1,2,3)---(4-1)$

那么对第j列的n个数据对象所实施的中心化变换为：

${x^{,}}_{ij}=x_{ij}-{\bar{x}}_j,(i=1,2...n),(j=1,2,3)---(4-2)$

经过上式的变换后，每列变量的均值为0，即每列变量的取值都有相同的基点， 然后将中心变换后的结果重新生成为新的矩阵T,:

$T^{,}=\left({x^{,}}_{ij}\right)=\left(\begin{array}{ccc}{x^{,}}_{11}& {x^{,}}_{12}& {x^{,}}_{13}\\ {x^{,}}_{21}& {x^{,}}_{22}& {x^{,}}_{23}\\ \begin{array}{c}·\\ ·\\ ·\end{array}& \begin{array}{c}·\\ ·\\ ·\end{array}& \begin{array}{c}·\\ ·\\ ·\end{array}\\ {x^{,}}_{n1}& {x^{,}}_{n2}& {x^{,}}_{n3}\end{array}\right),(i=1,2,...n),(j=1,2,3).$

②数据标准化

在此采用最大最小归一法对数据进行标准化处理，即将矩阵中的某一数据除以 该列数据绝对值的最大值，具体见下式：

${x_{ij}}^{,,}=\frac{{x_{ij}}^{,}}{{\max }_i\left|{x_{ij}}^{,}\right|}---(4-3)$

最终生成新的矩阵：

$T^{,,}=\left({x^{,,}}_{ij}\right)=\left(\begin{array}{ccc}{x^{,,}}_{11}& {x^{,,}}_{12}& {x^{,,}}_{13}\\ {x^{,,}}_{21}& {x^{,,}}_{22}& {x^{,,}}_{23}\\ \begin{array}{c}·\\ ·\\ ·\end{array}& \begin{array}{c}·\\ ·\\ ·\end{array}& \begin{array}{c}·\\ ·\\ ·\end{array}\\ {x^{,,}}_{n1}& {x^{,,}}_{n2}& {x^{,,}}_{n3}\end{array}\right)$

③特征生成

上述矩阵T”属于低维矩阵，为描述简便、形象，将T”的各个行向量投影于三维 直角坐标系中形成n个点，并以(x”_i1,x”_i2,x”_i3)作为这n个点的坐标，记为(x_i,y_i,z_i)；如 果有坐标重复的点即T”内的行向量数值相同，因为其ID不同，予以保留；

④法向量的聚类分析

DBSCAN聚类分析法是一种基于高密度联通区域的聚类算法，本方法运用并修 改了该聚类算法，从而达到提取结构面信息的目的；

首先，取ε＝3，MinPts＝5，T”中两点q₁(x_q1,y_q1,z_q2),q₂(x_q2,y_q2,z_q2)，其中，参数ε 为领域半径；以一点为圆心，ε为半径的圆，称为该点的ε领域；参数MinPts为最 小核心对象数；如果某点的ε领域内包含的点的数量大于等于MinPts，则该点称为 核心对象；具体工程中参数ε、MinPts的值需要调整；然后，计算点q₁到点q₂的距 离d，若d≤3,则将点q₂归于点q₁的ε领域中，

$d=\sqrt{{(x_{q1}-x_{q2})}^2+{(y_{q1}-y_{q2})}^2+{(z_{q1}-z_{q2})}^2}---(4-4)$

之后，输入T”中的所有点的坐标，计算并找出能够归于点q₁的ε领域内的所有点；

假设得到q₁的ε领域内包含的m个点为{q₂,q₃,…q_m+1}，则q₁的ε领域内包含的点 数为m；

如果m≥5，那么将q₁归为核心对象集合A₁内的点；如果0＜m＜5，则将q₁归为边 缘对象集合B₁内的点；如果m＝0,则将q₁归为噪音；

再对q₁的领域内的其他点{q₂,q₃,…q_m+1}按如上方法进行计算并判别是否为核心对 象，将点q₁的领域内所有被归为核心对象的点归于集合A₁，所有被归为边缘对象的点 归于集合B₁，删除被归为噪音的点；

重复输入T”内的所有点坐标，计算出每个点的核心对象集合{A₁,A₂…A_n}、边缘 对象集合{B₁,B₂…B_n}，其中噪音的核心对象集合与边缘对象集合都记为0，再将这 两个大集合按如下方法分别进行合并；

取{A₁,A₂…A_n}内的两个集合A₁、A₂，若集合即有重复点，则将A₂内不与 A₁相交的样本并入A₁中形成簇A₁’，并将A₁’作为新的合并对象参与下一次合并判断；若 集合即没有重复点，则将A₁赋值到簇A₁’中，A₂赋值到簇A₂’中，并将A₁’、A₂’作为 新的合并对象参与下一次合并判断，然后取{A₁,A₂…A_n}中A₃分别与A₁’、A₂’进行上述合 并判别，如果且则将A₃赋值到簇A₃’中，并将A₃’作为新的合并对 象参与下一次合并判断。如此重复，直至{A₁,A₂…A_n}中所有样本都进行完合并判别， 最终得到新的集合{A₁’,A₂’…A_v’}；

边缘对象集合{B₁,B₂…B_n}的合并也按上述原则进行合并，得到{B₁’,B₂’…B_w’}。

{A₁’,A₂’…A_v’}和{B₁’,B₂’…B_w’}则是聚类的最终结果，再利用聚类结果中各元素的ID 编号就能提取出对应法向量的φ_t、θ_t、r值，也能标记出对应的三角面片；

2)聚类结果处理

①三角面片合并

该步将除无效三角面片，弥补DBSCAN算法的不足，并将聚类后的结果提取、 标示，通过合并三角面构建出岩体结构面。

Ⅰ、通过DBSCAN的聚类结果利用ID标记出各个簇中对应的三角面片，删除 重复的三角面片，并对各个簇进行标号。

Ⅱ、在簇A_h’(A_h’为{A₁’,A₂’…A_v’}中任意一样本)中任取一三角面片i，若其邻近存 在同属于A_h’的其他三角面片，则将i保留；否则，delete。

Ⅲ、经过Ⅱ生成的新簇A_h’中，任取两三角面片ii、iii，若ii、iii能通过同属于A_h’ 中的任意多个三角面片相连通，则将ii、iii归入结构面X₁；否则ii归入X₁，iii归入 X₂。如此，将A_h’中的三角面片进行结构面分组；

Ⅳ、分别将{A₁’,A₂’…A_v’}中对应的三角面片进行Ⅱ、Ⅲ操作，得到X₁,X₂…X_s；

Ⅴ、分别将{B₁’,B₂’…B_w’}中对应三角面片进行Ⅱ、Ⅲ操作，得到Y₁,Y₂…Y_r。

最终得到的X₁,X₂…X_s，Y₁,Y₂…Y_r就是构成岩体结构面的三角面片，将这些三角 面片借由ID标示出来后即能得到岩体的各个结构面；

②剔除开挖面

实际工程中，岩体结构面的产状信息中包含了大量的开挖面信息，而开挖面并 不属于岩体的结构面。由于开挖面产状具有重复性、单一性的特征，只需要用人工 剔除的方法在RiscanPro中就能将其识别并将开挖面的三角面片删除；

③结构面产状计算

结构面的产状信息包括结构面的倾角与倾向，只要对X₁,X₂…X_s中三角面的参数 进行计算，就能得到结构面的产状信息，具体算法如下：

假设X₁对应H个三角面，计算X₁中三角面的法向量坐标(a_Td,b_Td,c_Td)、法向量 φ_T、θ_T值(φ_Td、θ_Td)(d与X₁中对应三角面的ID相同)的算数平均值：

$a_{X1}=\frac{1}{H}{\mathrm{\Sigma a}}_{Td}$

$b_{X1}=\frac{1}{H}{\mathrm{\Sigma b}}_{Td}---(4-5)$

$c_{X1}=\frac{1}{H}{\mathrm{\Sigma c}}_{Td}$

${\phi }_{X1}=\frac{1}{H}{\mathrm{\Sigma \phi }}_{Td}$

${\theta }_{X1}=\frac{1}{H}{\mathrm{\Sigma \theta }}_{Td}$

则β_X1＝|90°-θ_X1|(0°≤β_X1≤90°)即为结构面X₁的倾角；

因倾向角度范围为0°～360°，为精确计算倾向方向，需要利用X₁的法向量与Y、 Z轴的余弦值来判断X₁的法向量指向的卦限。具体判别方法如下：

假设X₁的法向量为直线L为X₁与坐标系中X-Y面的交线，向量为的 水平投影向量，根据倾向方向的定义，为结构面X₁倾向；α_X1、γ_X1分别为与X、 Z轴的夹角；k_X1为与Y轴的夹角。则有：

${\mathrm{cosk}}_{X1}=\frac{b_{X1}}{\sqrt{{a_{X1}}^2+{b_{X1}}^2}}$

${\mathrm{cos\alpha }}_{X1}=\frac{a_{X1}}{\sqrt{{a_{X1}}^2+{b_{X1}}^2+{c_{X1}}^2}}$

${\mathrm{cos\gamma }}_{X1}=\frac{c_{X1}}{\sqrt{{a_{X1}}^2+{b_{X1}}^2+{c_{X1}}^2}}$

计算出的即为结构面X₁的倾向；

重复上述步骤，带入剩余结构面的参数，就能算出所有结构面的倾角与倾向；

④补全结构面

在采集、处理等一系列操作中可能造成少量点云、三角面片的缺失，影响最终 展示效果，将得到的三角面片导入到Geomagic、Pointcloud等软件中进行补全，最终 得到的各簇所对应的轮廓面即可视为岩体结构面；

第五步：影像数据特征提取

1)图像光学矫正

根据数码相机成像数学模型，矫正参数分为数码相机的内部参数和外部参数， 式5-1为数码相机成像数学模型公式，

$Z_c\left(\begin{array}{c}U\\ V\\ 1\end{array}\right)=\left(\begin{array}{ccc}\frac{1}{d_x}& 0& u_0\\ 0& \frac{1}{d_y}& v_0\\ 0& 0& 1\end{array}\right)\left(\begin{array}{cccc}f& 0& 0& 0\\ 0& f& 0& 0\\ 0& 0& 1& 0\end{array}\right)\left(\begin{array}{cc}R& t\\ 0^T& 1\end{array}\right)\left(\begin{array}{c}{X}_W\\ Y_W\\ Z_W\\ 1\end{array}\right)---(5-1)$

令：

$M_1=\left(\begin{array}{ccc}\frac{1}{d_x}& 0& u_0\\ 0& \frac{1}{d_y}& v_0\\ 0& 0& 1\end{array}\right)M_2=\left(\begin{array}{cccc}f& 0& 0& 0\\ 0& f& 0& 0\\ 0& 0& 1& 0\end{array}\right)M_3=\left(\begin{array}{cc}R& t\\ 0^T& 1\end{array}\right)$

其中，Z_c为照片成像平面相对于镜头的距离；U与V为照片影像在成像平面内 像素坐标系下的坐标；X_W、Y_W、Z_W为现实物体在全局世界坐标系下的坐标；M₁、 M₂为数码相机的内部参数，内部参数只与相机内部结构和镜头有关，可由相机及镜 头的出厂说明书与技术规格表中获得；M₃为数码相机的外部参数，外部参数可由相 机标定实验经相机成像数学模型公式反演获得；

当得到矫正参数后，即可按照公式5-1对影像图片进行矫正，一般可采用成熟商 业软件或者三维激光扫描仪附带成套软件来进行。

2)图像灰度化

数码相机所拍摄数字图像的初始状态是彩色图像，首先需要对其转成灰度图像。 灰度图像是指只含有亮度信息的数字图像，且亮度值变化连续。将高分辨率数码图 像转换为灰度图并且数字化，实际上就是将图像转化为一个灰度值矩阵F(M,N)。 即表明图像大小为M×N个像素由M×N阶的矩阵表达，矩阵中每一个值表达为像素 单元的灰度值。式5-2为灰度值矩阵的表达形式。

3)轮廓线提取

轮廓线提取的方法是采用计算机图像处理方法中的图像分割法。本发明提出了 一种适合岩石节理轮廓线提取的混合全局和局部阀值法。

第一步利用岩石节理区的灰度值一般为区域局部最小值的特点，采用局部阀值 法，经过一系列阀值测试，可以得到局部阀值法中的min(x,y)为以(x,y)为中心的7×7 像素格网邻域内的局部最小灰度值，这样可以减小节理区的噪声干扰；令a(x,y)为 以(x,y)为中心的70×70像素格网邻域内的平均灰度值。首先，通过局部阀值法找出原 始图像中满足条件min(x,y)≥a(x,y)的像素格网点。第二步令满足条件像素格网点的 f(x,y)＝0形成中间过渡图像。第三步对中间过渡图像采用全局阀值法的Otsu法来确 定分割阈值T_O,进而进行图像分割来提取岩石节理轮廓线。Otsu法是根据统计理论 来寻找阈值的，Otsu法的最佳阀值是由背景图像与目标图像的类间方差最大值来确 定的。设图像中像素点的总和为N，灰度级l上的像素点总数为n₁，N与n₁的关系如 式5-3表示，图像直方图像素点的概率分布p_l符合式5-4，Otsu法的最佳阈值T_O最终 由式5-5求出。

$N={\Sigma }_{l=0}^{l-1}{n}_l---(5-3)$

$p_l=\frac{n_l}{N},{\Sigma }_{l=0}^{l-1}{p}_l=1---(5-4)$

$T_O=\underset{0<l<1}{\mathrm{argmax}}\left\{{\sigma }_B^2\left(l\right)\right\}---(5-5)$

式中，为类间方差。

4)图像去噪，

图像噪声对后续的在岩石节理骨架及拓扑关系提取会产生很大的干扰，所以应 对岩石节理轮廓线图像去噪。

具体方法如下：

Ⅰ、对岩石节理轮廓线图像进行二值化处理；

Ⅱ、消除岩石节理轮廓线图像中的岩石节理区里的黑斑或者岩石区里的白斑。 运用数字图像形态学中膨胀运算算子。原图像A被结构元素B膨胀可定义为：将结 构元素B的反射平移X个像素后仍与A有交点的所有的点x组成。即运算公式为：

本发明使用线结构元素对节理图像进行膨胀运算。

Ⅲ、膨胀处理后的节理区具有良好的连通性，运用Matlab软件中的Bwareaopen 函数，根据连通区大小来过滤掉图像中不需要的小面积部分；

Ⅳ、在第Ⅲ步处理后节理轮廓线的边缘会有一些不规则的细小毛刺，对后续 节理骨架提取产生影响，所以需要对不规则细小毛刺进行边缘光滑。本文运用中值 滤波法进行光滑处理。中值滤波原理是：给定的D个数值{a₁,a₂...a_D}按大小有序排列， 当D为奇数时，位于中间位置的那个数值被称为这D个数值的中值；当D为偶数时， 为于中间位置的那两个数值的平均值为这D个数值的中值，记作med{a₁,a₂...a_D}，邻 域窗口内所有像素点的灰度中值作为窗口中被滤波的像素点的灰度值。即：图像为 [x(I,J)]_M×N的矩阵，领域窗口为A_D,中值滤波后像素点x(I,J)的图像输出y(I,J)记为：

$y(I,J)={\mathrm{med}}_{A_d(I,J)}\left[x(I,J)\right]---(5-7)$

本文使用5×5的正方形窗口对图像进行中值滤波，

5)节理迹线骨架提取

节理轮廓线不能直观地表达岩石节理骨架特征及其扑拓结构，需要进一步提取 节理迹线骨架，提取工作分为两个步骤：①、图像细化；②、去除断枝；

①、图像细化就是将二值图像中的像素点在保持原有形状与连通性的基础上进 行一层一层地像素点剥离，直到残余图像骨架，图像骨架即为图像的中轴线，图像 细化后图像剩余像素宽度为l；

剥离条件如下：

★判断二值图中的点8个邻域内是否满足白点个数为2到6个；

★判断二值图中的点8个邻域内是否满足白点都连续；

★判断二值图中的点是否满足上，左，右不全为白点；

★判断二值图中的点是否满足上，左，下不全为白点；

②、去除断枝就是在图像细化后依旧存在毛刺断枝需要剔除。去除断枝需要对 节理骨架线的端点与交叉点进行识别；

端点识别：在二值图像中令黑色像素点值为0；白色像素点值为1，采用3×3的 正方形检测窗口迭代遍历整个图像，检查所有像素点值为1的像素点，检验其邻域中 8个相邻点中像素点值为1的点个数，若个数为1，则检验得到该像素点为端点；

交叉点识别：交叉点存在三线交叉或者四线叉点。第一步采用3×3的正方形检 测窗口检查遍历像素点值为0的像素点邻域中8个相邻点中像素点值为1的点个数为3 或4；第二步扩大检测窗口为5×5正方形，对3×3窗口最外侧的像素点值为1的点个数 进行统计，若与3×3窗口内侧的像素点值为1的点个数相同，则检验得到该像素点为 交叉点；

去除断枝：应用端点和交叉点将图像分割成各自单独的连通区域，计算每个连 通区域的大小，即连通域内像素点个数N_i，设置连通域内像素点个数阈值N_t，若 N_i＜N_t，就将其当作断枝，予以删除；

6)节理迹线标记

节理迹线信息统计需要对每一条岩石节理迹线进行标识，方法如下：

Ⅰ、搜索图像像素点集的连通性，并以不同颜色标记，利用Matlab软件的Bwlabel 函数寻找图像中的每一个连通对象，并且根据寻找到的顺序用不同的整数值标注，

Ⅱ、Bwlabel函数返回值是double类型的标记矩阵，因而无法实现图像显示，接 下来利用Matlab软件的label2rgb函数，将标记矩阵中指定每一个连通区的整数值定 义为颜色矩阵，以RGB值图像显示；

第六步：产状信息获取即影像部分

根据迹线标记与数据配准，获得每一条标记迹线起始坐标点与任一中间点的坐 标，假设迹线的起始坐标点与任一中间点的坐标分别为E(x_E,y_E,z_E)，F(x_F,y_F,z_F)， G(x_G,y_G,z_G)，三点可确定一个三角形，可求出对应的三个内角θ_E，θ_F，θ_G：

记：

$\overrightarrow{EF}=\overrightarrow{m}=(x_F-x_E,y_F-y_E,z_F-z_E)=(a_G,b_G,c_G)$

$\overrightarrow{EG}=\overrightarrow{n}=(x_G-x_E,y_G-y_E,z_G-z_E)=(a_F,b_F,c_F)$

$\overrightarrow{FG}=\overrightarrow{p}=(x_G-x_F,y_G-y_F,z_G-z_F)=(a_E,b_E,c_E)$

由${\mathrm{cos\theta }}_E=\frac{\overrightarrow{m}·\overrightarrow{n}}{\left|\overrightarrow{m}\right|\left|\overrightarrow{n}\right|}$得：

${\theta }_E=\arccos \frac{\overrightarrow{m}·\overrightarrow{n}}{\left|\overrightarrow{m}\right|\left|\overrightarrow{n}\right|}=\arccos \frac{a_G{a}_F+b_G{b}_F+c_G{c}_F}{\sqrt{a_G^2+b_G^2+c_G^2}\sqrt{a_F^2+b_F^2+c_F^2}}---(6-1)$

再运用公式6-1同理可求出θ_F，θ_G

记θ_J＝max{θ_E，θ_F，θ_G}

①如果θ_J≤120°，则此迹线上三点所构三角形可视为出露面，运用式6-2，即可 求此出露面法向量

$\left(\begin{array}{c}{a}_E{a}_J+b_E{b}_J+c_E{c}_J=0\\ a_F{a}_J+b_F{b}_J+c_F{c}_J=0\\ a_G{a}_J+b_G{b}_J+c_G{c}_J=0\end{array}\right)---(6-2)$

法向量与X-Y平面法向量{0，0}的交角，即出露面的倾角，表示为

因倾向角度范围为0°～360°，为精确计算倾向方向，需要利用法向向量与Y， Z轴的余弦值来判断结构面的法向向量指向的卦限；具体判别方法如下：假设平面P 为一结构面，直线L为该结构面与水平面的交线，向量为平面P的法向向量的水 平投影向量，根据倾向方向的定义，为结构面的倾向；α_J，γ_J分别为法向向量与X，Z轴的角；k_J为向量与Y轴的夹角。则有

${\mathrm{cosk}}_J=\frac{b_J}{\sqrt{{a_J}^2+{b_J}^2}}$

${\mathrm{cos\alpha }}_J=\frac{a_J}{\sqrt{{a_J}^2+{b_J}^2+{c_J}^2}}$

${\mathrm{cos\gamma }}_J=\frac{c_J}{\sqrt{{a_J}^2+{b_J}^2+{c_J}^2}}$

②如果θ＞120°，则将此迹线视为普通迹线；

利用E、G点坐标即可求出线EG方程：z＝Qx+Py，任取线上一点(a_L,b_L,c_L)

线EG与坐标平面X-Y的夹角，即为迹线倾角，表示为

因倾向角度范围为0°～360°，为精确计算倾向方向，需要利用向量与Y，Z 轴的余弦值来判断结构面的法向向量指向的卦限；具体判别方法如下：假设向量向量为的水平投影向量，根据倾向方向的定义，为迹线EG的倾向；α_L，γ_L分别 为向量与X，Z轴的角；k_L为向量与Y轴的夹角，则有

${\mathrm{cosk}}_L=\frac{b_L}{\sqrt{{a_L}^2+{b_L}^2}}$

${\mathrm{cos\alpha }}_L=\frac{a_L}{\sqrt{{a_L}^2+{b_L}^2+{c_L}^2}}$

${\mathrm{cos\gamma }}_L=\frac{c_L}{\sqrt{{a_L}^2+{b_L}^2+{c_L}^2}}$

第七步：激光点云与数字影像融合

1)近似拟合特征平面

前面的工作中已获得各簇三角面片，根据ID提取出三角面的顶点坐标，假设X₁中某三角形三个顶点坐标为(x_R1,y_R1,z_R1)，(x_R2,y_R2,z_R2)，(x_R3,y_R3,z_R3)，那么可根据下 式求得该三角形的形心坐标(x_R,y_R,z_R)用形心表征此三角面：

$x_R=\frac{x_{R1}+x_{R2}+x_{R3}}{3}$

$y_R=\frac{y_{R1}+y_{R2}+y_{R3}}{3}---(7-1)$

$z_R=\frac{z_{R1}+z_{R2}+z_{R3}}{3}$

计算出同簇三角面片的形心坐标后则可通过最小二乘拟合的方法拟合出近似的 特征平面，空间平面的方程通常形如：A_Rx+B_Ry+C_Rz+1＝0，用同簇的H个三角面 形心坐标拟合出这个平面可以表示成以下矩阵形式：

$\left(\begin{array}{ccc}{x}_{R1}& y_{R1}& z_{R1}\\ \begin{array}{c}·\\ ·\\ ·\end{array}& \begin{array}{c}·\\ ·\\ ·\end{array}& \begin{array}{c}·\\ ·\\ ·\end{array}\\ x_{RH}& y_{RH}& z_{RH}\end{array}\right)\left(\begin{array}{c}{A}_R\\ B_R\\ C_R\end{array}\right)=\left(\begin{array}{c}-1\\ -1\\ -1\end{array}\right)$

经过运算可化成以下形式：

$\left(\begin{array}{ccc}{\mathrm{\Sigma x}}_g^2& {\mathrm{\Sigma x}}_g{y}_g& {\mathrm{\Sigma x}}_g{z}_g\\ {\mathrm{\Sigma x}}_g{y}_g& {\mathrm{\Sigma y}}_g^2& {\mathrm{\Sigma y}}_g{z}_g\\ {\mathrm{\Sigma x}}_g{z}_g& {\mathrm{\Sigma y}}_g{z}_g& {\mathrm{\Sigma z}}_g^2\end{array}\right)\left(\begin{array}{c}{A}_R\\ B_R\\ C_R\end{array}\right)=\left(\begin{array}{c}-\Sigma x_g\\ -\Sigma y_g\\ -{\mathrm{\Sigma z}}_g\end{array}\right)(g=1,\mathrm{2...}H)$

所以得到：$\left(\begin{array}{c}{A}_R\\ B_R\\ C_R\end{array}\right)={\left(\begin{array}{ccc}{\mathrm{\Sigma x}}_g^2& {\mathrm{\Sigma x}}_g{y}_g& {\mathrm{\Sigma x}}_g{z}_g\\ {\mathrm{\Sigma x}}_g{y}_g& {\mathrm{\Sigma y}}_g^2& {\mathrm{\Sigma y}}_g{z}_g\\ {\mathrm{\Sigma x}}_g{z}_g& {\mathrm{\Sigma y}}_g{z}_g& {\mathrm{\Sigma z}}_g^2\end{array}\right)}^{-1}\left(\begin{array}{c}-\Sigma x_g\\ -{\mathrm{\Sigma y}}_g\\ -{\mathrm{\Sigma z}}_g\end{array}\right)(g=1,\mathrm{2...}H),$系数A_R,B_R,C_R即可求得，平面方程即可得出；

2)迹线与结构面归并

从点云数据中，提取出迹线所对应的各个点的坐标，逐点求出点到各个簇所拟 合出的近似平面的距离，取最小距离即能描述该迹线与各簇结构面的位置关系，

点到平面距离的算法：$D=\frac{Ax+By+Cz+1}{\sqrt{A^2+B^2+C^2}}---(7-2)$

每个簇周围都存在与其距离最近的1条或多条迹线，可将这1条或多条迹线与该 簇归并入同一个组，那么该组的产状信息即为对应区域的产状信息。若存在某迹线 周围不存在相近的结构面，则将该迹线单独列为一组，用以表征所在区域的产状信 息；若存在某结构面周围不存在相近的迹线，则将该结构面单独列为一组，用以表 征所在区域的产状信息。