CN105180890A

CN105180890A - Rock mass structural plane attitude measuring method integrating laser point cloud and digital image

Info

Publication number: CN105180890A
Application number: CN201510452792.0A
Authority: CN
Inventors: 张鹏; 戴静; 李俊才; 蒋立辰
Original assignee: Nanjing Tech University
Current assignee: Nanjing Tech University
Priority date: 2015-07-28
Filing date: 2015-07-28
Publication date: 2015-12-23
Anticipated expiration: 2035-07-28
Also published as: CN105180890B

Abstract

A fully automatic digital measurement method for the occurrence of rock mass discontinuities that combines laser point clouds and digital images. The present invention uses photogrammetry technology to extract the contour line of the track by using the mixed global and local threshold method, extracts the joint track through skeleton extraction and track mark, and combines other image processing methods to realize the occurrence of the rock joint track Finally, according to the distance between the trace and the structural surface, the positional relationship between the trace and the structural surface is judged, and the trace and structural surface are merged and grouped, and the occurrence of the rock mass structural surface is represented by each group of occurrences. The invention integrates the measurement method of laser point cloud and digital image, fully utilizes the fast, efficient, non-contact and terrain-ignoring properties of these two measurement methods, improves the accuracy of measurement at the same time, and has a high degree of automation, universal It has the advantages of strong adaptability and comprehensive expression.

Description

Measuring method of rock mass discontinuity plane occurrence by fusion of laser point cloud and digital image

技术领域technical field

本发明涉及一种融合激光点云和数字影像的岩体结构产状全自动数字化测量方法，属于工程地质勘察技术领域。The invention relates to a fully automatic digital measurement method for the occurrence of rock mass structure by combining laser point clouds and digital images, and belongs to the technical field of engineering geological survey.

背景技术Background technique

岩体结构面产状测量方法归纳起来主要有四类：There are four main types of methods for measuring the occurrence of rock mass discontinuities:

①、测线法、窗口统计法即通过皮尺和罗盘人工现场逐一测量结构面产状信息。该方法在面对实际工程中岩体结构面分布广、数量多、随机性强的特点，低效、费力、耗时使之难以满足工程需求，而且遇到人迹罕至陡崖险坡受现场条件的影响，测量方法便难以实施，从而导致测量数据的全面性与代表性也受到限制。①. Line survey method and window statistics method are to manually measure the occurrence information of structural planes one by one by measuring tape and compass on site. In the face of the characteristics of wide distribution, large number and strong randomness of rock mass structural planes in actual engineering, this method is inefficient, laborious and time-consuming, making it difficult to meet the engineering needs, and it is difficult to meet the requirements of the site conditions when encountering inaccessible steep cliffs and dangerous slopes. Therefore, the measurement method is difficult to implement, which leads to the limitation of the comprehensiveness and representativeness of the measurement data.

②、钻孔定向取芯法(Rosengren,1968)、孔内照相法(Eoek&Pentz,1968)、数字全景钻孔摄像法(王川婴,2001)即通过钻孔取芯后对岩芯中所含结构面产状的量测或通过孔内摄像对孔壁内侧结构面纹理进行数字图像量测来获取钻孔深度范围内相交的结构面产状信息。该方法可以获得深部岩体结构面产状信息，但是孔径范围测量尺寸使得该方法对于长大结构面的测量精度较低，同时钻孔破岩与孔位垂直度也会影响到测量精度。②. Borehole directional coring method (Rosengren, 1968), in-hole photography method (Eoek&Pentz, 1968), digital panoramic borehole camera method (Wang Chuanying, 2001), that is, to analyze the structural planes contained in the core after drilling and coring. Occurrence measurement or digital image measurement of the structure surface texture inside the hole wall by in-hole camera to obtain the occurrence information of the intersecting structure surface within the drilling depth range. This method can obtain the occurrence information of deep rock mass structural planes, but the measurement accuracy of the method for long structural planes is low due to the measurement of the aperture range, and the verticality between the rock breaking and the hole position will also affect the measurement accuracy.

③、摄影测量法(Ross-Brown&At-kinson,1972；Hagan,1978)即基于数字图像与摄影测量的基本原理，应用计算机三维成像技术、影像匹配、图像插值、模式识别等多学科理论与方法相融合，可以获取地物表面的相对坐标空间几何信息，通过计算这些几何信息可以得到暴露地表岩体结构面产状信息。该方法采用非接触的测量手段，可以不受现场地形条件的制约，但是该方法精度受摄影图像质量与双目摄像间距控制，同时仍然需要一定的人工干预，单条岩体结构面的工作效率比精测线法有所提高，而面对实际工程中众多结构面大量人工干预仍然是无法接受的。③. Photogrammetry (Ross-Brown & At-kinson, 1972; Hagan, 1978) is based on the basic principles of digital images and photogrammetry, using computer 3D imaging technology, image matching, image interpolation, pattern recognition and other multidisciplinary theories and methods. Fusion can obtain the geometric information of the relative coordinate space of the surface of the ground object, and by calculating the geometric information, the occurrence information of the structural plane of the exposed surface rock mass can be obtained. This method adopts non-contact measurement means, which can not be restricted by the terrain conditions on site, but the accuracy of this method is controlled by the quality of photographic images and the distance between binocular cameras. At the same time, a certain amount of manual intervention is still required. The precise line method has been improved, but it is still unacceptable to deal with a large number of manual interventions on many structural surfaces in actual engineering.

④、地面三维激光扫描法(TerrestrialLaserScanning，简称TLS)，即通过对地物目标的逐点激光扫描测距可快速获得描述地物表面相对几何坐标信息的高精度点云数据(激光数据)，通过一定人工干预与数据截取可以计算出那些暴露在岩壁表面上的岩体结构面产状。该方法取消了摄影测量法中相对繁琐的摄影测量算法，提高了测量速度与精度，但是仍然需要人工干预在海量点云数据中手动截取描述暴露岩体结构面产状的点云数据，而且该方法对于未暴露岩体结构面或以迹线形式显现的岩体结构面产状是无法测量。④ Terrestrial Laser Scanning (TLS for short), that is, through the point-by-point laser scanning ranging of the ground object, the high-precision point cloud data (laser data) that describes the relative geometric coordinate information of the surface of the ground object can be quickly obtained. Certain manual intervention and data interception can calculate the occurrence of rock mass structural planes exposed on the rock wall surface. This method cancels the relatively cumbersome photogrammetry algorithm in the photogrammetry method, and improves the measurement speed and accuracy. However, manual intervention is still required to manually intercept the point cloud data describing the occurrence of the exposed rock mass structural plane from the massive point cloud data. The method cannot measure the occurrence of unexposed rock mass structural planes or rock mass structural planes that appear in the form of traces.

综上所述，面对实际工程中岩体结构面数量多、分布广和随机性强等特点，目前技术瓶颈在于如何兼顾非接触式(取消场地制约)、全自动(取消人工干预)、测量速度快(算法简洁)、普适性(适合多种赋存环境岩体结构面)。To sum up, facing the characteristics of large number of rock mass structural surfaces, wide distribution and strong randomness in actual engineering, the current technical bottleneck lies in how to take into account non-contact (cancellation of site constraints), fully automatic (cancellation of manual intervention), measurement Fast (simple algorithm), universal (suitable for rock mass structural planes in various occurrence environments).

发明内容Contents of the invention

技术问题：本发明目的在于借助地面三维激光扫描技术的快速精准空间测量手段和计算机视觉图像识别技术的目标信息自动化筛取，将含有精确几何信息激光点云数据与目标信息自动化筛取后的数字图像数据进行数据融合与几何配准，提出一种基于TLS三维激光扫描仪的全自动、非接触式岩体结构产状数字化测量方法，即融合激光点云和数字影像的岩体结构面产状测量方法。Technical problem: The purpose of the present invention is to use the rapid and accurate spatial measurement means of ground three-dimensional laser scanning technology and the automatic screening of target information of computer vision image recognition technology to automatically screen the laser point cloud data containing accurate geometric information and target information. Data fusion and geometric registration of image data are carried out, and a fully automatic, non-contact digital measurement method for the occurrence of rock mass structure based on TLS 3D laser scanner is proposed, that is, the occurrence of rock mass structure plane which is fused with laser point cloud and digital image Measurement methods.

技术方案：主要针对野外岩体结构面的数量多、分布广和随机性强等特征，传统岩体结构面产状测量方法工作效率低下，甚至遇到人迹罕至陡崖险坡根本无法实施。本发明借助于地面激光激光的精准空间测量技术和数字影像的计算机视觉识别技术，提出了一种全自动、非接触式的岩体结构面产状数字化测量方法。Technical solution: mainly for the characteristics of large number, wide distribution and strong randomness of rock mass structural planes in the field, the traditional method of measuring the occurrence of rock mass structural planes is inefficient, and it is even impossible to implement when encountering inaccessible and steep cliffs and dangerous slopes. The present invention proposes a fully automatic, non-contact digital measurement method for the occurrence of rock mass structural planes by means of the precise spatial measurement technology of ground lasers and the computer vision recognition technology of digital images.

融合激光点云和数字影像的岩体结构面产状全自动数字化测量方法，其测量内容包括了岩体结构面产状的点云数据提取、岩体节理迹线产状的图像数据与点云数据相融合的提取，其表述方式采用迹线与结构面匹配分组、共同表述的方法。A fully automatic digital measurement method for the occurrence of rock mass structural planes that combines laser point clouds and digital images. The measurement content includes point cloud data extraction of rock mass structural plane occurrences, image data and point cloud For the extraction of data fusion, the expression method adopts the method of matching grouping and common expression between traces and structural surfaces.

在点云数据的提取中，先用滤波法进行去噪，再对点云数据进行三角剖分从而建立初步的岩体三维模型，之后充分发挥点云数据精确度高的特点，通过对三角面、三角面法向量的计算构造出特征矩阵并采用基于高密度联通区域的DBSCAN聚类分析方法对结构面进行初步分组，随后借由三角面的连通性判别进行结构面的二次分组，最后计算出结构面的产状；岩石节理迹线需从图像数据中提取出来，先对数字影像进行灰度化处理，再通过图像混合阈值分割法提取出岩石的轮廓线，经过像素剥离和去除断肢提取出迹线的骨架信息，最后将迹线标记出来并进行产状计算；结构面和迹线的产状测量完成后，将迹线与结构面合理分组，通过计算出迹线和结构面之间的位置关系并进行判别归并，用归并后的最终结果表征岩体结构面的产状信息。In the extraction of point cloud data, the filter method is used to denoise first, and then the point cloud data is triangulated to establish a preliminary three-dimensional model of rock mass. , The calculation of the normal vector of the triangular surface constructs the characteristic matrix, and uses the DBSCAN cluster analysis method based on the high-density connected area to group the structural surfaces initially, and then uses the connectivity of the triangular surfaces to carry out the secondary grouping of the structural surfaces, and finally calculates The occurrence of structural planes; the rock joint traces need to be extracted from the image data. Firstly, the digital image is processed in gray scale, and then the contour line of the rock is extracted by the image mixing threshold segmentation method. After pixel stripping and removal of severed limbs Extract the skeleton information of the trace, and finally mark the trace and calculate the occurrence; The positional relationship among them is identified and merged, and the final result after merging is used to characterize the occurrence information of rock mass discontinuities.

