Rock mass structural plane identification and correlation analysis method and system based on spatial constraint

Abstract

The invention belongs to the technical field of rock mass structural surface recognition, and particularly discloses a rock mass structural surface recognition and correlation analysis method and system based on spatial constraint, and the method comprises the steps: obtaining the point cloud data of a continuous tunnel face of a to-be-measured rock mass, calculating the point cloud normal of each tunnel face, and displaying the overall spatial distribution of the normal through a line drawing joint pole map; density distribution of the point cloud data is obtained based on a kernel density estimation algorithm, and the main direction of the point cloud data is determined in combination with geological priori knowledge; introducing a Fisher probability distribution model, and optimizing lumped parameters of the Fisher probability distribution model to obtain a final structural surface grouping result; performing intra-group clustering on the structural surfaces by using a DBSCAN density clustering algorithm, and extracting each independent structural surface in each group of structural surfaces; and carrying out space correlation analysis on the independent structural surface on the continuous tunnel face based on the three-dimensional space constraint. According to the method, the spatial ductility and occurrence information of the structural surface can be completely restored.

Description

Technical Field

[0001] The invention relates to the technical field of rock mass structural surface identification, and in particular to a rock mass structural surface identification and correlation analysis method and system based on spatial constraints. Background Art

[0002] The statements in this section merely provide background information related to the present invention and do not necessarily constitute prior art.

[0003] Rock mass structural surfaces refer to various types of cracks, layers, joints and other surfaces formed in the rock mass due to the action of internal and external forces. Their distribution, morphology and mutual relationship have an important impact on the design, construction and safety assessment of underground projects. The identification and analysis of structural surfaces is a basic work in geotechnical engineering, which directly affects the accurate assessment of key parameters such as rock mass stability, permeability and deformation characteristics. With the continuous development of modern geotechnical engineering, especially the increasing demand for deep underground projects, tunnels and karst projects, it is particularly important to identify and analyze rock mass structural surfaces efficiently and accurately.

[0004] Traditional rock structural surface identification methods mainly rely on field observation, manual drawing and two-dimensional data analysis, which usually have problems such as low efficiency, poor accuracy and poor adaptability to complex geological bodies. In order to improve the accuracy and operability of rock structural surface identification, in recent years, with the development of three-dimensional imaging technology and computer graphics, rock structural surface identification methods based on three-dimensional spatial data have received widespread attention. The application of three-dimensional laser scanning technology, drone aerial photography technology and geological radar and other equipment has made three-dimensional modeling of rock structural surfaces possible, greatly improving the accuracy and efficiency of data acquisition. In addition, with the help of modern image processing and deep learning technology, the automatic identification of rock structural surfaces has gradually become a research hotspot.

[0005] However, although the existing three-dimensional recognition methods have improved the accuracy of structural surface recognition to a certain extent, there are still many challenges in dealing with complex geological conditions and structural surface correlation analysis. For example, how to effectively use three-dimensional spatial data to accurately determine the relationship between structural surfaces, how to establish reasonable judgment criteria under different rock mass conditions, and how to eliminate the influence of external disturbance factors on the analysis results have not been completely solved. The rock mass structural surface recognition and correlation analysis method based on spatial constraints, especially setting the coplanarity judgment criteria and identifying the internal relationship between structural surfaces through spatial correlation analysis, provides a new research direction for solving these problems.

[0006] At present, the research on structural plane correlation analysis has gradually expanded from single identification to multi-dimensional and multi-angle analysis. Researchers have introduced spatial constraints and used geometric principles to judge the coplanarity and interaction of structural planes, thereby achieving more accurate correlation analysis of structural planes on continuous tunnel faces. However, existing research still has problems such as insufficient data processing capabilities, low algorithm efficiency, and poor adaptability in complex geological environments. Summary of the invention

[0007] In order to solve the above problems, the present invention proposes a rock structure surface identification and association analysis method and system based on spatial constraints, combines the Fisher probability distribution model and the maximum likelihood estimation algorithm to accurately classify all point clouds, sets the coplanarity judgment criteria through three-dimensional spatial constraints, and performs spatial association analysis on the structural surfaces on the continuous tunnel surface.

[0008] In some embodiments, the following technical solutions are adopted:

[0009] A rock mass structural plane identification and correlation analysis method based on spatial constraints, comprising:

[0010] Obtain the point cloud data of the continuous tunnel faces of the rock mass to be tested, calculate the point cloud normal of each tunnel face, and display the overall spatial distribution of the normal by drawing the joint pole diagram;

[0011] The density distribution of point cloud data is obtained based on the kernel density estimation algorithm, and the main direction of point cloud data is determined in combination with geological prior knowledge;

[0012] The Fisher probability distribution model is introduced, the initial concentration parameters are preset, and the probability of each central direction is preliminarily grouped according to the points to be classified; the concentration parameters of the Fisher probability distribution model are optimized to obtain the final structural surface grouping results;

[0013] Use DBSCAN density clustering algorithm to cluster the structural surface groups and extract each independent structural surface in each group of structural surfaces;

[0014] Spatial correlation analysis of independent structural surfaces on the continuous tunnel face is carried out based on three-dimensional spatial constraints.

[0015] Furthermore, the density distribution of point cloud data is obtained based on the kernel density estimation algorithm, specifically:

[0016] For each grid point cloud data (x, y), its density estimate Specifically:

[0017]

[0018] Among them, h x and h yis the smoothing parameter in the kernel density estimation, (x i ,y i ) is a given two-dimensional point cloud dataset, i = 1, 2, …, n; n is the number of sample points.

