A rock mass structure plane rapid automatic identification method and device based on three-dimensional point cloud

Abstract

The application provides a kind of rock mass structure plane fast automatic identification method and device based on three-dimensional point cloud, it is related to image recognition, rock mass structure plane identification technical field, comprising: S1, obtains the three-dimensional point cloud model of the engineering rock mass of survey area;S2, obtains the super voxel model by the super voxel segmentation algorithm and the super voxel normal vector estimation RANSAC algorithm to three-dimensional subsampling;S3, obtains the structure plane model by the structure plane identification of super voxel model to the double criterion constraint region growing algorithm;S4, respectively calculates the occurrence of structure plane model structure plane.This application effectively reduces the seed point number that region growing algorithm needs iteration, improves processing efficiency, can accurately identify structure plane, according to the structure plane occurrence analysis generated by this application, can obtain more accurate and more real rock mass condition, solves the problem of traditional measurement means in efficiency, accuracy and safety.

Description

Technical Field

[0001] This invention relates to the fields of image recognition and rock mass structure surface recognition, and in particular to a method and apparatus for rapid automatic recognition of rock mass structure surfaces based on three-dimensional point clouds. Background Technology

[0002] Rock mass structural planes are key factors influencing rock mass deformation and failure. Accurately identifying the orientation of rock mass structural planes is crucial for preventing potential geological hazards. Currently widely used methods for measuring structural planes (traditional geological compasses, inclinometers, etc.) have significant limitations, specifically: traditional methods are time-consuming and inefficient; the representativeness of measurement data samples is insufficient, with available samples often limited to the accessible area of ​​the rock mass, and sample selection relying on the subjective experience of the operators; measurement operations involve safety risks, as operators must work close to the free face when the rock mass is unstable, facing potential risks such as rock spalling, water and mud inrushes, and even localized collapses. These technical deficiencies severely restrict the timeliness and accuracy of geological information collection, making it difficult to meet the needs of modern intelligent geological surveys.

[0003] Compared to traditional surveying methods, non-contact surveying methods such as digital photogrammetry and 3D laser scanning offer significant advantages in efficiency and cost-effectiveness. These technologies enable rapid acquisition of rock mass information, constructing detailed 3D point cloud models of the rock mass with millimeter-level precision, and allowing for more comprehensive and objective geological analysis. For geological survey personnel, non-contact on-site data collection methods and highly efficient acquisition processes greatly reduce the risks associated with surveying operations; simultaneously, the collected data can be persistently stored, allowing for the reconstruction of the rock mass scene at the time of acquisition.

[0004] Non-contact measurement typically employs unsupervised methods for structural surface identification. Among these, region growing is a segmentation algorithm that boasts advantages such as simple principles, easy implementation, and detailed segmentation. However, in traditional region growing algorithms, the clustering and segmentation of point set data usually relies on a single fixed threshold to control the growth process, such as the distance threshold between adjacent points and the angle threshold between normal vectors. In complex geological scenarios, rock surfaces are often uneven rather than ideally smooth structural surfaces, and they also face problems such as sparse point clouds, noise interference, or varying scales of discontinuous surfaces. Traditional region growing algorithms are prone to insufficient or excessive region growth, severely affecting the accuracy and completeness of identification.

[0005] Furthermore, downsampling is a common step in preprocessing large point cloud models. Among them, the voxelization mesh method is the most common downsampling method. This method first creates a voxel mesh (3D bounding box in space) covering the input point cloud, and then selects the centroid of all points in each voxel as the sampled points, thereby reducing the number of points in the point cloud and making the point cloud more uniform. However, the voxelization method is prone to point cloud resolution loss and detail loss, resulting in cross-surface merging and boundary blurring, which in turn affects the accuracy of structural surface recognition.

[0006] Therefore, in order to address the above problems, there is an urgent need to propose a method and device for rapid and automatic identification of rock mass structural surfaces based on three-dimensional point clouds to solve the aforementioned technical issues. Summary of the Invention

[0007] To address the aforementioned technical problems in related technologies, this invention proposes a rapid and automatic identification method and device for rock mass structural surfaces based on three-dimensional point clouds. The aim is to overcome the shortcomings of traditional measurement methods in terms of efficiency, accuracy, and safety, and to introduce a new, more efficient, accurate, and safer rock mass structural surface identification technology to promote the development of engineering rock mass geological surveys towards intelligence and informatization.

[0008] In a first aspect, the present invention provides a method for rapid and automatic identification of rock mass structural surfaces based on three-dimensional point clouds, comprising the following steps:

[0009] S1. Obtain a three-dimensional point cloud model of the engineering rock mass in the survey area;

[0010] S2. Adaptive clustering of the 3D point cloud model is performed using the supervoxel segmentation algorithm to obtain a list of supervoxel blocks and the point set to which the supervoxel blocks belong. The RANSAC algorithm for supervoxel normal vector estimation is used to perform plane fitting on each supervoxel block in the list of supervoxel blocks to obtain its corresponding fitting plane, forming a list of supervoxel block fitting planes. A supervoxel model is constructed based on the list of supervoxel blocks, the point set to which the supervoxel blocks belong, and the list of supervoxel block fitting planes.

[0011] S3. Obtain the structural surface model by identifying the structural surface of the supervoxel model using the dual-criteria constrained region growing algorithm;

[0012] S4. Calculate the attitude of each structural surface in the structural surface model.

[0013] Specifically, the supervoxel segmentation algorithm in step S2 includes:

[0014] S21. Determine the neighborhood point set of each point in the 3D point cloud model according to the nearest neighbor algorithm, and calculate the similarity metric between each point in the 3D point cloud model and its neighborhood points respectively.

