Rapid identification method of rock mass structure based on laser point cloud based on RANSAC algorithm

Abstract

The present invention provides a method for rapid identification of rock mass structure by laser point cloud based on RANSAC algorithm, comprising step 1: using a three-dimensional laser scanner to scan the rock mass structure surface to obtain the point cloud data of the rock mass structure surface; step 2: pre-processing the above point cloud data to obtain the pre-processed point cloud data, establishing a data index for the pre-processed point cloud, using the nearest neighbor method to solve the normal vector for each point, and obtaining the point cloud containing the normal vector; step 3: using the RANSAC algorithm to extract the structure surface point cloud from the point cloud containing the normal vector, and obtaining the optimized structure surface point set; step 4: using the regional growing algorithm to identify the structure surface according to the optimized structure surface point set and calculate the inclination and inclination of the structure surface. The jointed rock mass model and joint inclination and inclination analysis generated according to the present invention can obtain a more accurate and realistic rock mass condition, and provide a more reliable basis for the engineering site.

Description

Technical Field

[0001] The invention relates to the field of geotechnical engineering, and in particular to a laser point cloud rock mass structure rapid identification method based on a RANSAC algorithm. Background Art

[0002] Rock mass refers to a geological body with discontinuity, heterogeneity and anisotropy, which is composed of various types of rocks with weak structural surfaces within a certain engineering scope. Rock mass is formed in a long geological history and has a certain structure and texture. Rock mass is composed of various rocks, and has been destroyed and transformed by various internal and external geological forces such as structural changes, weathering and unloading during its formation. Therefore, rock mass is often cut by various structural surfaces, making it a discontinuous medium with multiple fractures.

[0003] 3D laser scanning technology is a technology that uses the principle of laser ranging to quickly reconstruct the 3D model of the measured target by recording the 3D coordinates, reflectivity, texture and other information of a large number of dense points on the surface of the measured object. Since the 3D laser scanning system can densely obtain a large number of data points of the target object, compared with the traditional single-point measurement, 3D laser scanning technology is also called a revolutionary technological breakthrough from single-point measurement to surface measurement. Scanning the surface of the rock mass using 3D laser scanning can obtain point cloud data of the rock mass surface, reconstruct the 3D grid model of the rock mass based on the point cloud data, and accurately identify the structural surface of the rock mass surface. The plane fitting technology of point cloud mainly includes least squares fitting and random sampling consensus method (RANSAC). Among them, the traditional least squares fitting only considers the error in the observation vector and ignores the error in the coefficient matrix. The accuracy of the fitted plane is poor and it is not robust. The random sampling consensus method (RANSAC) is an algorithm that determines the initial model by random sampling and iterates multiple times to seek the best parameter estimation based on the distance threshold. It can fit a plane model with a specified regularity when the data contains a large number of noise points.

[0004] The existing rock mass structural surface recognition technology is mainly based on image recognition, which obtains the excavation surface image and slices it, and uses the pre-trained rock mass structural surface recognition model for recognition. However, this method has extremely high requirements on the quality of the acquired image, cannot accurately eliminate redundant noise points, and has low structural surface recognition accuracy. Summary of the invention

[0005] The technical problem to be solved by the present invention is to provide a laser point cloud rock structure rapid identification method based on RANSAC algorithm in view of the shortcomings of the prior art, which can be used for identification under construction conditions with higher identification accuracy, thus providing more accurate support for on-site geological analysis.

[0006] The laser point cloud rock mass structure rapid identification method based on RANSAC algorithm includes the following steps:

[0007] Step 1: Use a 3D laser scanner to scan the rock mass structure surface to obtain point cloud data of the rock mass structure surface;

[0008] Step 2: Preprocess the above point cloud data to obtain preprocessed point cloud data, use the kd-tree method to establish a point cloud data index for the preprocessed point cloud, and use the nearest neighbor method to solve the normal vector for each point cloud to obtain a point cloud containing the normal vector;

[0009] The pre-processing steps include:

[0010] Step 2.1: Use CloudCompare to extract the ROI area of ​​the original rock structure surface point cloud and remove interference objects including underground pipelines and vehicles;

[0011] Step 2.2: Traverse the point cloud and calculate the average distance between each point cloud and its n neighboring points; then, calculate the standard deviation of the average distance between the point cloud and its neighboring points s and mean m , distance threshold d max The calculation formula is d max =μ+α×σ ,in α is the scale factor, i.e. the number of domain points; finally, traverse the point cloud again and remove points with an average distance greater than d max point.