[0019] Furthermore, the main direction of the point cloud data is determined by combining geological prior knowledge, specifically:

[0020] After obtaining the density estimation value of each sample point, the filter window is used to detect the sample points whose density estimation value is higher than the set value;

[0021] For the detected sample points, first filter out the sample points whose density estimation values ​​are less than the minimum threshold, and then sort the remaining sample points from large to small according to the density estimation values. If the difference between the normal vectors of two adjacent sample points is less than the set angle value, the sample point with a large density estimation value is retained; finally, the filtered sample points are obtained as the main direction of the point cloud data.

[0022] Furthermore, the Fisher probability distribution model is introduced, the initial concentration parameters are preset, and the probability of each central direction is preliminarily grouped according to the points to be classified; specifically:

[0023] Set the initial concentration parameter k, calculate the probability that all point cloud data belongs to each main direction, and complete the preliminary classification according to the probability value;

[0024] The probability f(x; μ, k) of the point cloud data belonging to each main direction is calculated as:

[0025]

[0026] Where x is the position of the point cloud normal on the unit sphere, u is the main direction of the distribution, k is the concentration parameter, and θ is the normal angle between x and u.

[0027] Furthermore, the centralized parameters of the Fisher probability distribution model are optimized as follows:

[0028] For a given point cloud data normal vector set X = {x1, x2, ..., x n}, construct the log-likelihood function;

[0029] Substitute the probability density function of the Fisher probability distribution into the log-likelihood function and simplify it to obtain a simplified log-likelihood function; convert the simplified log-likelihood function into a negative log-likelihood function; and use the negative log-likelihood function as the loss function;

[0030] Based on the initial concentration parameter k, the sum of the loss functions of each group of structural surfaces is calculated;

[0031] The BFGS algorithm is used to update the k value of each group of structural surfaces, and the sum of the loss functions of each group of structural surfaces is recalculated; this process is repeated until the change in the sum of the loss functions of each group of structural surfaces is less than the set threshold, and the optimal concentration parameter k is obtained.

[0032] Furthermore, the DBSCAN density clustering algorithm is used to perform clustering within the structural face group, specifically:

[0033] ①: Randomly select a data point p that has not been checked in the point cloud data set D within the same set of structural surfaces, search its neighborhood, and if the number of objects contained is not less than MinPts, establish a new cluster C and add all the points in it to the candidate set N;

[0034] ②: For all data points q in the candidate set N that have not been processed, search their neighborhoods. If they contain at least MinPts objects, add these objects to N; if the data point q is not classified into any cluster, add the data point q to cluster C;

[0035] ③: Repeat step ② and continue to check the unprocessed objects in N until the candidate set N is empty;

[0036] ④: Repeat steps ① to ③. When an object is not included in any cluster, it is regarded as noise data until all objects are classified into a cluster or marked as noise.

[0037] ⑤: After the algorithm is completed, the clustering result cluster C is output, which is the independent structural surface finally extracted.

[0038] Furthermore, the spatial correlation analysis of the independent structural surfaces on the continuous tunnel face is carried out based on the three-dimensional spatial constraints, specifically:

[0039] After identifying the independent structural planes on each tunnel face, each structural plane is represented by its normal vector and a center point;

[0040] The coplanarity between the structural faces is determined based on the criteria set based on the three-dimensional space constraints, specifically:

[0041] Select any structural surface P from the first tunnel face, and calculate the angle θ between the normal vectors of all structural surfaces on the second tunnel face and P;

[0042] All structural surfaces on the second tunnel face whose normal vector angle θ is less than the set angle threshold are spatially adjusted so that their normal vectors are parallel to each other;

[0043] Calculate the distance between P and each candidate structural surface adjusted on the second face. If the distance is less than the set distance threshold, the structural surface P is spatially coplanar with the candidate structural surface. Repeat this process to obtain all structural surfaces on the second face that are spatially coplanar with the structural surface P.

[0044] According to the above process, all structural surfaces on the first tunnel face are traversed to obtain all coplanar structural surfaces on the two tunnel faces.

[0045] In other embodiments, the following technical solutions are adopted:

[0046] A rock mass structural surface recognition and correlation analysis system based on spatial constraints, comprising:

[0047] The data acquisition module is used to obtain the point cloud data of the continuous tunnel faces of the rock mass to be tested, calculate the point cloud normal of each tunnel face, and display the overall spatial distribution of the normal by drawing the joint pole diagram;

[0048] The point cloud main direction determination module is used to obtain the density distribution of point cloud data based on the kernel density estimation algorithm and determine the main direction of point cloud data in combination with geological prior knowledge;

[0049] The structural surface grouping module is used to introduce the Fisher probability distribution model, preset the initial concentration parameters, and preliminarily group the probability of each central direction according to the points to be classified; optimize the concentration parameters of the Fisher probability distribution model to obtain the final structural surface grouping results;

[0050] The module for clustering within structural surface groups is used to cluster structural surface groups using the DBSCAN density clustering algorithm and extract each independent structural surface in each group of structural surfaces;

[0051] The spatial correlation analysis module is used to perform spatial correlation analysis on independent structural surfaces on a continuous tunnel face based on three-dimensional spatial constraints.

[0052] In other embodiments, the following technical solutions are adopted:

[0053] A terminal device comprises a processor and a memory, wherein the processor is used to implement instructions; the memory is used to store a plurality of instructions, wherein the instructions are suitable for being loaded by the processor and executing the above-mentioned rock structure surface identification and association analysis method based on spatial constraints.

