Three-dimensional rock mass structural plane trace extraction method based on voxel morphology

Abstract

The invention relates to the technical field of structural plane trace extraction, and discloses a voxel morphology-based three-dimensional rock mass structural plane trace extraction method, which comprises the following steps of: performing point cloud normal vector estimation on rock mass three-dimensional point cloud by adopting robust random Hough transform, calculating a three-dimensional curved surface feature value of the point cloud through a normal vector neighborhood difference method, and extracting the three-dimensional curved surface feature value of the point cloud according to the three-dimensional curved surface feature value. Separating trace points from plane points by adopting a clustering method; carrying out voxel grid division on the point cloud space according to the point density, and carrying out expansion and corrosion morphological preprocessing on a voxel grid; a branch structure connection method is adopted to vectorize the voxel skeleton topology so as to complete connection and penetration of trace paths, and redundant branches and single rings are removed at the same time. According to the invention, the trace extraction effect and precision are improved.

Description

Technical Field

[0001] This invention relates to the field of structural surface trace extraction technology, specifically to a method for extracting three-dimensional rock mass structural surface traces based on voxel morphology. Background Technology

[0002] In natural rock masses, whether on the surface or inside, the distribution of structural planes is highly random and uncertain. While cutting the rock mass to create a complex topological structure, they also define the mechanically weak interfaces within the rock mass. Therefore, identifying structural planes and obtaining their corresponding geological information elements is a crucial foundation for rock mass stability analysis and prediction, surrounding rock classification, dynamic support, and other engineering projects.

[0003] Traces, as one of the structural features, exist as intersection lines at the junction of the structural surface and the excavation face, exhibiting relatively obvious linear characteristics. They possess the property of directly characterizing the geometric dimensions and development scale of the structural surface. Therefore, the geometrical digital depiction of traces is an important component of rock mass exploration. However, for traditional manual contact measurement, in addition to issues of safety, efficiency, and range limitations, the subjectivity of human experience is also a major influencing factor in the logging results. Therefore, after abandoning the original manual measurement methods, in order to advance the research on automated trace extraction towards a direction that balances efficiency and accuracy, scholars have conducted a great deal of research in this area. Automated trace extraction was first applied in the field of simple line segment detection and pixel feature separation in two-dimensional image processing. Subsequently, various edge detection methods such as Canny, Roberts, and Sobel operators have gradually demonstrated good extraction results.

[0004] However, it should be noted that the aforementioned two-dimensional image trace extraction methods are only suitable for segmenting and extracting pixels with obvious features and simple structures. Their computational power is somewhat insufficient for images with poor quality or complex traces. Therefore, some researchers have combined edge extraction and region segmentation methods to achieve more refined extraction of rock fracture images. They have also combined artificial neural networks with noise and pseudo-boundary pixel filtering of trace boundary lines. Other researchers have used multi-angle processing to sequentially digitize rock fracture images, achieving semi-automatic detection and recognition of fracture edges and extraction of edge coordinates. With the increasing sophistication of rock mass image recognition technology, more research is focused on ensuring the robustness and noise resistance of the extraction results. Image resolution is improved through image interpolation methods, and improved Hough transform and Livewire methods achieve clearer extraction for fracture images of different shapes. Furthermore, some researchers have proposed a CNN-based FraSegNet deep learning image segmentation model for automatic segmentation and extraction of rock fractures and traces, demonstrating impressive noise reduction capabilities.

[0005] It is worth mentioning that the information acquisition range of two-dimensional images is limited, and it has high requirements for dust-free environments and lighting conditions. Meanwhile, the continuously developing three-dimensional model construction technology has long since led modern research towards methods for extracting geological information elements at a higher dimension. Therefore, in recent years, many scholars have begun to work on the extraction of traces from three-dimensional point cloud models of rock masses. Initially, this was done through manual interaction, and without automation, its computational cost was even slightly higher than that of two-dimensional extraction. Later, some scholars, based on the automatic identification of individual structural surfaces, used the two farthest points as traces, but this definition differs somewhat from reality. To more closely approximate the geometric shape of real traces, two main methods are currently used: the first is based on the definition of the essence, treating the intersection lines between adjacent structural surfaces as traces after obtaining the single surface of the structure, i.e., extracting traces through the common points of adjacent surfaces; the second is based on three-dimensional spatial features, separating feature points in the point cloud model with significant changes in principal curvature values ​​as trace points. However, considering that the extraction effect of the first method depends on the accuracy of the structural surface division, and the relationship between common points on various surfaces is also difficult to resolve, it presents certain challenges. Therefore, the second type of application is more widespread. Currently, some methods rely on the minimum cost path method combined with RGB values ​​and gradient cost functions to connect point clouds to form traces. A more accurate method is to construct a Digital Surface Model (DSM) from the point cloud model, calculate feature points in the triangular mesh based on tensor voting theory, and then achieve refined trace extraction through a series of operations such as feature point grouping, region growing and connecting, and redundant line segment removal. In addition, with further optimization of tensor computation, recognition efficiency and trace smoothness have also been greatly improved. Summary of the Invention