[0012] Step 3: extracting the structural surface point cloud from the point cloud data containing the normal vector by using the RANSAC method, establishing a structural surface point cloud set, optimizing the parameters of the structural surface point cloud by using a nonlinear model method for the point clouds in the structural surface point cloud set, and obtaining an optimized structural surface point cloud data set;

[0013] Step 4: Use the region growing algorithm to identify the structural surface based on the optimized structural surface point set and calculate the structural surface inclination and dip angle.

[0014] Preferably, the step of establishing a point cloud data index using the kd-tree method in step 2 includes:

[0015] Step 2.3: Sort the preprocessed point cloud on the three coordinate axes respectively, calculate the standard deviation of the preprocessed point cloud data on the three dimensions respectively and sort them according to the size of the standard deviation, and take the dimension with the largest standard deviation as the first sorting coordinate axis. The order of the standard deviation size is the replacement order of the sorting coordinate axis;

[0016] Step 2.4: Calculate the median of all point cloud data on the three sorting coordinate axes, and put the median into the current partition node as the root node of the binary tree, partition the point cloud less than or equal to the median into the left subtree of the binary tree, and partition the point cloud greater than the median into the right subtree of the binary tree;

[0017] Step 2.5: After completing the division, change the sorting coordinate axis according to the order in step 2.3;

[0018] Step 2.6: Repeat steps 2.4 to 2.5 until the kd-tree index of the point cloud data is established.

[0019] Preferably, step 3 is as follows:

[0020] Step 3.1: Randomly extract 4 non-collinear sample data from the point cloud data containing the normal vector, calculate the transformation matrix H, denoted as model M, and the calculation formula of matrix H is: , where (x, y) represents the corner position of the target image, (x', y') represents the corner position of the scene image, and s is the scale parameter;

[0021] Step 3.2: Calculate the projection error between all point cloud data containing normal vectors and the model M. If it is less than a threshold, add it to the internal point set Q. The threshold is determined according to the size of the data noise.

[0022] Step 3.3: When the number of elements in the interior point set Q is greater than the optimal interior point set Q-best, update Q-best=Q and update the maximum number of iterations. K , ,in p is the confidence level, w is the ratio of the elements contained in the internal point set Q to the elements in all point cloud data containing normal vectors; Step 3.4: If the number of iterations is greater than K, exit; otherwise, increase the number of iterations by 1 and repeat steps 3.1 to 3.3.

[0023] Further preferably, in step 3.1, the scale parameter is taken as the standard deviation of the point cloud data.

[0024] Further preferably, the threshold in step 3.2 is taken as 0.01.

[0025] Further preferably, the confidence in step 3.3 p Take 0.995.

[0026] Preferably, step 4 is as follows:

[0027] Step 4.1: Randomly select a point in the point cloud area as the initial point;

[0028] Step 4.2: Calculate the angle of the normal vector of the initial point and the angle of the normal vector of all points in the neighborhood around the initial point;

[0029] Step 4.3: If the absolute value of the difference between the normal vector angles of the points in the neighborhood around the initial point and the initial point is less than the flatness constraint threshold e , then merge the point into the point cloud area where the initial point is located;

[0030] Step 4.4: Find a new initial point among the points merged into the point cloud region until no new points are merged into the point cloud region, thereby obtaining a merged set of all point cloud regions;

[0031] Step 4.5: Fit the normal of the plane where the structural surface is located based on all points in each merged set, identify the structural surface and calculate the inclination and dip angle of the structural surface.

[0032] Further preferably, in step 4.3 e The value is 0.1. The beneficial effects of the present invention are:

[0033] The present invention proposes a laser point cloud rock structure rapid identification method based on RANSAC algorithm, uses a three-dimensional laser scanner to scan the rock structure surface, obtains point cloud data of the rock structure surface, pre-processes the acquired point cloud data to obtain pre-processed point cloud data, establishes a data index for the pre-processed point cloud, uses the nearest neighbor method to solve the normal vector for each point, obtains the point cloud containing the normal vector, uses the RANSAC algorithm to extract the structural surface point cloud from the point cloud containing the normal vector, obtains the optimized structural surface point set, and uses the regional growing algorithm to identify the structural surface according to the optimized structural surface point set and calculate the structural surface inclination and dip angle.

[0034] The jointed rock mass model and joint inclination angle analysis generated by the present invention can more accurately eliminate redundant noise points, obtain more accurate and realistic rock mass conditions, and achieve the purpose of improving the structural surface recognition precision and recognition accuracy. In addition, the rock mass structural surface recognition method can be identified under construction conditions and does not require on-site environmental quality, providing accurate support for on-site geological analysis. BRIEF DESCRIPTION OF THE DRAWINGS

