Rock mass structural surface intelligent identification method based on three-dimensional point cloud

Abstract

The invention provides a rock mass structural plane intelligent identification method based on three-dimensional point clouds, and belongs to the technical field of rock mass structural plane intelligent identification, and the method comprises the steps: selecting a target collection point, constructing an initial point cloud set, calculating an initial normal vector, setting a spherical accumulator, and rotating the point clouds in unstructured point clouds formed by a rock mass structural plane; obtaining an accurate normal vector; according to an improved region growing algorithm, obtaining a better number of structural surfaces, and obtaining a single structural surface of a plurality of rock mass point clouds; and traversing a single structural surface of each rock mass point cloud, replacing a distance function with an angle error quadratic sum according to an improved K-means + + clustering method, determining an optimal clustering number, and completing the intelligent identification method for the rock mass structural surface. According to the method, the high precision of the normal vector at the edge junction of the structural plane is realized under the normal vector smoothing phenomenon, the high standard of K-means clustering iteration and the high accuracy of the optimal cluster number of selected clustering are realized, and the problems of density anisotropy and the marginalization effect of the point cloud structural plane are solved.

Description

Technical Field

[0001] The invention belongs to the technical field of intelligent identification of rock mass structural surfaces, and in particular relates to an intelligent identification method of rock mass structural surfaces based on three-dimensional point clouds. Background Art

[0002] The rock mass structural surface has the morphological characteristics of random distribution and anisotropy, which cuts the rock into intricate unequal structures in the natural state. The structural surface usually refers to one or more discontinuous surfaces in the rock mass or rock formation, including cracks, faults, joints, etc. Because it defines the properties of the weak surface in the rock mass and contains a large number of geometric parameters such as inclination, dip, spacing, etc., most rock slope stability analysis relies on the information and characterization of the structural surface. Not only that, the precision of the structural surface data is also the basis of large-scale projects such as tunnel support, surrounding rock classification, underground seepage, and mine blasting.

[0003] At present, most of the characteristic information of rock mass structural surfaces is obtained by analyzing and calculating point cloud data obtained after being collected and processed by three-dimensional laser scanning and photogrammetry technology. Obviously, the quality of point cloud data directly determines the results of subsequent analysis and calculation, and the point cloud normal vector is one of the important evaluation criteria. Many algorithms and software also rely on point cloud normal vectors as basic data, such as point-based rendering, surface reconstruction, 3D segmented plane reconstruction and shape extraction. Therefore, it can be seen that the normal vector, as a basic three-dimensional feature, is an important prerequisite for three-dimensional space calculation, and its accuracy has a great influence on the calculation results. In the prior art, the regression method, namely PCA (principal component analysis), is widely used. The prior art proposes effective calculation regression fitting to estimate the normal vector of the approximate tangent plane. In addition, the prior art also has special plane fitting methods such as spheres, quadratic surfaces and MLS (moving least squares method). Such methods are efficient and robust to noise, but their defects are that in the case of highly curved surfaces (such as sharp edges or complex geometric shapes), the normal vector is distorted, and the smoothing of sharp edges causes distortion. Sharp and deflected features; In order to solve the above problems, the existing technology calculates the three-dimensional geometric structure based on the Voronoi diagram (Voronoi diagram) and retains the sharp features, but the noise resistance is weak and it is sensitive to the boundary effect of the point cloud. The existing technology improves and proposes Voronoi-PCA (Voronoi-Principal Component Analysis) to optimize the calculation results by adjusting the smoothness during fitting; In recent years, the robust noise calculation and its corresponding scale parameterization are combined with a similar RANSAC (random sampling consistency) method to solve the noise and sharp features well, but there is still a deficiency that the edge density change cannot be solved well. In response to the above problems, the existing technology proposes a robust random Hough transform method. This method is based on the principle of numerical statistics. By constructing a triple in the three-dimensional point cloud data set and combining a certain confidence interval to select the corresponding normal vector in the "random" normal distribution peak, it is not only highly robust to noise, but also can obtain the corresponding accurate normal vector when calculating the sharp edge area of ​​the rock structure surface.

[0004] After obtaining the normal vector, the eigenvalue can be used to further divide the cluster to obtain the dominant group structure surface. There are also corresponding clustering methods for the normal vector data type. Fuzzy K-means clustering was proposed earlier. It directly clusters the direction information of each face of the grid model, iterates the Fisher-K value to eliminate abnormal points, and selects the cluster size with the largest K value. In addition, the prior art further improves the method, using the Firefly algorithm to determine the optimal position of the initial clustering center of the fuzzy K-means algorithm, and applies it to the three-dimensional point cloud; the prior art proposes a spectral clustering algorithm, which converts the point cloud normal vector into a K-dimensional space, where the point cloud coordinates are obtained by the eigenvector of the similarity matrix normalized by the normal vector, and then the original point cloud is matched with the corresponding transformed point cloud for clustering; the prior art proposes a recheck clustering method to estimate the average direction of the structural surface, in the first round, the point cloud is divided and the maximum and minimum clusters are defined, and then the PCA (principal component analysis) is used to calculate the normal vector of each cluster, and finally the dominant group is identified by K-means, where the optimal K value is determined by the MDL (minimum description length) statistic; the prior art proposes a method for automatically clustering the mesh surface into the structural surface using an improved K-means clustering algorithm, based on finding the initial cluster using sample density and calculating the average contour effective index to estimate the optimal value of k. Finally, the plane parameters are determined by the RANSAC (random sampling consensus) random sample consensus. Different from the traditional method of using clustering to segment structural surfaces, a region growing algorithm has been derived in recent years to obtain all structural surfaces from 3D point clouds. This method expands the region by selecting seed points from the point cloud and then adding neighborhood points that meet certain growth criteria. The advantage of this method is that it can completely extract discontinuous structural surfaces, but the computational efficiency will be greatly reduced when the number of point clouds exceeds 1 million. To address this problem, the existing technology proposes two solutions: gridding data and optimizing the region growing algorithm to a certain extent. Among them, gridding converts the original point cloud into grid data with evenly spaced rows and columns, where each point can be accessed using row and column indexes to save data retrieval and processing time in the point cloud. Secondly, more effective growth criteria are used as expansion conditions, and only four neighborhood points are used to estimate the curvature. Compared with conventional region growing, the calculation time is reduced by 97.66% and the occurrence of the automatically segmented structural surface is closer to on-site manual measurement.