[0006] To address the aforementioned shortcomings in existing technologies, this invention provides a method for extracting traces of three-dimensional rock mass structure surfaces based on voxel morphology, thereby improving the effectiveness and accuracy of trace extraction and providing better technical support for intelligent rock mass logging.

[0007] To achieve the above-mentioned objectives, the technical solution adopted by this invention is as follows: A method for extracting three-dimensional rock mass structural surface traces based on voxel morphology includes the following steps: Robust stochastic Hough transform is used to estimate the normal vector of the three-dimensional point cloud of the rock mass. The three-dimensional surface feature value of the point cloud is calculated by the normal vector neighborhood difference method, and the clustering method is used to separate the trace points and the planar points. The point cloud space is divided into voxel meshes according to the point density, and the voxel meshes are subjected to morphological preprocessing of dilation and erosion. A branch structure connection method is used to topologically vectorize the voxel skeleton to complete the connection and continuity of the trace path, while removing redundant branches and single loops.

[0008] Furthermore, the point cloud normal vector estimation for the three-dimensional point cloud of the rock mass using robust stochastic Hough transform includes: A point is randomly selected as the center point, and a set number of neighboring points are searched to form a set of points to be calculated. Two other points in the set are randomly selected to form multiple nonlinear triples. At the same time, a spherical accumulator is constructed to statistically analyze the distribution of the normal vector values ​​of all current triples. The normal vectors are converted into dip angles and then accumulated by voting according to the corresponding grid range. Clustering is performed in the grid peak area by a set angle threshold to obtain the final normal vector.

[0009] Furthermore, the calculation of the 3D surface eigenvalues ​​of the point cloud using the normal vector neighborhood difference method includes: The angle between the target center point and its neighboring points is calculated based on the point cloud normal vectors, specifically as follows:

[0010] in, The angle between the target center point and its nearest neighboring points. The target center point x i Normal vector Neighboring points x j Normal vector; The 3D surface feature values ​​of the point cloud are calculated based on the angle between the target center point and its neighboring points, specifically as follows:

[0011] in, The three-dimensional surface feature value at the target center point. N Let be the number of neighboring points.

[0012] Furthermore, the clustering method is used to separate trace points from planar points, including: Using the K-means++ clustering method, the objective summation function is set as follows:

[0013] in, Let the objective sum function be... For sample points x and the center point of its cluster μ i The difference between the neighborhood values ​​of the normal vector. n This represents the current number of clusters. C i For the first i A cluster; The minimum value of the target summation function is calculated by iterating the center cluster position under different cluster numbers. The corresponding line graph is plotted and the elbow rule is used to select the larger cluster number as the ideal cluster number.

[0014] Furthermore, dividing the point cloud space into voxel meshes according to point density includes: We select twice the average point cloud spacing as the voxel mesh edge length, initialize the voxel index containing the point cloud, and use radius filtering to remove abnormal voxels.

[0015] Furthermore, the morphological preprocessing of the voxel mesh, including dilation and erosion, includes: Two erosion operations are performed. The second erosion operation involves sliding the structuring element sequentially across the foreground. If the structuring element is completely contained within the target region, all voxels in the structuring element except the central voxel are added to the initial set to be deleted. This set is then traversed sequentially, and if the element itself is also completely contained within the target region, it is removed from the initial set. The final set to be deleted is the result of this filtering process.

[0016] Furthermore, the voxel skeleton topology vectorization using a branched connection method includes: The voxel skeleton is classified into endpoint voxels, simple voxels, and core voxels. Among them, a foreground voxel that is adjacent to a single voxel only by a point, line, or surface is an endpoint voxel. If the Euler characteristic of the current foreground voxel and the background voxel does not change after deleting a foreground voxel, it is called a simple voxel. Other voxels are called core voxels. The voxels are divided into multiple directional subsets according to the 26-adjacency method. Each time, only voxels in a specific direction are deleted. Then, each foreground voxel is traversed. If the foreground voxel is a simple voxel and not an endpoint voxel, then the foreground voxel is deleted and marked as a background voxel. This continues until the voxel skeleton converges.

[0017] Furthermore, the Euler characteristic is specifically:

[0018] in, Let be an Euler characteristic. The number of connected components. The number of 1-dimensional holes, The number of 2D closed cavities.

[0019] Furthermore, the voxel skeleton topology vectorization using a branched connection method includes: Count the 26 adjacencies of each voxel; if the number of adjacencies is 0, it is classified as an isolated voxel; if the number of adjacencies is 1, it is classified as an endpoint voxel; if the number of adjacencies is 2, it is classified as a connected voxel; if the number of adjacencies is ≥3, it is classified as a branch voxel. Construct the parent-child node relationships of voxels; starting from each endpoint voxel, calculate the topological distance between each endpoint, and take the pair of endpoint voxels with the largest sum of path branch voxels and hop count as the skeleton body; starting from other endpoint voxels as child nodes and the nearest branch voxel as parent node, take any end of the skeleton body as the root child node and the other end as the terminal parent node; after determining the global direction, abstract the connection of the remaining branch nodes into graph theory relationships, and then use the breadth-first search method to perform a global search starting from the root child node, and record the parent-child relationships of each voxel; For each voxel, only its parent node and the connection points of the first and second branches are recorded, and topological vectorization is performed in a binary form.