[0015] S22. Set estimation parameters The initial value is the median of the similarity metric between each point and its neighbors in all three-dimensional point cloud models;

[0016] S23. Perform voxel downsampling on the 3D point cloud model according to the super-voxel resolution R to obtain the downsampled point cloud model. The number of points in the downsampled point cloud model is DK.

[0017] S24. Initialize each point in the 3D point cloud model as an independent supervoxel block, select the point as the representative point of each supervoxel block and initialize the point set to which the supervoxel block belongs; initialize the list of adjacent supervoxel blocks of each supervoxel block as a sequence of supervoxel blocks constructed from the neighborhood points of each point in the 3D point cloud model.

[0018] S25. Traverse all current supervoxes. For each supervoxe, traverse the list of adjacent supervoxes. If the following condition is met, merge the adjacent supervoxes into the current supervoxe, merge the set of points belonging to the adjacent supervoxes into the set of points belonging to the current supervoxe, and merge the list of adjacent supervoxes:

[0019] ;

[0020] in, The number of points contained in the adjacent hypervoxel block; and These represent the representative points of the current superpixel and its adjacent superpixels, respectively. express and Similarity measure between two points;

[0021] S26. Estimate the parameters Double the number of supervoxels and continue with step S25 until the number of supervoxels is DK, finally obtaining a list of supervoxels of number DK and their corresponding point sets.

[0022] Specifically, the similarity metric between each point in the 3D point cloud model and its neighboring points is calculated using the following formula:

[0023] ;

[0024] Where p and N represent points in the 3D point cloud model and their neighboring points, respectively; A measure of similarity between points p and N; Let p be the normal vector of point p; This represents the normal vector of point N; R represents the Euclidean distance between point p and point N; R is the supervoxel resolution.

[0025] Specifically, the RANSAC algorithm for estimating the hypervoxel normal vector includes:

[0026] S27. Randomly and without repetition select a point from the set of points belonging to the supervoxel block, and use the spatial coordinates and normal vector of this point to solve for and determine the candidate plane:

[0027] S28. Calculate the interior ratio of the candidate plane using the interior ratio formula. The formula for the interior point ratio is as follows:

[0028] ;

[0029] Where point p is a point belonging to the point set s of the supervoxel block; is the number of points in the point set s to which the supervoxel belongs; n is the normal vector of the candidate plane; d is the intercept of the candidate plane; It is the interior point threshold; indicator function. This means that the value is 1 when the condition is met and 0 when the condition is not met.

[0030] S29. Repeat steps S27-S28 until the number of repetitions reaches the preset maximum number or there are no points to choose from within the super-void block. At this point, select the candidate plane with the largest interior point ratio as the fitting plane of the super-void block. That is, the plane equation of the candidate plane with the largest interior point ratio is the fitting plane equation of the super-void block. At the same time, the normal vector of the candidate plane with the largest interior point ratio is used as the fitting plane normal vector of the super-void block.

[0031] Specifically, in step S27, the candidate plane is determined according to the following formula:

[0032] ;

[0033] Where n is the normal vector of a point in the point set to which the supervoxel belongs; d is the intercept calculated from the spatial coordinates of that point; and X represents the coordinate variable of a point belonging to that plane.

[0034] Specifically, the dual-criteria constrained region growth algorithm includes:

[0035] S31. Determine the neighborhood supervoxel block set of each supervoxel block in the supervoxel block list of the supervoxel model according to the nearest neighbor algorithm, wherein the position of the supervoxel block is determined by calculating the average coordinate of each point contained in the supervoxel block.

[0036] S32. Randomly select a supervoxel block from the list of supervoxel blocks in the supervoxel model as the current structural surface. Select the supervoxel block as the seed point. Use the fitting plane of the supervoxel block as the fitting plane of the current structural surface. Use the fitting plane equation of the supervoxel block as the fitting plane equation of the current structural surface. Calculate the interior point ratio of the current fitting plane according to the interior point ratio formula based on the point set to which the current structural surface belongs and the fitting plane equation of the current structural surface. Use the interior point ratio of the current fitting plane as the interior point ratio of the fitting plane of the current structural surface.

[0037] S33. Iterate through all neighboring supervoids in the set of neighboring supervoids of the seed point, and include the neighboring supervoids that satisfy the dual criterion constraints and the set of points to which the neighboring supervoids belong into the current structural surface:

[0038] S34. If the dual criterion constraint in step S33 is satisfied, then update the fitting plane equation of the current structural surface. Calculate the interior point ratio of the fitting plane of the super-voxel block in the neighboring area of ​​the seed point based on the interior point ratio formula using the fitting plane equation of the current structural surface's point set and the seed point's neighborhood super-voxel block. If the interior point ratio of the fitting plane of the seed point's neighborhood super-voxel block is greater than the interior point ratio of the current structural surface's fitting plane, use it as the interior point ratio of the current structural surface's fitting plane, and update the fitting plane equation of the structural surface to the fitting plane equation of the seed point's neighborhood super-voxel block.

[0039] S35. Continue to select non-repeating supervoxels from the current structural plane as seed points, and repeat steps S33-S34 until there are no more supervoxels to select in the current structural plane. At this point, the division of the structural plane is complete.

[0040] S36. Continue to randomly select a supervoxel block that is not included in the current structural plane from the list of supervoxel blocks in the supervoxel model, and use it as the new current structural plane. Select the supervoxel block as the seed point and repeat steps S32-S35 until all supervoxel blocks are divided into independent structural planes.

[0041] Specifically, the dual-criteria constraint conditions described in step S33 are as follows:

[0042] ;

[0043] in, The normal vector of the fitted plane for the seed point; The normal vector of the fitting plane of the hypervoxel block in the neighborhood of the seed point; This is the normal vector of the fitting plane for the current structural surface; The norm of a vector; The first judgment threshold; This is the second judgment threshold.