[0005] In summary, the quality of normal vectors is a prerequisite for ideal structural surface recognition, and a clustering method that is robust to the anisotropy of point cloud data sets is also required. Therefore, the calculation of normal vectors urgently needs to overcome the difficulty of smoothing the sharp edges of rock masses, so as to achieve the computing power that can accurately obtain normal vectors when the natural occurrence of discontinuous structural surfaces in point clouds changes sharply; and for large-scale point cloud normal vector data sets, clustering and segmentation can achieve better single discontinuous structural surface extraction effects. Summary of the invention

[0006] In view of the above-mentioned deficiencies in the prior art, the present invention provides an intelligent identification method for rock structural surfaces based on three-dimensional point clouds, which solves the problem of normal vector smoothing at the edge junction of two or more structural surfaces, the inability of existing algorithms to well solve the density anisotropy and point cloud structural surface marginalization effect, and the K-means clustering problem of non-standard clustering iteration and inaccurate selection of the optimal number of clusters.

[0007] In order to achieve the above purpose, the technical solution adopted by the present invention is: a method for intelligent identification of rock mass structural surface based on three-dimensional point cloud, comprising the following steps:

[0008] S1. Filter out the noise points and outliers in the unstructured point cloud formed by the rock mass structural surface, randomly select the target acquisition point, search for the nearest point of the target acquisition point to construct an initial point cloud set, construct a triplet based on the initial point cloud set, and calculate the initial normal vector based on the triplet, set a gridded spherical accumulator, rotate the point cloud, and perform cluster average calculation on the initial normal vector to obtain the precise normal vector corresponding to the target acquisition point;

[0009] S2. Using the improved region growing algorithm, the seed point is calculated according to the precise normal vector of the corresponding target acquisition point, and the similarity region is defined by the similarity criterion. The optimal number of structural planes is obtained by setting the gradient threshold of the normal vector of the initial growth condition and drawing the number of structural planes, and a single structural plane of several rock point clouds is obtained;

[0010] S3. Traverse the individual structural surfaces of each rock mass point cloud, calculate the analytical normal vector and use it as the feature representative. Based on the feature representative, use the improved K-means++ clustering method to replace the distance function with the sum of squared angle errors to obtain the optimal number of clusters, determine the division of the dominant group structural surface, and complete the intelligent identification method of the rock mass structural surface.

[0011] The beneficial effects of the present invention are as follows: the intelligent identification method of rock mass structural surfaces constructed by the present invention changes the process from multiple target acquisition points directly to group structural surfaces, and adopts segmentation from multiple target acquisition points to a single structural surface, and then clustering identification from a single structural surface to a group structural surface, thereby avoiding the defect of normal vector distortion caused by sharp edges and density anisotropy of three-dimensional point cloud structural surfaces. The present invention improves the accuracy and efficiency of identifying rock mass structural surfaces, and has good applicability in complex and changeable environments.

[0012] Furthermore, the S1 comprises the following steps:

[0013] S101, using the KD tree index to search for the nearest points of the target collection point, and using the searched nearest points as the initial point cloud set;

[0014] S102, balancing different density weights of each point cloud in the initial point cloud set by weighting the point cloud density;

[0015] S103, randomly selecting three nonlinear points from the initial point cloud set to form a triplet, randomly selecting two three-dimensional vectors from the constructed triplet, and obtaining an initial normal vector by vector cross multiplication calculation;

[0016] S104, setting a gridded spherical accumulator to simulate the voting process, making the area of ​​each grid in the spherical accumulator the same, and setting an independent correlation distribution between the number of spherical accumulators and the number of triplets;

[0017] S105, according to the setting result of S104, using the probability empirical formula, determine the relationship between the actual value and the empirical value, and quote the Hoefding inequality into the probability empirical formula to calculate the minimum triple number inequality;

[0018] S106. According to the central limit theorem, a random variable of a process of determining the larger one between two grids with the most votes in the spherical accumulator through a confidence interval is obtained, and a disjoint inequality is used in combination with a minimum triplet number inequality, and in response to the confidence interval of the grid with the most votes in the spherical accumulator and the confidence interval of the grid with the second most votes not intersecting, triplet construction is stopped, wherein, in the triplet construction process, an initial normal vector in the spherical accumulator forms a normal distribution and generates a peak area, and the peak of the peak area is made to be located in the corresponding grid of the spherical accumulator by rotating the point cloud;

[0019] S107, performing cluster average calculation on the initial normal vectors of the peak area within a certain threshold range to obtain the precise normal vectors corresponding to the target acquisition points.

[0020] The beneficial effects of the above further scheme are: the present invention adopts a method of weighted processing of point cloud density to balance the different density weights of each point cloud, solves the anisotropy of point cloud data density, ensures that the normal distribution formed by the initial normal vector in random selection has the correct peak, and uses a rotated point cloud to reduce the error probability and improve the accuracy of the calculated normal vector.

[0021] Furthermore, the expression of the independent correlation distribution between the number of spherical accumulators and the number of triples is as follows:

[0022] (X m,t ) m∈{1,…,M}，t∈{1,…,T}