[0020] Furthermore, removing redundant branches and single loops includes: Starting from the terminal voxel of any branch, for the voxel sequence, calculate the normalized local direction vector starting from the terminal voxel; Calculate the angle between two consecutive direction vectors and the corresponding angle function value based on the normalized local direction vector; The geometric curvature of the branch voxel is determined by accumulating all the included angle function values ​​in the branch voxel, and branch voxels that exceed the geometric curvature threshold are identified as redundant branches. Record the adjacency relationships of all branch points; traverse each branch point pair, and if there are two paths that only contain the intermediate point connecting the point pair, it is determined to be a single cycle, and all intermediate points of the longer path are deleted while updating the node relationships.

[0021] The present invention has the following beneficial effects: This invention addresses the problems of high risk, low efficiency, and poor accuracy associated with traditional methods, as well as the insufficient dimensionality of geological information elements that can be acquired. It proposes a voxel morphology-based method for extracting 3D rock mass structural surface traces. First, a robust stochastic Hough transform is used to calculate the point cloud normal vectors. Then, the surface eigenvalues ​​are calculated using the neighborhood difference method based on the normal vectors, and an improved K-means++ clustering method is used to separate and extract structural surface trace points from planar points. Next, a cubic mesh is set based on the point density distribution to voxelize the spatial location of the point cloud. Further, the existing voxel mesh is preprocessed by expansion and contraction based on voxel morphology. Finally, a branch structure connection method is used to topologically vectorize the voxel skeleton and remove redundant branches and single-loop parts of the skeleton, thereby improving the effectiveness and accuracy of trace extraction. Attached Figure Description

[0022] Figure 1 This is a schematic diagram of a method for extracting surface traces of three-dimensional rock mass structure based on voxel morphology. Figure 2 Schematic diagram of a spherical accumulator; Figure 3A diagram illustrating the elbow rule; Figure 4 A schematic diagram of the voxelized point cloud space; Figure 5 This is a schematic diagram of index voxel adjacency; Figure 6 A schematic diagram of the structural unit and its index; Figure 7 A diagram illustrating voxel classification; Figure 8 This is a schematic diagram of graph theory relations; Figure 9 This is a schematic diagram of the slope site photos; Figure 10 This is a schematic diagram of the original point cloud data; Figure 11 This is a diagram illustrating the optimal number of clusters. Figure 12 This is a diagram illustrating the clustering and merging results; Figure 13 For detailed views of the index area; Figure 14 This is a schematic diagram of the structural unit; Figure 15 This is a schematic diagram of morphological restoration. Figure 16 This is a schematic diagram of voxel skeleton extraction and topological vector connection. Figure 17 This is a schematic diagram illustrating the removal of redundant branches and single loops. Figure 18 A visual representation of the overall extraction results; Figure 19 This is a diagram showing the comparison of trace extraction results. Detailed Implementation

[0023] The specific embodiments of the present invention are described below to enable those skilled in the art to understand the present invention. However, it should be understood that the present invention is not limited to the scope of the specific embodiments. For those skilled in the art, various changes are obvious as long as they are within the spirit and scope of the present invention as defined and determined by the appended claims. All inventions utilizing the concept of the present invention are protected.

[0024] This invention proposes a method for extracting traces of 3D rock mass structural surfaces based on voxel morphology. The core idea of ​​this method is as follows: First, a robust stochastic Hough transform is used to estimate the normal vectors of the point cloud. Then, the 3D surface feature values ​​of the point cloud are calculated using the normal vector neighborhood difference method, and an improved K-means++ clustering method is combined to perform preliminary separation and extraction of planar points and trace points. Next, the point cloud space is divided into voxel meshes of appropriate size according to the point density, and the existing voxel meshes undergo morphological preprocessing involving dilation and erosion. Finally, to find the simplest topological voxel paths, a branch structure connection method is used to topologically vectorize the voxel skeleton to complete the connection and continuity of the trace paths, while removing redundant parts.

[0025] like Figure 1 As shown, an embodiment of the present invention provides a method for extracting three-dimensional rock mass structure surface traces based on voxel morphology, comprising the following steps S1 to S3: S1. Robust stochastic Hough transform is used to estimate the normal vector of the three-dimensional point cloud of the rock mass. The three-dimensional surface feature value of the point cloud is calculated by the normal vector neighborhood difference method, and the clustering method is used to separate the trace points and the plane points. In an optional embodiment of the present invention, for rock mass structures with complex geometric surfaces, especially where sharp edges abruptly change direction, the quality requirements for the normal vector are higher. Therefore, this embodiment employs robust stochastic Hough transform to calculate the point cloud normal vector. First, a point is randomly selected as the center point, and a given number of neighboring points are searched to form a set of points to be calculated. Then, two other points in the set are arbitrarily selected to form several nonlinear triples. The result of the vector cross product is used as the normal vector of this triple. Simultaneously, a "spherical accumulator" is constructed, such as... Figure 2 As shown, this is used to statistically analyze the range distribution of normal vector values ​​for all current triplet pairs, converting the normal vectors into dip angles. Representative tendencies, (Representative inclination) and then accumulate votes based on the corresponding grid range.