[0044] Specifically, step S4 includes: extracting structural surfaces from the structural surface list of the structural surface model, and calculating the attitude of each structural surface according to the following formula; the attitude includes dip. and tilt angle :

[0045] ;

[0046] ;

[0047] in, These are the components of the normal vector n of the fitting plane of the structural surface in the X, Y, and Z directions, respectively.

[0048] Specifically, the structure of the three-dimensional point cloud model includes point cloud spatial coordinates and point cloud normal vectors.

[0049] Secondly, the present invention provides a rapid automatic identification device for rock mass structural surfaces based on three-dimensional point clouds, and a rapid automatic identification method for rock mass structural surfaces based on three-dimensional point clouds as described in the first aspect above, comprising the following units:

[0050] The point cloud model generation unit is used to obtain a three-dimensional point cloud model of the engineering rock mass in the survey area;

[0051] The supervoxel model construction unit is used to adaptively cluster the 3D point cloud model through the supervoxel segmentation algorithm to obtain a list of supervoxel blocks and the point set to which the supervoxel blocks belong. Then, the supervoxel normal vector estimation RANSAC algorithm is used to perform plane fitting on each supervoxel block in the list of supervoxel blocks to obtain its corresponding fitting plane, forming a list of supervoxel block fitting planes. The supervoxel model is constructed based on the list of supervoxel blocks, the point set to which the supervoxel blocks belong, and the list of supervoxel block fitting planes.

[0052] The structural surface model construction unit is used to identify structural surfaces in the supervoxel model and obtain the structural surface model by using a dual-criteria constrained region growth algorithm.

[0053] The attitude analysis unit is used to calculate the attitude of each structural surface in the structural surface model.

[0054] This invention provides a method and apparatus for rapid and automatic identification of rock mass structural surfaces based on 3D point clouds. By first segmenting the 3D point cloud of the rock mass into hypervoxel blocks with planar characteristics, the number of seed points required for subsequent region growing algorithms is effectively reduced, greatly improving processing efficiency. Simultaneously, an improved dual-criteria constrained region growing algorithm is used to accurately identify structural surfaces. Analysis of the attitude (dip angle) of the structural surfaces generated in this embodiment yields more accurate and realistic rock mass conditions, providing a more reliable basis for engineering site assessments. This invention overcomes the shortcomings of traditional measurement methods in terms of efficiency, accuracy, and safety, and promotes the development of engineering rock mass geological surveys towards intelligence and informatization.

[0055] Furthermore, in this embodiment, multiple high-precision normal vector estimations are performed (including normal vector estimation of the original point cloud, planar normal vector estimation of the hypervoxel, and incremental planar normal vector estimation of the structural surface), so that the fitted planar normal vector of the extracted structural surface is always close to the normal vector of the original point cloud. Therefore, it is beneficial to maintain the accuracy of the attitude calculation and can realize the processing of large-scale rock mass point clouds and structural surface identification.

[0056] Furthermore, the RANSAC algorithm for estimating supervoxel normal vectors in this invention can effectively reduce the computational cost of least squares fitting and improve the estimation efficiency of supervoxel blocks; and the normal vectors are directly derived from points in the supervoxel, avoiding the spurious plane problem that is prone to occur in the traditional RANSAC algorithm.

[0057] Furthermore, traditional region growing algorithms rely on a single threshold for growth and are sensitive to threshold changes. They are prone to oversegmentation or undersegmentation in sparse or undulating point cloud regions, and have low accuracy in structural surface recognition. The dual-criteria constrained region growing algorithm proposed in this invention can simultaneously consider the local geometric planarity of the point cloud and the overall planarity of the structural surface, which greatly improves the problem of oversegmentation or undersegmentation that traditional algorithms are prone to. Attached Figure Description

[0058] To more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the drawings used in the embodiments will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0059] Figure 1 A schematic diagram of a method for rapid and automatic identification of rock mass structure surfaces based on three-dimensional point clouds, provided in an embodiment of the present invention;

[0060] Figure 2 Real-life images of the engineering rock mass in the survey area provided in this embodiment of the invention;

[0061] Figure 3 The correspondence between attitude (dipping, dip angle) and color provided in the embodiments of the present invention;

[0062] Figure 4 The three-dimensional point cloud model and its corresponding orientation cloud map are provided for embodiments of the present invention; Figure 4 (a) in the figure is a three-dimensional point cloud model; Figure 4 (b) in the figure is the attitude cloud map corresponding to the three-dimensional point cloud model;

[0063] Figure 5 The superbody pixel model and its corresponding orientation cloud map provided in the embodiments of the present invention; Figure 5 (a) in the figure represents the supervoxel model; Figure 5 (b) in the figure is the orientation cloud map corresponding to the supervoxel model;

[0064] Figure 6 The structural surface model and its corresponding attitude cloud map are provided for embodiments of the present invention; Figure 6 (a) in the diagram represents the structural surface model; Figure 6 (b) in the figure is the attitude contour map corresponding to the structural surface model;

[0065] Figure 7 This is a schematic diagram of a rapid automatic identification device for rock mass structure surfaces based on three-dimensional point clouds, provided as an embodiment of the present invention. Detailed Implementation

[0066] The present invention will be explained in detail through the following embodiments. The purpose of this invention is to protect all technical improvements within its scope. In the description of this invention, it should be understood that the terms "first" and "second" are used for descriptive purposes only and should not be construed as indicating or implying relative importance or implicitly specifying the number of indicated technical features. Therefore, a feature defined as "first" or "second" may explicitly or implicitly include one or more of that feature. In the description of this invention, "a plurality of" means two or more, unless otherwise explicitly specified.