[0023] Where m represents the mth spherical accumulator, t represents the tth triplet, and X m,t represents the independent correlation distribution between variables m and t, M represents the number of spherical accumulators, and T represents the number of triplets;

[0024] The probability empirical formula is as follows:

[0025]

[0026]

[0027] Among them, P m Indicates the actual value, represents the empirical value, α represents the confidence interval, and α∈[0,1], T min represents the minimum number of triples, δ represents the distance between the actual value and the empirical value, and δ∈[0,1], ln(·) represents the logarithmic function.

[0028] Furthermore, the expression of the random variable of the process of determining the larger one is as follows:

[0029]

[0030] Among them, Q represents the random variable of the process of determining the larger one, and its distribution converges to the standard normal random variable N(0,1), p m Indicates the actual value, represents the experience value, T represents the number of triplets;

[0031] The expression of the disjoint inequality is as follows:

[0032]

[0033] Among them, m1 represents the grid with the most votes, and m2 represents the grid with the second most votes. Indicates the experience value of the grid with the most votes, Indicates the experience value of the grid with the second most votes.

[0034] The beneficial effect of the above further scheme is that the present invention adopts the probability empirical formula and the central limit theorem, which saves unnecessary computing power in the process of constructing triples and reduces the computing time of the intelligent recognition process.

[0035] Furthermore, the S2 comprises the following steps:

[0036] S201, setting a normal vector angle threshold and a curvature magnitude threshold between two points, and adding a distance threshold and an average curvature threshold of the neighborhood of the current target acquisition point as auxiliary items to define a similarity criterion;

[0037] S202, calculating the point cloud curvature using a normal vector difference method according to the precise normal vector of the target acquisition point, and selecting seed points in order by sorting the point cloud curvature from small to large;

[0038] S203, traversing the seed points in sequence, constructing an index sequence according to the seed points, adding the points in the initial point cloud set that meet the similarity criterion to the index sequence, defining the same index sequence as a similarity region, marking the processed seed points, and packaging and coloring each similarity region until all points in the unstructured point cloud are processed;

[0039] S204, calculating the number of structural surfaces according to the packaged and colored similarity regions, and drawing a structural surface number curve according to the set normal vector angle threshold between two points, curvature magnitude threshold, distance threshold, and average curvature threshold of the neighborhood of the current target acquisition point, determining the inflection point of the curve and combining visual interpretation to achieve a compromise segmentation effect;

[0040] S205. Setting the gradient threshold of the normal vector of the initial growth condition, observing the change in the number of structural planes as the curvature increases under different normal vectors under a compromise segmentation effect, obtaining the optimal number of structural planes completed by the improved regional growing algorithm, and obtaining a single structural plane of several rock point clouds.

[0041] The beneficial effects of the above further scheme are as follows: the present invention improves the region growing algorithm, sets the normal vector angle threshold between two points and the curvature size threshold as the main growing conditions, the distance threshold and the average curvature threshold of the neighborhood of the current target acquisition point as the secondary growing conditions, and determines the optimal threshold range through the growth condition and the structural surface quantity curve, thereby achieving a segmentation effect that is a compromise between over-segmentation and under-segmentation while achieving fast calculation.

[0042] Furthermore, the S3 comprises the following steps:

[0043] S301, traversing the individual structural planes of each rock mass point cloud, calculating the analytical normal vector of each individual structural plane using the principal component analysis method, and taking the analytical normal vector corresponding to each individual structural plane as a characteristic representative of the corresponding individual structural plane;

[0044] S302, using the improved K-means++ clustering method, taking the feature representatives of each single structural surface as a data set, using the squared sum of angle errors of the normal vector to replace the distance function, and using it as a point cloud normal vector clustering judgment criterion, wherein the angle judgment criterion of the squared sum of angle errors adopts the square of the cosine of the angle between two normal vectors;

[0045] S303, according to the point cloud normal vector clustering evaluation criteria, draw an angle error square sum curve, and according to the angle error square sum curve, obtain the optimal number of clusters completed by the improved K-means++ clustering method;

[0046] S304. Combining the number of optimal structural surfaces completed by the improved region growing algorithm and the optimal number of clusters completed by the improved K-means++ clustering method, the number of surfaces completed by region growing and clustering is the same, the division of the dominant group structural surfaces is determined, and the intelligent identification method of rock mass structural surfaces is completed.

[0047] The beneficial effects of the above further scheme are as follows: the present invention improves and perfects the K-means++ clustering method according to the characteristics of point cloud normal vector data, changes the distance into the normal vector angle to calculate the total data error and iterates, which is used as the point cloud normal vector clustering evaluation criterion, and combines the curve to enable the algorithm to automatically find the inflection point to determine the optimal number of clusters, thereby realizing intelligent recognition, optimizing the problem of the influence of manual subjective selection in K-means clustering, and realizing the automation of the entire algorithm process. BRIEF DESCRIPTION OF THE DRAWINGS

[0048] Figure 1 The present invention is a flow chart of the method.

[0049] Figure 2 This is a schematic diagram of the improved region growing algorithm based on random Hough transform of the present invention.

[0050] Figure 3 A normal vector graph is calculated by randomly constructing triples for this embodiment.

[0051] Figure 4 Schematic diagram of the gridded spherical accumulator of this embodiment.

[0052] Figure 5 This is a diagram of the model region growth process in this embodiment.

[0053] Figure 6 Graph showing the over-segmentation and under-segmentation of the model in this embodiment.

[0054] Figure 7 This is the completed diagram of the improved region growing algorithm for the model of this embodiment.

[0055] Figure 8 This is a screening diagram for the “better” threshold range of the model in this embodiment.

[0056] Fig. 9 This is the normal vector diagram of the region to be clustered in the model of this embodiment.

