Mine weathering zone range identification method based on unmanned aerial vehicle point cloud data

Abstract

The present application provides a mine weathering zone range identification method based on unmanned aerial vehicle point cloud data, comprising: obtaining unmanned aerial vehicle point cloud data of a mine slope after blasting operation, calculating local geometric roughness of each data point, constructing an initial candidate weathering point cloud set of data points greater than a fragmentation threshold, and obtaining a blast hole connecting line strike vector and a step free surface strike vector, performing morphological feature analysis on the candidate set, identifying a three-dimensional crack trace network and a spatial direction vector of each crack trace, traversing the three-dimensional crack trace network with a sliding window of a preset size, calculating the spatial included angle between the crack trace direction and the blast hole connecting line strike and the step free surface strike, calculating the directional aggregation index, when the directional aggregation index is greater than a blasting control fragmentation threshold, determining that the data points in the window are pseudo-weathering points and removing them, clustering the remaining point cloud, and demarcating the true weathering zone range. The present application uses crack direction features to distinguish between natural weathering and blasting fragmentation, improving identification accuracy.

Description

Technical Field

[0001] This invention belongs to the field of image data processing and mining engineering geology technology, specifically a method for identifying the range of weathered zones in mines based on UAV point cloud data. Background Technology

[0002] Accurate delineation of the weathered zone is crucial for evaluating slope stability, designing slope angles, calculating stripping ratios, and developing reinforcement plans in open-pit mines. The degree of weathering directly determines the values ​​of the rock mass's physical and mechanical parameters. Deviations in boundary determination can lead to distorted engineering parameters and disruptions to production plans, or even inaccurate slope stability assessments, threatening production safety.

[0003] Currently, weathering zone identification mainly relies on two types of technologies. The first type is drilling and geophysical methods, including the wave velocity ratio method of deep-hole sonic logging, borehole coring and logging, and resistivity methods. These methods directly obtain deep parameters of the rock mass using boreholes as sampling units, offering high single-point accuracy and reliability. However, they suffer from low sampling density, high cost, and long cycle times, making it difficult to depict the continuous distribution of weathering zones on slope surfaces and meet the needs of rapid identification in large-area mining areas. The second type is UAV remote sensing methods, including the analysis of geometric features such as roughness, slope, and curvature of airborne LiDAR point clouds, multispectral / hyperspectral image alteration mineral mapping, and machine learning point cloud semantic segmentation. These methods can quickly acquire large-area, high-resolution three-dimensional information of the surface and have become an important direction for rapid weathering zone identification.

[0004] However, in the application of open-pit mine production slopes, the two methods are contradictory: traditional methods such as deep hole wave velocity ratio reflect the internal integrity of the rock mass. On the slope of a mining area that has just undergone bench blasting, the rock mass below about 1 meter from the ground surface is less disturbed by the blasting, and the test results are still hard and intact unweathered bedrock. However, the blasting creates a damaged area near the ground surface, and the surface of the bench is covered with blasting fissures and loose gravel. The geometric indicators such as roughness and volumetric fragmentation rate extracted from the point cloud are significantly increased. Based on this, the identification algorithm misclassifies the blasting damaged area as a heavily weathered zone, resulting in false weathering zone misclassification. Moreover, it is often distributed in a continuous pattern along the bench slope, directly contaminating the weathering zone boundary fitting results, leading to contradictions between the identification results and the actual borehole measurement data.

[0005] Existing point cloud-based weathering zone identification methods only utilize the magnitude of surface geometric features, such as roughness and other scalar indices, without considering the spatial directionality of fracture morphology. Natural weathering exhibits isotropic and randomly dispersed fracture directions; while blasting fracture, controlled by the arrangement of boreholes and the free surface, exhibits significant anisotropy, typically parallel to the borehole or free surface direction. Current technologies lack the means to distinguish between natural weathering and blasting fracture using fracture direction features, and to remove blasting-damaged areas before fitting weathering boundaries. Therefore, the blasting collapse effect becomes a significant source of error affecting the accuracy of weathering zone identification.

[0006] Therefore, distinguishing between fracturing caused by natural weathering and fracturing caused by human blasting, and eliminating the interference of blasting collapse effects on weathering zone identification, is a technical problem that urgently needs to be solved in this field. To this end, this invention provides a method for identifying the range of weathering zones in mines based on UAV point cloud data. Summary of the Invention

[0007] In order to overcome the shortcomings of the prior art, at least one technical problem raised in the background art is solved.

[0008] The technical solution adopted by this invention to solve its technical problem is: This invention provides a method for identifying the extent of weathered zones in mines based on UAV point cloud data, including: After blasting operations, UAV point cloud data of the mine slope is obtained, the local geometric roughness of each data point in the UAV point cloud data is calculated, data points with local geometric roughness greater than the fracture threshold are obtained, and an initial candidate weathering point cloud set is constructed. At the same time, the blast hole connection vector and the step free surface directional vector of the mine slope blasting operation are obtained. Morphological feature analysis was performed on the initial candidate weathering point cloud set to identify the three-dimensional fracture trace network of the mine slope and obtain the spatial direction vector of each fracture trace in the three-dimensional fracture trace network. The three-dimensional fracture trace network is traversed through a sliding window of a preset size. The spatial angle between the spatial direction vector of all fracture traces in the current sliding window and the direction vector of the borehole connection line and the direction vector of the free surface of the step is calculated. Based on the spatial angle, the directional aggregation index of the fracture traces in the current sliding window is calculated. Determine whether the directional aggregation index is greater than the blasting control fracture threshold. If it is, determine that the data point in the current sliding window is a pseudo weathering point and remove the pseudo weathering point from the initial candidate weathering point cloud set to obtain the real weathering point cloud data set. Spatial clustering is performed on data points in the real weathering point cloud dataset to delineate the actual range of the mine weathering zone.

[0009] As a further improvement of the present invention, the specific process for obtaining UAV point cloud data of the mine slope after the blasting operation is as follows: A drone equipped with an airborne lidar device was used to scan the mine slope and obtain point cloud data of the mine slope containing only the exposed rock surface.