[0067] To make the technical problems, technical solutions and advantages of the present invention clearer, a detailed description will be given below in conjunction with the accompanying drawings and specific embodiments.

[0068] Example 1

[0069] refer to Figure 1 This embodiment provides a method for rapid and automatic identification of rock mass structural surfaces based on three-dimensional point clouds, including the following steps:

[0070] S1. Obtain a three-dimensional point cloud model of the engineering rock mass in the survey area;

[0071] Specifically, a three-dimensional point cloud model of the engineering rock mass in the survey area is acquired through digital photogrammetry or three-dimensional laser scanning;

[0072] The digital photogrammetry specifically involves taking multiple, multi-view photos of the engineering rock mass in the survey area using devices with high-definition video capabilities, such as mobile phones, cameras, and drones. Then, digital photogrammetry software (such as COLMAP, Agisoft Metashape, etc.) is used to perform three-dimensional reconstruction of the acquired images, and a three-dimensional point cloud model of the engineering rock mass in the survey area is established and exported.

[0073] The three-dimensional laser scanning specifically involves scanning and collecting data on the engineering rock mass in the survey area using a handheld or base station-type laser scanner, and then exporting a three-dimensional point cloud model of the engineering rock mass in the survey area using the laser scanner's processing software.

[0074] The structure of the three-dimensional point cloud model includes point cloud spatial coordinates and point cloud normal vectors;

[0075] It is understandable that a 3D point cloud model is a collection of a large number of discrete 3D points. Each point contains its point cloud spatial coordinates in 3D space. Estimating the point cloud normal vector of the point cloud data is a mature existing technology. For example, principal component analysis (PCA) can be used to estimate the point cloud data, or point cloud processing software such as Agisoft Metashape and Cloudcompare can be used to estimate the point cloud data. All of these methods can obtain the point cloud normal vector, that is, the normal vector of each point in the point cloud data.

[0076] In this embodiment, a base station-type laser scanner is used to scan and collect data on the engineering rock mass in the survey area, such as... Figure 2 As shown. For ease of explanation, a correspondence between attitude (dip, dip angle) or normal vector and color has been established, as follows: Figure 3 As shown, this relationship allows visualization of the model's attitude or normal vector distribution. The 3D point cloud model established in this embodiment is as follows: Figure 4 As shown in (a) above, the 3D point cloud model corresponds to the attitude cloud map as follows: Figure 4 As shown in (b);

[0077] S2. Adaptive clustering of the 3D point cloud model is performed using the supervoxel segmentation algorithm to obtain a list of supervoxel blocks and the point set to which the supervoxel blocks belong. The RANSAC algorithm for supervoxel normal vector estimation is used to perform plane fitting on each supervoxel block in the list of supervoxel blocks to obtain its corresponding fitting plane, forming a list of supervoxel block fitting planes. A supervoxel model is constructed based on the list of supervoxel blocks, the point set to which the supervoxel blocks belong, and the list of supervoxel block fitting planes.

[0078] The hypervoxel segmentation algorithm specifically includes:

[0079] S21. Determine the neighborhood point set of each point in the 3D point cloud model using the nearest neighbor algorithm, and calculate the similarity metric between each point in the 3D point cloud model and its neighborhood points according to the following formula:

[0080] ;

[0081] Where p and N represent points in the 3D point cloud model and their neighboring points, respectively, both of which can be represented by 3D coordinates; A measure of similarity between points p and N; Let p represent the normal vector. This represents the normal vector of point N; R represents the Euclidean distance between point p and point N; R is the supervoxel resolution, used to control the size of the generated supervoxel block, and its value is determined by the minimum structural surface size that needs to be identified; in this embodiment, R is taken as 0.1m;

[0082] The nearest neighbor algorithm, or K-Nearest Neighbor (KNN) classification algorithm, is one of the simplest methods in data mining classification techniques. K-Nearest Neighbor means that each sample can be represented by its K nearest neighbors. The nearest neighbor algorithm is a method for classifying each record in a dataset; this is existing technology and will not be elaborated upon here.

[0083] S22. Set estimation parameters The initial value is the median of the similarity metric between each point and its neighbors in all three-dimensional point cloud models;

[0084] S23. Perform voxel downsampling on the 3D point cloud model according to the super-voxel resolution R to obtain the downsampled point cloud model. The number of points in the downsampled point cloud model is DK.

[0085] Voxel downsampling is a commonly used data compression and downsampling method in point cloud algorithms. Its principle is to divide the point cloud data into small voxels (3D pixels), and then retain only one point in each voxel as a representative.

[0086] The super-voxel resolution R is used to control the size of the voxels. The larger R is, the larger the spatial range covered by the voxels and the fewer the number of points DK after downsampling; conversely, the smaller R is, the more details are preserved and the larger DK is.

[0087] S24. Initialize each point in the 3D point cloud model as an independent supervoxel, select the point as the representative point of each supervoxel, and initialize the point set to which the supervoxel belongs. At this time, each supervoxel contains only a unique point, and the point set to which the supervoxel belongs contains the representative point. Initialize the list of adjacent supervoxels of each supervoxel as a sequence of supervoxels constructed from the neighborhood points of each point in the 3D point cloud model.

[0088] It is understandable that the set of points to which a supervox belongs is the set of points contained in that supervox.

[0089] S25. Traverse all current supervoxes. For each supervoxe, traverse the list of adjacent supervoxes. If the following condition is met, merge the adjacent supervoxes into the current supervoxe, merge the set of points belonging to the adjacent supervoxes into the set of points belonging to the current supervoxe, and merge the list of adjacent supervoxes:

[0090] ;

[0091] in, The number of points contained in the adjacent hypervoxel block; and These represent the representative points of the current super-void block and the adjacent super-void blocks, respectively. express and Similarity measure between two points;