[0054] In other embodiments, the following technical solutions are adopted:

[0055] A computer-readable storage medium stores a plurality of instructions, wherein the instructions are suitable for being loaded by a processor of a terminal device and executing the above-mentioned rock mass structural surface identification and association analysis method based on spatial constraints.

[0056] Compared with the prior art, the present invention has the following beneficial effects:

[0057] (1) The present invention displays the overall spatial distribution of point cloud normals by drawing joint pole diagrams, and determines the number of groups and central directions based on the kernel density estimation algorithm and geological prior knowledge; introduces the Fisher probability distribution model, and performs preliminary classification of the probability of each central direction according to the points to be classified; and iteratively updates the concentration parameters through the maximum likelihood estimation algorithm to obtain the optimal concentration parameters for each group of structural surfaces, thereby achieving refined identification of structural surfaces.

[0058] (2) After obtaining the refined structural surface, the present invention performs correlation analysis on the continuous structural surfaces between different faces based on three-dimensional spatial constraints, thereby improving the accurate judgment of the spatial distribution of large structural surfaces. Compared with traditional methods, the present invention combines the structural surface normal vector and spatial geometric characteristics (such as the normal vector angle and the spatial distance) for accurate judgment, which can overcome the problem of structural surface fragmentation caused by the step-by-step excavation of the face, thereby completely restoring the spatial ductility and occurrence information of the structural surface. This not only helps to grasp the macroscopic distribution law of rock structural surfaces, but also provides high-precision data support for the construction of structural surface network models in underground engineering. In addition, the identification and association of continuous structural surfaces lays the foundation for subsequent engineering analyses such as geological stability analysis, rock mass classification, and potential sliding surface determination, greatly improving the engineering application value and reliability of rock structural surface identification.

[0059] Other features and advantages of additional aspects of the present invention will be given in part in the following description, and in part will become obvious from the following description, or will be learned through the practice of the present invention. BRIEF DESCRIPTION OF THE DRAWINGS

[0060] Figure 1 It is a flow chart of a rock mass structural surface identification and correlation analysis method based on spatial constraints in an embodiment of the present invention;

[0061] Fig. 2 (a)-(b) are respectively a joint extreme point diagram and a density extreme value diagram in an embodiment of the present invention;

[0062] Figure 3 Schematic diagram of the main direction after screening in an embodiment of the present invention;

[0063] Figure 4 Schematic diagram of data distribution of different centralized parameters K in an embodiment of the present invention;

[0064] Figure 5 Schematic diagram of probability classification in an embodiment of the present invention;

[0065] Figure 6 Schematic diagram of the DBSCAN clustering algorithm in an embodiment of the present invention;

[0066] Figure 7 This is a schematic diagram of spatial plane fitting in an embodiment of the present invention;

[0067] Figure 8 It is a schematic diagram of the spatial extensibility of the structural surface in an embodiment of the present invention;

[0068] Fig. 9 A schematic diagram of structural surface space adjustment in an embodiment of the present invention;

[0069] Fig.10 A schematic diagram of spatial coplanarity in an embodiment of the present invention;

[0070] Fig.11 It is a diagram of the extreme points of the rock mass joints at the tunnel face in the embodiment of the present invention;

[0071] Figure 12 (a)-(b) are the extreme value diagrams of the rock joints at the tunnel face before and after the density extreme value point screening, respectively;

[0072] Figure 13 (a)-(b) are the thermal images of rock mass point cloud probability distribution before and after the centralized parameter update;

[0073] Fig.14 It is a schematic diagram of the classification structure between structural surface groups in an embodiment of the present invention;

[0074] FIG. 15( a ) is a schematic diagram of the first group of structural surfaces before and after clustering in an embodiment of the present invention;

[0075] FIG15( b ) is a schematic diagram of the second group of structural surfaces before and after clustering in an embodiment of the present invention;

[0076] FIG15( c ) is a schematic diagram of the third group of structural surfaces before and after clustering in an embodiment of the present invention;

[0077] FIG15( d ) is a schematic diagram of the fourth group of structural surfaces before and after clustering in an embodiment of the present invention;

[0078] FIG. 15( e ) is a schematic diagram of the fifth group of structural surfaces before and after clustering in an embodiment of the present invention. DETAILED DESCRIPTION

[0079] It should be noted that the following detailed descriptions are illustrative and are intended to provide further explanation of the present application. Unless otherwise specified, all technical and scientific terms used in the present invention have the same meanings as those commonly understood by those skilled in the art to which the present application belongs.

[0080] It should be noted that the terms used herein are only for describing specific embodiments and are not intended to limit the exemplary embodiments according to the present application. As used herein, unless the context clearly indicates otherwise, the singular form is also intended to include the plural form. In addition, it should be understood that when the terms "comprise" and / or "include" are used in this specification, it indicates the presence of features, steps, operations, devices, components and / or combinations thereof.

[0081] Embodiment 1

[0082] In one or more embodiments, a rock mass structural surface identification and correlation analysis method based on spatial constraints is disclosed, combined with Figure 1 , specifically including the following process:

[0083] S101: Obtain continuous face point cloud data of the rock mass to be tested, calculate the point cloud normal, and display the overall spatial distribution of the normal by drawing a joint pole diagram.

[0084] In this embodiment, the point cloud data is collected by laser scanning, photogrammetry and other methods during the construction of a tunnel or underground project, and reflects the three-dimensional spatial distribution information of the rock structure surface (such as joints, cracks, etc.) on the tunnel face.