[0010] As a further improvement of the present invention, the specific process for calculating the local geometric roughness of each data point in the UAV point cloud data is as follows: For any data point in the UAV point cloud data, a local spatial sphere is constructed using the roughness analysis radius; Obtain all adjacent data points within the local spatial sphere, and use the least squares method to perform local plane fitting on all adjacent data points to obtain the best fitting plane; Calculate the orthogonal projection distance from each adjacent data point within the local spatial sphere to the best-fit plane; The standard deviation or root mean square error of all orthogonal projected distances are used as the local geometric roughness of the data points.

[0011] As a further improvement of the present invention, the specific process of constructing the initial candidate weathering point cloud set is as follows: The local geometric roughness of each data point in the UAV point cloud data is obtained, and data points with local geometric roughness less than or equal to the fracture threshold are identified as unweathered intact bedrock points and removed. Data points with local geometric roughness greater than the fracture threshold are extracted and summarized to construct an initial candidate weathering point cloud set.

[0012] As a further improvement of the present invention, the specific process of simultaneously obtaining the blast hole connection vector and the bench free surface traverse vector of the mine slope blasting operation is as follows: The coordinates of the borehole openings in the same row are obtained using historical point cloud data before blasting. Based on principal component analysis, spatial straight line fitting is performed on the borehole coordinates of the same row of boreholes. The eigenvector corresponding to the largest eigenvalue is extracted and normalized. The normalized eigenvector is then used as the direction vector of the borehole connection. Based on UAV point cloud data, an edge detection algorithm is used to obtain the macroscopic contour trajectory of the top or bottom slope of the current step. Principal component analysis is used to perform spatial linear fitting on the macroscopic contour trajectory. The feature vector corresponding to the largest feature value is extracted and normalized. The normalized feature vector is then used as the free surface orientation vector of the step.

[0013] As a further improvement of the present invention, the specific process of identifying the three-dimensional fracture trace network of the mine slope is as follows: Obtain the neighboring data points of any data point in the initial candidate weathering point cloud set, and calculate the local principal curvature of the data point; Set a morphological curvature threshold, and obtain data points whose local principal curvature is greater than the morphological curvature threshold as crack feature points; A density-based 3D region growing algorithm is adopted. The connectivity search radius is set, and the fracture feature points with spatial distances smaller than the connectivity search radius are iteratively connected topologically to construct the initial fracture connectivity domain. For any data point within the initial fracture connectivity domain, a spatial shrinkage algorithm is used for iterative processing until the shrinkage reaches the maximum number of iterations to form a one-dimensional point skeleton. The nearest neighbor connection method is used to connect the one-dimensional point skeleton sequentially according to the spatial topology, refining the initial fracture connectivity domain into a three-dimensional spatial curve with a width of one pixel. All the refined three-dimensional spatial curves and their intersections together constitute a three-dimensional fracture trace network of the mine slope.

[0014] As a further improvement of the present invention, the specific process for obtaining the spatial direction vector of each fracture trace in the three-dimensional fracture trace network is as follows: Each discrete three-dimensional data point constituting a one-dimensional point-like skeleton is denoted as a skeleton data point. The topological connectivity of each skeleton data point is obtained. Skeleton data points with a topological connectivity greater than or equal to 3 are marked as intersections of three-dimensional spatial curves, and skeleton data points with a topological connectivity equal to 1 are marked as endpoints. Using the intersections of the three-dimensional spatial curves and the endpoints as breakpoints, the three-dimensional fracture trace network is topologically broken to obtain multiple initial fracture segments. Calculate the local curvature of each skeleton data point on each initial fracture segment. If the local curvature of a skeleton data point is greater than the data point mutation threshold, it is taken as the curvature mutation point. At the curvature mutation point, the initial fracture segment is broken by curvature to obtain multiple fracture traces. For each fracture trace, principal component analysis is used to perform spatial linear fitting on the three-dimensional spatial coordinates of all skeleton data points constituting the fracture trace. After extracting the eigenvector corresponding to the largest eigenvalue and normalizing it, the spatial direction vector of the fracture trace is obtained.

[0015] As a further improvement of the present invention, the specific process for calculating the spatial angle between the spatial direction vector of all fracture traces within the current sliding window and the direction vector of the borehole connection line and the direction vector of the free surface of the step is as follows: Create a sliding window and set its preset size; According to the preset spatial step size, the sliding window is translated and traversed in the three-dimensional space where the three-dimensional fracture trace network is located; At each traversal position, extract all crack traces that fall completely or partially within the current sliding window; For each fracture trace within the current sliding window, the spatial direction vector of the fracture trace is extracted. Using the three-dimensional vector dot product formula and the inverse cosine function, the first spatial angle between the spatial direction vector of the fracture trace and the direction vector of the borehole connection line, and the second spatial angle between the spatial direction vector of the fracture trace and the direction vector of the free surface of the step are calculated respectively.

[0016] As a further improvement of the present invention, the specific process of obtaining the real weathering point cloud data set is as follows: For each crack trace within the current sliding window, obtain the first spatial angle and the second spatial angle corresponding to the crack trace; Compare the size of the first spatial angle and the second spatial angle and select the minimum angle value as the main control angle for blasting; Set a parallel tolerance threshold, and count the number of fracture traces within the current sliding window whose blasting control angle is less than or equal to the parallel tolerance threshold. These are recorded as the number of directional fracture traces. The ratio of the number of directional fracture traces to the total number of fracture traces within the current sliding window is used as the directional aggregation index of fracture traces within the current sliding window. Compare the directional aggregation index calculated within the current sliding window with the blasting control fracture threshold. If the directional aggregation index is greater than the blasting control rupture threshold, then the data point in the current sliding window is determined to be a pseudo-weathering point; False weathering points are removed from the initial candidate weathering point cloud set, and all remaining data points after the removal operation are obtained to form the real weathering point cloud data set.