由于通过点云数据提取迹线产状是较为困难的，因此本发明融合图像数据和点云数据进行迹线产状的提取。岩石节理迹线需从图像数据中提取出来，本发明先对数字影像进行灰度化处理，再通过图像混合阈值分割法提取出岩石的轮廓线，经过像素剥离和去除断肢提取出迹线的骨架信息，最后将迹线标记出来并进行产状计算。Since it is difficult to extract the occurrence of the trace through the point cloud data, the present invention fuses the image data and the point cloud data to extract the occurrence of the trace. The rock joint trace needs to be extracted from the image data. The present invention first performs gray-scale processing on the digital image, and then extracts the contour line of the rock through the image mixing threshold segmentation method. Skeleton information, and finally mark the trace and calculate the occurrence.

结构面和迹线的产状测量完成后，为更全面、准确地表述岩体结构面产状，本发明将迹线与结构面合理分组，通过计算出迹线和结构面之间的位置关系并进行判别归并，用归并后的最终结果表征岩体结构面的产状信息。After the occurrence measurement of structural surfaces and traces is completed, in order to describe the occurrence of rock mass structural surfaces more comprehensively and accurately, the present invention groups traces and structural surfaces reasonably, and calculates the positional relationship between traces and structural surfaces And carry out discrimination and merging, and use the final result after merging to characterize the occurrence information of rock mass structural planes.

融合激光点云和数字影像的岩体结构面产状全自动数字化测量过程主要包括如下步骤：The process of fully automatic digital measurement of rock mass discontinuity plane occurrence by fusing laser point cloud and digital image mainly includes the following steps:

第一步：采集现场数据；The first step: collecting field data;

第二步：将点云数据和影像数据进行匹配、融合；Step 2: Match and fuse point cloud data and image data;

第三步：对点云数据进行三角剖分，计算出剖分后三角形法相量的参数，对法向量进行聚类分析；Step 3: Triangulate the point cloud data, calculate the parameters of the triangular normal phasor after division, and perform cluster analysis on the normal vector;

第四步：将聚类分析后的结果进行合并，对合并后的数据进行产状计算；Step 4: Merge the results after cluster analysis, and calculate the occurrence of the merged data;

第五步：对影像数据进行矫正、灰度化，然后提取出节理的轮廓线和骨架并在点云中标记出骨架；Step 5: Correct and grayscale the image data, then extract the outline and skeleton of the joint and mark the skeleton in the point cloud;

第六步：计算迹线对应点云任意三点构成三角形的内角，对迹线进行产状判别并进行产状计算；Step 6: Calculate the inner angle of the triangle formed by any three points corresponding to the point cloud of the trace, and then judge the occurrence of the trace and calculate the occurrence;

第七步：近似拟合提取出的结构面并计算结构面与各迹线的距离，然后按照结构面和迹线的位置关系将结构面和迹线进行分组；Step 7: Approximately fit the extracted structural surface and calculate the distance between the structural surface and each trace, and then group the structural surface and trace according to the positional relationship between the structural surface and the trace;

第八步：测量工作完成。Step 8: The measurement work is completed.

有益效果：本发明主要针对野外岩体结构面的数量多、分布广和随机性强等特征，传统岩体结构面产状测量方法工作效率低下，甚至遇到人迹罕至陡崖险坡根本无法实施。本发明借助于地面激光激光的精准空间测量技术和数字影像的计算机视觉识别技术，提出了一种全自动、非接触式的岩体结构面产状数字化测量方法。Beneficial effects: the present invention is mainly aimed at the characteristics of large number, wide distribution and strong randomness of rock mass structural planes in the field. The traditional method for measuring the occurrence of rock mass structural planes is inefficient, and even cannot be implemented when encountering inaccessible steep cliffs and dangerous slopes. The present invention proposes a fully automatic, non-contact digital measurement method for the occurrence of rock mass structural planes by means of the precise spatial measurement technology of ground lasers and the computer vision recognition technology of digital images.

在激光点云数据处理的过程中，本发明采用基于密度的法向量聚类算法以自动化的形式提取出主要的结构面产状信息，并通过TIN中三角形的连通性判别将结构面进行分组。在此过程中，没有现有的K值聚类、模糊聚类等需要人工假定聚类结果组数的缺点，数据处理结果不受人工干预的影响，因此能自动化识别出各种形状的岩体结构面。In the process of laser point cloud data processing, the present invention uses a density-based normal vector clustering algorithm to automatically extract the occurrence information of the main structural planes, and groups the structural planes through the connectivity of triangles in the TIN. In this process, there is no disadvantage of manually assuming the number of clustering result groups in the existing K-value clustering, fuzzy clustering, etc., and the data processing results are not affected by manual intervention, so rock masses of various shapes can be automatically identified structural surface.

在数字影像数据处理的过程中，本发明采用混合全局和局部阈值法提取出迹线的轮廓线，并通过骨架提取、迹线标记提取出迹线的产状信息。经实践检验，混合全局和局部阈值法在岩体节理迹线的轮廓线提取方面能够达到较高精度，骨架提取所采用判定条件也完全适用于迹线处理，而随后的迹线标记更能同时发挥数字影像数据信息丰富与激光点云数据精度高的优点。In the process of digital image data processing, the present invention extracts the outline of the trace by using the mixed global and local threshold method, and extracts the occurrence information of the trace by extracting the skeleton and marking the trace. It has been verified by practice that the mixed global and local threshold method can achieve high accuracy in the contour extraction of rock mass joint traces, and the judgment conditions used in skeleton extraction are also fully applicable to trace processing, and the subsequent trace marks can be more accurate at the same time. Take advantage of the rich information of digital image data and the high precision of laser point cloud data.

综上所述，本发明在产状测量过程中无人工干预，采集数据后结合一系列数字处理方法能获取高精度的岩体结构面产状信息，是一种全自动、非接触式的岩体结构面产状数字化测量方法。In summary, the present invention has no manual intervention in the process of occurrence measurement, and can obtain high-precision rock mass discontinuity information by combining a series of digital processing methods after collecting data. It is a fully automatic, non-contact rock mass A digital measurement method for the occurrence of body structural planes.

附图说明Description of drawings

图1为本测量方法的技术方案流程图。Figure 1 is a flow chart of the technical solution of the measurement method.

具体实施方式Detailed ways

本发明的融合激光点云和数字影像的岩体结构面产状全自动数字化测量方法，主要包括现场数据采集、数据融合与配准、点云数据特征提取、点云部分的产状信息获取、影像特征提取、影像部分的产状信息获取、激光点云与数字影像融合，具体步骤如下：The fully automatic digital measurement method for the occurrence of rock mass discontinuities fused with laser point clouds and digital images of the present invention mainly includes on-site data collection, data fusion and registration, feature extraction of point cloud data, acquisition of occurrence information of point cloud parts, The image feature extraction, the acquisition of the occurrence information of the image part, and the fusion of laser point cloud and digital image, the specific steps are as follows:

第一步：现场数据采集Step 1: Field Data Collection

1)仪器现场安装1) On-site installation of the instrument

根据测量区域位置与大小、岩体结构特征以及复杂工程环境对目标对象可能造成光线遮挡与反射的影响，确定TLS三维激光扫描仪架设测点的站点数量、测站点具体位置以及仪器架设的高度、扫描角度，并且按照上述要求组装和架设TLS三维激光扫描仪，若需要对不同站点的扫描数据进行拼接，还需要确定合理的控制点布设位置，以便对不同站点的点云数据进行配准；According to the position and size of the measurement area, the structural characteristics of the rock mass, and the possible impact of light occlusion and reflection on the target object by the complex engineering environment, determine the number of stations for the TLS 3D laser scanner to set up the measurement points, the specific location of the measurement points, and the height of the instrument. Scanning angle, and assemble and erect the TLS 3D laser scanner according to the above requirements. If it is necessary to splicing the scanned data of different sites, it is also necessary to determine a reasonable control point layout position in order to register the point cloud data of different sites;

2)现场扫描测量2) On-site scanning measurement

根据拟测量区的岩体结构完整性，选择TLS三维激光扫描仪的合适扫描分辨率、校正参数、扫描角度。根据仪器架设位置距测区的距离与光线，确定高分辨率数码相机的曝光度、光圈、快门时间、焦距等参数，以期保证能够获得高质量数字图像数据；According to the structural integrity of the rock mass in the area to be measured, select the appropriate scanning resolution, calibration parameters, and scanning angle of the TLS 3D laser scanner. Determine the exposure, aperture, shutter time, focal length and other parameters of the high-resolution digital camera according to the distance and light from the instrument installation position to the measurement area, in order to ensure that high-quality digital image data can be obtained;

设置好参数后，用扫描仪对需要扫描的区域进行激光扫描，从而获取点云数据，目标物体距离扫描仪越远，目标物体表面扫描点的间隔就越大，根据需要适当地调整扫描精度；After setting the parameters, use the scanner to scan the area to be scanned by laser to obtain point cloud data. The farther the target object is from the scanner, the larger the interval between the scanning points on the surface of the target object, and adjust the scanning accuracy as needed;

在激光扫描的同时，利用置于扫描仪顶部的高分辨率数码相机进行岩体结构特征的图像采集，从而获取图像数据；At the same time as laser scanning, a high-resolution digital camera placed on the top of the scanner is used to collect images of rock mass structural features to obtain image data;

第二步：数据配准与融合Step 2: Data Registration and Fusion

1)相机内部参数获取1) Acquisition of camera internal parameters

三维激光扫描仪的自带数码相机一般均为固定焦距相机，其焦距可以直接从镜头标识中获得，数码相机的内部参数只与相机内部结构和镜头有关，可由相机及镜头的出厂说明书与技术规格表中获得；The built-in digital camera of the 3D laser scanner is generally a fixed-focus camera, and its focal length can be obtained directly from the lens mark. The internal parameters of the digital camera are only related to the internal structure of the camera and the lens, which can be determined by the factory manual and technical specifications of the camera and lens. obtained from the table;