[0057] Fig.10 This is a curve chart for finding the optimal number of clusters for the SASE inflection point model in this embodiment.

[0058] Fig.11 The improved k-means++ clustering algorithm is used to complete the graph for this embodiment model. DETAILED DESCRIPTION

[0059] The specific implementation modes of the present invention are described below so that those skilled in the art can understand the present invention. However, it should be clear that the present invention is not limited to the scope of the specific implementation modes. For those of ordinary skill in the art, as long as various changes are within the spirit and scope of the present invention as defined and determined by the attached claims, these changes are obvious, and all inventions and creations utilizing the concept of the present invention are protected.

[0060] Before describing this embodiment, the following terms are explained:

[0061] PCA: principal component analysis;

[0062] RANSAC: Random Sampling Consensus;

[0063] SSE: Sum of Squared Errors;

[0064] SASE: Sum Angle of Squared Errors.

[0065] Example

[0066] like Figure 1 As shown, the present invention provides a method for intelligently identifying rock mass structural surfaces based on three-dimensional point clouds, and the implementation method thereof is as follows:

[0067] S1. Filter out the noise points and outliers in the unstructured point cloud formed by the rock mass structural surface, randomly select the target acquisition point, search for the nearest point of the target acquisition point to construct the initial point cloud set, construct the triplet according to the initial point cloud set, and calculate the initial normal vector according to the triplet, set the gridded spherical accumulator, rotate the point cloud, and perform cluster average calculation on the initial normal vector to obtain the precise normal vector corresponding to the target acquisition point. The implementation method is as follows:

[0068] S101, using the KD tree index to search for the nearest points of the target collection point, and using the searched nearest points as the initial point cloud set;

[0069] S102, balancing different density weights of each point cloud in the initial point cloud set by weighting the point cloud density;

[0070] S103, randomly selecting three nonlinear points from the initial point cloud set to form a triplet, randomly selecting two three-dimensional vectors from the constructed triplet, and obtaining an initial normal vector by vector cross multiplication calculation;

[0071] S104, setting a gridded spherical accumulator to simulate the voting process, making the area of ​​each grid in the spherical accumulator the same, and setting an independent correlation distribution between the number of spherical accumulators and the number of triplets;

[0072] S105, according to the setting result of S104, using the probability empirical formula, determine the relationship between the actual value and the empirical value, and quote the Hoefding inequality into the probability empirical formula to calculate the minimum triple number inequality;

[0073] S106. According to the central limit theorem, a random variable of a process of determining the larger one between two grids with the most votes in the spherical accumulator through a confidence interval is obtained, and a disjoint inequality is used in combination with a minimum triplet number inequality, and in response to the confidence interval of the grid with the most votes in the spherical accumulator and the confidence interval of the grid with the second most votes not intersecting, triplet construction is stopped, wherein, in the triplet construction process, an initial normal vector in the spherical accumulator forms a normal distribution and generates a peak area, and the peak of the peak area is made to be located in the corresponding grid of the spherical accumulator by rotating the point cloud;

[0074] S107, performing cluster average calculation on the initial normal vectors of the peak area within a certain threshold range to obtain the precise normal vectors corresponding to the target acquisition points.

[0075] In this embodiment, in the unstructured point cloud formed by the rock mass structural surface, pre-processing is first performed to filter out noise points and outliers in the unstructured point cloud, and target collection points are randomly selected, and target collection points are continuously randomly selected and processed.

[0076] In this embodiment, the problems of difficulty in smoothing the sharp edges of rock masses, difficulty in accurately obtaining the normal vector when the natural attitude of discontinuous structural surfaces in point clouds changes drastically, and inability to achieve good extraction effects on single discontinuous structural surfaces when clustering and segmenting large-scale point cloud normal vector data sets are solved.

[0077] like Figure 2 As shown, in this embodiment, the point cloud normal vector is calculated based on robust random Hough transform to avoid the defect of normal vector distortion caused by sharp edges and density anisotropy of the three-dimensional point cloud structure surface, and the structure surface is segmented by an improved regional growing algorithm. Finally, the improved k-means++ clustering is used to obtain the dominant structure surface with the normal vector as the single-surface representative eigenvalue for several existing structure surfaces. The present invention refines the direct clustering process from points to dominant groups, and realizes the division of the dominant group structure surface by first segmenting several single surfaces and then using a method that is more suitable for the characteristics of the data to be clustered.

[0078] In this embodiment, an initial point cloud set is established during calculation for the unstructured point cloud formed by the rock mass structural surface: a target acquisition point is randomly selected and a certain amount of nearest neighboring points are searched through the KD tree index as the initial point cloud set. At this time, three nonlinear points are randomly selected from the initial point cloud set to form a triplet, and two three-dimensional vectors A and B are selected from them, and the initial normal vector is obtained by vector cross multiplication calculation;

[0079] like Figure 3As shown, the initial normal vector is perpendicular to the plane where vectors A and B are located;

[0080] The calculation formula of the vector cross product is as follows:

[0081]

[0082] Where A represents a three-dimensional vector A=(a x ,a y ,a z ), B represents a three-dimensional vector B = (b x ,b y ,b z ), a x 、a y 、a z , b x , b y , b z Both represent vector values.

[0083] In this embodiment, the initial point cloud set will derive two calculation assumptions: the first is a common form, that is, the current target point P is far away from the sharp edge of the rock mass. Obviously, the normal vector obtained by the triple constructed by the point cloud set p to be processed under this relatively ideal regular plane is the normal vector of the current target point; the second is a special form, that is, the current target point P is close to or at the mutation point of two point cloud planes, and the drawn triplet may be on both sides of the edge. However, with its "random" normal characteristics, the highest probability is still the corresponding correct normal vector.