[0085] The point cloud normal is mainly used to represent the local surface direction of each point in the point cloud data. The calculation of the point cloud normal usually relies on neighborhood estimation and least squares plane fitting, which can be obtained according to the existing technology.

[0086] The point cloud normal can reflect local morphological information such as surface curvature and refraction, thereby describing the geometric characteristics of the point cloud; by using the normal direction, a pole diagram is generated to display the overall spatial distribution of the structural surface and realize the drawing of the joint pole diagram.

[0087] The joint pole diagram is a commonly used tool in rock mass structure analysis. It is drawn by projecting the normal poles of all point clouds. The radial lines represent the inclination (0° to 360°), and the concentric circles represent the inclination (0° to 90° from the center to the circumference). Figure 2(a) shows an example of a joint pole diagram.

[0088] S102: Obtain the density distribution of the point cloud data based on a kernel density estimation algorithm, and determine the main direction of the point cloud data in combination with geological prior knowledge.

[0089] Normally, the highest concentration point is selected as the representation of the central direction. In this embodiment, the extreme point with higher density is calculated by the kernel density estimation algorithm, and then the orientation of the main dominant structural surface is automatically identified as the corresponding independent Fisher center in combination with geological prior knowledge.

[0090] Kernel density estimation estimates the overall density by expanding each data point into a kernel function and summing all kernel functions. For a given two-dimensional point cloud dataset {(x1,y1),(x2,y2),…,(x n ,y n )}, its kernel density estimate can be expressed as:

[0091]

[0092] in, is the density estimate at point (x, y), n is the number of sample data points, and h x and h y is the bandwidth (smoothing parameter), which determines the degree of smoothing in the x and y directions respectively; K(·,·) is the two-dimensional kernel function.

[0093] The two-dimensional kernel function K(x, y) is a non-negative function whose integral is 1. The two-dimensional kernel function selected in this embodiment is a Gaussian kernel function, which has smoothness and robustness to noise. The specific formula is as follows:

[0094]

[0095] Bandwidth x and h y It is a smoothing parameter in kernel density estimation. A bandwidth that is too small will cause the estimation result to be too volatile (overfitting), and a bandwidth that is too large will cause the estimation result to be too smooth (underfitting). This embodiment automatically calculates the appropriate bandwidth size based on the empirical rule method.

[0096] For each grid point (x, y), substitute formula (2) into formula (1) to calculate the density estimate have to:

[0097]

[0098] In addition, when searching for density extreme sample points, the size of the filter window will affect the detection results of the local maximum. A larger filter window value will result in fewer and smoother local maxima being detected, while a smaller filter window value will result in more local maxima being detected. In order to accurately and comprehensively find all density extreme points, this embodiment selects a smaller filter window, and uses a set of randomly generated data to illustrate the method, and the calculation results are shown in Figure 2(b).

[0099] This embodiment combines geological prior knowledge to set the following criteria to screen sample points:

[0100] ①The density of sample points must be greater than the minimum threshold;

[0101] ② If the normal vector difference between two sample points is less than 30°, select the sample point with higher density.

[0102] The specific steps of sample point screening in this embodiment are as follows:

[0103] (1) Perform minimum value filtering on all sample points according to criterion ①;

[0104] (2) The remaining points are sorted from large to small according to density, and some sample points are filtered according to criterion ②.

[0105] After obtaining the sample points with high density, the normal vector direction of the sample points (normal vector n) is taken as the main direction. The main direction refers to the main spatial distribution direction of the rock structure surface (such as joints and cracks), which reflects the regularity and dominant orientation of the rock structure surface in three-dimensional space.

[0106] Figure 3 A schematic diagram of the main directions after screening is given.

[0107] S103: Introduce the Fisher probability distribution model, preset the initial concentration parameters, and preliminarily group the probability of each central direction according to the points to be classified; optimize the concentration parameters of the Fisher probability distribution model to obtain the final structural surface grouping result.

[0108] In this embodiment, Fisher probability distribution is a commonly used model for directional data, and is particularly suitable for data on a unit sphere. Its probability density function is given by the following formula:

[0109]

[0110] Among them, x is the position of the point cloud normal on the unit sphere, μ is the center direction of the distribution, k is the concentration parameter of the distribution, which is used to control the concentration of data points around u; θ is the normal angle between x and u. The normal angle is the angle between two normal vectors, which is used to measure their relative orientation in space. It reflects the degree of difference in the spatial direction of two surfaces (or planes).

[0111] Different concentration parameters k represent different degrees of dispersion, such as Figure 4 As shown in the figure, the smaller k is, the more dispersed the data is around the central direction. When performing preliminary grouping, we first assume that the concentration parameter k of each group of structural surfaces is 10 based on empirical values, calculate the probability that all points belong to each central direction, and complete the preliminary grouping based on the probability value. The probability classification results are shown in the figure. Figure 5 shown.

[0112] Assuming that the occurrence of each group of structural planes follows the Fisher distribution and the main directions of each group are known, it is necessary to find the best centralized parameter to improve the fitting quality of the data. The specific process of optimizing the centralized parameter of the Fisher probability distribution model in this embodiment is as follows:

[0113] S1031: For a given sample set X = {x1, x2, ..., x n}, construct a log-likelihood function; in this embodiment, the given sample set refers to the set of normal vectors in the point cloud data, specifically, the direction of the normal vector corresponding to each sample point; these normal vectors are calculated by local plane fitting at each point in the point cloud data, reflecting the spatial normal direction of the local surface of each point.