2)相机外部参数获取2) Acquisition of camera external parameters

利用三维扫描仪系统采集数据时，相机是随着扫描仪转动而转动的，因此在进行相机外部参数确定时需要考虑相机相对于初始坐标系的旋转矩阵，扫描仪系统会涉及到四个坐标系：相机坐标系CMCS、扫描仪坐标系SOCS、工程项目坐标系PRCS、世界坐标系GLCS；三维扫描仪厂商提供了坐标旋转、平移相关矩阵参数，主要有①Mounting矩阵，相机坐标系与扫描仪坐标系之间的坐标转换参数矩阵；②COP矩阵，相机拍摄瞬间相机坐标系和初始相机坐标系的旋转矩阵；③SOP矩阵，各扫描站点坐标系相当于工程坐标系的旋转平移矩阵；When using a 3D scanner system to collect data, the camera rotates as the scanner rotates. Therefore, when determining the external parameters of the camera, it is necessary to consider the rotation matrix of the camera relative to the initial coordinate system. The scanner system involves four coordinate systems. : Camera coordinate system CMCS, scanner coordinate system SOCS, engineering project coordinate system PRCS, world coordinate system GLCS; 3D scanner manufacturers provide coordinate rotation and translation related matrix parameters, mainly including ①Mounting matrix, camera coordinate system and scanner coordinate system The coordinate transformation parameter matrix between them; ②COP matrix, the rotation matrix of the camera coordinate system and the initial camera coordinate system at the moment of camera shooting; ③SOP matrix, the rotation and translation matrix of each scanning station coordinate system equivalent to the engineering coordinate system;

3)坐标系统的相互转换3) Mutual conversion of coordinate systems

①世界坐标系与相机坐标系① World coordinate system and camera coordinate system

世界坐标P_w(X_w,Y_w,Z_w)，相机坐标P_u(X_u,Y_u,Z_u)，两者转换关系：World coordinates P _w (X _w , Y _w , Z _w ), camera coordinates P _u (X _u , Y _u , Z _u ), the conversion relationship between them:

$[\begin{matrix} {X x}_{u u} \\ {Y Y}_{u u} \\ {Z Z}_{u u} \\ 11 \end{matrix}] = = [\begin{matrix} R R & t t \\ 00^{T T} & 11 \end{matrix}] [\begin{matrix} {X x}_{w w} \\ {Y Y}_{w w} \\ {Z Z}_{w w} \\ 11 \end{matrix}] = = [\begin{matrix} {r r}_{1111} & {r r}_{1212} & {r r}_{1313} & {t t}_{11} \\ {r r}_{21 twenty one} & {r r}_{22 twenty two} & {r r}_{23 twenty three} & {t t}_{22} \\ {r r}_{3131} & {r r}_{3232} & {r r}_{3333} & {t t}_{33} \\ 00 & 00 & 00 & 11 \end{matrix}] [\begin{matrix} {X x}_{w w} \\ {Y Y}_{w w} \\ {Z Z}_{w w} \\ 11 \end{matrix}] = = M m \cdot &Center Dot; {COP COP}^{- - 11} \cdot &Center Dot; {SOP SOPs}^{- - 11} [\begin{matrix} {X x}_{w w} \\ {Y Y}_{w w} \\ {Z Z}_{w w} \\ 11 \end{matrix}] - - - - - - ((22 - - 11))$

其中，旋转矩阵t是世界坐标系原点在相机坐标系中的位置坐标，旋转矩阵R为正交旋转矩阵，满足：Among them, the rotation matrix t is the position coordinate of the origin of the world coordinate system in the camera coordinate system, and the rotation matrix R is an orthogonal rotation matrix, which satisfies:

$\{\begin{matrix} {r r}_{1111}^{22} + + {r r}_{1212}^{22} + + {r r}_{1313}^{22} = = 11 \\ {r r}_{21 twenty one}^{22} + + {r r}_{22 twenty two}^{22} + + {r r}_{23 twenty three}^{22} = = 11 \\ {r r}_{3131}^{22} + + {r r}_{3232}^{22} + + {r r}_{3333}^{22} = = 11 \end{matrix} - - - - - - ((22 - - 22))$

②相机坐标系与图像坐标系②Camera coordinate system and image coordinate system

相机坐标P_u(X_u,Y_u,Z_u)，图像坐标为(x,y)，两者转换关系：The camera coordinates P _u (X _u , Y _u , Z _u ), the image coordinates are (x, y), and the conversion relationship between them is:

${Z Z}_{u u} [\begin{matrix} x x \\ y the y \\ 11 \end{matrix}] = = [\begin{matrix} f f & 00 & 00 & 00 \\ 00 & f f & 00 & 00 \\ 00 & 00 & 11 & 00 \end{matrix}] [\begin{matrix} {X x}_{u u} \\ {Y Y}_{u u} \\ {Z Z}_{u u} \\ 11 \end{matrix}] - - - - - - ((22 - - 33))$

其中，f为相机焦距。Among them, f is the focal length of the camera.

③世界坐标系与图像坐标系③World coordinate system and image coordinate system

世界坐标P_w(X_w，Y_w，Z_w)，图像坐标为(x，y)，两者转换关系：The world coordinate P _w (X _w , Y _w , Z _w ), the image coordinate is (x, y), the conversion relationship between the two:

${Z Z}_{u u} [\begin{matrix} x x \\ y the y \\ 11 \end{matrix}] = = [\begin{matrix} f f & 00 & 00 & 00 \\ 00 & f f & 00 & 00 \\ 00 & 00 & 11 & 00 \end{matrix}] [\begin{matrix} {X x}_{u u} \\ {Y Y}_{u u} \\ {Z Z}_{u u} \\ 11 \end{matrix}] = = [\begin{matrix} f f & 00 & 00 & 00 \\ 00 & f f & 00 & 00 \\ 00 & 00 & 11 & 00 \end{matrix}] M m \cdot \cdot {COP COP}^{- - 11} \cdot \cdot {SOP SOPs}^{- - 11} [\begin{matrix} {X x}_{w w} \\ {Y Y}_{w w} \\ {Z Z}_{w w} \\ 11 \end{matrix}] - - - - - - ((22 - - 44))$

第三步：点云特征信息提取The third step: point cloud feature information extraction

1)点云数据去噪1) Point cloud data denoising

本发明采用一般的中值滤波法对点云数据进行去噪处理，中值滤波法的原理就是用点云目标点的邻域内诸点坐标的中间值来替代点云目标点的坐标值，The present invention adopts general median filtering method to carry out denoising processing to point cloud data, and the principle of median filtering method is exactly to replace the coordinate value of point cloud target point with the intermediate value of all point coordinates in the neighborhood of point cloud target point,

2)三角剖分2) Triangulation

三角剖分是拓扑学中最基本的方法，它能将零散的点云剖分为无数曲边三角形，在此将用Delaunay三角剖分法进行三角剖分，运用Matlab中的delaunay函数实现三维点云的Delaunay剖分，delaunay函数输出的数据是完成划分后各个三角形的编号和顶点坐标，delaunay函数在Matlab中的用法如下：Triangulation is the most basic method in topology. It can divide the scattered point cloud into countless curved triangles. Here we will use the Delaunay triangulation method for triangulation, and use the delaunay function in Matlab to realize the three-dimensional point The Delaunay subdivision of the cloud, the data output by the delaunay function is the number and vertex coordinates of each triangle after the division is completed. The usage of the delaunay function in Matlab is as follows:

Ⅰ、输入所有点的x,y,z坐标；Ⅰ. Input the x, y, z coordinates of all points;

Ⅱ、运用delaunay函数进行三角划分，得到各三角形的ID编号和顶点坐标：Ⅱ. Use the delaunay function to perform triangulation, and obtain the ID number and vertex coordinates of each triangle:

Tri＝delaunay(x,y)；Tri=delaunay(x,y);

Ⅲ、调用z坐标，运用trimesh或者trisurf指令画出三维曲面图。Ⅲ. Call the z coordinate, and use the trimesh or trisurf command to draw a three-dimensional surface map.

3)三角形的法向向量计算3) Calculation of the normal vector of the triangle

利用剖分后三角形的三个顶点的坐标，计算出三角形的法向量坐标和三角形的形心坐标，并进一步计算出法向量与坐标Z轴、X-Z面的夹角φ_T、θ_T,形心距坐标原点的距离r_T，Using the coordinates of the three vertices of the triangle after division, calculate the coordinates of the normal vector of the triangle and the coordinates of the centroid of the triangle, and further calculate the angles φ _T , θ _T between the normal vector and the coordinate Z axis, XZ plane, and the centroid the distance r _T from the coordinate origin,

设Matlab中算得的某三角形的三个顶点为Let the three vertices of a triangle calculated in Matlab be

A(x_T1,y_T1,z_T1),B(x_T2,y_T2,z_T2)，C(x_T3,y_T3,z_T3)，则A(x _T1 ,y _T1 ,z _T1 ), B(x _T2 ,y _T2 ,z _T2 ), C(x _T3 ,y _T3 ,z _T3 ), then

$\overset{&RightArrow; &Right Arrow;}{A A B B} = = (({x x}_{T T 22} - - {x x}_{T T 11},, {y the y}_{T T 22} - - {y the y}_{T T 11},, {z z}_{T T 22} - - {z z}_{T T 11})),, \overset{&RightArrow; &Right Arrow;}{A A C C} = = (({x x}_{T T 33} - - {x x}_{T T 11},, {y the y}_{T T 33} - - {y the y}_{T T 11},, {z z}_{T T 33} - - {z z}_{T T 11})),,$

三角形的形心坐标即为则形心O距坐标原点的距离 $r_{T} = \sqrt{{(Σ_{k = 1}^{3} x_{T k} / 3)}^{2} + {(Σ_{k = 1}^{3} y_{T k} / 3)}^{2} + {(Σ_{k = 1}^{3} z_{T k} / 3)}^{2}} - - - (3 - 1)$ The coordinates of the centroid of the triangle are Then the distance from the centroid O to the origin of coordinates $r_{T} = \sqrt{{(Σ_{k = 1}^{3} x_{T k} / 3)}^{2} + {(Σ_{k = 1}^{3} {the y}_{T k} / 3)}^{2} + {(Σ_{k = 1}^{3} z_{T k} / 3)}^{2}} - - - (3 - 1)$

设ΔABC的法向量为通过构建如下方程组可得求法向量坐标：Let the normal vector of ΔABC be The coordinates of the normal vector can be obtained by constructing the following equations:

$\{\begin{matrix} (({x x}_{T T 22} - - {x x}_{T T 11})) {a a}_{T T} + + (({y the y}_{T T 22} - - {y the y}_{T T 11})) {b b}_{T T} + + (({z z}_{T T 22} - - {z z}_{T T 11})) {c c}_{T T} = = 00 \\ (({x x}_{T T 22} - - {x x}_{T T 33})) {a a}_{T T} + + (({y the y}_{T T 22} - - {y the y}_{T T 33})) {b b}_{T T} + + (({z z}_{T T 22} - - {z z}_{T T 33})) {c c}_{T T} = = 00 \\ (({x x}_{T T 33} - - {x x}_{T T 11})) {a a}_{T T} + + (({y the y}_{T T 33} - - {y the y}_{T T 11})) {b b}_{T T} + + (({z z}_{T T 33} - - {z z}_{T T 11})) {c c}_{T T} = = 00 \end{matrix} - - - - - - ((33 - - 22))$

随后利用三角函数关系求得法向量的φ_T、θ_T值：Then use the trigonometric function relationship to find the normal vector φ _T , θ _T values of :

${θ θ}_{T T} = = arccos arccos \frac{{c c}_{T T}}{\sqrt{{a a}_{T T}^{22} + + {b b}_{T T}^{22} + + {c c}_{T T}^{22}}}$

第四步：产状信息获取即点云部分Step 4: Occurrence information acquisition is the point cloud part

1)法向量聚类分析1) Normal vector cluster analysis

①数据中心化①Data centralization

数据中心化是聚类分析前的必要准备工作，本发明中数据中心化的方法，是用单个三角形的φ_T、θ_T、r_T值减去所有三角形φ_T、θ_T、r_T值的平均值；Data centralization is a necessary preparatory work before cluster analysis. The method of data centralization in the present invention is to subtract all triangle φ _T , θ _T , r _T values from the values of φ _T , θ _T , r _T of a single triangle. average value;

点云数据总共被剖分为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：The point cloud data is divided into n triangles in total, then the coordinates (a _Ti ,b _Ti ,1) (i=1,2,3...n) and φ _T , θ _T , r _T values φ _Ti , θ _Ti , r _Ti (i=1,2,3....n); use the n normal vectors φ _T , θ _T , r _T values as follows n×3 The order matrix T is expressed, and the ID is preserved during data processing:

$T = (x_{i j}) = [\begin{matrix} x_{11} & x_{12} & x_{13} \\ x_{21} & x_{22} & x_{23} \\ \begin{matrix} \cdot \\ \cdot \\ \cdot \end{matrix} & \begin{matrix} \cdot \\ \cdot \\ \cdot \end{matrix} & \begin{matrix} \cdot \\ \cdot \\ \cdot \end{matrix} \\ x_{n 1} & x_{n 2} & x_{n 3} \end{matrix}],$ φ_Ti＝x_i1,θ_Ti＝x_i2，r_T＝(x_i3),(i＝1,2,…n),(j＝1,2,3) $T = (x_{i j}) = [\begin{matrix} x_{11} & x_{12} & x_{13} \\ x_{twenty one} & x_{twenty two} & x_{twenty three} \\ \begin{matrix} \cdot \\ &Center Dot; \\ \cdot \end{matrix} & \begin{matrix} &Center Dot; \\ &Center Dot; \\ \cdot \end{matrix} & \begin{matrix} \cdot \\ &Center Dot; \\ &Center Dot; \end{matrix} \\ x_{no 1} & x_{no 2} & x_{no 3} \end{matrix}],$ φ _Ti =x _i1 , θ _Ti =x _i2 , r _T =(x _i3 ),(i=1,2,...n),(j=1,2,3)

记第j列的平均值为：Record the average value of column j as:

$\overset{&OverBar; &OverBar;}{{x x}_{j j}} = = \frac{11}{n no} {Σ Σ}_{i i = = 11}^{n no} {x x}_{i i j j} ((i i = = 11,, 22,, ... ... n no)),, ((j j = = 11,, 22,, 33)) - - - - - - ((44 - - 11))$

那么对第j列的n个数据对象所实施的中心化变换为：Then the centralization transformation implemented on the n data objects in the jth column is:

${x x}^{,,}_{i i j j} = = {x x}_{i i j j} - - {\overset{&OverBar; &OverBar;}{x x}}_{j j},, ((i i = = 11,, 22 ... ... n no)),, ((j j = = 11,, 22,, 33)) - - - - - - ((44 - - 22))$

经过上式的变换后，每列变量的均值为0，即每列变量的取值都有相同的基点，然后将中心变换后的结果重新生成为新的矩阵T,:After the transformation of the above formula, the mean value of each column variable is 0, that is, the value of each column variable has the same base point, and then the result of the center transformation is regenerated into a new matrix T,:

${T T}^{,,} = = (({x x}^{,,}_{i i j j})) = = [\begin{matrix} {x x}^{,,}_{1111} & {x x}^{,,}_{1212} & {x x}^{,,}_{1313} \\ {x x}^{,,}_{21 twenty one} & {x x}^{,,}_{22 twenty two} & {x x}^{,,}_{23 twenty three} \\ \begin{matrix} \cdot &Center Dot; \\ \cdot \cdot \\ \cdot \cdot \end{matrix} & \begin{matrix} \cdot \cdot \\ \cdot &Center Dot; \\ \cdot \cdot \end{matrix} & \begin{matrix} \cdot &Center Dot; \\ \cdot &Center Dot; \\ \cdot &Center Dot; \end{matrix} \\ {x x}^{,,}_{n no 11} & {x x}^{,,}_{n no 22} & {x x}^{,,}_{n no 33} \end{matrix}],, ((i i = = 11,, 22,, ... ... n no)),, ((j j = = 11,, 22,, 33)) . .$

②数据标准化②Data standardization

在此采用最大最小归一法对数据进行标准化处理，即将矩阵中的某一数据除以该列数据绝对值的最大值，具体见下式：Here, the maximum and minimum normalization method is used to standardize the data, that is, to divide a certain data in the matrix by the maximum value of the absolute value of the column data, see the following formula for details:

${x x}_{i i j j}^{,,,,} = = \frac{{x x}_{i i j j}^{,,}}{{max max}_{i i} | | {x x}_{i i j j}^{,,} | |} - - - - - - ((44 - - 33))$

最终生成新的矩阵：Finally a new matrix is generated:

${T T}^{,,,,} = = (({x x}^{,,,,}_{i i j j})) = = [\begin{matrix} {x x}^{,,,,}_{1111} & {x x}^{,,,,}_{1212} & {x x}^{,,,,}_{1313} \\ {x x}^{,,,,}_{21 twenty one} & {x x}^{,,,,}_{22 twenty two} & {x x}^{,,,,}_{23 twenty three} \\ \begin{matrix} \cdot &Center Dot; \\ \cdot &Center Dot; \\ \cdot &Center Dot; \end{matrix} & \begin{matrix} \cdot &Center Dot; \\ \cdot &Center Dot; \\ \cdot &Center Dot; \end{matrix} & \begin{matrix} \cdot &Center Dot; \\ \cdot &Center Dot; \\ \cdot &Center Dot; \end{matrix} \\ {x x}^{,,,,}_{n no 11} & {x x}^{,,,,}_{n no 22} & {x x}^{,,,,}_{n no 33} \end{matrix}]$

③特征生成③ Feature generation