[0084] like Figure 4 As shown, a gridded spherical accumulator is set to simulate the "voting" process (assuming that θ represents the inclination and φ represents the inclination angle). In order to ensure the convenience of subsequent statistics and the result of averaging after subsequent normal vector clustering, the area of ​​each grid in the accumulator is made the same, and the relationship between the number of slices along the z-axis of the accumulator and the number of grids is defined, as shown in Table 1 below. Table 1 is a table showing the relationship between the number of slices along the z-axis of the accumulator and the number of grids.

[0085] Table 1

[0086]

[0087]

[0088] In this process, the more triplets there are, the more the authenticity of the normal distribution peak can be guaranteed. However, in the constructed initial point cloud set, the number of triplets that can be generated is approximately the cube of the number of points in the initial point cloud set. In order to save unnecessary computing power, the drawing of triplets can be stopped when a certain confidence level is reached based on the empirical distribution in the accumulator. The number of spherical accumulators can be set to M, the number of triplets to T, the independent correlation distribution between the number of spherical accumulators and the number of triplets can be set, the confidence value α∈[0,1] is set between the actual value and the empirical value, the distance value δ∈[0,1] between the actual value and the empirical value is set, and the relationship between the actual value and the empirical value is determined by the probability empirical formula.

[0089] The expression for the independent correlation distribution between the number of spherical accumulators and the number of triples is as follows:

[0090] (X m,t ) m∈{1,…,M}，t∈{1,…,T}

[0091] Where m represents the mth spherical accumulator, t represents the tth triplet, and X m,t represents the independent correlation distribution between variables m and t, M represents the number of spherical accumulators, and T represents the number of triplets;

[0092] The probability empirical formula is as follows:

[0093]

[0094] Among them, P m Indicates the actual value, represents the empirical value, α represents the confidence interval, and α∈[0,1], δ represents the distance between the actual value and the empirical value, and δ∈[0,1];

[0095] Then, we quoted the Hoefding inequality and substituted it into the probability empirical formula to calculate the inequality of the minimum number of triples;

[0096]

[0097] T min represents the minimum number of triplets, and ln(·) represents the logarithmic function.

[0098] In this embodiment, the unstructured point cloud data obtained by scanning the natural rock mass will have its actual quality reduced due to uncertain factors such as the scanning distance of the laser scanner equipment and the drastic changes in the rock mass structural surface. For example, point cloud voids, noise outliers, etc. will be generated. Most algorithms are affected by the point cloud density when calculating sampling, so the anisotropy of the density of the obtained point cloud data often leads to large errors.

[0099] In this embodiment, in the random selection of triples, according to the characteristics of data statistics, the normal distribution peak is always on the side of the target point, but this is obviously only true when the point cloud density is the same. When several point clouds with obvious density differences intersect each other, it is easy to know that more triplets will be established in the high-density area, resulting in a certain amplitude deviation of the target point normal vector. Therefore, a method of weighted processing of density is adopted to balance the weights of different densities of each point cloud, so that the correct normal distribution peak appears in the random selection.

[0100] In this embodiment, as the number of triplets reaches a certain number, the normal vectors in the spherical accumulator will form a normal distribution and produce a peak area. When the triplets continue to "vote" in the same grid of the spherical accumulator, it can be determined that the current grid is the ideal correct interval, so the subsequent "voting" is of little significance to a certain extent. Therefore, in order to save unnecessary computing power, the larger of the two grids with the most votes is determined by the confidence interval and will not change, and then the selection of triplets is stopped. According to the central limit theorem, the random variable Q of this process can be defined. The random variable Q converges to the standard normal random variable N(0,1) in distribution, and r∈[0,1] is set as the confidence level. For example, using a confidence level of r=95%, at this confidence level, according to the non-intersection inequality, the generation of triplets can be stopped when the grid with the most votes and the grid with the second most votes do not intersect.

[0101] The expression of the random variable of the larger determination process is as follows:

[0102]

[0103] Among them, Q represents the random variable of the process of determining the larger one, and its distribution converges to the standard normal random variable N(0,1), p m Indicates the actual value, represents the experience value, T represents the number of triplets;

[0104] The expression of the disjoint inequality is as follows:

[0105]

[0106] Among them, m1 represents the grid with the most votes, and m2 represents the grid with the second most votes. Indicates the experience value of the grid with the most votes, Indicates the experience value of the grid with the second most votes.

[0107] Since the accumulator is gridded, there will be a grid discretization effect, that is, the normal distribution peak mentioned above is not necessarily located inside the grid but falls on two or more grid boundaries, resulting in large errors in the next step of statistical clustering. Therefore, in order to make the target peak in the corresponding grid for clustering calculation, the spherical accumulator is generally rotated to simply and quickly handle this problem. And from the empirical value, it can be known that the rotation angle is adjusted to 5° and the number of rotations is set to 5 times to offset the error caused by the discretization effect, but rotating the accumulator will lead to tedious calculations in Python programming. From the implementation point of view, rotating the point cloud is more effective than rotating the accumulator itself.

[0108] In this embodiment, a method of rotating point cloud is used to make the peak value of the peak area fall within the grid of the corresponding spherical accumulator, and finally cluster average calculation is performed on the normal vector of the peak area within a certain threshold range to obtain the precise normal vector of the corresponding target acquisition point.

[0109] S2. Using the improved region growing algorithm, the seed point is calculated according to the precise normal vector of the corresponding target acquisition point, and the similarity region is defined by the similarity criterion. The optimal number of structural surfaces is obtained by setting the gradient threshold of the normal vector of the initial growth condition and drawing the number of structural surfaces, and a single structural surface of several rock point clouds is obtained. The implementation method is as follows:

[0110] S201, setting a normal vector angle threshold and a curvature magnitude threshold between two points, and adding a distance threshold and an average curvature threshold of the neighborhood of the current target acquisition point as auxiliary items to define a similarity criterion;