[0114] For a given sample set X = {x1, x2, ..., x n}, assuming that the samples are independent and identically distributed, the likelihood function represents the joint probability of all these sample data, that is, the product form of the likelihood function is as follows:

[0115]

[0116] The original likelihood function is the product of multiple probabilities. When the number of samples is large, these products may be very small, and direct calculation may lead to numerical underflow or calculation precision problems. Usually, taking the logarithm of the function and converting the product into a sum can effectively avoid such problems. The log-likelihood function L can be expressed as:

[0117]

[0118] S1032: Substitute the probability density function of the Fisher probability distribution into the log-likelihood function and simplify it to obtain a simplified log-likelihood function; convert the simplified log-likelihood function into a negative log-likelihood function; and use the negative log-likelihood function as a loss function.

[0119] Specifically, substituting the probability density function of the Fisher distribution, the log-likelihood function is:

[0120]

[0121] It can be further simplified to:

[0122]

[0123] Among them, θ i Represents the angle between the normal vector of the ith data point and the main direction.

[0124] In order to improve the speed and stability of calculation, this embodiment converts the maximization of the log-likelihood function into the minimization of the negative log-likelihood function problem. The negative log-likelihood function can be expressed as:

[0125]

[0126] S1033: Based on the initial lumped parameter k, calculate the sum of the loss functions of each group of structural surfaces - L(k).

[0127] S1034: Use the BFGS algorithm (Broyden-Fletcher-Goldfarb-Shanno algorithm) to update the k value of each group of structural surfaces and recalculate the sum of the loss functions of each group of structural surfaces; repeat this process until the change in the sum of the loss functions of each group of structural surfaces is less than the set threshold, and the optimal concentration parameter k is obtained.

[0128] Since the number of rock point clouds is usually hundreds of thousands or even millions, an efficient algorithm is needed to update the parameters. The BFGS algorithm has a fast convergence speed and good numerical stability. It performs well in dealing with large-scale optimization problems and is an ideal choice for solving such optimization problems.

[0129] The specific steps of iterative updating of the centralized parameter K in this embodiment are as follows:

[0130] ①After the preliminary classification is completed, the total loss function -L(k) is calculated, which is the sum of the loss values ​​of each group of structural surfaces;

[0131] ② Use the BFGS algorithm to update the K value of each group of structural surfaces and simultaneously update the -L(k) value;

[0132] ③ Repeat step ②. When the total loss function change value is less than the set threshold, end the update and output the optimal K value.

[0133] Based on the optimized optimal concentration parameter K value, the final structural surface grouping result is obtained through the Fisher probability distribution model.

[0134] S104: using the DBSCAN density clustering algorithm to perform clustering within the structural surface group, and extracting each independent structural surface in each group of structural surfaces;

[0135] The core idea of ​​the DBSCAN algorithm is to start from a core point and continuously expand to the density-reachable area, so as to obtain a maximized area containing core points and boundary points. Any two points in the area are densely connected. This method can find clusters of any shape in spatial data. Intuitively speaking, DBSCAN can find all dense areas in the sample points and treat them as clusters one by one. The DBSCAN algorithm mainly contains two parameters: search radius (eps) and minimum number of included points (MinPts). If the number of sample points within the search radius of a given object P is greater than or equal to MinPts, the object P is called a core point; for non-core point B, if B is within the search radius of any core point P, then sample B is called a boundary point; if B is not within the search radius of any core point P, then sample B is called a noise point.

[0136] Combination Figure 6 ,The specific process of DBSCAN algorithm clustering is as follows:

[0137] ① Randomly select an object p that has not been checked in the total data set D (a data point randomly selected from the point cloud data set, representing the starting point of the clustering process), search its neighborhood, and if the number of objects contained is not less than MinPts, establish a new cluster C and add all its points to the candidate set N;

[0138] The total data set D refers to a set of point cloud data within a certain set of structural surfaces, that is, the same set of point cloud data after Fisher distribution classification; these data points represent a part of the point cloud in the same direction and may belong to the same structural surface.

[0139] MinPts is a parameter in the DBSCAN algorithm, which indicates the minimum number of neighborhood points of a core point. MinPts is used to determine whether a point in a point cloud is a core point, that is, the neighborhood around the point must contain at least MinPts points to be considered as part of a structural surface. The setting of MinPts depends on the density of the point cloud data and can be adjusted according to the point cloud sampling accuracy and noise conditions of the actual project.

[0140] All the points here refer to all points in the neighborhood of a core point p whose Euclidean distance from the core point p is less than the set search radius ∈, such as: the core point p itself, and other points in the neighborhood of the core point (including core points and boundary points).

[0141] ② For all objects q in the candidate set N that have not been processed (the data points currently being checked when searching the neighborhood are used to determine whether they belong to the current cluster), search their neighborhoods. If they contain at least MinPts objects, add these objects to N; if q does not belong to any cluster, add q to cluster C;

[0142] ③ Repeat step ② and continue to check the unprocessed objects in N until the candidate set N is empty;

[0143] ④ Repeat steps ① to ③. When an object is not included in any cluster, it is regarded as noise data until all objects are classified into a cluster or marked as noise.

[0144] ⑤After the algorithm is completed, the clustering result cluster C is output, thus obtaining the final extracted independent structural surface.

[0145] After the point cloud of each independent structural surface is extracted, the final recognition result is obtained by plane fitting. This embodiment uses the random sampling consensus algorithm for spatial plane fitting (RANSAC). The basic principle of this algorithm is to randomly select a part of the data points from the data set in an iterative manner, calculate the model parameters, and verify whether the remaining data points conform to the model. Through multiple iterations, the model with the most inliers (data points that conform to the model) is selected as the optimal model.