[0026] This process is set One accumulator grid, If there are three triples, then the two triples satisfy the following condition: The independent correlation distribution means that as the number of triples increases, the normal vector on the accumulator grid will gradually exhibit a normal distribution. However, the number of triples required to construct a triple is approximately the cube of the number of points in the set. To save computational resources, we introduce an empirical probability formula and Hoeffding's inequality to obtain the minimum number of triples inequality, expressed as:

[0027]

[0028]

[0029] in, Represents probability. Represents the empirical distribution value. Represents the actual distribution value. This represents the difference between the empirical distribution and the actual distribution. Indicates the confidence level. This represents the minimum number of triplets. Finally, clustering is performed in the grid peak region using a set angle threshold to obtain the final normal vector.

[0030] After calculating the point cloud normal vectors, to quantify the three-dimensional spatial surface features of the point cloud and to quickly separate planar points from trace points, this embodiment uses the neighborhood difference method for normal vectors of non-smooth features instead of the commonly used principal curvature values. First, the angle between the target center point and a neighboring point is calculated. , is represented as:

[0031] in, Represents the normal vector of the target center point. This represents the normal vector of the nearest point.

[0032] Then calculate the difference value of its normal vector neighborhood. , is represented as:

[0033] in, Indicates the number of neighboring points. The angle between the normal vectors in Equation 4. This yields the angle for each point cloud. After the values ​​are calculated, the improved K-means++ clustering method is used to separate planar points from trace points. It is known that K-means++ optimizes the selection of initial centroids compared to the original K-means, i.e., selecting points far from the already selected centroids based on probability distribution. However, its iterative stage still determines the final clusters by minimizing the WCSS (Sum of Squared Errors within Clusters), and the number of clusters needs to be pre-set, making it significantly susceptible to subjective influence. The original WCSS calculation method was:

[0034] in, Indicates the current number of clusters. Represents each data point, Indicates the first Clusters, This indicates the center point of the current cluster.

[0035] Therefore, this embodiment sets a target summation function based on this. Find the function value for each cluster size:

[0036] in, Represents sample points and the center point of its cluster The difference in the neighborhood of the normal vector. Then, the position of the central cluster is calculated iteratively. Find the minimum value for different numbers of clusters, plot the corresponding line graphs, and combine this with the elbow rule: at the ideal number of clusters... Below and correspond The value will show a sharp drop, with a clear inflection point on the line chart. The magnitude of the drop is defined as follows:

[0037] in, Indicates the amplitude threshold. For example... Figure 3 As shown, the number of clusters corresponding to the previous and next two clusters If the difference is greater than a threshold, the latter is determined to be an inflection point. Point The conditions are met, but The condition is not met. Finally, the point with the larger cluster number is selected. point As the ideal number of clusters.

[0038] S2. Divide the point cloud space into voxel meshes according to the point density, and perform morphological preprocessing of the voxel meshes by dilation and erosion. In an optional embodiment of the present invention, point clouds, as a type of three-dimensional data with large data volume, no topological structure, and anisotropic density, require a certain degree of organization and encoding to give them a regular spatial structure for subsequent applications. Voxels, as a type of cubic mesh in three-dimensional space, are similar to pixels in a two-dimensional image, but with added... Information along the axis. Compared to other data structures such as triangular meshes, voxels can not only clearly define the spatial location and topological structure of discrete point clouds, but can also be regarded as three-dimensional units with fixed volume, coordinates and features, making them a convenient method for compressing and storing point cloud information elements.

[0039] Unlike vector graphics, voxels are discrete, and their resolution depends on the size of the grid. For example... Figure 4As shown, when the voxel mesh resolution is too high, although it can better represent the spatial morphology of the overall trace points, holes are prone to appear between adjacent points, and it will cause a large computational burden, which is manifested as over-voxelization. On the other hand, when the resolution is too low, although it can achieve fast calculation, it will lose a lot of geometric details and cannot perform effective morphological operations, which is manifested as under-voxelization. Therefore, in order to achieve the effect of preserving structural surface information elements and avoiding empty voxels from destroying connectivity, this embodiment selects twice the average point cloud spacing as the voxel mesh edge length, initializes the voxel index containing the point cloud, and uses radius filtering to remove abnormal voxels.

[0040] Morphology is an image processing technique based on set theory and topology, used to analyze geometric structures. This section uses morphology to optimize the topological space of trace point voxel meshes.

[0041] First, in the voxel mesh, voxels containing trace points are marked as "foreground" and assigned an index of 1; conversely, voxels without trace points are marked as "background" and assigned an index of 0. The overall result after these assignments is called a "binary 3D image".

[0042] For each voxel, there are 6 adjacencies based on the voxel's 6 faces, 12 adjacencies based on the voxel's 12 edges, and 8 adjacencies based on the voxel's 8 vertices. There are also various other adjacency methods; for ease of classification and explanation, voxel mesh index assignment is defined here. For example... Figure 5 As shown, the indices of voxels adjacent in three ways—point, edge, and face—relative to the central voxel are respectively... , , .