[0035] Figure 1 This is a flow chart of the rapid identification method of rock structure from laser point cloud based on RANSAC algorithm. DETAILED DESCRIPTION

[0036] In order to facilitate the understanding of the present invention, the present invention is fully and in detail described below in conjunction with the accompanying drawings.

[0037] like Figure 1 As shown, a laser point cloud rock structure rapid identification method based on RANSAC algorithm includes:

[0038] Step 1: Use a 3D laser scanner to scan the rock structure surface to obtain point cloud data of the rock structure surface.

[0039] Step 2: Preprocess the above point cloud data to obtain preprocessed point cloud data, use the kd-tree method to establish a point cloud data index for the preprocessed point cloud, and use the nearest neighbor method to solve the normal vector for each point cloud to obtain a point cloud containing the normal vector.

[0040] Step 2.1: Use CloudCompare to extract the ROI area of ​​the original rock structure surface point cloud and remove interference objects including underground pipelines and vehicles.

[0041] Step 2.2: Traverse the point cloud and calculate the average distance between each point cloud and its n neighboring points; then, calculate the standard deviation of the average distance between the point cloud and its neighboring points s and mean m , distance threshold d max The calculation formula is d max =μ+α×σ ,in s is the scale factor, i.e. the number of domain points; finally, traverse the point cloud again and remove points with an average distance greater than d max point.

[0042] Step 2.3: Sort the preprocessed point cloud on the three coordinate axes respectively, calculate the standard deviation of the preprocessed point cloud data on the three dimensions respectively and sort them according to the size of the standard deviation, and take the dimension with the largest standard deviation as the first sorting coordinate axis. The order of the standard deviation size is the replacement order of the sorting coordinate axis;

[0043] Step 2.4: Calculate the median of all point cloud data on the three sorting coordinate axes, and put the median into the current partition node as the root node of the binary tree, partition the point cloud less than or equal to the median into the left subtree of the binary tree, and partition the point cloud greater than the median into the right subtree of the binary tree;

[0044] Step 2.5: After completing the division, change the sorting coordinate axis according to the order in step 2.3;

[0045] Step 2.6: Repeat steps 2.4 to 2.5 until the kd-tree index of the point cloud data is established.

[0046] Step 3: Use the RANSAC method to extract the structural surface point cloud from the point cloud data containing the normal vector, establish a structural surface point cloud set, and use a nonlinear model method to optimize the parameters of the structural surface point cloud in the point cloud set to obtain an optimized structural surface point cloud data set.

[0047] Step 3.1: Randomly extract 4 non-collinear sample data from the above point cloud data containing normal vectors, calculate the transformation matrix H, denoted as model M, and the calculation formula of matrix H is: , where (x, y) represents the corner position of the target image, (x', y') represents the corner position of the scene image, and s is the scale parameter, which is generally taken as the standard deviation of the point cloud data.

[0048] Step 3.2: Calculate the projection error between all point cloud data containing normal vectors and the model M. If it is less than the threshold, add it to the internal point set Q. The threshold is determined according to the size of the data noise and is generally taken as 0.01.

[0049] Step 3.3: When the number of elements in the interior point set Q is greater than the optimal interior point set Q-best, update Q-best=Q and update the maximum number of iterations. K , ,in p is the confidence level, which is generally taken as 0.995. w It is the ratio of the elements contained in the internal point set Q to the elements in all point cloud data containing normal vectors.

[0050] Step 3.4: If the number of iterations is greater than K , then exit; otherwise, increase the number of iterations by 1 and repeat steps 3.1 to 3.3.

[0051] Step 4: Use the region growing algorithm to identify the structural surface based on the optimized structural surface point set and calculate the structural surface inclination and dip angle.

[0052] Step 4.1: Randomly select a point in the point cloud area as the initial point;

[0053] Step 4.2: Calculate the angle of the normal vector of the initial point and the angle of the normal vector of all points in the neighborhood around the initial point;

[0054] Step 4.3: If the absolute value of the difference between the normal vector angles of the points in the neighborhood around the initial point and the initial point is less than the flatness constraint threshold e , then merge the point into the point cloud area where the initial point is located. e The general value is 0.1.

[0055] Step 4.4: Find a new initial point among the points merged into the point cloud region until no new points are merged into the point cloud region, thereby obtaining a merged set of all point cloud regions;

[0056] Step 4.5: Fit the normal of the plane where the structural surface is located based on all points in each merged set, identify the structural surface and calculate the inclination and dip angle of the structural surface.