[0017] As a further improvement of the present invention, the specific process of delineating the actual range of the mine weathering zone is as follows: A density-based spatial clustering algorithm is used to perform three-dimensional spatial clustering on data points in the real weathering point cloud dataset. Data points that meet the conditions of clustering neighborhood radius and minimum number of points contained in the core point are aggregated into high-density weathering clusters, and isolated data points are filtered out. The three-dimensional spatial coordinates of all data points constituting the high-density weathering clusters are projected onto a two-dimensional projection surface parallel to the macroscopic slope surface of the mine slope. The concave hull extraction algorithm is used to extract the outer contour boundary lines of each high-density weathering cluster on the two-dimensional projection surface. The area enclosed by the outer contour boundary of all high-density weathering clusters is defined as the actual mining weathering zone.

[0018] The beneficial effects of this invention are as follows: This invention addresses the problem of false weathering misjudgment caused by blasting collapse effect. It utilizes the spatial difference between the random dispersion of natural weathering fissures and the control of blasting fracture direction by the borehole connection line and the free surface of the step. By using a directional aggregation index, it identifies and eliminates false weathering points caused by blasting damage, thereby eliminating the interference of blasting damage areas on weathering zone identification and improving the accuracy of weathering zone range identification. Furthermore, by using local geometric roughness pre-screening and fissure trace skeletonization, it reduces the data scale of subsequent processing and improves identification efficiency.

[0019] This invention enables rapid, non-contact, large-area weathering zone identification based on UAV point cloud data, with low cost and short cycle. By performing spatial clustering and concave hull boundary extraction on the real weathering point cloud after removing false weathering points, continuous and regular weathering zone boundaries can be obtained, providing a reliable basis for open-pit mine slope stability evaluation, bench slope angle design, and stripping ratio calculation. Attached Figure Description

[0020] The invention will now be further described with reference to the accompanying drawings.

[0021] Figure 1 This is a flowchart of the steps in the mine weathering zone range identification method based on UAV point cloud data of the present invention; Figure 2 This is a logic judgment diagram of the mine weathering zone range identification method based on UAV point cloud data of the present invention. Detailed Implementation

[0022] To make the technical means, creative features, objectives and effects of this invention easier to understand, the invention will be further described below in conjunction with specific embodiments.

[0023] Example 1: As Figure 1 and Figure 2 As shown in the embodiment of the present invention, the method for identifying the range of weathered zones in mines based on UAV point cloud data includes: Step 1: Obtain UAV point cloud data of the mine slope after blasting operation, calculate the local geometric roughness of each data point in the UAV point cloud data, obtain data points with local geometric roughness greater than the fracture threshold, construct an initial candidate weathering point cloud set, and at the same time obtain the blast hole connection vector and the step free surface orientation vector of the mine slope blasting operation. In some embodiments, the specific process of obtaining UAV point cloud data of the mine slope after blasting operations in step one is as follows: After the blasting operation on the mine slope is completed and the blasting dust has settled naturally, a drone equipped with an airborne lidar device is used to perform slope-following flight scanning of the target mine slope to obtain the original three-dimensional point cloud data. The original three-dimensional point cloud data was processed by mathematical morphology filtering algorithm to remove non-ground noise such as vegetation, construction machinery and suspended dust, and obtain UAV point cloud data of the mine slope after blasting operation containing only the exposed rock surface.

[0024] In other embodiments, the specific process of using a mathematical morphological filtering algorithm to filter the original three-dimensional point cloud data, removing non-ground noise such as vegetation, construction machinery, and suspended dust, and obtaining UAV point cloud data containing only the exposed rock surface of the mine slope after blasting is as follows: Let the original 3D point cloud data be denoted as set. ,in The total number of original data points. Represents the three-dimensional spatial coordinates of the i-th data point. The elevation value of the data point; Set a radius as The filtering window, denoted as structuring element. ; For any data point in the original 3D point cloud data The horizontal coordinates of the data points Centered on the structural element, obtain the elements falling into the structural element. For all locally adjacent points within the range, perform morphological erosion operation, calculate the lowest elevation value within the filter window, and denot it as the erosion elevation. The calculation formula is: In the formula, For falling into structural elements The horizontal coordinate of the j-th adjacent point within the range, This corresponds to the original elevation value; The physical meaning of the corrosion operation is to forcibly pull the elevation reference plane within the window down to the local lowest point, thereby eliminating the noise from the vegetation canopy or construction machinery. Based on the corrosion elevation Perform morphological dilation operation to calculate the opening elevation. The calculation formula is: In the formula, For falling into structural elements The horizontal coordinate of the kth neighboring point within the range. The corrosion elevation corresponding to the adjacent points; The physical meaning of the expansion operation is to restore the original elevation profile of the real terrain while keeping isolated bulge noise in a state of elimination. Calculate data points Original elevation value With the opening operation elevation elevation difference ,Right now: ; Set elevation difference threshold If the elevation difference Then determine the data point Non-ground point noise is excluded; if the elevation difference Then determine the data point Represent the actual ground points and preserve them; Traverse the set After all data points have been judged and eliminated as described above, all the remaining ground points constitute the UAV point cloud data of the mine slope after blasting operations, which only includes the exposed rock surface.

[0025] It should be noted that the radius of the filtering window is specifically determined by the spatial resolution of the UAV point cloud data and the physical scale of the noise of the ground objects to be filtered out, for example, it is set to 0.5 meters to 2.0 meters; the specific setting of the elevation difference threshold usually depends on the macroscopic topographic undulation of the target mine slope and the estimated height of non-ground interference objects, for example, it is set to 0.3 meters to 1.5 meters. This invention does not limit the radius of the filtering window or the elevation difference threshold.

[0026] In some embodiments, the specific process of calculating the local geometric roughness of each data point in the UAV point cloud data in step one is as follows: For any data point in the UAV point cloud data, a local spatial sphere is constructed with the data point as the center and a preset roughness analysis radius. In other embodiments, the specific process of constructing a local spatial sphere with the data point as the center and a preset roughness analysis radius for any data point in the UAV point cloud data is as follows: For any target data point in the UAV point cloud data With the target data points Centered on the sphere, with a preset roughness analysis radius. Construct a local spatial sphere; Traverse the UAV point cloud data and calculate the relationship between the remaining data points and the target data point. The Euclidean distance is used to extract the roughness analysis radius if and only if the Euclidean distance is less than or equal to the specified radius. All data points constitute the set of adjacent data points. Where K is the total number of adjacent data points within the local spatial sphere. This represents the three-dimensional spatial coordinates of k adjacent data points.