[0092] S26. Estimate the parameters Double the number of supervoxels and continue with step S25 until the number of supervoxels is DK, and finally obtain a list of supervoxels with a number of DK and the corresponding point set to which the supervoxels belong.

[0093] After obtaining the list of supervoxels and the corresponding point set to which they belong, the RANSAC algorithm for estimating supervoxel normal vectors is used to perform plane fitting on each supervoxel in the list.

[0094] The RANSAC algorithm for estimating the hypervoxel normal vector specifically includes:

[0095] S27. Randomly and without repetition select a point from the set of points belonging to the supervoxel block. Use the spatial coordinates and normal vector of this point to solve for the candidate plane based on the following plane equation:

[0096] ;

[0097] Where n is the normal vector of the point; d is the intercept calculated from the spatial coordinates of the point; and X represents the coordinate variable of the point belonging to the plane.

[0098] It is understandable that the normal vector of a plane is contained within the plane's equation. In the equation, n represents the normal vector of the candidate plane, and d represents the intercept of the candidate plane.

[0099] S28. Calculate the interior ratio of the candidate plane using the interior ratio formula. The formula for the interior point ratio is as follows:

[0100] ;

[0101] In the formula, point p is a point belonging to the point set s of the supervoxel block; is the number of points in the point set s to which the supervoxel belongs; n is the normal vector of the candidate plane; d is the intercept of the candidate plane; This is the interior point threshold, which is set to 0.001m in this embodiment; indicator function. The value is 1 when the condition is met and 0 when the condition is not met. The inlier rate of the plane refers to the proportion of points in all data points that are less than a given inlier threshold, and is used to measure the quality of the plane fitting.

[0102] S29. Repeat steps S27-S28 until the number of repetitions reaches the preset maximum number or there are no points to choose from within the super-void block. At this point, select the candidate plane with the largest interior point ratio as the fitting plane of the super-void block. That is, the plane equation of the candidate plane with the largest interior point ratio is the fitting plane equation of the super-void block. At the same time, the normal vector of the candidate plane with the largest interior point ratio is used as the fitting plane normal vector of the super-void block.

[0103] The preset maximum number of times is selected according to actual needs;

[0104] RANSAC (Random Sample Consensus) is an efficient method for estimating mathematical model parameters from a sample set containing outliers. It was proposed by Fischler and Bolles in 1981. The algorithm assumes that the data contains interior points that can be described by the model and outliers that deviate from the model. It initializes the model by randomly sampling a minimum set of samples and iteratively filters the interior point set to optimize the parameters based on a set threshold.

[0105] The RANSAC algorithm for estimating the supervoxel normal vector in this embodiment is an improved form of the traditional RANSAC algorithm. The difference between the traditional RANSAC and the traditional method is that the algorithm randomly selects one point in the supervoxel block as a candidate plane each time without repetition; the traditional method performs least squares fitting of the plane by randomly selecting three points in the supervoxel block each time without repetition.

[0106] The advantages of the RANSAC algorithm for estimating supervoxel normal vectors in this embodiment are:

[0107] 1. It reduces the computational cost of least squares fitting and improves the estimation efficiency of supervoid blocks;

[0108] 2. The normal vector is directly derived from the points in the supervoxel, avoiding the false plane problem that is prone to occur in traditional algorithms.

[0109] The RANSAC algorithm is used to estimate the fitting planes for each supervoxel block using supervoxel normal vector estimation. These fitting planes are then listed to obtain a supervoxel block fitting plane list. Finally, the supervoxel block list, the point set to which the supervoxel block belongs, and the supervoxel block fitting plane list together constitute the supervoxel model. In this embodiment, the supervoxel model generated according to the above steps is as follows: Figure 5 As shown in (a) in the figure; the orientation contour map corresponding to the supervoxel block model is as follows. Figure 5 As shown in (b);

[0110] It is understandable that there is a one-to-one correspondence between super-body blocks and super-body block fitting planes. Therefore, there is a correspondence between the super-body blocks in the list of super-body blocks and the fitting planes in the list of super-body block fitting planes. Both can be found through the super-body blocks. That is, the corresponding super-body block fitting plane and the fitting plane normal vector of the super-body block can be obtained through the super-body block.

[0111] Step S2 is essentially a 3D downsampling process, which includes two key steps: First, the 3D point cloud model is processed by a supervoxel segmentation algorithm to adaptively cluster and obtain DK supervoxel blocks, thereby obtaining a list of supervoxel blocks. At this point, a basic supervoxel model has been obtained, because the supervoxel model is composed of supervoxel blocks. Each supervoxel block contains a certain number of 3D points of the original point cloud model (i.e., the point set to which the supervoxel block belongs). Then, the RANSAC algorithm for supervoxel normal vector estimation is used to perform plane fitting on each supervoxel block to obtain the fitting plane of each supervoxel block. This is also an extension of the basic supervoxel model attributes.

[0112] S3. Obtain the structural surface model by identifying the structural surface of the supervoxel model using the dual-criteria constrained region growing algorithm;

[0113] The dual-criteria constrained region growth algorithm specifically includes:

[0114] S31. Determine the neighborhood supervoxel block set of each supervoxel block in the supervoxel block list of the supervoxel model according to the nearest neighbor algorithm, wherein the position of the supervoxel block is determined by calculating the average coordinate of each point contained in the supervoxel block.

[0115] S32. Randomly select a supervoxel block from the list of supervoxel blocks in the supervoxel model as the current structural surface. Select the supervoxel block as the seed point. Use the fitting plane of the supervoxel block as the fitting plane of the current structural surface. Use the fitting plane equation of the supervoxel block as the fitting plane equation of the current structural surface. Calculate the interior point ratio of the current fitting plane according to the interior point ratio formula based on the point set to which the current structural surface belongs and the fitting plane equation of the current structural surface. Use the interior point ratio of the current fitting plane as the interior point ratio of the fitting plane of the current structural surface.