[0146] The plane fitting process is as follows Figure 7 As shown, the blue point cloud represents the internal points, and the red point cloud represents the noise points. After the plane fitting is completed, the plane parameters A, B, C, and D are obtained. The plane equation can be expressed as follows:

[0147] A·x+B·y+C·z+D=0(10)

[0148] Among them, x, y, and z represent the coordinates of a point in three-dimensional space.

[0149] S105: Perform spatial correlation analysis on independent structural surfaces on the continuous tunnel face based on three-dimensional spatial constraints.

[0150] In rock mass engineering, the distribution of structural planes usually has local regularity, and the entire underground rock mass can be regarded as a huge crack network. During the tunnel excavation process, the cracks on the tunnel face are gradually exposed. For larger structural planes, there is a high probability that they will run through several consecutive tunnel faces, such as Figure 8 In view of this situation, this embodiment performs correlation analysis between continuous structural surfaces based on three-dimensional space constraints and extracts geometric parameters such as the occurrence of the structural surfaces.

[0151] Independent structural planes on each tunnel face have been identified above. Each structural plane can be represented by its normal vector and a center point. Suppose there are two structural planes distributed on different tunnel faces:

[0152] The equation of the first structural surface is A·x+B·y+C·z+D=0, and the normal vector is n1=(n 1x ,n 1y ,n 1z ), the center point P1=(x1,y1,z1); the equation of the second structural surface is A·x+B·y+C·z+D=0, and the normal vector is n2=(n 2x ,n 2y ,n 2z ), center point P2 = (x2, y2, z2).

[0153] Two criteria are set based on three-dimensional space constraints to determine the coplanarity between structural faces:

[0154] ① Angle threshold θ: the normal vector angle between the two structural surfaces is less than the angle threshold;

[0155] ②Distance threshold d: The distance between the two structural surfaces is less than the distance threshold.

[0156] The specific process of determining the coplanarity between structural surfaces is as follows:

[0157] (1) Select any structural surface P from the first tunnel face and calculate the angles between the normal vectors of all structural surfaces on the second tunnel face and P. The angle θ between the normal vectors can be calculated using the following formula:

[0158]

[0159] Among them, n1·n2 represents the dot product of the vector, and |n1| and |n2| represent the modulus of the vector respectively. The specific calculation is as follows:

[0160]

[0161] (2) All structural surfaces with normal vector angles less than the angle threshold are spatially adjusted so that their normal vectors are parallel to each other, which is convenient for calculating the spacing, such as Fig. 9 shown.

[0162] (3) Take the adjusted structural surfaces on the second face as candidate structural surfaces and calculate the distance between P and the candidate structural surfaces. The calculation formula is as follows:

[0163]

[0164] (4) If the distance D between the candidate structural surface and P is less than the distance threshold, the two structural surfaces are considered to be the same structural surface in three-dimensional space, such as Fig.10 shown.

[0165] (5) Repeat steps (1) to (4), traverse all structural surfaces of the first tunnel face, and search for structural surfaces that meet the coplanarity condition in the second tunnel face until all structural surfaces have been searched.

[0166] At the same time, the geometric parameters of the structural surface such as the strike and spacing can also be extracted, where the strike is the dip angle data of the structural surface; the spacing is the distance between each group of structural surfaces in each tunnel face, which can be calculated using formula (13).

[0167] The following is a specific case analysis using the method of this embodiment.

[0168] Project the point cloud normal calculation results to the upper unit hemisphere and convert them into dip angle data, and draw the joint pole diagram, such as Fig.11 shown.

[0169] By calculating the density based on the kernel density estimation algorithm, due to the complex rock structure of the face and the use of a small filter window value, multiple density extreme points will be detected, as shown in Figure 12(a). The density extreme points are screened based on the screening criteria to obtain the final main direction result, as shown in Figure 12(b).

[0170] The face data is initially classified according to the main direction using Fisher probability distribution, and the probability distribution heat map of some points is shown in Figure 13(a). As can be seen from the figure, when faced with complex rock mass data classification, although the probability values ​​of most points have obvious differences, there are still some points that are difficult to directly determine the belonging through probability values. For example, the probability value of data point 2 belonging to the third group is 0.39, and the probability value of belonging to the fourth group is 0.43. Since the normal calculation may still have inaccurate errors, it is not possible to absolutely determine that data point 2 belongs to the fourth group at this time. It is necessary to find the concentrated parameters that best match the actual rock mass structure of the project to complete the final classification and provide accurate parameters for the subsequent construction of the fracture network model. As can be seen from Figure 13(b), after iterative updating of the concentrated parameter K, the probability value of data point 2 belonging to the fourth group increases to 0.58, and the probability of belonging to the third group decreases to 0.31. Data point 2 can be more certainly classified as the fourth group. In addition, other data points can also be more clearly grouped by the probability values ​​calculated by the optimal parameters; the final grouping results are shown in Fig.14 shown.

[0171] After the inter-group classification is completed, each independent structural surface is extracted by the DBSCAN clustering algorithm, sorted according to the number of point clouds contained in each structural surface, and visualized by different colors. In addition, since some unnatural structural surfaces may be generated during the blasting of the face rock mass, these blasting surfaces are often scattered and contain a small number of point clouds. Therefore, the clusters with point clouds less than the set threshold are set to blue and directly deleted. From the first group to the fifth group, the independent structural surface extraction results of each group are shown in Figure 15 (a)-(e).