[0027] All adjacent data points within the local spatial sphere are obtained, and the least squares method is used to perform local microplane fitting on all adjacent data points to obtain the best fitting plane. Calculate the orthogonal projection distance from each adjacent data point within the local spatial sphere to the best-fit plane; use the standard deviation or root mean square error of all the orthogonal projection distances as the local geometric roughness of the data points.

[0028] In other embodiments, the specific process of obtaining all adjacent data points within the local spatial sphere, performing local microplane fitting on all adjacent data points using the least squares method to obtain the best-fit plane, calculating the orthogonal projection distance from each adjacent data point within the local spatial sphere to the best-fit plane, and using the standard deviation or root mean square error of all the orthogonal projection distances as the local geometric roughness of the data points is as follows: Calculate the set of adjacent data points 3D geometric center coordinates The calculation formula is: Based on the principle of total least squares, the set of adjacent data points is constructed. of covariance matrix : For the covariance matrix Perform eigenvalue decomposition and extract the eigenvector corresponding to the smallest eigenvalue as the unit normal vector of the best-fit plane. ,in ; Therefore, the spatial analytical equation of the best-fit plane is obtained as follows: Calculate the set of adjacent data points Each adjacent data point Orthogonal projection distance to the best-fit plane : Calculate all the orthogonal projection distances. The root mean square error is used as the target data point. Local geometric roughness The calculation formula is: , It should be noted that the roughness analysis radius is usually set to 0.2 meters to 0.5 meters. The range of 0.2 meters to 0.5 meters can effectively filter out the low-frequency undulation interference of the macroscopic step terrain of the mine slope, while accurately capturing the microscopic fracture texture of the rock surface caused by weathering or blasting.

[0029] In some embodiments, the specific process of acquiring data points with local geometric roughness greater than the fracturing threshold and constructing an initial candidate weathering point cloud set in step one is as follows: Historical point cloud roughness statistics are obtained for known unweathered intact bedrock sample areas within the target mining area, and the 95th percentile of the historical point cloud roughness statistics is used as the fracturing threshold. The local geometric roughness of each data point in the UAV point cloud data is obtained, and the data points whose local geometric roughness is less than or equal to the fracture threshold are identified as unweathered intact bedrock points and are removed. Data points with local geometric roughness greater than the fracture threshold are extracted and summarized to construct the initial candidate weathering point cloud set.

[0030] It should be understood that the surface of intact bedrock is relatively flat and has extremely low roughness; when the local geometric roughness is greater than the fracture threshold, it indicates that the surface of the rock mass has been fractured. However, at this time, it is impossible to distinguish whether the weathering is due to natural geological weathering or blasting. Therefore, it is uniformly included in the initial candidate weathering point cloud set as the basic data for subsequent false detection.

[0031] In some embodiments, the specific process of simultaneously obtaining the blast hole alignment vector and the bench free surface alignment vector in step one is as follows: The coordinates of the borehole openings in the same row are obtained using historical point cloud data before blasting. Based on principal component analysis, spatial straight line fitting is performed on the borehole coordinates of the same row of boreholes. The eigenvector corresponding to the largest eigenvalue is extracted, and after normalization, the normalized eigenvector is used as the direction vector of the borehole connection.

[0032] In other embodiments, the specific process of performing spatial straight-line fitting on the borehole coordinates of the same row of boreholes, and using the fitted spatial straight-line direction as the direction vector of the borehole connection line, is as follows: The obtained borehole coordinates of the same row of blast holes are set as a set. ,in This represents the total number of boreholes in the same row. Represents the three-dimensional spatial coordinates of the m-th borehole; Calculate the set 3D geometric center coordinates The calculation formula is as follows: Based on the principle of principal component analysis (PCA), the set is constructed. of covariance matrix : For the covariance matrix Eigenvalue decomposition is performed to extract the unit eigenvector corresponding to the largest eigenvalue. The normalized unit eigenvector corresponding to the largest eigenvalue is then used as the borehole alignment vector. .

[0033] Based on the UAV point cloud data, an edge detection algorithm is used to obtain the macroscopic contour trajectory of the top or bottom line of the current step; principal component analysis is used to perform spatial linear fitting on the macroscopic contour trajectory, and the feature vector corresponding to the largest feature value is extracted and normalized. The normalized feature vector is then used as the free surface orientation vector of the step.

[0034] In other embodiments, the specific process of obtaining the macroscopic contour trajectory of the top or bottom slope of the current step based on the UAV point cloud data, and performing linear fitting on the macroscopic contour trajectory, using the fitted straight line direction as the direction vector of the free surface of the step, is as follows: An edge detection algorithm is used to extract discrete trajectory points on the top or bottom line of the current step from the UAV point cloud data, and these points are denoted as a set. Where S is the total number of trajectory points. Represents the three-dimensional spatial coordinates of the s-th trajectory point; Calculate the set 3D geometric center coordinates The calculation formula is: , Construct the set of covariance matrix : , For the covariance matrix Eigenvalue decomposition is performed to extract the unit eigenvector corresponding to the largest eigenvalue. The unit eigenvector corresponding to the largest eigenvalue is then normalized, and the normalized unit eigenvector corresponding to the largest eigenvalue is used as the orientation vector of the free surface of the step. .

[0035] The feature vectors are normalized to obtain unit direction vectors with a magnitude of 1, which represent direction but do not contain length information. The purpose of normalization is to eliminate the interference of vector magnitude, coordinate scale, and the number of data points involved in the fitting on subsequent spatial angle calculations, so that the direction vectors obtained from different borehole rows and different step regions have a unified measurement benchmark, thereby ensuring the accuracy of spatial angle threshold judgment and direction consistency comparison.

[0036] It should be understood that in mine slope blasting operations, the superposition and reflection of explosive stress waves are mainly controlled by the arrangement direction of the blast holes and the direction of the free face. The damage cracks generated by the blast will inevitably develop in a directional manner along the direction vector of the line connecting the blast holes or the direction vector of the free face of the step. These two vectors provide an absolute spatial reference for distinguishing between isotropic natural weathering cracks and anisotropic blast damage cracks.