[0116] The interior point ratio formula is the same as the interior point ratio formula in step S28.

[0117] S33. Iterate through all neighboring supervoids in the set of neighboring supervoids of the seed point, and include the neighboring supervoids that satisfy the dual criterion constraints and the set of points to which the neighboring supervoids belong into the current structural surface:

[0118] The dual-criteria constraints are as follows:

[0119] ;

[0120] In the formula, The normal vector of the fitted plane for the seed point; The normal vector of the fitting plane of the hypervoxel block in the neighborhood of the seed point; This is the normal vector of the fitting plane for the current structural surface; The first judgment threshold; This is the second judgment threshold; Represents the norm of a vector; first threshold judgment. With the second judgment threshold The value is selected according to actual needs. In this embodiment, Take 0.94, Take 0.87.

[0121] S34. If the dual criterion constraint in step S33 is satisfied, then consider updating the fitting plane equation of the current structural surface. Calculate the interior point ratio of the fitting plane of the super-voxel block in the neighboring area of ​​the seed point based on the interior point ratio formula, using the fitting plane equation of the current structural surface's point set and the seed point's neighborhood super-voxel block. If the interior point ratio of the fitting plane of the seed point's neighborhood super-voxel block is greater than the interior point ratio of the current structural surface's fitting plane, use it as the interior point ratio of the current structural surface's fitting plane, and update the fitting plane equation of the structural surface to the fitting plane equation of the seed point's neighborhood super-voxel block.

[0122] As is understandable, a structural plane is essentially a set of multiple super-voids, so the set of points to which a structural plane belongs is the set of points to which the super-voids in the structural plane belong; the fitting plane equation of a structural plane is actually selected from the fitting plane equations of the super-voids to which the structural plane belongs.

[0123] S35. Continue to select non-repeating supervoxels from the current structural plane as seed points, and repeat steps S33-S34 until there are no more supervoxels to select in the current structural plane. At this point, the division of the structural plane is complete.

[0124] It is understandable that step S35 requires traversing all super-void blocks of the current structural plane without repetition, in order to avoid repeating the calculation or traversal steps for the same super-void blocks, until there are no super-void blocks to select in the current structural plane. At this point, the division of the current structural plane is completed, and the final structural plane will correspond to a fitted plane equation for that structural plane.

[0125] It is understandable that steps S32 to S35 involve region growing of a single structural surface, i.e., identifying and extracting a single structural surface from the supervoxel model;

[0126] S36. Continue to randomly select a supervoxel block that is not included in the current structural plane from the list of supervoxel blocks in the supervoxel model, and use it as the new current structural plane. Select the supervoxel block as the seed point and repeat steps S32-S35 until all supervoxel blocks are divided into independent structural planes.

[0127] The structural surfaces and their corresponding fitted plane equations obtained in the above steps are stored in lists. Finally, the list of structural surfaces and the list of fitted plane equations together form the structural surface model. The structural surface model generated in this embodiment is as follows: Figure 6 As shown in (a) in the figure; the orientation contour map corresponding to the structural surface model is as follows. Figure 6 As shown in (b);

[0128] It is understandable that step S36 is an extension of the single-structure surface region growth in steps S32 to S35, that is, step S36 realizes the identification and extraction of all structural surfaces of the supervoxel model.

[0129] Traditional region growing algorithms rely on a single threshold for growth and are sensitive to threshold changes. They are prone to oversegmentation or undersegmentation in sparse or undulating point cloud regions, and have low accuracy in structural surface recognition. The dual-criteria constrained region growing algorithm proposed in this embodiment can simultaneously consider the local geometric planarity of the point cloud and the overall planarity of the structural surface, which greatly improves the problem of oversegmentation or undersegmentation that traditional algorithms are prone to.

[0130] S4. Calculate the attitude of each structural surface in the structural surface model;

[0131] Specifically, this includes: extracting structural surfaces from the structural surface list of the structural surface model, and calculating the attitude of each structural surface according to the following formula; the attitude includes dip direction. and tilt angle :

[0132] ;

[0133] ;

[0134] in, These are the components of the normal vector n of the fitting plane of the structural surface in the X, Y, and Z directions, respectively.

[0135] The calculated orientation results of some structural planes in this embodiment are shown in Table 1.

[0136] Table 1 Structural surface orientation results

[0137]

[0138] This embodiment provides a rapid and automatic identification method for rock mass structural surfaces based on 3D point clouds. By first segmenting the 3D point cloud of the rock mass into hypervoxel blocks with planar characteristics, the number of seed points required for subsequent region growing algorithms is effectively reduced, greatly improving processing efficiency. At the same time, an improved dual-criteria constrained region growing algorithm is used to accurately identify structural surfaces. Based on the analysis of the attitude (dip angle) of the structural surfaces generated in this embodiment, a more accurate and realistic rock mass condition can be obtained, providing a more reliable basis for engineering sites. This method compensates for the shortcomings of traditional measurement methods in terms of efficiency, accuracy, and safety, and promotes the development of engineering rock mass geological surveys towards intelligence and informatization.

[0139] Furthermore, this embodiment performs multiple high-precision normal vector estimations (including normal vector estimation of the original point cloud, planar normal vector estimation of the hypervoxel block, and incremental planar normal vector estimation of the structural surface), ensuring that the fitted planar normal vector of the extracted structural surface is always close to the normal vector of the original point cloud. This helps maintain the accuracy of the attitude calculation and enables the processing of large-scale rock mass point clouds and structural surface identification.

[0140] Example 2