[0172] After the intra-group clustering is completed, the plane parameters of each independent structural surface are extracted by the RANSAC fitting algorithm, and converted into the occurrence information according to the spatial plane parameter extraction. In order to verify the accuracy and applicability of the identification method, the Compass geological toolbox in the point cloud processing software Cloudcompare is used to manually extract the occurrence of some structural surfaces. Compass is often used to interpret and analyze virtual outcrop models, and can intuitively measure and calculate the relevant information of the structural surface. Table 1 is a comparison table of the results of the method of this embodiment and manual extraction. It can be seen from the table that most of the results are very consistent with the actual occurrence. The average errors of the dip direction and dip angle of the 15 structural surfaces are 1.27° and 1.6°, respectively. Among them, the maximum deviation of the dip angle is 5°, and the maximum deviation of the dip angle is 6°. The large deviations of these results may be due to other complex properties such as surface roughness. In general, the occurrence information extracted after identification can meet the accuracy requirements of actual engineering applications.

[0173] Table 1 Comparison of manual extraction of occurrence and extraction results of this method

[0174]

[0175] Embodiment 2

[0176] In one or more embodiments, a rock mass structural surface identification and correlation analysis system based on spatial constraints is disclosed, comprising:

[0177] The data acquisition module is used to obtain the point cloud data of the continuous tunnel faces of the rock mass to be tested, calculate the point cloud normal of each tunnel face, and display the overall spatial distribution of the normal by drawing the joint pole diagram;

[0178] The point cloud main direction determination module is used to obtain the density distribution of point cloud data based on the kernel density estimation algorithm and determine the main direction of point cloud data in combination with geological prior knowledge;

[0179] The structural surface grouping module is used to introduce the Fisher probability distribution model, preset the initial concentration parameters, and preliminarily group the probability of each central direction according to the points to be classified; optimize the concentration parameters of the Fisher probability distribution model to obtain the final structural surface grouping results;

[0180] The module for clustering within structural surface groups is used to cluster structural surface groups using the DBSCAN density clustering algorithm and extract each independent structural surface in each group of structural surfaces;

[0181] The spatial correlation analysis module is used to perform spatial correlation analysis on independent structural surfaces on a continuous tunnel face based on three-dimensional spatial constraints.

[0182] The specific implementation of the above modules is the same as that in Example 1 and will not be described in detail.

[0183] Embodiment 3

[0184] In one or more embodiments, a terminal device is disclosed, which includes a processor and a memory, the processor is used to implement instructions; the memory is used to store multiple instructions, and the instructions are suitable for being loaded by the processor and executing the rock structure surface identification and association analysis method based on spatial constraints described in Example 1.

[0185] It should be understood that in this embodiment, the processor may be a central processing unit CPU, and the processor may also be other general-purpose processors, digital signal processors DSP, application-specific integrated circuits ASIC, off-the-shelf programmable gate arrays FPGA or other programmable logic devices, discrete gate or transistor logic devices, discrete hardware components, etc. The general-purpose processor may be a microprocessor or the processor may also be any conventional processor, etc.

[0186] The memory may include a read-only memory and a random access memory, and provide instructions and data to the processor. A portion of the memory may also include a non-volatile random access memory. For example, the memory may also store information about the device type.

[0187] In the implementation process, each step of the above method can be completed by an integrated logic circuit of hardware in a processor or an instruction in the form of software.

[0188] Embodiment 4

[0189] In one or more embodiments, a computer-readable storage medium is disclosed, in which a plurality of instructions are stored, wherein the instructions are suitable for being loaded by a processor of a terminal device and executing the rock structure surface identification and association analysis method based on spatial constraints described in Example 1.

[0190] Although the above describes the specific implementation mode of the present invention in conjunction with the accompanying drawings, it is not intended to limit the scope of protection of the present invention. Those skilled in the art should understand that various modifications or variations that can be made by those skilled in the art on the basis of the technical solution of the present invention without creative work are still within the scope of protection of the present invention.

Claims

1. A method for identifying and correlating rock mass structural surfaces based on spatial constraints, characterized in that: include: Obtain the point cloud data of the continuous tunnel faces of the rock mass to be tested, calculate the point cloud normal of each tunnel face, and display the overall spatial distribution of the normal by drawing the joint pole diagram; The density distribution of point cloud data is obtained based on the kernel density estimation algorithm, and the main direction of point cloud data is determined in combination with geological prior knowledge; The Fisher probability distribution model is introduced, the initial concentration parameters are preset, and the probability of each central direction is preliminarily grouped according to the points to be classified; the concentration parameters of the Fisher probability distribution model are optimized to obtain the final structural surface grouping results; Use DBSCAN density clustering algorithm to cluster the structural surface groups and extract each independent structural surface in each group of structural surfaces; Spatial correlation analysis of independent structural surfaces on the continuous tunnel face is carried out based on three-dimensional spatial constraints.

2. A method for identifying and analyzing rock mass structural surfaces based on spatial constraints according to claim 1, characterized in that: The density distribution of point cloud data is obtained based on the kernel density estimation algorithm, specifically: For each grid point cloud data (x, y), its density estimate Specifically: Among them, h x and h y is the smoothing parameter in the kernel density estimation, (x i ,y i ) is a given two-dimensional point cloud dataset, i = 1, 2, …, n; n is the number of sample points.