[0037] Step 2: Perform morphological feature analysis on the initial candidate weathering point cloud set to identify the three-dimensional fracture trace network of the mine slope and obtain the spatial direction vector of each fracture trace in the three-dimensional fracture trace network. In some embodiments, the specific process of performing morphological feature analysis on the initial candidate weathering point cloud set to identify the three-dimensional fracture trace network of the mine slope in step two is as follows: Traverse each data point in the initial candidate weathering point cloud set, obtain the neighboring data points of any data point, and calculate the local principal curvature of the data point; A morphological curvature threshold is set, and data points with local principal curvature greater than the morphological curvature threshold are extracted as crack feature points; It should be noted that the specific setting of the morphological curvature threshold is usually determined based on the spatial resolution of the UAV point cloud data and the expected crack size of the target crack. For example, it can be set to 2 to 4 times the average curvature of the local neighborhood point cloud, or to 0.15 to 0.30 of the normalized surface curvature value. This invention does not limit this. The physical significance of the morphological curvature threshold is to identify micro-crack areas on the rock surface that exhibit obvious depressions, broken lines, or V-shaped grooves, thereby filtering out ordinary surface undulations.

[0038] A density-based three-dimensional region growing algorithm is adopted. A connectivity search radius is set, and the fracture feature points with a spatial distance smaller than the connectivity search radius are iteratively connected topologically to construct an initial fracture connectivity domain with a certain width. In other embodiments, the specific process of using a density-based three-dimensional region growing algorithm to topologically connect spatially close fracture feature points and construct an initial fracture connected domain of a certain width is as follows: Set the connectivity search radius ; Randomly select an unprocessed seed point from the set of fracture feature points, and use the seed point as the center of a sphere, with the connected search radius as the boundary. Perform a spherical neighborhood search to obtain all crack feature points that fall within the spherical neighborhood; All the acquired crack feature points are assigned to the connected cluster where the current seed point is located, and these newly added crack feature points are used as new seed points. The spherical neighborhood search and assignment operation is iteratively executed until no new crack feature points are added to the current connected cluster. Repeat the above seed point selection and iterative search process until all data points in the fracture feature point set have been processed, resulting in several independent connected clusters. Each connected cluster is used as an initial fracture connected domain with a certain width.

[0039] It should be noted that the connectivity search radius is usually set based on the spatial resolution of the UAV point cloud data and the expected discontinuity spacing of the micro-cracks on the rock surface. For example, it is set to 2 to 5 times the average point spacing of the point cloud. The specific value can be set according to the actual situation, and this invention does not limit it.

[0040] For any data point within the initial fracture connected domain, a spatial shrinkage algorithm is used for iterative processing until the number of shrinkage processes reaches a preset maximum number of iterations to converge and form a one-dimensional point skeleton. The nearest neighbor connection method is used to connect the converged one-dimensional point skeleton sequentially according to spatial topological relationships, thereby refining the initial fracture connected domain with a certain width into a three-dimensional spatial curve with a single pixel width. All the refined three-dimensional spatial curves and their intersections together constitute the three-dimensional fracture trace network of the mine slope.

[0041] In other embodiments, the process of iteratively processing any data point within the initial fracture connected domain using a spatial shrinkage algorithm until the number of shrinkage processes reaches a preset maximum number of iterations to converge and form a one-dimensional point skeleton, and then using the nearest neighbor connection method to sequentially connect the converged one-dimensional point skeleton according to spatial topological relationships, thereby refining the connected domain with a certain width into a three-dimensional spatial curve with a single pixel width, and the specific process of forming the three-dimensional fracture trace network of the mine slope by all the refined three-dimensional spatial curves and their intersections, is as follows: The Laplace shrinkage algorithm is used as the spatial shrinkage algorithm to iteratively shrink the data points in the initial fracture connected domain. For any data point within the initial fracture connected domain In the t-th iteration, its position update formula is: , In the formula, and These are the three-dimensional spatial coordinates of the data point at the t-th and t+1-th iterations, respectively. The set of local neighborhood points of the data point. To shrink the step size weight, Spatial distance weights; Through multiple iterations, the data points within the initial fracture connectivity domain shrink toward their local geometric center until they converge to form a one-dimensional point-like skeleton. The nearest neighbor method is used to connect the converged one-dimensional point skeleton sequentially according to the spatial topology to generate a three-dimensional spatial curve with a single pixel width. All the three-dimensional spatial curves and their spatial intersection nodes are obtained, which together constitute the three-dimensional fracture trace network of the mine slope.

[0042] In some embodiments, the specific process of obtaining the spatial direction vector of each fracture trace in the three-dimensional fracture trace network in step two is as follows: Each discrete three-dimensional data point constituting the one-dimensional point skeleton is denoted as a skeleton data point. The topological connectivity of each skeleton data point is obtained. Skeleton data points with a topological connectivity greater than or equal to 3 are marked as intersections of the three-dimensional spatial curves, and skeleton data points with a topological connectivity equal to 1 are marked as endpoints. The three-dimensional fracture trace network is topologically broken using the intersections of the three-dimensional spatial curves and the endpoints as breakpoints to obtain multiple initial fracture segments. Calculate the local curvature of each skeleton data point on each initial fracture segment. If the local curvature of a skeleton data point is greater than the set data point mutation threshold, then the skeleton data point is taken as the curvature mutation point, and the initial fracture segment is broken by curvature at the curvature mutation point. It should be noted that the physical meaning of the data point mutation threshold is to identify the spatial location where the fracture direction is drastically deflected, so as to ensure that the fracture segment obtained after the break meets the flatness requirements of the subsequent principal component analysis (PCA) spatial line fitting. Its specific value can be set according to the actual situation, and the present invention does not limit it. After the above topology breaking and curvature breaking processes, multiple independent and approximately straight crack traces are obtained.