[0043] like Figure 6 As shown, different structural units and corresponding indices can be set according to the existing voxel mesh size. Taking 26 adjacencies as an example (i.e., 3×3×3 outer circumference), for each voxel index... The adjacency index is .in ,and (Hereinafter referred to as 26 adjacency).

[0044] In morphology, dilation and erosion are based on the core idea of ​​locally modifying existing unit elements using specific structuring elements to extract or enhance geometric features. The purpose of dilation is to expand the target region (foreground) in a binary 3D image outwards, filling edges and gaps in connecting voxels. The principle is as follows: a structuring element is slid across the foreground; if the structuring element intersects with at least one voxel in the foreground, the central voxel of that structuring element is designated as the foreground. The mathematical expression for dilation can be defined as:

[0045] The purpose of erosion is to shrink the foreground portion of a binary 3D image inward, removing edges and isolated points. Its principle is as follows: as the structuring element slides sequentially across the foreground, when the structuring element is completely contained within the target region, the central voxel is retained while the other voxels are set as background. The mathematical expression for erosion can be defined as:

[0046] in, Represents a binary 3D image. Represents a structural element. Indicates the position of the voxel. Indicates will Move to .

[0047] In summary, the setting of the structuring element is the decisive factor in the expansion and erosion operation. If the structuring element is too large, it will cause the foreground area to merge and stick together excessively during expansion, resulting in the loss of geometric details. It will also cause the deletion of effective structures or even the destruction of connectivity during erosion. On the other hand, if the structuring element is too small, it will cause the broken areas to be unable to be effectively connected during expansion, and may even produce more burrs. At the same time, it will amplify the noise effect during erosion.

[0048] However, based on the aforementioned properties and the special characteristics of trace geometry, the structural units of expansion and erosion often need to be different, and the execution order should be expansion first followed by erosion. This also results in the latter's spatial shape being less than or equal to the former. Therefore, this embodiment adopts a two-stage erosion method to minimize noise generation while preserving connectivity to the greatest extent. The improved erosion method is as follows: As the structural elements slide sequentially on the foreground, when a structural element is completely contained within the target region (i.e., completely surrounded), all voxels in the structural element except the central voxel are added to the initial set to be deleted. This set is then traversed sequentially, and if it is also completely surrounded, it is removed from the initial set. The final set after filtering is the final set to be deleted. This method can be defined as:

[0049] in, This represents the voxel currently being traversed. Indicates the current structural element. Represents the original set of indexes.

[0050] S3. The voxel skeleton topology is vectorized using a branch structure connection method to complete the connection and connection of trace paths, while removing redundant branches and single loops.

[0051] In an optional embodiment of the present invention, after the above morphological processing, the broken voxels at the trace points and the void portions can be connected while preserving the overall geometric properties. However, to extract the voxel skeleton, it is necessary to further incorporate topological principles to achieve the simplest topological connection.

[0052] The steps for voxel skeleton extraction are as follows: 1. First, classify the key voxels into three types: endpoint voxels, simple voxels, and core voxels. An endpoint voxel is a foreground voxel that is adjacent to a single voxel only by a point, line, or surface. A simple voxel is one in which the Euler characteristic of the current foreground and background voxels does not change after deleting a foreground voxel. All other voxels are called core voxels.

[0053] Among them, Euler characteristic for:

[0054] in, Indicates the number of connected components. Indicates the number of 1-dimensional holes (such as tori). This indicates the number of 2D closed cavities (such as spherical cavities). Number of connected components: refers to the total number of mutually unconnected but internally connected sub-parts in a binary 3D image (using 26-adjacent connectivity).

[0055] like Figure 7 As shown, taking the current voxel group as an example, the whole can be regarded as two tori connected along one edge, which is still a connected space. Therefore The value is 1. Simultaneously, each torus contributes two one-dimensional toruses, resulting in... The value is 4. The internal cavity is connected... It is also 1. Therefore, the Euler characteristic of the voxel group is 1. As can be seen from the figure, after classifying into three key voxels, there are exactly three simple voxels that meet the conditions, excluding the endpoint voxels.

[0056] 2. Divide the voxels into multiple directional subsets based on 26-adjacency. Each time, delete only voxels in a specific direction. Then iterate through each foreground voxel. ,like It is a simple voxel and If it is not an endpoint voxel, then... Delete and mark as background voxels.

[0057] 3. Iterate through step 2 until the voxel skeleton converges, which means that after deleting any voxel, the Euler characteristic of the foreground and background will change (i.e., simple voxels no longer exist).

[0058] After obtaining the minimal topological skeleton of the trace point voxels, the next step is to vectorize it using a branching structure method to fit the real trace, which is used to record the topological connectivity of each voxel.

[0059] The basic steps are as follows: 1. Count the 26 adjacencies of each voxel. If the number of adjacencies is 0, it is classified as an isolated voxel; if the number of adjacencies is 1, it is classified as an endpoint voxel; if the number of adjacencies is 2, it is classified as a connected voxel; if the number of adjacencies is ≥3, it is classified as a branch voxel.