上述矩阵T”属于低维矩阵，为描述简便、形象，将T”的各个行向量投影于三维直角坐标系中形成n个点，并以(x”_i1,x”_i2,x”_i3)作为这n个点的坐标，记为(x_i,y_i,z_i)；如果有坐标重复的点即T”内的行向量数值相同，因为其ID不同，予以保留；The above-mentioned matrix T" is a low-dimensional matrix. For the convenience of description and visualization, each row vector of T" is projected into a three-dimensional Cartesian coordinate system to form n points, and (x" _i1 , x" _i2 , x" _i3 ) are used as The coordinates of these n points are recorded as ( _xi , y _i , _zi ); if there are points with repeated coordinates, that is, the row vectors in T” have the same value, because their IDs are different, they are reserved;

④法向量的聚类分析④ Cluster analysis of normal vectors

DBSCAN聚类分析法是一种基于高密度联通区域的聚类算法，本方法运用并修改了该聚类算法，从而达到提取结构面信息的目的；The DBSCAN clustering analysis method is a clustering algorithm based on high-density connected areas. This method uses and modifies the clustering algorithm to achieve the purpose of extracting structural plane information;

首先，取ε＝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₁的ε领域中，First, take ε=3, MinPts=5, two points q ₁ (x _q1 , y _q1 , z _q2 ), q ₂ (x _q2 , y _q2 , z _q2 ) in T”, where the parameter ε is the domain radius; A circle with a point as the center and ε as the radius is called the ε field of the point; the parameter MinPts is the minimum number of core objects; if the number of points contained in the ε field of a point is greater than or equal to MinPts, the point is called the core object ; The values of parameters ε and MinPts need to be adjusted in specific projects; then, calculate the distance d from point q ₁ to point q ₂ , if d≤3, then assign point q ₂ to the ε domain of point q ₁ ,

$d d = = \sqrt{{(({x x}_{q q 11} - - {x x}_{q q 22}))}^{22} + + {(({y the y}_{q q 11} - - {y the y}_{q q 22}))}^{22} + + {(({z z}_{q q 11} - - {z z}_{q q 22}))}^{22}} - - - - - - ((44 - - 44))$

之后，输入T”中的所有点的坐标，计算并找出能够归于点q₁的ε领域内的所有点；Afterwards, input the coordinates of all points in T", calculate and find out all points in the ε field that can be attributed to point _q1 ;

假设得到q₁的ε领域内包含的m个点为{q₂,q₃,…q_m+1}，则q₁的ε领域内包含的点数为m；Assuming that the m points contained in the ε domain of q ₁ are {q ₂ , q ₃ ,…q _m+1 }, then the number of points contained in the ε domain of q ₁ is m;

如果m≥5，那么将q₁归为核心对象集合A₁内的点；如果0＜m＜5，则将q₁归为边缘对象集合B₁内的点；如果m＝0,则将q₁归为噪音；If m≥5, then classify q ₁ as a point in core object set A ₁ ; if 0<m<5, classify q ₁ as a point in edge object set B ₁ ; if m=0, classify q ₁ is classified as noise;

再对q₁的领域内的其他点{q₂,q₃,…q_m+1}按如上方法进行计算并判别是否为核心对象，将点q₁的领域内所有被归为核心对象的点归于集合A₁，所有被归为边缘对象的点归于集合B₁，删除被归为噪音的点；Then calculate other points {q ₂ , q ₃ ,…q _m+1 } in the field of q ₁ according to the above method and judge whether they are core objects, and classify all points in the field of point q ₁ that are classified as core objects Belong to the set A ₁ , all points classified as edge objects belong to the set B ₁ , delete the points classified as noise;

重复输入T”内的所有点坐标，计算出每个点的核心对象集合{A₁,A₂…A_n}、边缘对象集合{B₁,B₂…B_n}，其中噪音的核心对象集合与边缘对象集合都记为0，再将这两个大集合按如下方法分别进行合并；Repeatedly input the coordinates of all points in T”, and calculate the core object set {A ₁ ,A ₂ …A _n }, the edge object set {B ₁ ,B ₂ …B _n } of each point, and the noise core object set and the edge object collection are recorded as 0, and then these two large collections are merged as follows;

取{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’}；Take two sets A ₁ , A ₂ in {A ₁ ,A ₂ …A _n }, if the set That is, if there are duplicate points, the samples in A ₂ that do not intersect with A ₁ will be merged into A ₁ to form cluster A ₁ ', and A ₁ ' will be used as a new merging object to participate in the next merging judgment; if the set That is, if there is no duplicate point, assign A ₁ to cluster A ₁ ', assign A ₂ to cluster A ₂ ', and use A ₁ ' and A ₂ ' as new merge objects to participate in the next merge judgment, and then take { In A ₁ , A ₂ ... A _n }, A ₃ and A ₁ ' and A ₂ ' are respectively subjected to the above combined discrimination, if and Then assign A ₃ to cluster A ₃ ', and use A ₃ ' as a new merge object to participate in the next merge judgment. Repeat this until all the samples in {A ₁ ,A ₂ ...A _n } are completed and discriminated, and finally a new set {A ₁ ',A ₂ '...A _v '} is obtained;

边缘对象集合{B₁,B₂…B_n}的合并也按上述原则进行合并，得到{B₁’,B₂’…B_w’}。The merging of the edge object set {B ₁ , B ₂ ...B _n } is also merged according to the above principles to obtain {B ₁ ', B ₂ '...B _w '}.

{A₁’,A₂’…A_v’}和{B₁’,B₂’…B_w’}则是聚类的最终结果，再利用聚类结果中各元素的ID编号就能提取出对应法向量的φ_t、θ_t、r值，也能标记出对应的三角面片；{A ₁ ',A ₂ '…A _v '} and {B ₁ ',B ₂ '…B _w '} are the final results of clustering, and then use the ID number of each element in the clustering results to extract Corresponding to the φ _t , θ _t , and r values of the normal vector, the corresponding triangular surface can also be marked;

2)聚类结果处理2) Clustering result processing

①三角面片合并①Merge triangle faces

该步将除无效三角面片，弥补DBSCAN算法的不足，并将聚类后的结果提取、标示，通过合并三角面构建出岩体结构面。In this step, invalid triangular faces are removed to make up for the shortcomings of the DBSCAN algorithm, and the clustered results are extracted and marked, and the rock mass structural plane is constructed by merging the triangular faces.

Ⅰ、通过DBSCAN的聚类结果利用ID标记出各个簇中对应的三角面片，删除重复的三角面片，并对各个簇进行标号。Ⅰ. Through the clustering results of DBSCAN, use ID to mark the corresponding triangles in each cluster, delete repeated triangles, and label each cluster.

Ⅱ、在簇A_h’(A_h’为{A₁’,A₂’…A_v’}中任意一样本)中任取一三角面片i，若其邻近存在同属于A_h’的其他三角面片，则将i保留；否则，delete。Ⅱ. In the cluster A _h ' (A _h ' is any sample in {A ₁ ', A ₂ '…A _v '}), randomly select a triangular facet i, if there is other triangular facet i that belongs to A _h ' nearby triangular face, then keep i; otherwise, delete.

Ⅲ、经过Ⅱ生成的新簇A_h’中，任取两三角面片ii、iii，若ii、iii能通过同属于A_h’中的任意多个三角面片相连通，则将ii、iii归入结构面X₁；否则ii归入X₁，iii归入X₂。如此，将A_h’中的三角面片进行结构面分组；Ⅲ. In the new cluster A _h ' generated by Ⅱ, randomly select two triangular faces ii and iii. If ii and iii can be connected through any number of triangular faces belonging to A _h ', then ii and iii belong to the structure plane X ₁ ; otherwise, ii belongs to X ₁ , and iii belongs to X ₂ . In this way, the triangular faces in A _h ' are grouped into structural faces;

Ⅳ、分别将{A₁’,A₂’…A_v’}中对应的三角面片进行Ⅱ、Ⅲ操作，得到X₁,X₂…X_s；Ⅳ. Perform operations II and III on the corresponding triangular faces in {A ₁ ', A ₂ '…A _v '} to obtain X ₁ , X ₂ ...X _s ;

Ⅴ、分别将{B₁’,B₂’…B_w’}中对应三角面片进行Ⅱ、Ⅲ操作，得到Y₁,Y₂…Y_r。Ⅴ. Perform operations II and III on the corresponding triangle faces in {B ₁ ', B ₂ '…B _w '} respectively to obtain Y ₁ , Y ₂ ...Y _r .

最终得到的X₁,X₂…X_s，Y₁,Y₂…Y_r就是构成岩体结构面的三角面片，将这些三角面片借由ID标示出来后即能得到岩体的各个结构面；The final obtained X ₁ , X ₂ ...X _s , Y ₁ , Y ₂ ...Y _r are the triangular faces that constitute the structural plane of the rock mass. After marking these triangular faces with ID, the various structures of the rock mass can be obtained noodle;

②剔除开挖面② Eliminate the excavation surface

实际工程中，岩体结构面的产状信息中包含了大量的开挖面信息，而开挖面并不属于岩体的结构面。由于开挖面产状具有重复性、单一性的特征，只需要用人工剔除的方法在RiscanPro中就能将其识别并将开挖面的三角面片删除；In actual engineering, the occurrence information of the rock mass structural plane contains a large amount of excavation face information, but the excavation face does not belong to the structural plane of the rock mass. Since the occurrence of the excavation surface has the characteristics of repetition and singleness, it only needs to be manually eliminated in RiscanPro to identify it and delete the triangular faces of the excavation surface;

③结构面产状计算③Occurrence calculation of structural plane

结构面的产状信息包括结构面的倾角与倾向，只要对X₁,X₂…X_s中三角面的参数进行计算，就能得到结构面的产状信息，具体算法如下：The occurrence information of the structural plane includes the inclination and tendency of the structural plane. As long as the parameters of the triangular plane in X ₁ , X ₂ ...X _s are calculated, the occurrence information of the structural plane can be obtained. The specific algorithm is as follows:

假设X₁对应H个三角面，计算X₁中三角面的法向量坐标(a_Td,b_Td,c_Td)、法向量φ_T、θ_T值(φ_Td、θ_Td)(d与X₁中对应三角面的ID相同)的算数平均值：Assuming that X ₁ corresponds to H triangular surfaces, calculate the normal vector coordinates (a _Td , b _Td , c _Td ), normal vector φ _T , θ _T values (φ _Td , θ _Td ) of the triangular surfaces in X ₁ (d and X ₁ The arithmetic mean of the IDs of the corresponding triangle faces in the same):

${a a}_{X x 11} = = \frac{11}{H h} {Σa Σa}_{T T d d}$

${b b}_{X x 11} = = \frac{11}{H h} {Σb Σb}_{T T d d} - - - - - - ((44 - - 55))$

${c c}_{X x 11} = = \frac{11}{H h} {Σc Σc}_{T T d d}$

${φ φ}_{X x 11} = = \frac{11}{H h} {Σφ Σφ}_{T T d d}$

${θ θ}_{X x 11} = = \frac{11}{H h} {Σθ Σθ}_{T T d d}$

则β_X1＝|90°-θ_X1|(0°≤β_X1≤90°)即为结构面X₁的倾角；Then β _X1 = |90°-θ _X1 |(0°≤β _X1 ≤90°) is the inclination angle _of structural plane X1;

因倾向角度范围为0°～360°，为精确计算倾向方向，需要利用X₁的法向量与Y、Z轴的余弦值来判断X₁的法向量指向的卦限。具体判别方法如下：Since the inclination angle ranges from 0° to 360°, in order to accurately calculate the inclination direction, it is necessary to use the normal vector of X ₁ and the cosine values of the Y and Z axes to determine the hexagram limit to which the normal vector of X ₁ points. The specific identification method is as follows:

假设X₁的法向量为直线L为X₁与坐标系中X-Y面的交线，向量为的水平投影向量，根据倾向方向的定义，为结构面X₁倾向；α_X1、γ_X1分别为与X、Z轴的夹角；k_X1为与Y轴的夹角。则有：Suppose the normal vector _of X1 is The straight line L is the intersection line between X ₁ and the XY plane in the coordinate system, the vector for The horizontal projection vector of , according to the definition of the dip direction, is the orientation of structural plane X ₁ ; α _X1 and γ _X1 are respectively The included angle with the X and Z axes; k _X1 is The angle with the Y axis. Then there are:

${cosk cosk}_{X x 11} = = \frac{{b b}_{X x 11}}{\sqrt{{a a}_{X x 11}^{22} + + {b b}_{X x 11}^{22}}}$

${cosα cosα}_{X x 11} = = \frac{{a a}_{X x 11}}{\sqrt{{a a}_{X x 11}^{22} + + {b b}_{X x 11}^{22} + + {c c}_{X x 11}^{22}}}$

${cosγ cosγ}_{X x 11} = = \frac{{c c}_{X x 11}}{\sqrt{{a a}_{X x 11}^{22} + + {b b}_{X x 11}^{22} + + {c c}_{X x 11}^{22}}}$

计算出的即为结构面X₁的倾向；Calculated That is the tendency of structural plane X ₁ ;

重复上述步骤，带入剩余结构面的参数，就能算出所有结构面的倾角与倾向；Repeat the above steps and bring in the parameters of the remaining structural surfaces to calculate the inclination and inclination of all structural surfaces;

④补全结构面④ Completing the structural plane

在采集、处理等一系列操作中可能造成少量点云、三角面片的缺失，影响最终展示效果，将得到的三角面片导入到Geomagic、Pointcloud等软件中进行补全，最终得到的各簇所对应的轮廓面即可视为岩体结构面；In a series of operations such as collection and processing, a small number of point clouds and triangular patches may be missing, which will affect the final display effect. Import the obtained triangular patches into Geomagic, Pointcloud and other software for completion, and finally obtain all clusters. The corresponding profile surface can be regarded as the rock mass structural surface;

第五步：影像数据特征提取Step 5: Image data feature extraction

1)图像光学矫正1) Image optical correction

根据数码相机成像数学模型，矫正参数分为数码相机的内部参数和外部参数，式5-1为数码相机成像数学模型公式，According to the digital camera imaging mathematical model, the correction parameters are divided into internal parameters and external parameters of the digital camera. Equation 5-1 is the digital camera imaging mathematical model formula,

${Z Z}_{c c} [\begin{matrix} U u \\ V V \\ 11 \end{matrix}] = = [\begin{matrix} \frac{11}{{d d}_{x x}} & 00 & {u u}_{00} \\ 00 & \frac{11}{{d d}_{y the y}} & {v v}_{00} \\ 00 & 00 & 11 \end{matrix}] [\begin{matrix} f f & 00 & 00 & 00 \\ 00 & f f & 00 & 00 \\ 00 & 00 & 11 & 00 \end{matrix}] [\begin{matrix} R R & t t \\ 00^{T T} & 11 \end{matrix}] [\begin{matrix} {X x}_{W W} \\ {Y Y}_{W W} \\ {Z Z}_{W W} \\ 11 \end{matrix}] - - - - - - ((55 - - 11))$

令：make:

${M m}_{11} = = [\begin{matrix} \frac{11}{{d d}_{x x}} & 00 & {u u}_{00} \\ 00 & \frac{11}{{d d}_{y the y}} & {v v}_{00} \\ 00 & 00 & 11 \end{matrix}] {M m}_{22} = = [\begin{matrix} f f & 00 & 00 & 00 \\ 00 & f f & 00 & 00 \\ 00 & 00 & 11 & 00 \end{matrix}] {M m}_{33} = = [\begin{matrix} R R & t t \\ 00^{T T} & 11 \end{matrix}]$

其中，Z_c为照片成像平面相对于镜头的距离；U与V为照片影像在成像平面内像素坐标系下的坐标；X_W、Y_W、Z_W为现实物体在全局世界坐标系下的坐标；M₁、M₂为数码相机的内部参数，内部参数只与相机内部结构和镜头有关，可由相机及镜头的出厂说明书与技术规格表中获得；M₃为数码相机的外部参数，外部参数可由相机标定实验经相机成像数学模型公式反演获得；Among them, Z _c is the distance between the photo imaging plane and the lens; U and V are the coordinates of the photo image in the pixel coordinate system in the imaging plane; X _W , Y _W , Z _W are the coordinates of the real object in the global world coordinate system ; M ₁ and M ₂ are the internal parameters of the digital camera, which are only related to the internal structure of the camera and the lens, and can be obtained from the factory manuals and technical specification sheets of the camera and lens; M ₃ is the external parameter of the digital camera, and the external parameters can be obtained from The camera calibration experiment is obtained through the inversion of the camera imaging mathematical model formula;

当得到矫正参数后，即可按照公式5-1对影像图片进行矫正，一般可采用成熟商业软件或者三维激光扫描仪附带成套软件来进行。After the correction parameters are obtained, the image can be corrected according to formula 5-1. Generally, mature commercial software or a complete set of software attached to the 3D laser scanner can be used for correction.

2)图像灰度化2) Image grayscale