[0043] For each of the aforementioned crack traces, principal component analysis (PCA) is used to perform spatial linear fitting on the three-dimensional spatial coordinates of all the skeleton data points constituting the crack trace. Specifically, this involves constructing the covariance matrix of the crack trace coordinates and performing eigenvalue decomposition. The eigenvector corresponding to the largest eigenvalue is extracted, and after normalization, it is used as the spatial direction vector of the crack trace.

[0044] In other embodiments, the specific process of using principal component analysis to perform spatial linear fitting on the three-dimensional spatial coordinates of all skeleton data points, extracting the eigenvector corresponding to the largest eigenvalue, and then normalizing it to serve as the spatial direction vector of the crack trace is as follows: Let the set be the set of all the skeleton data points that constitute a certain fracture trace. ,in This represents the total number of skeleton data points on the crack trace. Represents the three-dimensional spatial coordinates of the u-th skeleton data point; Calculate the set 3D geometric center coordinates The calculation formula is: Construct the set of covariance matrix : For the covariance matrix Perform eigenvalue decomposition to obtain three eigenvalues ​​and their corresponding eigenvectors; Extract the eigenvector corresponding to the largest eigenvalue. Normalize it to obtain a unit vector. : , the unit vector This serves as the spatial direction vector of the aforementioned crack trace.

[0045] It should be noted that the physical meaning of the spatial direction vector is to accurately characterize the actual extension direction of each micro-crack on the rock surface in three-dimensional space; extracting the spatial direction vector is the core premise of the entire algorithm for eliminating falsehoods and retaining truth, and it provides a vector basis for subsequent judgment on whether weathering is controlled by blasting stress or by natural weathering.

[0046] Step 3: Traverse the three-dimensional fracture trace network with a sliding window of a preset size, calculate the spatial angle between the spatial direction vector of all fracture traces in the current sliding window and the direction vector of the borehole connection line and the direction vector of the free surface of the step, and calculate the orientation aggregation index of the fracture traces in the current sliding window based on the spatial angle. In some embodiments, the specific process of step three, which involves traversing the three-dimensional fracture trace network using a sliding window of a preset size and calculating the spatial angle between the spatial direction vector of all fracture traces within the current sliding window and the orientation vector of the borehole connection line and the orientation vector of the free surface of the step, is as follows: Construct a sliding window and set the preset size of the sliding window (i.e., the side length of the cube); it should be noted that the preset size is usually set to a three-dimensional space that can cover at least 3 to 5 steps of blast hole spacing, for example, the side length is set to 3 meters to 10 meters, to ensure that the window contains a sufficient number of fracture samples for statistical analysis. This invention does not limit this. According to the preset spatial step size (i.e., the distance the sliding window translates in space each time), the sliding window is translated and traversed in the three-dimensional space where the three-dimensional fracture trace network is located. At each traversal position, extract all crack traces that fall completely or partially within the current sliding window; For each fracture trace within the current sliding window, extract the spatial direction vector of the fracture trace calculated in step two. Using the three-dimensional vector dot product formula and the inverse cosine function, calculate the first spatial angle between the spatial direction vector of the fracture trace and the direction vector of the borehole connection line, and the second spatial angle between the spatial direction vector of the fracture trace and the direction vector of the free surface of the step, thereby obtaining the set of spatial angles of all fracture traces within the current sliding window.

[0047] In other embodiments, the specific process of extracting the spatial direction vector of the fracture trace calculated in step two for each fracture trace within the current sliding window, and using the three-dimensional vector dot product formula and the inverse cosine function to calculate the first spatial angle between the spatial direction vector and the borehole connection vector, and the second spatial angle between the spatial direction vector and the step free surface orientation vector, is as follows: Suppose there are a total of Extract the first crack trace. Crack traces Spatial direction vector ; Obtain the direction vector of the borehole connection determined in step one. Vector of the free surface of the step ; The spatial direction vector is calculated using the three-dimensional vector dot product formula and the inverse cosine function. Vector of the line connecting the boreholes First spatial angle : Calculate the spatial direction vector With respect to the direction vector of the free surface of the step The second spatial angle : Iterate through all the cubes within the current sliding window. By drawing a fracture trace, the first spatial angle set is obtained. and the set of angles in the second space Together, they constitute the set of spatial angles of all crack traces within the current sliding window.

[0048] In some embodiments, the specific process of calculating the directional aggregation index of the crack traces within the current sliding window based on the spatial angle in step three is as follows: For each crack trace within the current sliding window, obtain the first spatial angle and the second spatial angle corresponding to the crack trace; Since the rock fracture surface caused by the blasting tensile stress wave is usually parallel to the borehole connection line or the step free surface, the size of the first spatial angle and the second spatial angle are compared, and the minimum angle value between the two is selected as the blasting main control deflection angle of the fracture trace. Set a parallel tolerance threshold, which is usually set to 10 to 15 degrees, to compensate for the slight directional deviation caused by the heterogeneity of the rock mass in actual blasting. The number of fracture traces within the current sliding cube-shaped moving window whose blasting control angle is less than or equal to the parallel tolerance threshold is counted and recorded as the number of directional fracture traces. Calculate the ratio of the number of directional fracture traces to the total number of fracture traces in the current sliding window, and use the ratio (the ratio value is between 0 and 1) as the directional aggregation index of fracture traces in the current sliding window.

[0049] It should be understood that the closer the directional aggregation index is to 1, the more concentrated the direction of the fractures within the window is on the blasting surface, exhibiting strong anisotropy; the closer the directional aggregation index is to 0, the more randomly the fractures within the window are dispersed, exhibiting isotropic natural geological weathering characteristics.

[0050] Step 4: Determine whether the directional aggregation index is greater than the blasting control fracture threshold. If it is greater, determine that the data points in the current sliding window are pseudo-weathering points and remove them from the initial candidate weathering point cloud set to obtain the real weathering point cloud data set. In some embodiments, in step four, the specific process of determining whether the directional aggregation index is greater than the blasting control fracture threshold, and if it is greater, then determining that the data point in the current sliding window is a pseudo-weathering point is as follows: Historical fracture statistics of fresh bench areas (i.e., blast damage areas) in the target mine slope that are known to have been affected by blasting but have not undergone geological weathering are obtained. The directional aggregation index of the fresh bench areas is calculated, and the lower limit extreme value or 5% quantile of its statistical distribution is extracted as the blasting control fracture threshold. Compare the directional aggregation index calculated within the current sliding window with the blasting control fracture threshold. If the directional aggregation index is greater than the blasting-controlled fracture threshold, it indicates that the fracture network in the current sliding window is highly parallel to the borehole connection line or the free surface of the step in the spatial direction, and the weathering mechanism is mainly dominated by the directional tearing of the blasting stress wave, rather than the natural geological weathering. Based on this judgment logic, all data points within the space covered by the current sliding window are marked as fracture points controlled by the blasting structural surface, and blasting damage attribute labels are added to them, thereby determining that the data points in the current sliding window are pseudo-weathering points.