[0060] 2. Construct parent-child node relationships for voxels. Starting from each endpoint voxel, calculate the topological distance (hop count) between each endpoint. Take the pair of endpoint voxels with the highest sum of path branch voxels and hop counts as the main skeleton. Starting from the nearest branch voxel as the parent node, take any one end of the main skeleton as the root child node and the other end as the terminal parent node. After determining the global direction, as follows... Figure 8 As shown, the connection of the remaining branch nodes can be abstracted into graph theory relations, and then BFS (Breadth-First Search) can be used to perform a global search starting from the root child node, and the parent-child relationship of each voxel can be recorded.

[0061] 3. To facilitate recording parent-child relationships, each voxel only records its parent node and the connection points of the first and second branches, and topologically vectorizes it in a binary form.

[0062] Redundant branches and loops can occur during skeleton extraction due to various reasons, such as uneven point density, noise, and over-dilation. Therefore, to ensure that the extraction results closely resemble reality, it is necessary to perform post-processing to remove redundant branches and loops.

[0063] In the previous step of establishing the parent-child node relationship, the situation from other endpoint voxels to their respective nearest branch voxels was already recorded. Redundant branches usually appear in these situations and manifest as geometrically irregular "spurs." Here, the cumulative angle change is calculated to quantify the degree of geometric curvature of the branches. First, starting from the endpoint voxel of any branch... Starting from this point, the branch voxel to which this branch belongs is defined as Then for this voxel sequence ,from Start calculating the normalized local direction vector :

[0064]

[0065] Calculate the angle between two consecutive direction vectors:

[0066] It can be known The range of values ​​is Therefore, through the included angle function measure Size:

[0067] Finally, the geometric curvature of the branch is determined by summing up all the included angle function values ​​in the branch:

[0068] And set This is the geometric bending threshold for the branch; if it exceeds the threshold, it is considered a redundant branch.

[0069] Redundant single cycles typically occur on paths formed by a single intermediate point. Therefore, in addition to obtaining the global relationship through BFS, the adjacency relationships of all branch points are also recorded. For each pair of branch points, if two paths connecting this pair contain only the intermediate point, it is determined to be a single cycle, and all intermediate points of the longer path are deleted while updating the node relationships.

[0070] like Figure 9 , Figure 10 As shown, this embodiment uses an open-source dataset from RockBench as an example. The location is a highway slope in Colorado, USA. The scanning device is an Optech Ilris 3D laser scanner, and the total number of point clouds is 1,515,722.

[0071] First, data preprocessing was performed, including cropping, denoising, and downsampling of the point cloud data, resulting in a final point cloud of 1,085,373 points. Second, the point cloud normal vectors were calculated based on a robust stochastic Hough transform, with a target point neighborhood search size of 500, a minimum triplet count of 700, and a grid clustering threshold of [value missing]. After obtaining the normal vector, the neighborhood difference value of the normal vector of each point is calculated, and clustering is performed using the improved K-means++ clustering method proposed in this embodiment. like Figure 11 As shown. An amplitude threshold of 0.6 is set, and the output inflection point is locked when the number of clusters equals 3. To ensure sufficient trace points are selected, the cluster number is rounded down by one. Finally, the cluster with the smallest average normal vector neighborhood difference among the four types of points is selected as the plane point, and the other three types are merged as trace points. The classification and color assignment results are shown below. Figure 12 As shown.

[0072] After selecting three times the average point spacing as the voxel side length, to clearly visualize the trace extraction process, this embodiment first divides the entire slope into several blocks of 10×10×10 along the three axes. The block with index 435 is used as an example. For example... Figure 13As shown, the selected region has a relatively well-developed structural surface and a certain number of broken surfaces, which can be used as a verification object for accuracy and robustness.

[0073] In the morphological pretreatment stage, this embodiment employs one expansion followed by two etching processes. The first etching is normal etching, and the second is modified etching. The structural units used in the three processes are as follows: Figure 14 As shown.

[0074] like Figure 15 As shown, after one expansion and two contractions, the morphological repair not only improved the voxel connectivity in areas with low point density, but also maintained the smoothness of the structural edges and the details of the geometric shape through uniform contraction.

[0075] After morphological processing, the simplest topological skeleton was further extracted. After classifying three types of keypoints, simple voxels were iteratively removed, and the voxel skeletons were topologically vectorized and connected. The result is shown below. Figure 16 As shown.

[0076] Finally, the cumulative angle change was calculated, redundant branches and single loops in the region were removed, and a selection of comparison results were shown below. Figure 17 As shown.

[0077] To verify the accuracy, the traces proposed by this method were compared with those obtained manually, and the accuracy and error were evaluated using several parameters. A comparison of global effect visualization and specified region extraction was also provided. Figure 18 , Figure 19 As shown in Table 1 (where white lines represent the results extracted by this method, and red lines represent the results of manual delineation).

[0078] Table 1 Comparison of Trajectory Extraction Data