[0141] refer to Figure 7 This embodiment provides a rapid automatic identification device for rock mass structural surfaces based on three-dimensional point clouds, which is based on the rapid automatic identification method for rock mass structural surfaces based on three-dimensional point clouds described in Embodiment 1, and includes the following units:

[0142] The point cloud model generation unit is used to obtain a three-dimensional point cloud model of the engineering rock mass in the survey area;

[0143] The supervoxel model construction unit is used to adaptively cluster the 3D point cloud model through the supervoxel segmentation algorithm to obtain a list of supervoxel blocks and the point set to which the supervoxel blocks belong. Then, the supervoxel normal vector estimation RANSAC algorithm is used to perform plane fitting on each supervoxel block in the list of supervoxel blocks to obtain its corresponding fitting plane, forming a list of supervoxel block fitting planes. The supervoxel model is constructed based on the list of supervoxel blocks, the point set to which the supervoxel blocks belong, and the list of supervoxel block fitting planes.

[0144] The structural surface model construction unit is used to identify structural surfaces in the supervoxel model and obtain the structural surface model by using a dual-criteria constrained region growth algorithm.

[0145] The attitude analysis unit is used to calculate the attitude of the structural surfaces of the structural surface model.

[0146] This specification is described with reference to flowchart illustrations and / or block diagrams of methods, apparatus (systems), and computer program products according to embodiments of this specification. It will be understood that each block of the flowchart illustrations and / or block diagrams, and combinations of blocks in the flowchart illustrations and / or block diagrams, can be implemented by computer program instructions. These computer program instructions can be provided to a processor of a general-purpose computer, special-purpose computer, embedded processor, or other programmable data processing apparatus to produce a machine, such that the instructions, which execute via the processor of the computer or other programmable data processing apparatus, create a machine for implementing the flowchart illustrations and / or block diagrams. Figure 1 A process, multiple processes, and / or boxes Figure 1 Devices that specify the functions in one or more boxes.

[0147] These computer program instructions may also be stored in a computer-readable storage medium that can direct a computer or other programmable data processing device to function in a particular manner, such that the instructions stored in the computer-readable storage medium produce an article of manufacture including an instruction device, which is implemented in a process Figure 1 One or more processes and / or boxes Figure 1 The function specified in one or more boxes.

[0148] These computer program instructions may also be loaded onto a computer or other programmable data processing equipment to cause a series of operational steps to be performed on the computer or other programmable equipment to produce a computer-implemented process, thereby providing instructions that execute on the computer or other programmable equipment for implementing the process. Figure 1 One or more processes and / or boxes Figure 1 The steps of the function specified in one or more boxes.

[0149] The parts of this invention not described in detail are prior art. It will be apparent to those skilled in the art that this invention is not limited to the details of the above exemplary embodiments, and that the invention can be implemented in other specific forms without departing from the spirit or essential characteristics of the invention. Therefore, the embodiments should be regarded as exemplary and non-limiting in all respects, and are intended to encompass all changes falling within the meaning and scope of equivalents within this invention.

Claims

1. A method for quickly and automatically identifying rock mass structural planes based on a three-dimensional point cloud, characterized in that, The method comprises the following steps: S1, obtaining a three-dimensional point cloud model of the engineering rock mass in the survey area; S2, performing adaptive clustering on the three-dimensional point cloud model by using a super voxel segmentation algorithm to obtain a super voxel block list and a point set to which the super voxel blocks belong, and performing plane fitting on each super voxel block in the super voxel block list by using a super voxel normal vector estimation RANSAC algorithm to obtain a corresponding fitting plane, thereby forming a super voxel block fitting plane list, and constructing a super voxel model according to the super voxel block list, the point set to which the super voxel blocks belong and the super voxel block fitting plane list; S3, identifying the structural surface model by using a double-criterion constraint region growing algorithm on the super voxel model; The double-criterion constraint region growing algorithm specifically comprises: S31, determining a neighborhood super voxel block set of each super voxel block in the super voxel block list of the super voxel model according to a neighborhood algorithm, wherein the position of the super voxel block is determined according to the average value of the coordinates of the points contained in the super voxel block; S32, randomly selecting a super voxel block from the super voxel block list of the super voxel model as a current structural surface, selecting the super voxel block as a seed point, taking the fitting plane of the super voxel block as the fitting plane of the current structural surface, taking the fitting plane equation of the super voxel block as the fitting plane equation of the current structural surface, and calculating the inlier rate of the current fitting plane according to an inlier rate formula based on the point set to which the current structural surface belongs and the fitting plane equation of the current structural surface, and taking the inlier rate of the current fitting plane as the inlier rate of the fitting plane of the current structural surface; S33, sequentially traversing all neighborhood super voxel blocks in the neighborhood super voxel block set of the seed point, and including the neighborhood super voxel blocks and the point set to which the neighborhood super voxel blocks belong in the current structural surface, which meet the double-criterion constraint condition; S34, if the double-criterion constraint condition is met in step S33, updating the fitting plane equation of the current structural surface, calculating the inlier rate of the fitting plane of the seed point neighborhood super voxel block according to an inlier rate formula based on the point set to which the current structural surface belongs and the fitting plane equation of the seed point neighborhood super voxel block, and if the inlier rate of the fitting plane of the seed point neighborhood super voxel block is greater than the inlier rate of the fitting plane of the current structural surface, taking the inlier rate of the fitting plane of the seed point neighborhood super voxel block as the inlier rate of the fitting plane of the current structural surface, and updating the fitting plane equation of the structural surface to the fitting plane equation of the seed point neighborhood super voxel block; S35, continuing to select non-repeated super voxel blocks from the current structural surface as seed points, and repeating steps S33-S34 until there is no super voxel block available for selection in the current structural surface, at which time the division of the structural surface is completed; S36, continuing to randomly select super voxel blocks that are not included in the current structural surface from the super voxel block list of the super voxel model as new current structural surfaces, and selecting the super voxel blocks as seed points, and repeating steps S32-S35 until all super voxel blocks are divided into independent structural surfaces; S4, calculating the occurrence of the structural surfaces of the structural surface model respectively.