数码相机所拍摄数字图像的初始状态是彩色图像，首先需要对其转成灰度图像。灰度图像是指只含有亮度信息的数字图像，且亮度值变化连续。将高分辨率数码图像转换为灰度图并且数字化，实际上就是将图像转化为一个灰度值矩阵F(M,N)。即表明图像大小为M×N个像素由M×N阶的矩阵表达，矩阵中每一个值表达为像素单元的灰度值。式5-2为灰度值矩阵的表达形式。The initial state of the digital image taken by the digital camera is a color image, which needs to be converted into a grayscale image first. A grayscale image refers to a digital image that only contains brightness information, and the brightness value changes continuously. Converting a high-resolution digital image into a grayscale image and digitizing it is actually converting the image into a grayscale value matrix F(M,N). That is to say, the size of the image is M×N pixels, which are expressed by a matrix of M×N order, and each value in the matrix is expressed as the gray value of the pixel unit. Equation 5-2 is the expression form of the gray value matrix.

3)轮廓线提取3) Contour extraction

轮廓线提取的方法是采用计算机图像处理方法中的图像分割法。本发明提出了一种适合岩石节理轮廓线提取的混合全局和局部阀值法。The method of contour line extraction is to adopt the image segmentation method in the computer image processing method. The invention proposes a hybrid global and local threshold value method suitable for extracting rock joint outlines.

第一步利用岩石节理区的灰度值一般为区域局部最小值的特点，采用局部阀值法，经过一系列阀值测试，可以得到局部阀值法中的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求出。The first step is to use the characteristic that the gray value of the rock joint area is generally the local minimum value, and use the local threshold method. After a series of threshold tests, the min(x, y) in the local threshold method can be obtained as ( x, y) as the local minimum gray value in the 7×7 pixel grid neighborhood, which can reduce the noise interference in the joint area; let a(x, y) be 70 with (x, y) as the center The average gray value within a grid neighborhood of ×70 pixels. First, the pixel grid points satisfying the condition min(x,y)≥a(x,y) in the original image are found through the local threshold method. In the second step, f(x,y)=0 of pixel grid points satisfying the condition forms an intermediate transition image. In the third step, the Otsu method of the global threshold method is used for the intermediate transition image to determine the segmentation threshold T _O , and then image segmentation is performed to extract rock joint contour lines. The Otsu method is based on the statistical theory to find the threshold, and the optimal threshold of the Otsu method is determined by the maximum variance between the background image and the target image. Assuming that the sum of pixels in the image is N, the total number of pixels on gray level l is n ₁ , the relationship between N and n ₁ is shown in Equation 5-3, and the probability distribution p _l of image histogram pixels conforms to Equation 5- 4. The optimal threshold T _O of the Otsu method is finally obtained from formula 5-5.

$N N = = {Σ Σ}_{l l = = 00}^{l l - - 11} {n no}_{l l} - - - - - - ((55 - - 33))$

${p p}_{l l} = = \frac{{n no}_{l l}}{N N},, {Σ Σ}_{l l = = 00}^{l l - - 11} {p p}_{l l} = = 11 - - - - - - ((55 - - 44))$

${T T}_{O o} = = \underset{00 < < l l < < 11}{argmax argmax} {{{σ σ}_{B B}^{22} ((l l))}} - - - - - - ((55 - - 55))$

式中，为类间方差。In the formula, is the between-class variance.

4)图像去噪，4) image denoising,

图像噪声对后续的在岩石节理骨架及拓扑关系提取会产生很大的干扰，所以应对岩石节理轮廓线图像去噪。Image noise will cause great interference to the subsequent extraction of rock joint skeleton and topological relationship, so the rock joint contour image should be denoised.

具体方法如下：The specific method is as follows:

Ⅰ、对岩石节理轮廓线图像进行二值化处理；Ⅰ. Binarize the rock joint contour line image;

Ⅱ、消除岩石节理轮廓线图像中的岩石节理区里的黑斑或者岩石区里的白斑。运用数字图像形态学中膨胀运算算子。原图像A被结构元素B膨胀可定义为：将结构元素B的反射平移X个像素后仍与A有交点的所有的点x组成。即运算公式为：Ⅱ. Eliminate black spots in the rock joint area or white spots in the rock area in the rock joint outline image. Use the expansion operator in digital image morphology. The dilation of the original image A by the structural element B can be defined as: after the reflection of the structural element B is translated by X pixels, it is composed of all points x that still have intersections with A. That is, the operation formula is:

本发明使用线结构元素对节理图像进行膨胀运算。The present invention uses line structure elements to perform expansion operations on joint images.

Ⅲ、膨胀处理后的节理区具有良好的连通性，运用Matlab软件中的Bwareaopen函数，根据连通区大小来过滤掉图像中不需要的小面积部分；Ⅲ. The joint area after expansion processing has good connectivity. Use the Bwareaopen function in Matlab software to filter out unnecessary small areas in the image according to the size of the connected area;

Ⅳ、在第Ⅲ步处理后节理轮廓线的边缘会有一些不规则的细小毛刺，对后续节理骨架提取产生影响，所以需要对不规则细小毛刺进行边缘光滑。本文运用中值滤波法进行光滑处理。中值滤波原理是：给定的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)记为：Ⅳ. There will be some irregular small burrs on the edge of the joint outline after the treatment in step Ⅲ, which will affect the subsequent joint skeleton extraction, so it is necessary to smooth the edges of the irregular and small burrs. In this paper, the median filtering method is used for smoothing. The principle of median filtering is: given D values {a ₁ , a ₂ ... a _D } are arranged in order according to size, when D is an odd number, the value in the middle position is called the value of the D values Median value; when D is an even number, the average value of the two values in the middle is the median value of the D values, which is recorded as med{a ₁ ,a ₂ ...a _D }, within the neighborhood window The gray median value of all pixels is used as the gray value of the filtered pixel in the window. That is: the image is a matrix of [x(I,J)] _M×N , the domain window is _AD , and the image output y(I,J) of pixel x(I,J) after median filtering is recorded as:

$y the y ((I I,, J J)) = = {med med}_{{A A}_{d d} ((I I,, J J))} [[x x ((I I,, J J))]] - - - - - - ((55 - - 77))$

本文使用5×5的正方形窗口对图像进行中值滤波，In this paper, a 5×5 square window is used to perform median filtering on the image.

5)节理迹线骨架提取5) Joint trace skeleton extraction

节理轮廓线不能直观地表达岩石节理骨架特征及其扑拓结构，需要进一步提取节理迹线骨架，提取工作分为两个步骤：①、图像细化；②、去除断枝；Joint contours cannot intuitively express the characteristics of the rock joint skeleton and its topological structure, and it is necessary to further extract the joint trace skeleton. The extraction work is divided into two steps: ①, image refinement; ②, removal of broken branches;

①、图像细化就是将二值图像中的像素点在保持原有形状与连通性的基础上进行一层一层地像素点剥离，直到残余图像骨架，图像骨架即为图像的中轴线，图像细化后图像剩余像素宽度为l；①. Image thinning is to peel off the pixels in the binary image layer by layer on the basis of maintaining the original shape and connectivity until the residual image skeleton. The image skeleton is the central axis of the image. The remaining pixel width of the image after thinning is l;