[0079] As shown in the charts above, the maximum and minimum error rates differ significantly in example regions 1 and 2, while region 3 is more stable. This is because the structural surfaces in region 3 are of similar size, and the overall trace length distribution is more regular, resulting in smaller identification deviations. However, the varying degrees of development of structural surfaces in regions 1 and 2 lead to larger deviations. Furthermore, it can be observed that this method ignores a few small structural surface traces or merges them into a single line, due to connectivity issues during voxel expansion, leading to similar misjudgments. Analysis reveals that the main reason for the errors in this method is that the clustering and separation of trace points slightly absorbs some planar points, resulting in fewer remaining planar points in excessively small structural surfaces. Consequently, after voxel expansion and erosion, traces that should form "topological one-dimensional loops" from small structural surfaces are mistakenly extracted as points. Clearly, reducing the voxel resolution can ensure the identification of such structural surface traces, but this would also significantly reduce the connectivity of other voxels. Therefore, balancing overall identification performance and accuracy should be the dominant factor in voxel size selection, and the impact of small structural surfaces is within an acceptable range in terms of accuracy and error parameters.

[0080] In summary, the three-dimensional traces obtained by the method in this embodiment through a series of feature extraction steps may be affected to some extent by relatively special small structural surfaces. However, through visual interpretation and combined with various evaluation indicators, it can be seen that the three-dimensional traces have a good fit with the real traces in complex three-dimensional point cloud scenes. Compared with traditional methods, it can provide richer digital information support for rock masses.

[0081] Current methods and techniques for depicting rock mass structural traces are relatively underdeveloped, and some methods are limited by the dimensionality of geological information elements acquired. To address these issues, this study uses voxel morphology to extract target trace points from a point cloud divided by voxel mesh. Through analysis and data comparison, the following conclusions are drawn: 1. This embodiment obtains the surface features of the point cloud based on robust stochastic Hough transform and normal vector neighborhood difference method, and then uses an improved K-means++ clustering method to initially separate trace points. Simultaneously, it combines morphology and topology to perform voxel mesh dilation and erosion, minimal topological skeleton extraction, topological vectorization of trace paths, and removal of redundant branches and single loops, thus realizing trace extraction of three-dimensional rock mass structural surfaces.

[0082] 2. This embodiment proposes a method for dividing the space of trace point clouds using voxels, which solves the problem that the spatial discreteness and density heterogeneity of point clouds make it difficult to quantify the connection information between a large number of points. Essentially, this method compresses the information of the target point cloud, preserving not only the geometric details of the original data but also achieving dimensionality reduction and ordered organization of the point cloud data through a regularized grid.

[0083] 3. This embodiment selected a highway slope cutting project with well-developed structural surfaces to verify the feasibility of the method. By comparing the extraction and manual delineation results of three example areas, and combining the quantitative identification results with multiple evaluation indicators, the method has good accuracy and stability.

[0084] This invention is described with reference to flowchart illustrations and / or block diagrams of methods, apparatus (systems), and computer program products according to embodiments of the invention. 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, generate instructions for implementing the flowchart illustrations and / or block diagrams. Figure 1 One or more processes and / or boxes Figure 1 A device that provides the functions specified in one or more boxes.

[0085] 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 instruction means, which are implemented in a process Figure 1 One or more processes and / or boxes Figure 1 The function specified in one or more boxes.

[0086] 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.

[0087] Specific embodiments have been used to illustrate the principles and implementation methods of this invention. The descriptions of the embodiments above are only for the purpose of helping to understand the method and core ideas of this invention. At the same time, for those skilled in the art, there will be changes in the specific implementation methods and application scope based on the ideas of this invention. Therefore, the content of this specification should not be construed as a limitation of this invention.

[0088] Those skilled in the art will recognize that the embodiments described herein are intended to help the reader understand the principles of the invention, and should be understood that the scope of protection of the invention is not limited to such specific statements and embodiments. Those skilled in the art can make various other specific modifications and combinations based on the technical teachings disclosed in this invention without departing from the spirit of the invention, and these modifications and combinations are still within the scope of protection of this invention.

Claims

1. A method for extracting three-dimensional rock mass structural surface traces based on voxel morphology, characterized in that, Includes the following steps: Robust stochastic Hough transform is used to estimate the normal vector of the three-dimensional point cloud of the rock mass. The three-dimensional surface feature value of the point cloud is calculated by the normal vector neighborhood difference method, and the clustering method is used to separate the trace points and the planar points. The point cloud space is divided into voxel meshes according to the point density, and the voxel meshes are subjected to morphological preprocessing of dilation and erosion. A branch structure connection method is used to topologically vectorize the voxel skeleton to complete the connection and continuity of the trace path, while removing redundant branches and single loops.

2. The method for extracting three-dimensional rock mass structure surface traces based on voxel morphology according to claim 1, characterized in that, Point cloud normal vector estimation for 3D point clouds of rock masses using robust stochastic Hough transform includes: A point is randomly selected as the center point, and a set number of neighboring points are searched to form a set of points to be calculated. Two other points in the set are randomly selected to form multiple nonlinear triples. At the same time, a spherical accumulator is constructed to statistically analyze the distribution of the normal vector values ​​of all current triples. The normal vectors are converted into dip angles and then accumulated by voting according to the corresponding grid range. Clustering is performed in the grid peak region by setting an angle threshold to obtain the final normal vector.