[0051] In some embodiments, the specific process of removing false weathering points from the initial candidate weathering point cloud set in step four to obtain the real weathering point cloud data set is as follows: Traverse the attribute labels of all data points in the initial candidate weathering point cloud set; The pseudo-weathering points with blast damage attribute labels are removed from the initial candidate weathering point cloud set by using set difference operation or Boolean elimination algorithm; All remaining data points after the removal operation are obtained. These remaining data points satisfy the condition of high local geometric roughness and eliminate the interference of directional blasting fracturing, thus forming the real weathering point cloud data set formed solely by natural geological weathering.

[0052] Step 5: Perform spatial clustering on the data points in the real weathering point cloud dataset to delineate the actual range of the mine weathering zone.

[0053] In some embodiments, the specific process of spatially clustering the data points in the real weathering point cloud dataset to delineate the actual range of the mine weathering zone in step five is as follows: A density-based spatial clustering algorithm is used to perform three-dimensional spatial clustering on the data points in the real weathered point cloud dataset; Set the cluster neighborhood radius and the minimum number of points contained in the core point. The cluster neighborhood radius is usually set according to the minimum continuous physical span of the weathered zone in the mining geological specifications, for example, it is set to 1 meter to 3 meters. Traverse the real weathering point cloud data set, aggregate data points whose spatial distance is less than the clustering neighborhood radius and whose local density reaches the minimum number of points contained in the core point into the same high-density weathering cluster, and filter out isolated data points outside the high-density weathering cluster as random geological noise. For each high-density weathering cluster generated by clustering, the three-dimensional spatial coordinates of all data points constituting the high-density weathering cluster are extracted, and the three-dimensional spatial coordinates are projected onto a two-dimensional projection plane parallel to the macroscopic slope surface of the mine slope. The concave hull extraction algorithm is used to obtain the outer contour boundary lines of each high-density weathering cluster on the two-dimensional projection surface; The macroscopic area encompassed by the outer contour boundary lines of all high-density weathering clusters is determined and delineated as the actual range of the mine weathering zone.

[0054] It should be understood that by eliminating the interference of scattered and isolated natural weathering fragments through density clustering, discrete micro data points are transformed into the main body of the weathering zone with macroscopic geological continuity, and finally the boundary range of the weathering zone can be directly used for mine safety production and slope support design.

[0055] The foregoing has shown and described the basic principles, main features, and advantages of the present invention. Those skilled in the art should understand that the present invention is not limited to the above embodiments. The embodiments and descriptions in the specification are merely illustrative of the principles of the invention. Various changes and modifications can be made to the invention without departing from its spirit and scope, and all such changes and modifications fall within the scope of the present invention as claimed. The scope of protection of the present invention is defined by the appended claims and their equivalents.

Claims

1. A method for identifying the range of weathered zones in mines based on UAV point cloud data, characterized in that: include: After blasting operations, UAV point cloud data of the mine slope is obtained, the local geometric roughness of each data point in the UAV point cloud data is calculated, data points with local geometric roughness greater than the fracture threshold are obtained, and an initial candidate weathering point cloud set is constructed. At the same time, the blast hole connection vector and the step free surface directional vector of the mine slope blasting operation are obtained. Morphological feature analysis was performed on the initial candidate weathering point cloud set to identify the three-dimensional fracture trace network of the mine slope and obtain the spatial direction vector of each fracture trace in the three-dimensional fracture trace network. The three-dimensional fracture trace network is traversed through a sliding window of a preset size. The spatial angle between the spatial direction vector of all fracture traces in the current sliding window and the direction vector of the borehole connection line and the direction vector of the free surface of the step is calculated. Based on the spatial angle, the directional aggregation index of the fracture traces in the current sliding window is calculated. Determine whether the directional aggregation index is greater than the blasting control fracture threshold. If it is, determine that the data point in the current sliding window is a pseudo weathering point and remove the pseudo weathering point from the initial candidate weathering point cloud set to obtain the real weathering point cloud data set. Spatial clustering is performed on data points in the real weathering point cloud dataset to delineate the actual range of the mine weathering zone.

2. The method for identifying the range of weathered zones in mines based on UAV point cloud data according to claim 1, characterized in that, The specific process for obtaining UAV point cloud data of the mine slope after the blasting operation is as follows: A drone equipped with an airborne lidar device was used to scan the mine slope and obtain point cloud data of the mine slope containing only the exposed rock surface.

3. The method for identifying the range of weathered zones in mines based on UAV point cloud data according to claim 1, characterized in that, The specific process for calculating the local geometric roughness of each data point in the UAV point cloud data is as follows: For any data point in the UAV point cloud data, a local spatial sphere is constructed using the roughness analysis radius; Obtain all adjacent data points within the local spatial sphere, and use the least squares method to perform local plane fitting on all adjacent data points to obtain the best fitting plane; Calculate the orthogonal projection distance from each adjacent data point within the local spatial sphere to the best-fit plane; The standard deviation or root mean square error of all orthogonal projected distances are used as the local geometric roughness of the data points.

4. The method for identifying the range of weathered zones in mines based on UAV point cloud data according to claim 1, characterized in that, The specific process for constructing the initial candidate weathering point cloud set is as follows: The local geometric roughness of each data point in the UAV point cloud data is obtained, and data points with local geometric roughness less than or equal to the fracture threshold are identified as unweathered intact bedrock points and removed. Data points with local geometric roughness greater than the fracture threshold are extracted and summarized to construct an initial candidate weathering point cloud set.