[0111] S202, calculating the point cloud curvature using a normal vector difference method according to the precise normal vector of the target acquisition point, and selecting seed points in order by sorting the point cloud curvature from small to large;

[0112] S203, traversing the seed points in sequence, constructing an index sequence according to the seed points, adding the points in the initial point cloud set that meet the similarity criterion to the index sequence, defining the same index sequence as a similarity region, marking the processed seed points, and packaging and coloring each similarity region until all points in the unstructured point cloud are processed;

[0113] S204, calculating the number of structural surfaces according to the packaged and colored similarity regions, and drawing a structural surface number curve according to the set normal vector angle threshold between two points, curvature magnitude threshold, distance threshold, and average curvature threshold of the neighborhood of the current target acquisition point, determining the inflection point of the curve and combining visual interpretation to achieve a compromise segmentation effect;

[0114] S205. Setting the gradient threshold of the normal vector of the initial growth condition, observing the change in the number of structural planes as the curvature increases under different normal vectors under a compromise segmentation effect, obtaining the optimal number of structural planes completed by the improved regional growing algorithm, and obtaining a single structural plane of several rock point clouds.

[0115] Since the region growing method was originally used for two-dimensional image segmentation, its goal is to divide the image into regions with certain similar features. It is a pixel-based region segmentation method that starts from one or more seed points and adds neighboring pixels to the seed point related sequence according to certain similarity criteria, and so on until all pixels are partitioned.

[0116] In this embodiment, the region growing algorithm is improved and applied to the three-dimensional point cloud without topological structure. By setting the normal vector angle threshold and the curvature size threshold between two points, and adding the distance threshold as an auxiliary item with the current neighborhood average curvature threshold, the similarity criterion is defined together, aiming to divide the point cloud into rock structure surfaces that are practical and easy to identify through several approximate criteria. According to the obtained accurate normal vector of the corresponding target acquisition point, the normal vector difference method is used here to calculate the point cloud curvature: that is, the average normal vector angle difference between the target point and its neighboring points is calculated as the corresponding curvature. And because of the characteristics of the curvature value, the point cloud with the smallest curvature is in a more regular plane point cloud concentration, so the point cloud curvature is sorted from small to large to select seed points in turn, and the seed points are traversed in turn. The points in the current neighborhood that meet the set similarity criterion are added to the index sequence formed by the seed points, that is, they grow into the same area, which can be defined as the similarity area. Finally, the processed seed points are marked, and each similarity area is packaged and colored until all points in the structureless point cloud are processed and the calculation is terminated.

[0117] like Figure 5 As shown in the figure, a neighborhood set is constructed with the target point as the center. When the characteristic parameters of the neighborhood point cloud meet the set similarity criteria, that is, they grow into the same area and all points in the unstructured point cloud are processed, complete and accurate intelligent identification of the rock structure surface is achieved.

[0118] In this embodiment, the rock mass structural surface segmentation effects obtained under different growth conditions are quite different, mainly reflected in the number of structural surface divisions, and the final segmentation result to be clustered is preferably presented as a single natural structural surface that is divided, which is convenient for the final clustering of the dominant group structural surfaces. Different numbers of structural surfaces will be obtained for different growth condition settings, which can be roughly divided into two situations: over-segmentation and under-segmentation.

[0119] like Figure 6 As shown, Figure 6(a) indicates over-segmentation, and (b) indicates under-segmentation. It is easy to see that the number of structural surfaces is small in the case of under-segmentation, and most adjacent structural surfaces are divided into the same group. However, in fact, the occurrence of adjacent surfaces is still different because the threshold of the growth condition is too large, resulting in a strong inclusiveness of the spatial properties of the point cloud at the connection of adjacent surfaces, and the adjacent surfaces with small differences are merged. As a result, the segmented surface is easy to identify but the effect is rough, and the subsequent clustering is not very operable; on the contrary, in the case of over-segmentation, it is easy to see that the number of structural surfaces is large in this case. Except for the more regular plane area, the single structural surface that should not be segmented can be easily divided into several small pieces. The reason is that the threshold of the growth condition is set too small, resulting in the segmentation result on the slightly undulating structural surface showing strong fragmentation. As a result, there is a large space for the subsequent clustering material to operate, but a large amount of redundant calculations will be generated and the segmented surface is difficult to identify. Therefore, an ideal compromise needs to be achieved for the above segmentation effect, with both better single-face recognition and an appropriate total number of structural faces for subsequent clustering. Therefore, this embodiment sets multiple groups of thresholds and calculates the number of structural faces to draw the corresponding curve graph, finds the inflection point of the curve and combines visual interpretation to achieve a compromise segmentation effect and obtain several independent rock point cloud structural faces.

[0120] In this embodiment, in the growth condition, the curvature value range is smaller than the normal vector value range, which is 0 to 1, and is more important in three-dimensional space;

[0121] like Figure 8 As shown in the figure, the normal vector threshold of the initial growth condition is set at 10°, 20°, and 30° as the gradient, and the change in the number of structural surfaces as the curvature increases under different normal vectors is observed. Figure 8 As shown in the figure, the number of faces in all three cases is decreasing, and the optimal growth situation, i.e., 10 faces, can be achieved in the curvature range of 0.11 to 0.22. However, in the current "better" range, the above-mentioned extreme segmentation of structural faces, over-segmentation and under-segmentation, appeared, and it was verified that the curvature factor has a greater relative influence. Figure 7 As shown, it is easy to know that the optimal number of structural faces of this experimental model should be 10 faces.

[0122] S3, traverse the single structural surface of each rock mass point cloud, calculate the analytical normal vector and use it as the feature representative, and according to the feature representative, use the improved K-means++ clustering method to replace the distance function with the square sum of angle errors to obtain the optimal number of clusters, determine the division of the dominant group structural surface, and complete the intelligent identification method of the rock mass structural surface. The implementation method is as follows:

[0123] S301, traversing the individual structural planes of each rock mass point cloud, calculating the analytical normal vector of each individual structural plane using the principal component analysis method, and taking the analytical normal vector corresponding to each individual structural plane as a characteristic representative of the corresponding individual structural plane;

[0124] S302, using the improved K-means++ clustering method, taking the feature representatives of each single structural surface as a data set, using the squared sum of angle errors of the normal vector to replace the distance function, and using it as a point cloud normal vector clustering judgment criterion, wherein the angle judgment criterion of the squared sum of angle errors adopts the square of the cosine of the angle between two normal vectors;

[0125] S303, according to the point cloud normal vector clustering evaluation criteria, draw an angle error square sum curve, and according to the angle error square sum curve, obtain the optimal number of clusters completed by the improved K-means++ clustering method;

[0126] S304. Combining the number of optimal structural surfaces completed by the improved region growing algorithm and the optimal number of clusters completed by the improved K-means++ clustering method, the number of surfaces completed by region growing and clustering is the same, the division of the dominant group structural surfaces is determined, and the intelligent identification method of rock mass structural surfaces is completed.

[0127] In this embodiment, after obtaining several independent rock mass point cloud structure surfaces, it is necessary to merge the structural surfaces with the same or similar occurrence to accurately depict the dominant group structural surface. Different from directly clustering multiple points, this embodiment clusters the single structural surface after extracting it through the region growing algorithm. Traverse each single structural surface and use PCA (principal component analysis) to calculate the normal vector of the single structural surface as the feature representative of the corresponding single structural surface, such as Fig. 9 shown.

[0128] Compared with the existing improved watershed clustering and DBSCAN clustering based on density clusters, hierarchical clustering and agglomerative clustering based on hierarchical structure, etc., all of them directly cluster large-scale point cloud sets according to a certain eigenvalue, but the target eigenvalue of this embodiment is different, so the clustering method needs to be changed.

[0129] In this embodiment, the K-means++ clustering method is improved. First, K-means++ is an optimization of the traditional K-means clustering algorithm, which is used to improve the method of selecting the initial center point, thereby improving the clustering quality and accelerating the convergence speed. The main idea of ​​K-means++ is to select the initial center point through a probability distribution function so that these center points are as dispersed as possible. First, randomly select a data point from the cluster set as the first center point, calculate the square of the distance between each data and the nearest selected center point and record it as D(x), and construct a probability distribution function according to the value of D(x). Select a data point as the new center point. The probability of selecting a point farther away is greater, and the probability of selecting data point x as the new center point is as follows:

[0130]

[0131] Among them, X represents the set of clusters and x represents the data point.

[0132] After continuous selection, K center points are selected. Then the standard K-means algorithm is run on these initial center points, that is, clustering is performed by assigning data points to the nearest center point and iteratively updating the center point position until the algorithm converges. In essence, it is based on improving the initial center point by continuously iterating the error criterion function SSE (Sum of Squared Errors). SSE is a commonly used evaluation index in the K-means algorithm, which is used to measure the difference between the samples in the cluster and its cluster center.

[0133] In this embodiment, the sum of squared errors is improved, and the normal vectors calculated from all structural surface areas are represented as data sets. Due to the particularity of normal vector data, the distance function is not applicable here for iterative judgment, and it is proposed to replace it with the sum of squared angle errors of the normal vectors (SASE) as the iteration criterion, where the angle judgment criterion uses the square of the cosine of the angle between the two normal vectors cosine 2 θ, θ represents the angle between the two normal vectors. Then, the SASE curve is plotted to observe that the inflection point when SASE slows down significantly with the change of the number of clusters is the optimal number of clusters.

[0134] In this embodiment, the number of dominant groups k of rock mass structural planes is usually between 1 and 10, so only SASE within this range is considered for analysis. Due to the particularity of the icosahedral model, Fig.10 As shown in the figure, it can be seen that when the number of clusters is 4 and 7 respectively, the square sum of angle errors SASE is still large, which means that when SASE is zero, that is, the optimal number of clusters K = 10 is the inflection point. Fig.11 As shown, Fig.11 Improved K-means++ clustering algorithm for the model to complete the graph.

[0135] In this embodiment, the number of structural surfaces completed by the improved region growing algorithm and the optimal number of clusters completed by the improved K-means++ clustering method are combined to obtain the same number of surfaces completed by region growing and clustering, determine the division of the structural surfaces of the dominant group, and complete the intelligent identification method of rock mass structural surfaces.

Claims

1. A method for intelligent identification of rock mass structural surfaces based on three-dimensional point clouds, characterized in that: The following steps are involved: S1. Filter out the noise points and outliers in the unstructured point cloud formed by the rock mass structural surface, randomly select the target acquisition point, search for the nearest point of the target acquisition point to construct an initial point cloud set, construct a triplet based on the initial point cloud set, and calculate the initial normal vector based on the triplet, set a gridded spherical accumulator, rotate the point cloud, and perform cluster average calculation on the initial normal vector to obtain the precise normal vector corresponding to the target acquisition point; S2. Using the improved region growing algorithm, the seed point is calculated according to the precise normal vector of the corresponding target acquisition point, and the similarity region is defined by the similarity criterion. The optimal number of structural planes is obtained by setting the gradient threshold of the normal vector of the initial growth condition and drawing the number of structural planes, and a single structural plane of several rock point clouds is obtained; S3. Traverse the individual structural surfaces of each rock mass point cloud, calculate the analytical normal vector and use it as the feature representative. Based on the feature representative, use the improved K-means++ clustering method to replace the distance function with the sum of squared angle errors to obtain the optimal number of clusters, determine the division of the dominant group structural surface, and complete the intelligent identification method of the rock mass structural surface.