3. The method for extracting three-dimensional rock mass structure surface traces based on voxel morphology according to claim 1, characterized in that, The calculation of 3D surface eigenvalues ​​of point clouds using the normal vector neighborhood difference method includes: The angle between the target center point and its neighboring points is calculated based on the point cloud normal vectors, specifically as follows: in, The angle between the target center point and its nearest neighboring points. The target center point x i Normal vector Neighboring points x j Normal vector; The 3D surface feature values ​​of the point cloud are calculated based on the angle between the target center point and its neighboring points, specifically as follows: in, The three-dimensional surface feature value at the target center point, N Let be the number of neighboring points.

4. The method for extracting three-dimensional rock mass structure surface traces based on voxel morphology according to claim 1, characterized in that, Clustering methods are used to separate trace points from planar points, including: Using the K-means++ clustering method, the objective summation function is set as follows: in, Let the objective sum function be... For sample points x and the center point of its cluster μ i The difference between the neighborhood values ​​of the normal vector. n This represents the current number of clusters. C i For the first i A cluster; The minimum value of the target summation function is calculated by iterating the center cluster position under different cluster numbers. The corresponding line graph is plotted and the elbow rule is used to select the larger cluster number as the ideal cluster number.

5. The method for extracting three-dimensional rock mass structure surface traces based on voxel morphology according to claim 1, characterized in that, Dividing the point cloud space into voxel meshes according to point density includes: We select twice the average point cloud spacing as the voxel mesh edge length, initialize the voxel index containing the point cloud, and use radius filtering to remove abnormal voxels.

6. The method for extracting three-dimensional rock mass structure surface traces based on voxel morphology according to claim 1, characterized in that, Morphological preprocessing of voxel meshes, including dilation and erosion, includes: Two erosion operations are performed. The second erosion operation involves sliding the structuring element sequentially across the foreground. If the structuring element is completely contained within the target region, all voxels in the structuring element except the central voxel are added to the initial set to be deleted. This set is then traversed sequentially, and if the element itself is also completely contained within the target region, it is removed from the initial set. The final set to be deleted is the result of this filtering process.

7. The method for extracting three-dimensional rock mass structure surface traces based on voxel morphology according to claim 1, characterized in that, The voxel skeleton topology vectorization using the branch structure connection method includes: The voxel skeleton is classified into endpoint voxels, simple voxels, and core voxels. Among them, a foreground voxel that is adjacent to a single voxel only by a point, line, or surface is an endpoint voxel. If the Euler characteristic of the current foreground voxel and the background voxel does not change after deleting a foreground voxel, it is called a simple voxel. Other voxels are called core voxels. The voxels are divided into multiple directional subsets according to the 26-adjacency method. Each time, only voxels in a specific direction are deleted. Then, each foreground voxel is traversed. If the foreground voxel is a simple voxel and not an endpoint voxel, then the foreground voxel is deleted and marked as a background voxel. This continues until the voxel skeleton converges.

8. The method for extracting three-dimensional rock mass structure surface traces based on voxel morphology according to claim 7, characterized in that, Euler's characteristic is specifically: in, Let be an Euler characteristic. The number of connected components. The number of 1-dimensional holes, The number of 2D closed cavities.

9. The method for extracting three-dimensional rock mass structure surface traces based on voxel morphology according to claim 7, characterized in that, The voxel skeleton topology vectorization using the branch structure connection method includes: Count the 26 adjacencies of each voxel; if the number of adjacencies is 0, it is classified as an isolated voxel; if the number of adjacencies is 1, it is classified as an endpoint voxel; if the number of adjacencies is 2, it is classified as a connected voxel; if the number of adjacencies is ≥3, it is classified as a branch voxel. Construct the parent-child node relationships of voxels; starting from each endpoint voxel, calculate the topological distance between each endpoint, and take the pair of endpoint voxels with the largest sum of path branch voxels and hop count as the skeleton body; starting from other endpoint voxels as child nodes and the nearest branch voxel as parent node, take any end of the skeleton body as the root child node and the other end as the terminal parent node; after determining the global direction, abstract the connection of the remaining branch nodes into graph theory relationships, and then use the breadth-first search method to perform a global search starting from the root child node, and record the parent-child relationships of each voxel; For each voxel, only its parent node and the connection points of the first and second branches are recorded, and topological vectorization is performed in a binary form.

10. The method for extracting three-dimensional rock mass structure surface traces based on voxel morphology according to claim 1, characterized in that, Removing redundant branches and single loops includes: Starting from the terminal voxel of any branch, for the voxel sequence, calculate the normalized local direction vector starting from the terminal voxel; Calculate the angle between two consecutive direction vectors and the corresponding angle function value based on the normalized local direction vector; The geometric curvature of the branch voxel is determined by accumulating all the included angle function values ​​in the branch voxel, and branch voxels that exceed the geometric curvature threshold are identified as redundant branches. Record the adjacency relationships of all branch points; traverse each branch point pair, and if there are two paths that only contain the intermediate point connecting the point pair, it is determined to be a single cycle, and all intermediate points of the longer path are deleted while updating the node relationships.

Images