剥离条件如下：The stripping conditions are as follows:

★判断二值图中的点8个邻域内是否满足白点个数为2到6个；★Judging whether the number of white points in the 8 neighborhoods of the points in the binary image is 2 to 6;

★判断二值图中的点8个邻域内是否满足白点都连续；★Judging whether the white points in the 8 neighborhoods of the points in the binary image are continuous;

★判断二值图中的点是否满足上，左，右不全为白点；★Judging whether the points in the binary image meet the requirements that the upper, left, and right are not all white points;

★判断二值图中的点是否满足上，左，下不全为白点；★Judging whether the points in the binary image satisfy the requirement that the top, left, and bottom are not all white points;

②、去除断枝就是在图像细化后依旧存在毛刺断枝需要剔除。去除断枝需要对节理骨架线的端点与交叉点进行识别；②. Removing broken branches means that after the image is thinned, there are still burrs and broken branches that need to be removed. Removing broken branches requires identifying the endpoints and intersections of joint skeleton lines;

端点识别：在二值图像中令黑色像素点值为0；白色像素点值为1，采用3×3的正方形检测窗口迭代遍历整个图像，检查所有像素点值为1的像素点，检验其邻域中8个相邻点中像素点值为1的点个数，若个数为1，则检验得到该像素点为端点；Endpoint recognition: In the binary image, the black pixel value is 0; the white pixel value is 1, and a 3×3 square detection window is used to iterate through the entire image, check all the pixel points with a pixel value of 1, and check its neighbors. The number of points whose pixel value is 1 among the 8 adjacent points in the domain, if the number is 1, then the pixel point is checked as the endpoint;

交叉点识别：交叉点存在三线交叉或者四线叉点。第一步采用3×3的正方形检测窗口检查遍历像素点值为0的像素点邻域中8个相邻点中像素点值为1的点个数为3或4；第二步扩大检测窗口为5×5正方形，对3×3窗口最外侧的像素点值为1的点个数进行统计，若与3×3窗口内侧的像素点值为1的点个数相同，则检验得到该像素点为交叉点；Intersection recognition: There is a three-line intersection or a four-line intersection at the intersection. The first step is to use a 3×3 square detection window to check that the number of points with a pixel value of 1 among the 8 adjacent points in the neighborhood of traversing pixels with a pixel value of 0 is 3 or 4; the second step is to expand the detection window It is a 5×5 square, count the number of points whose value is 1 on the outermost pixel of the 3×3 window, and if it is the same as the number of points whose value is 1 on the inner side of the 3×3 window, then the pixel is obtained by inspection The point is the intersection point;

去除断枝：应用端点和交叉点将图像分割成各自单独的连通区域，计算每个连通区域的大小，即连通域内像素点个数N_i，设置连通域内像素点个数阈值N_t，若N_i＜N_t，就将其当作断枝，予以删除；Remove broken branches: use endpoints and intersection points to divide the image into separate connected regions, calculate the size of each connected region, that is, the number of pixels in the connected region N _i , set the threshold N _t of the number of pixels in the connected region, if N _i < N _t , it will be regarded as a broken branch and deleted;

6)节理迹线标记6) Joint trace mark

节理迹线信息统计需要对每一条岩石节理迹线进行标识，方法如下：The joint track information statistics need to identify each rock joint track, the method is as follows:

Ⅰ、搜索图像像素点集的连通性，并以不同颜色标记，利用Matlab软件的Bwlabel函数寻找图像中的每一个连通对象，并且根据寻找到的顺序用不同的整数值标注，Ⅰ. Search for the connectivity of the image pixel point set and mark it with different colors. Use the Bwlabel function of Matlab software to find each connected object in the image, and mark it with different integer values according to the order found.

Ⅱ、Bwlabel函数返回值是double类型的标记矩阵，因而无法实现图像显示，接下来利用Matlab软件的label2rgb函数，将标记矩阵中指定每一个连通区的整数值定义为颜色矩阵，以RGB值图像显示；Ⅱ. The return value of the Bwlabel function is a label matrix of double type, so the image display cannot be realized. Next, use the label2rgb function of Matlab software to define the integer value specifying each connected region in the label matrix as a color matrix, and display it as an RGB value image ;

第六步：产状信息获取即影像部分Step 6: Occurrence information acquisition is the image part

根据迹线标记与数据配准，获得每一条标记迹线起始坐标点与任一中间点的坐标，假设迹线的起始坐标点与任一中间点的坐标分别为E(x_E,y_E,z_E)，F(x_F,y_F,z_F)，G(x_G,y_G,z_G)，三点可确定一个三角形，可求出对应的三个内角θ_E，θ_F，θ_G：According to the registration of trace marks and data, the coordinates of the starting coordinate point and any intermediate point of each marked trace are obtained, assuming that the coordinates of the starting coordinate point of the trace and any intermediate point are E(x _E ,y _E ,z _E ), F(x _F ,y _F ,z _F ), G(x _G ,y _G ,z _G ), three points can determine a triangle, and the corresponding three interior angles θ _E , θ _F can be obtained , θ _G :

记：remember:

$\overset{&RightArrow; &Right Arrow;}{E E. F f} = = \overset{&RightArrow; &Right Arrow;}{m m} = = (({x x}_{F f} - - {x x}_{E E.},, {y the y}_{F f} - - {y the y}_{E E.},, {z z}_{F f} - - {z z}_{E E.})) = = (({a a}_{G G},, {b b}_{G G},, {c c}_{G G}))$

$\overset{&RightArrow; &Right Arrow;}{E E. G G} = = \overset{&RightArrow; &Right Arrow;}{n no} = = (({x x}_{G G} - - {x x}_{E E.},, {y the y}_{G G} - - {y the y}_{E E.},, {z z}_{G G} - - {z z}_{E E.})) = = (({a a}_{F f},, {b b}_{F f},, {c c}_{F f}))$

$\overset{&RightArrow; &Right Arrow;}{F f G G} = = \overset{&RightArrow; &Right Arrow;}{p p} = = (({x x}_{G G} - - {x x}_{F f},, {y the y}_{G G} - - {y the y}_{F f},, {z z}_{G G} - - {z z}_{F f})) = = (({a a}_{E E.},, {b b}_{E E.},, {c c}_{E E.}))$

由 ${cosθ}_{E} = \frac{\overset{&RightArrow;}{m} \cdot \overset{&RightArrow;}{n}}{| \overset{&RightArrow;}{m} | | \overset{&RightArrow;}{n} |}$ 得：Depend on ${cosθ}_{E.} = \frac{\overset{&Right Arrow;}{m} &Center Dot; \overset{&Right Arrow;}{no}}{| \overset{&Right Arrow;}{m} | | \overset{&Right Arrow;}{no} |}$ have to:

${θ θ}_{E E.} = = arccos arccos \frac{\overset{&RightArrow; &Right Arrow;}{m m} \cdot &Center Dot; \overset{&RightArrow; &Right Arrow;}{n no}}{| | \overset{&RightArrow; &Right Arrow;}{m m} | | | | \overset{&RightArrow; &Right Arrow;}{n no} | |} = = arccos arccos \frac{{a a}_{G G} {a a}_{F f} + + {b b}_{G G} {b b}_{F f} + + {c c}_{G G} {c c}_{F f}}{\sqrt{{a a}_{G G}^{22} + + {b b}_{G G}^{22} + + {c c}_{G G}^{22}} \sqrt{{a a}_{F f}^{22} + + {b b}_{F f}^{22} + + {c c}_{F f}^{22}}} - - - - - - ((66 - - 11))$

再运用公式6-1同理可求出θ_F，θ_G Then use the formula 6-1 to obtain θ _F , θ _G in the same way

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

①如果θ_J≤120°，则此迹线上三点所构三角形可视为出露面，运用式6-2，即可求此出露面法向量 ① If θ _J ≤ 120°, the triangle formed by the three points on this trace can be regarded as the exposed surface, and the normal vector of the exposed surface can be obtained by using Equation 6-2

$\{\begin{matrix} {a a}_{E E.} {a a}_{J J} + + {b b}_{E E.} {b b}_{J J} + + {c c}_{E E.} {c c}_{J J} = = 00 \\ {a a}_{F f} {a a}_{J J} + + {b b}_{F f} {b b}_{J J} + + {c c}_{F f} {c c}_{J J} = = 00 \\ {a a}_{G G} {a a}_{J J} + + {b b}_{G G} {b b}_{J J} + + {c c}_{G G} {c c}_{J J} = = 00 \end{matrix} - - - - - - ((66 - - 22))$

法向量与X-Y平面法向量{0，0}的交角，即出露面的倾角，表示为normal vector The intersection angle with the XY plane normal vector {0, 0}, that is, the inclination angle of the exposed surface, is expressed as

因倾向角度范围为0°～360°，为精确计算倾向方向，需要利用法向向量与Y，Z轴的余弦值来判断结构面的法向向量指向的卦限；具体判别方法如下：假设平面P为一结构面，直线L为该结构面与水平面的交线，向量为平面P的法向向量的水平投影向量，根据倾向方向的定义，为结构面的倾向；α_J，γ_J分别为法向向量与X，Z轴的角；k_J为向量与Y轴的夹角。则有Since the inclination angle ranges from 0° to 360°, in order to accurately calculate the inclination direction, it is necessary to use the normal vector and the cosine value of the Y and Z axes to determine the hexagram limit pointed by the normal vector of the structural surface; the specific identification method is as follows: Assume that the plane P is a structural surface, and the straight line L is the intersection line between the structural surface and the horizontal plane, and the vector is the normal vector of the plane P The horizontal projection vector of , according to the definition of the dip direction, is the tendency of the structural surface; α _J , γ _J are the normal vectors Angle with X, Z axis; k _J is a vector The angle with the Y axis. then there is

${cosk cosk}_{J J} = = \frac{{b b}_{J J}}{\sqrt{{a a}_{J J}^{22} + + {b b}_{J J}^{22}}}$

${cosα cosα}_{J J} = = \frac{{a a}_{J J}}{\sqrt{{a a}_{J J}^{22} + + {b b}_{J J}^{22} + + {c c}_{J J}^{22}}}$

${cosγ cosγ}_{J J} = = \frac{{c c}_{J J}}{\sqrt{{a a}_{J J}^{22} + + {b b}_{J J}^{22} + + {c c}_{J J}^{22}}}$

②如果θ＞120°，则将此迹线视为普通迹线；②If θ>120°, this trace is regarded as a normal trace;

利用E、G点坐标即可求出线EG方程：z＝Qx+Py，任取线上一点(a_L,b_L,c_L)Use the coordinates of points E and G to obtain the line EG equation: z=Qx+Py, take any point on the line (a _L , b _L , c _L )

线EG与坐标平面X-Y的夹角，即为迹线倾角，表示为The angle between the line EG and the coordinate plane X-Y is the inclination angle of the trace, expressed as

因倾向角度范围为0°～360°，为精确计算倾向方向，需要利用向量与Y，Z轴的余弦值来判断结构面的法向向量指向的卦限；具体判别方法如下：假设向量向量为的水平投影向量，根据倾向方向的定义，为迹线EG的倾向；α_L，γ_L分别为向量与X，Z轴的角；k_L为向量与Y轴的夹角，则有Since the inclination angle ranges from 0° to 360°, in order to accurately calculate the inclination direction, it is necessary to use the vector and the cosine values of the Y and Z axes to determine the hexagram limit pointed to by the normal vector of the structural surface; the specific method of discrimination is as follows: Assume that the vector vector for The horizontal projection vector of , according to the definition of the dip direction, is the tendency of the trace EG; α _L , γ _L are the vectors Angle with X, Z axis; k _L is a vector The angle with the Y axis, then there is

${cosk cosk}_{L L} = = \frac{{b b}_{L L}}{\sqrt{{a a}_{L L}^{22} + + {b b}_{L L}^{22}}}$

${cosα cosα}_{L L} = = \frac{{a a}_{L L}}{\sqrt{{a a}_{L L}^{22} + + {b b}_{L L}^{22} + + {c c}_{L L}^{22}}}$

${cosγ cosγ}_{L L} = = \frac{{c c}_{L L}}{\sqrt{{a a}_{L L}^{22} + + {b b}_{L L}^{22} + + {c c}_{L L}^{22}}}$

第七步：激光点云与数字影像融合Step 7: Fusion of laser point cloud and digital image

1)近似拟合特征平面1) Approximately fitting the feature plane