Claims

1. A laser point cloud rock structure rapid identification method based on RANSAC algorithm, characterized by: The following steps are involved: Step 1: Use a 3D laser scanner to scan the rock mass structure surface to obtain point cloud data of the rock mass structure surface; Step 2: preprocessing the point cloud data to obtain preprocessed point cloud data; using the kd-tree method to establish a point cloud data index for the preprocessed point cloud data; Use the nearest neighbor method to solve the normal vector for each point cloud to obtain point cloud data containing the normal vector; The pre-processing comprises: Step 2.1: Use CloudCompare to extract the ROI area of ​​the original rock structure surface point cloud data and remove interference objects including underground pipelines and vehicles; Step 2.2: Traverse the point cloud data and calculate the average distance between each point cloud and its n neighboring points; calculate the standard deviation σ and mean μ of the average distance between the point cloud and its neighboring points, and calculate the distance threshold , =μ+α×σ, where α is the proportional coefficient; traverse the point cloud and remove points whose average distance is greater than dmax; The kd-tree method is used to establish a point cloud data index, including: Step 2.11: sorting the preprocessed point cloud data on the three coordinate axes respectively; including: calculating the standard deviation of the preprocessed point cloud data on the three dimensions respectively, and sorting them according to the size of the standard deviation, taking the dimension with the largest standard deviation as the first sorting coordinate axis, and the order of the standard deviation size is the replacement order of the sorting coordinate axis; Step 2.12: Calculate the median of all point cloud data on the three sorting coordinate axes, and put the median into the current partition node as the root node of the binary tree, partition the point cloud less than or equal to the median into the left subtree of the binary tree, and partition the point cloud greater than the median into the right subtree of the binary tree; Step 2.13: Change the sorting axis according to step 2.11; Step 2.14: Repeat steps 2.12 to 2.13 until the kd-tree index of the point cloud data is established; Step 3: extracting the structural surface point cloud from the point cloud data containing the normal vector by using the RANSAC method, and establishing a structural surface point cloud set; optimizing the parameters of the structural surface point cloud by using a nonlinear model method for the point cloud in the structural surface point cloud set, and obtaining an optimized structural surface point cloud data set; The step 3 comprises: Step 3.1: Randomly extract 4 non-collinear sample data from the point cloud data containing the normal vector, calculate the transformation matrix H, and record it as model M; the calculation formula of the matrix H is: ； Among them, (x, y) represents the corner position of the target image, (x', y') is the corner position of the scene image, and s is the scale parameter; Taking the scale parameter as the standard deviation of the point cloud data; Step 3.2: Calculate the projection error between all point cloud data containing normal vectors and the model M. If it is less than a preset threshold, add it to the internal point set Q; the preset threshold is determined according to the size of the data noise; Step 3.3: When the number of elements in the interior point set Q is greater than the optimal interior point set Q-best, update Q-best=Q; and update the maximum number of iterations K: ； Among them, p is the confidence, w is the ratio of the elements contained in the internal point set Q to the elements in all point cloud data containing normal vectors; Step 3.4: If the number of iterations is greater than K, exit; otherwise, increase the number of iterations by 1 and repeat steps 3.1 to 3.3; Step 4: Use the region growing algorithm to identify the structural surface based on the optimized structural surface point set, and calculate the structural surface inclination and inclination; The step 4 comprises: Step 4.1: Randomly select a point in the point cloud area as the initial point; Step 4.2: Calculate the angle of the normal vector of the initial point and the angle of the normal vector of all points in the neighborhood around the initial point; Step 4.3: If the absolute value of the difference between the angle between the normal vector of the initial point and the points in the neighborhood around the initial point is less than the flatness constraint threshold ε, the points in the neighborhood around the initial point are merged into the point cloud area where the initial point is located; Step 4.4: Find a new initial point in the merged point cloud area until no new points are merged into the point cloud area, and obtain a merged set of all point cloud areas; Step 4.5: Fit all points in each merged set to obtain the normal of the plane where the structural surface is located, identify the structural surface and calculate the inclination and dip angle of the structural surface.

2. The laser point cloud rock mass structure rapid identification method based on RANSAC algorithm according to claim 1 is characterized in that: In step 3.2, the threshold is 0.

01.

3. The laser point cloud rock mass structure rapid identification method based on RANSAC algorithm according to claim 1 is characterized in that: In step 3.3, the confidence level p is 0.

995.

4. The laser point cloud rock mass structure rapid identification method based on RANSAC algorithm according to claim 1 is characterized in that: In step 4.3, the ε value is 0.1.

Images