2. The method for intelligently identifying rock mass structural surfaces based on three-dimensional point clouds according to claim 1, characterized in that: The S1 comprises the following steps: S101, using the KD tree index to search for the nearest points of the target collection point, and using the searched nearest points as the initial point cloud set; S102, balancing different density weights of each point cloud in the initial point cloud set by weighting the point cloud density; S103, randomly selecting three nonlinear points from the initial point cloud set to form a triplet, randomly selecting two three-dimensional vectors from the constructed triplet, and obtaining an initial normal vector by vector cross multiplication calculation; S104, setting a gridded spherical accumulator to simulate the voting process, making the area of ​​each grid in the spherical accumulator the same, and setting an independent correlation distribution between the number of spherical accumulators and the number of triplets; S105, according to the setting result of S104, using the probability empirical formula, determine the relationship between the actual value and the empirical value, and quote the Hoefding inequality into the probability empirical formula to calculate the minimum triple number inequality; S106. According to the central limit theorem, a random variable of a process of determining the larger one between two grids with the most votes in the spherical accumulator through a confidence interval is obtained, and a disjoint inequality is used in combination with a minimum triplet number inequality, and in response to the confidence interval of the grid with the most votes in the spherical accumulator and the confidence interval of the grid with the second most votes not intersecting, triplet construction is stopped, wherein, in the triplet construction process, an initial normal vector in the spherical accumulator forms a normal distribution and generates a peak area, and the peak of the peak area is made to be located in the corresponding grid of the spherical accumulator by rotating the point cloud; S107, performing cluster average calculation on the initial normal vectors of the peak area within a certain threshold range to obtain the precise normal vectors corresponding to the target acquisition points.

3. The method for intelligently identifying rock mass structural surfaces based on three-dimensional point clouds according to claim 2 is characterized in that: The expression for the independent correlation distribution between the number of spherical accumulators and the number of triples is as follows: (X m,t ) m∈ {1,…,M},t∈{1,…,T} Where m represents the mth spherical accumulator, t represents the tth triplet, and X m,t represents the independent correlation distribution between variables m and t, M represents the number of spherical accumulators, and T represents the number of triplets; The probability empirical formula is as follows: Among them, P(·) represents the probability, P m Indicates the actual value, represents the empirical value, α represents the confidence interval, and α∈[0,1], T min represents the minimum number of triples, δ represents the distance between the actual value and the empirical value, and δ∈[0,1], ln(·) represents the logarithmic function.

4. The method for intelligently identifying rock mass structural surfaces based on three-dimensional point clouds according to claim 2, characterized in that: The expression of the random variable of the larger determination process is as follows: Among them, Q represents the random variable of the process of determining the larger one, and its distribution converges to the standard normal random variable N(0,1), p m Indicates the actual value, represents the experience value, T represents the number of triplets; The expression of the disjoint inequality is as follows: Among them, m1 represents the grid with the most votes, and m2 represents the grid with the second most votes. Indicates the experience value of the grid with the most votes, Indicates the experience value of the grid with the second most votes.

5. The method for intelligently identifying rock mass structural surfaces based on three-dimensional point clouds according to claim 1, characterized in that: The S2 comprises the following steps: S201, setting a normal vector angle threshold and a curvature magnitude threshold between two points, and adding a distance threshold and an average curvature threshold of the neighborhood of the current target acquisition point as auxiliary items to define a similarity criterion; S202, calculating the point cloud curvature using a normal vector difference method according to the precise normal vector of the target acquisition point, and selecting seed points in order by sorting the point cloud curvature from small to large; S203, traversing the seed points in sequence, constructing an index sequence according to the seed points, adding the points in the initial point cloud set that meet the similarity criterion to the index sequence, defining the same index sequence as a similarity region, marking the processed seed points, and packaging and coloring each similarity region until all points in the unstructured point cloud are processed; S204, calculating the number of structural surfaces according to the packaged and colored similarity regions, and drawing a structural surface number curve according to the set normal vector angle threshold between two points, curvature magnitude threshold, distance threshold, and average curvature threshold of the neighborhood of the current target acquisition point, determining the inflection point of the curve and combining visual interpretation to achieve a compromise segmentation effect; S205. Setting the gradient threshold of the normal vector of the initial growth condition, observing the change in the number of structural planes as the curvature increases under different normal vectors under a compromise segmentation effect, obtaining the optimal number of structural planes completed by the improved regional growing algorithm, and obtaining a single structural plane of several rock point clouds.

6. The method for intelligently identifying rock mass structural surfaces based on three-dimensional point clouds according to claim 1, characterized in that: The S3 comprises the following steps: S301, traversing the individual structural planes of each rock mass point cloud, calculating the analytical normal vector of each individual structural plane using the principal component analysis method, and taking the analytical normal vector corresponding to each individual structural plane as a characteristic representative of the corresponding individual structural plane; S302, using the improved K-means++ clustering method, taking the feature representatives of each single structural surface as a data set, using the squared sum of angle errors of the normal vector to replace the distance function, and using it as a point cloud normal vector clustering judgment criterion, wherein the angle judgment criterion of the squared sum of angle errors adopts the square of the cosine of the angle between two normal vectors; S303, according to the point cloud normal vector clustering evaluation criteria, draw an angle error square sum curve, and according to the angle error square sum curve, obtain the optimal number of clusters completed by the improved K-means++ clustering method; S304. Combining the number of optimal structural surfaces completed by the improved region growing algorithm and the optimal number of clusters completed by the improved K-means++ clustering method, the number of surfaces completed by region growing and clustering is the same, the division of the dominant group structural surfaces is determined, and the intelligent identification method of rock mass structural surfaces is completed.

Images