前面的工作中已获得各簇三角面片，根据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)用形心表征此三角面：Each cluster of triangular faces has been obtained in the previous work, and the vertex coordinates of the triangular faces are extracted according to the ID. Assume that the coordinates of the three vertices of a triangle in X ₁ are (x _R1 , y _R1 , z _R1 ), (x _R2 , y _R2 ,z _R2 ), (x _R3 ,y _R3 ,z _R3 ), then the centroid coordinates (x _R ,y _R ,z _R ) of the triangle can be obtained according to the following formula:

${x x}_{R R} = = \frac{{x x}_{R R 11} + + {x x}_{R R 22} + + {x x}_{R R 33}}{33}$

${y the y}_{R R} = = \frac{{y the y}_{R R 11} + + {y the y}_{R R 22} + + {y the y}_{R R 33}}{33} - - - - - - ((77 - - 11))$

${z z}_{R R} = = \frac{{z z}_{R R 11} + + {z z}_{R R 22} + + {z z}_{R R 33}}{33}$

计算出同簇三角面片的形心坐标后则可通过最小二乘拟合的方法拟合出近似的特征平面，空间平面的方程通常形如：A_Rx+B_Ry+C_Rz+1＝0，用同簇的H个三角面形心坐标拟合出这个平面可以表示成以下矩阵形式：After calculating the centroid coordinates of the triangular faces of the same cluster, the approximate feature plane can be fitted by the method of least squares fitting. The equation of the space plane is usually in the form: A _R x+B _R y+C _R z+ 1=0, the plane can be expressed in the following matrix form by fitting the centroid coordinates of H triangular faces of the same cluster:

$[\begin{matrix} {x x}_{R R 11} & {y the y}_{R R 11} & {z z}_{R R 11} \\ \begin{matrix} \cdot &Center Dot; \\ \cdot \cdot \\ \cdot &Center Dot; \end{matrix} & \begin{matrix} \cdot &Center Dot; \\ \cdot &Center Dot; \\ \cdot &Center Dot; \end{matrix} & \begin{matrix} \cdot &Center Dot; \\ \cdot &Center Dot; \\ \cdot &Center Dot; \end{matrix} \\ {x x}_{R R H h} & {y the y}_{R R H h} & {z z}_{R R H h} \end{matrix}] [\begin{matrix} {A A}_{R R} \\ {B B}_{R R} \\ {C C}_{R R} \end{matrix}] = = [\begin{matrix} - - 11 \\ - - 11 \\ - - 11 \end{matrix}]$

经过运算可化成以下形式：After calculation, it can be transformed into the following form:

$[\begin{matrix} {Σx Σx}_{g g}^{22} & {Σx Σx}_{g g} {y the y}_{g g} & {Σx Σx}_{g g} {z z}_{g g} \\ {Σx Σx}_{g g} {y the y}_{g g} & {Σy Σy}_{g g}^{22} & {Σy Σy}_{g g} {z z}_{g g} \\ {Σx Σx}_{g g} {z z}_{g g} & {Σy Σy}_{g g} {z z}_{g g} & {Σz Σz}_{g g}^{22} \end{matrix}] [\begin{matrix} {A A}_{R R} \\ {B B}_{R R} \\ {C C}_{R R} \end{matrix}] = = [\begin{matrix} - - Σ Σ {x x}_{g g} \\ - - Σ Σ {y the y}_{g g} \\ - - {Σz Σz}_{g g} \end{matrix}] ((g g = = 11,, 2... 2... H h))$

所以得到： $[\begin{matrix} A_{R} \\ B_{R} \\ C_{R} \end{matrix}] = {[\begin{matrix} {Σx}_{g}^{2} & {Σx}_{g} y_{g} & {Σx}_{g} z_{g} \\ {Σx}_{g} y_{g} & {Σy}_{g}^{2} & {Σy}_{g} z_{g} \\ {Σx}_{g} z_{g} & {Σy}_{g} z_{g} & {Σz}_{g}^{2} \end{matrix}]}^{- 1} [\begin{matrix} - Σ x_{g} \\ - {Σy}_{g} \\ - {Σz}_{g} \end{matrix}] (g = 1, 2... H),$ 系数A_R,B_R,C_R即可求得，平面方程即可得出；so get: $[\begin{matrix} A_{R} \\ B_{R} \\ C_{R} \end{matrix}] = {[\begin{matrix} {Σx}_{g}^{2} & {Σx}_{g} {the y}_{g} & {Σx}_{g} z_{g} \\ {Σx}_{g} {the y}_{g} & {Σy}_{g}^{2} & {Σy}_{g} z_{g} \\ {Σx}_{g} z_{g} & {Σy}_{g} z_{g} & {Σz}_{g}^{2} \end{matrix}]}^{- 1} [\begin{matrix} - Σ x_{g} \\ - {Σy}_{g} \\ - {Σz}_{g} \end{matrix}] (g = 1, 2... h),$ The coefficients A _R , B _R , and C _R can be obtained, and the plane equation can be obtained;

2)迹线与结构面归并2) Merging traces and structural planes

从点云数据中，提取出迹线所对应的各个点的坐标，逐点求出点到各个簇所拟合出的近似平面的距离，取最小距离即能描述该迹线与各簇结构面的位置关系，From the point cloud data, the coordinates of each point corresponding to the trace are extracted, and the distance from the point to the approximate plane fitted by each cluster is calculated point by point, and the minimum distance can be used to describe the trace and the structural surface of each cluster positional relationship,

点到平面距离的算法： $D = \frac{A x + B y + C z + 1}{\sqrt{A^{2} + B^{2} + C^{2}}} - - - (7 - 2)$ Algorithm for point-to-plane distance: $D. = \frac{A x + B the y + C z + 1}{\sqrt{A^{2} + B^{2} + C^{2}}} - - - (7 - 2)$

每个簇周围都存在与其距离最近的1条或多条迹线，可将这1条或多条迹线与该簇归并入同一个组，那么该组的产状信息即为对应区域的产状信息。若存在某迹线周围不存在相近的结构面，则将该迹线单独列为一组，用以表征所在区域的产状信息；若存在某结构面周围不存在相近的迹线，则将该结构面单独列为一组，用以表征所在区域的产状信息。There is one or more traces closest to each cluster, and these one or more traces can be combined with the cluster into the same group, then the occurrence information of this group is the occurrence information of the corresponding area. status information. If there is no similar structural surface around a certain trace, the trace will be listed as a group separately to represent the occurrence information of the area; if there is no similar structural surface around a certain structural surface, the Structural surfaces are listed as a group separately to represent the occurrence information of the region.

Claims

1. A method for measuring the occurrence of rock mass discontinuities that fuses laser point clouds and digital images, is characterized in that the measurement content of the measurement method includes point cloud data extraction of rock mass discontinuity plane occurrences, rock mass joint trace production The extraction of shape image data and point cloud data is fused, and its expression method adopts the matching grouping of traces and structural

the method described;

In the extraction of point cloud data, the filter method is used to denoise first, and then the point cloud data is triangulated to establish a preliminary three-dimensional model of rock mass. , The calculation of the normal vector of the triangular surface constructs the characteristic matrix, and uses the DBSCAN cluster analysis method based on the high-density connected area to group the structural surfaces initially, and then uses the connectivity of the triangular surfaces to distinguish the secondary grouping of the structural surfaces, and finally calculates The occurrence of structural planes; the rock joint traces need to be extracted from the image data. Firstly, the digital image is processed in gray scale, and then the contour line of the rock is extracted by the image mixing threshold segmentation method. After pixel stripping and removal of severed limbs Extract the skeleton information of the trace, and finally mark the trace and calculate the occurrence; The positional relationship among them is identified and merged, and the final result after merging is used to characterize the occurrence information of rock mass structural planes.

2. the method for measuring the occurrence of rock mass discontinuities that fuses laser point clouds and digital images as claimed in claim 1, is characterized in that the measuring method mainly comprises the steps:

The first step: collecting field data;

Step 2: Match and fuse point cloud data and image data;

Step 3: Triangulate the point cloud data, calculate the parameters of the triangular normal phasor after division, and perform cluster analysis on the normal vector;

Step 4: Merge the results after cluster analysis, and calculate the occurrence of the merged data;

Step 5: Correct and grayscale the image data, then extract the outline and skeleton of the joint and mark the skeleton in the point cloud;

Step 6: Calculate the inner angle of the triangle formed by any three points corresponding to the point cloud of the trace, and then judge the occurrence of the trace and calculate the occurrence;

Step 7: Approximately fit the extracted structural surface and calculate the distance between the structural surface and each trace, and then group the structural surface and trace according to the positional relationship between the structural surface and the trace;

Step 8: The measurement work is completed.