5. The method for identifying the range of weathered zones in mines based on UAV point cloud data according to claim 1, characterized in that, The specific process for simultaneously obtaining the borehole alignment vector and the bench free surface alignment vector for mine slope blasting operations is as follows: The coordinates of the borehole openings in the same row are obtained using historical point cloud data before blasting. Based on principal component analysis, spatial straight line fitting is performed on the borehole coordinates of the same row of boreholes. The eigenvector corresponding to the largest eigenvalue is extracted and normalized. The normalized eigenvector is then used as the direction vector of the borehole connection. Based on UAV point cloud data, an edge detection algorithm is used to obtain the macroscopic contour trajectory of the top or bottom slope of the current step. Principal component analysis is used to perform spatial linear fitting on the macroscopic contour trajectory. The feature vector corresponding to the largest feature value is extracted and normalized. The normalized feature vector is then used as the free surface orientation vector of the step.

6. The method for identifying the range of weathered zones in mines based on UAV point cloud data according to claim 1, characterized in that, The specific process for identifying the three-dimensional fracture trace network of the mine slope is as follows: Obtain the neighboring data points of any data point in the initial candidate weathering point cloud set, and calculate the local principal curvature of the data point; Set a morphological curvature threshold, and obtain data points whose local principal curvature is greater than the morphological curvature threshold as crack feature points; A density-based 3D region growing algorithm is adopted. The connectivity search radius is set, and the fracture feature points with spatial distances smaller than the connectivity search radius are iteratively connected topologically to construct the initial fracture connectivity domain. For any data point within the initial fracture connectivity domain, a spatial shrinkage algorithm is used for iterative processing until the shrinkage reaches the maximum number of iterations to form a one-dimensional point skeleton. The nearest neighbor connection method is used to connect the one-dimensional point skeleton sequentially according to the spatial topology, refining the initial fracture connectivity domain into a three-dimensional spatial curve with a width of one pixel. All the refined three-dimensional spatial curves and their intersections together constitute a three-dimensional fracture trace network of the mine slope.

7. The method for identifying the range of weathered zones in mines based on UAV point cloud data according to claim 1, characterized in that, The specific process for obtaining the spatial direction vector of each fracture trace in the three-dimensional fracture trace network is as follows: Each discrete three-dimensional data point constituting a one-dimensional point skeleton is denoted as a skeleton data point. The topological connectivity of each skeleton data point is obtained. Skeleton data points with a topological connectivity greater than or equal to 3 are marked as intersections of three-dimensional space curves, and skeleton data points with a topological connectivity equal to 1 are marked as endpoints. Using the intersection of the three-dimensional spatial curves and the endpoints as breakpoints, the three-dimensional fracture trace network is topologically broken to obtain multiple initial fracture segments; Calculate the local curvature of each skeleton data point on each initial fracture segment. If the local curvature of a skeleton data point is greater than the data point mutation threshold, it is taken as the curvature mutation point. At the curvature mutation point, the initial fracture segment is broken by curvature to obtain multiple fracture traces. For each fracture trace, principal component analysis is used to perform spatial linear fitting on the three-dimensional spatial coordinates of all skeleton data points constituting the fracture trace. After extracting the eigenvector corresponding to the largest eigenvalue and normalizing it, the spatial direction vector of the fracture trace is obtained.

8. The method for identifying the range of weathered zones in mines based on UAV point cloud data according to claim 1, characterized in that, The specific process for calculating the spatial angle between the spatial direction vectors of all fracture traces within the current sliding window and the spatial orientation vectors of the borehole connection line and the free surface orientation vector of the step is as follows: Create a sliding window and set its preset size; According to the preset spatial step size, the sliding window is translated and traversed in the three-dimensional space where the three-dimensional fracture trace network is located; At each traversal position, extract all crack traces that fall completely or partially within the current sliding window; For each fracture trace within the current sliding window, the spatial direction vector of the fracture trace is extracted. Using the three-dimensional vector dot product formula and the inverse cosine function, the first spatial angle between the spatial direction vector of the fracture trace and the direction vector of the borehole connection line, and the second spatial angle between the spatial direction vector of the fracture trace and the direction vector of the free surface of the step are calculated respectively.

9. The method for identifying the range of weathered zones in mines based on UAV point cloud data according to claim 1, characterized in that, The specific process for obtaining the actual weathering point cloud dataset is as follows: For each crack trace within the current sliding window, obtain the first spatial angle and the second spatial angle corresponding to the crack trace; Compare the size of the first spatial angle and the second spatial angle and select the minimum angle value as the main control angle for blasting; Set a parallel tolerance threshold, and count the number of fracture traces within the current sliding window whose blasting control angle is less than or equal to the parallel tolerance threshold. These are recorded as the number of directional fracture traces. The ratio of the number of directional fracture traces to the total number of fracture traces within the current sliding window is used as the directional aggregation index of fracture traces within the current sliding window. Compare the directional aggregation index calculated within the current sliding window with the blasting control fracture threshold. If the directional aggregation index is greater than the blasting control rupture threshold, then the data point in the current sliding window is determined to be a pseudo-weathering point; False weathering points are removed from the initial candidate weathering point cloud set, and all remaining data points after the removal operation are obtained to form the real weathering point cloud data set.

10. The method for identifying the range of weathered zones in mines based on UAV point cloud data according to claim 1, characterized in that, The specific process for delineating the actual extent of the mine weathering zone is as follows: A density-based spatial clustering algorithm is used to perform three-dimensional spatial clustering on data points in the real weathering point cloud dataset. Data points that meet the conditions of clustering neighborhood radius and minimum number of points contained in the core point are aggregated into high-density weathering clusters, and isolated data points are filtered out. The three-dimensional spatial coordinates of all data points constituting the high-density weathering clusters are projected onto a two-dimensional projection surface parallel to the macroscopic slope surface of the mine slope. The concave hull extraction algorithm is used to extract the outer contour boundary lines of each high-density weathering cluster on the two-dimensional projection surface. The area enclosed by the outer contour boundary of all high-density weathering clusters is defined as the actual mining weathering zone.

Images