3. A method for identifying and analyzing rock mass structural surfaces based on spatial constraints according to claim 1, characterized in that: Combined with geological prior knowledge, the main direction of the point cloud data is determined as follows: After obtaining the density estimation value of each sample point, the filter window is used to detect the sample points whose density estimation value is higher than the set value; For the detected sample points, first filter out the sample points whose density estimation values ​​are less than the minimum threshold, and then sort the remaining sample points from large to small according to the density estimation values. If the difference between the normal vectors of two adjacent sample points is less than the set angle value, the sample point with a large density estimation value is retained; finally, the filtered sample points are obtained as the main direction of the point cloud data.

4. A method for identifying and analyzing rock mass structural surfaces based on spatial constraints according to claim 1, characterized in that: The Fisher probability distribution model is introduced, the initial concentration parameters are preset, and the probability of each central direction is preliminarily grouped according to the points to be classified; specifically: Set the initial concentration parameter k, calculate the probability that all point cloud data belongs to each main direction, and complete the preliminary classification according to the probability value; The probability f(x; μ, k) of the point cloud data belonging to each main direction is calculated as: Where x is the position of the point cloud normal on the unit sphere, u is the main direction of the distribution, k is the concentration parameter, and θ is the normal angle between x and u.

5. The method for identifying and analyzing rock mass structural surfaces based on spatial constraints according to claim 1, characterized in that: The centralized parameters of the Fisher probability distribution model are optimized as follows: For a given point cloud data normal vector set X = {x1, x2, ..., x n }, construct the log-likelihood function; Substitute the probability density function of the Fisher probability distribution into the log-likelihood function and simplify it to obtain a simplified log-likelihood function; convert the simplified log-likelihood function into a negative log-likelihood function; and use the negative log-likelihood function as the loss function; Based on the initial concentration parameter k, the sum of the loss functions of each group of structural surfaces is calculated; The BFGS algorithm is used to update the k value of each group of structural surfaces, and the sum of the loss functions of each group of structural surfaces is recalculated; this process is repeated until the change in the sum of the loss functions of each group of structural surfaces is less than the set threshold, and the optimal concentration parameter k is obtained.

6. The method for identifying and analyzing rock mass structural surfaces based on spatial constraints according to claim 1, characterized in that: The DBSCAN density clustering algorithm is used to perform clustering within the structural face group, specifically: ①: Randomly select a data point p that has not been checked in the point cloud data set D within the same set of structural surfaces, search its neighborhood, and if the number of objects contained is not less than MinPts, establish a new cluster C and add all the points in it to the candidate set N; ②: For all the data points q in the candidate set N that have not been processed, search their neighborhoods. If they contain at least MinPts objects, add these objects to N; If the data point q does not belong to any cluster, then add the data point q to cluster C; ③: Repeat step ② and continue to check the unprocessed objects in N until the candidate set N is empty; ④: Repeat steps ① to ③. When an object is not included in any cluster, it is regarded as noise data until all objects are classified into a cluster or marked as noise. ⑤: After the algorithm is completed, the clustering result cluster C is output, which is the independent structural surface finally extracted.

7. The method for identifying and analyzing rock mass structural planes based on spatial constraints according to claim 1, characterized in that: Based on three-dimensional spatial constraints, the spatial correlation analysis of independent structural surfaces on the continuous tunnel face is carried out, specifically: After identifying the independent structural planes on each tunnel face, each structural plane is represented by its normal vector and a center point; The coplanarity between the structural faces is determined based on the criteria set based on the three-dimensional space constraints, specifically: Select any structural surface P from the first tunnel face, and calculate the angle θ between the normal vectors of all structural surfaces on the second tunnel face and P; All structural surfaces on the second tunnel face whose normal vector angle θ is less than the set angle threshold are spatially adjusted so that their normal vectors are parallel to each other; Calculate the distance between P and each candidate structural surface adjusted on the second face. If the distance is less than the set distance threshold, the structural surface P is spatially coplanar with the candidate structural surface. Repeat this process to obtain all structural surfaces on the second face that are spatially coplanar with the structural surface P. According to the above process, all structural surfaces on the first tunnel face are traversed to obtain all coplanar structural surfaces on the two tunnel faces.

8. A rock mass structural surface identification and correlation analysis system based on spatial constraints, characterized in that: include: The data acquisition module is used to obtain the point cloud data of the continuous tunnel faces of the rock mass to be tested, calculate the point cloud normal of each tunnel face, and display the overall spatial distribution of the normal by drawing the joint pole diagram; The point cloud main direction determination module is used to obtain the density distribution of point cloud data based on the kernel density estimation algorithm and determine the main direction of point cloud data in combination with geological prior knowledge; The structural surface grouping module is used to introduce the Fisher probability distribution model, preset the initial concentration parameters, and preliminarily group the probability of each central direction according to the points to be classified; optimize the concentration parameters of the Fisher probability distribution model to obtain the final structural surface grouping results; The module for clustering within structural surface groups is used to cluster structural surface groups using the DBSCAN density clustering algorithm and extract each independent structural surface in each group of structural surfaces; The spatial correlation analysis module is used to perform spatial correlation analysis on independent structural surfaces on a continuous tunnel face based on three-dimensional spatial constraints.

9. A terminal device, comprising a processor and a memory, wherein the processor is used to implement instructions; and the memory is used to store multiple instructions, characterized in that: The instructions are suitable for being loaded by a processor and executing the rock structure surface identification and association analysis method based on spatial constraints as described in any one of claims 1-7.

10. A computer-readable storage medium storing a plurality of instructions, characterized in that: The instructions are suitable for being loaded by a processor of a terminal device and executing the rock structure surface identification and association analysis method based on spatial constraints as described in any one of claims 1 to 7.

Images