2. The method of claim 1, wherein, The super voxel segmentation algorithm in step S2 specifically comprises: S21, determining a neighborhood point set of each point in the three-dimensional point cloud model according to a neighborhood algorithm, and calculating the similarity metric of each point in the three-dimensional point cloud model and its neighborhood points respectively; S22, setting an estimation parameter The initial value of the similarity metric value is the median of the similarity metric values of each point in all three-dimensional point cloud models and the neighborhood points. S23, voxel down-sampling the three-dimensional point cloud model according to the super voxel resolution R to obtain a down-sampled point cloud model, and counting the number of points of the down-sampled point cloud model as DK; S24, initializing each point in the three-dimensional point cloud model as an independent super voxel block, selecting the point as a representative point of each super voxel block and initializing a point set to which the super voxel block belongs; initializing an adjacent super voxel block list of each super voxel block as a sequence of super voxel blocks constructed by the neighborhood points of each point in the three-dimensional point cloud model; S25, traversing all current super voxel blocks, for each super voxel block, traversing the adjacent super voxel block list of the super voxel block, if the following formula is satisfied, merging the adjacent super voxel block into the super voxel block, merging the point set to which the adjacent super voxel block belongs into the point set to which the super voxel block belongs, and merging the adjacent super voxel block list: ； wherein, is the number of points comprised by the adjacent super voxel block; and represent the representative points of the current super voxel block and its adjacent super voxel block, respectively; represents and is the similarity measure of two points. S26, estimating parameters Doubling, continue to step S25 until the number of super voxel blocks is DK, and finally obtain the list of super voxel blocks with the number of DK and the corresponding point set to which the super voxel blocks belong.

3. The method of claim 2, wherein, The similarity of each point in the three-dimensional point cloud model and its neighborhood points is calculated according to the following formula: ； where p and N represent a point and its neighborhood points in the three-dimensional point cloud model, respectively; represents a similarity measure between points p and N; represents a normal vector of point p; represents a normal vector of point N; represents the Euclidean distance between points p and N; and R is the super voxel resolution.

4. The method of claim 1, wherein, The super voxel normal vector estimation RANSAC algorithm specifically includes: S27, randomly selecting a point in the point set to which the super voxel block belongs, using the spatial coordinates and normal vector of the point to solve, and determining a candidate plane; S28, calculating the inlier rate of the candidate plane according to an inlier rate formula The inlier rate formula is as follows: ； wherein point p is a point belonging to the point set s to which the super voxel block belongs; is the number of points in the point set s to which the super voxel block belongs; n is the normal vector of the candidate plane; d is the intercept of the candidate plane; is the inlier threshold; indicator function denotes 1 when the condition is satisfied, and 0 when the condition is not satisfied. S29, repeating steps S27-S28 until the number of repetitions reaches a preset maximum number or there is no point in the super voxel block, at this time, selecting the candidate plane with the largest inlier rate as the fitting plane of the super voxel block, that is, the plane equation of the candidate plane with the largest inlier rate is the plane equation of the fitting plane of the super voxel block; at the same time, the normal vector of the candidate plane with the largest inlier rate is used as the fitting plane normal vector of the super voxel block.

5. The method of claim 4, wherein, In step S27, the candidate plane is determined according to the following formula: ； Wherein, n is the normal vector of a point in the point set to which the super voxel block belongs; d is the intercept solved according to the spatial coordinates of the point; X represents the coordinate variable belonging to the point of the plane.

6. The method of claim 1, wherein, The double criterion constraint condition in step S33 is as follows: ； wherein is the fitted plane normal vector of the seed point; is the fitted plane normal vector of the seed point neighborhood super voxel patch; is the fitted plane normal vector of the current structure plane; denotes the norm of a vector; is a first decision threshold; is a second decision threshold.

7. The method of claim 1, wherein, The step S4 specifically comprises: extracting a structural plane in a structural plane list of the structural plane model, and calculating the occurrence of the structural plane according to the following formula respectively; the occurrence comprises a dip direction and a dip angle : ； ； wherein are the components of the normal vector n of the fitting plane of the structural surface in the X, Y, Z directions, respectively.

8. The method of claim 1, wherein, The structure of the three-dimensional point cloud model includes point cloud spatial coordinates and point cloud normal vectors.

9. A device for quickly and automatically identifying rock mass structural planes based on a three-dimensional point cloud, based on the method for quickly and automatically identifying rock mass structural planes based on a three-dimensional point cloud according to any one of claims 1-8, characterized in that, It includes the following units: A point cloud model generation unit is configured to obtain a three-dimensional point cloud model of an engineering rock mass in a survey area; A super voxel model construction unit is configured to perform adaptive clustering on the three-dimensional point cloud model by a super voxel segmentation algorithm to obtain a super voxel block list and a point set to which each super voxel block belongs, and perform plane fitting on each super voxel block in the super voxel block list by a super voxel normal vector estimation RANSAC algorithm to obtain a corresponding fitting plane, thereby forming a super voxel block fitting plane list, and construct a super voxel model according to the super voxel block list, the point set to which each super voxel block belongs, and the super voxel block fitting plane list; A structural surface model construction unit is configured to perform structural surface identification on the super voxel model by a double criterion constraint region growing algorithm to obtain a structural surface model; An occurrence analysis unit is configured to calculate the occurrence of each structural surface of the structural surface model.

Images