An algorithm for extracting structural surfaces of rock slopes

Abstract

This invention discloses an algorithm for extracting structural surfaces of rock slopes, aiming to solve the problems of poor targeting and insufficient accuracy of structural surface extraction in existing point cloud denoising techniques. The algorithm comprises five steps: 3D point cloud data acquisition, neighborhood point set search, 3D orthogonal anisotropic Gaussian filtering and smoothing, covariance matrix construction and eigenvalue decomposition, and structural surface attitude information solution. By setting independent smoothing parameters for the x, y, and z directions, directional denoising is achieved, balancing data smoothness and structural details. The normal vector is obtained through covariance decomposition, accurately calculating the dip and tilt angle of the structural surface, with the overall error controlled within 2%. The algorithm is adaptable to various slope types, including long and steep slopes, and can seamlessly integrate with existing surveying equipment. It features high automation, strong anti-interference capabilities, and can provide high-precision data support for mountain engineering surveys, geological disaster early warning, and support scheme design, with broad application prospects.

Description

Technical Field

[0001] This invention relates to rock slope engineering survey, specifically an algorithm for extracting structural surfaces of rock slopes. Background Technology

[0002] In recent years, the engineering layout of industries such as water conservancy, mining, and construction in mountainous areas has become increasingly dense, and the excavation of rock slopes has become an indispensable key link in such projects.

[0003] The stability of rock slopes directly affects the safety, progress, and investment efficiency of engineering projects, and is one of the core factors restricting the sustainable development of engineering in mountainous areas. Extensive engineering practice shows that the vast majority of rock slope instability and failure are closely related to the geometric characteristics of the rock mass's structural surfaces. Their dip direction, angle, and other geometric parameters directly determine the slope's anti-sliding stability, failure mode, and instability risk level. Therefore, accurately obtaining the attitude information (dip direction, angle) of the rock mass's structural surfaces is the prerequisite and foundation for conducting rock slope stability assessments, developing scientific support plans, and preventing geological disasters.

[0004] Currently, the techniques for acquiring geometric parameters of rock mass structural planes are mainly divided into two categories, and each type of technique has significant limitations:

[0005] Traditional manual surveying, represented by the compass method, requires workers to be physically present on the slope, closely contacting the rock surface and directly reading the dip and inclination of the structural surfaces using a compass. Its drawbacks are significant:

[0006] 1. High safety risks: For steep slopes, fractured slopes or geological disaster hazard points, manual close-range measurement is very likely to cause safety accidents such as falls and rock collapses, which pose a serious threat to the life safety of measurement personnel;

[0007] 2. Low measurement efficiency: Manual measurement requires point-by-point and surface-by-surface measurement. For large-area and complex rock slopes, the measurement cycle is long and the labor intensity is high, which makes it difficult to meet the efficiency requirements of large-scale engineering surveys.

[0008] 3. Poor measurement accuracy: Measurement results are significantly affected by factors such as human operating experience, visual judgment errors, and environmental interference (such as wind and light), resulting in large data dispersion and making it impossible to guarantee the consistency of measurement accuracy. It is especially difficult to accurately measure hidden structural surfaces.

[0009] With the development of modern surveying technologies such as 3D laser scanning (TLS) and unmanned aerial vehicle (UAV) photogrammetry (UAV-SfM), intelligent extraction technology based on high-density point cloud data has gradually become mainstream. This technology acquires 3D point cloud datasets of rock slopes using laser scanners or UAVs, and then uses intelligent post-processing algorithms to extract structural surface information. Its core advantages lie in non-contact measurement, high data acquisition efficiency, and wide coverage. However, in practical engineering applications, this technology still faces significant technical bottlenecks:

[0010] 1. Severe noise interference in point cloud data: In actual engineering scenarios, the surface of rock masses naturally has uneven roughness. At the same time, the point cloud data acquisition process is inevitably affected by multiple factors such as instrument measurement errors (such as laser ranging accuracy and scanner angular resolution), changes in light intensity (such as strong light reflection and shadow occlusion), and environmental interference (such as airflow disturbance and dust influence). This causes point cloud data that should be on the same plane to exhibit irregular "jumping" phenomena in three-dimensional space, that is, there are a lot of random noise signals in the point cloud data.

[0011] 2. Existing filtering algorithms have poor noise reduction performance: Most existing intelligent post-processing algorithms use traditional homogenization filtering techniques (such as ordinary Gaussian filtering and median filtering) for noise processing. The essence of these filtering algorithms is to smooth the point cloud data in three-dimensional space as a whole, failing to distinguish the differences in noise distribution in the three orthogonal directions of x, y, and z, making it difficult to specifically weaken noise signals in different directions. For example, the surface roughness of rock mass mainly affects the point cloud data in the z-direction (elevation direction), while instrument measurement errors may produce more significant interference in the x and y directions (planar directions). Traditional filtering algorithms cannot adjust the smoothing intensity according to the noise characteristics in different directions, resulting in limited noise reduction effects—either the noise is not sufficiently removed, or excessive smoothing leads to the loss of structural surface details.

[0012] 3. Insufficient accuracy in structural surface extraction: Due to the failure to effectively weaken noise signals in the point cloud data, subsequent structural surface fitting and attitude calculation are severely interfered with. Existing algorithms construct planar models based on noisy point cloud data that deviate from the actual structural surfaces, resulting in large errors in the calculation of dip and tilt angles (typically above 5%-8%). This fails to meet the requirements of high-precision engineering surveys, thereby affecting engineers' judgment of slope stability and potentially leading to unreasonable support scheme designs and hidden safety hazards in the project.

[0013] In summary, existing technologies cannot simultaneously meet the requirements of "safety, efficiency, and high accuracy" for extracting structural surfaces of rock slopes. In particular, intelligent extraction technology based on 3D point clouds suffers from a core bottleneck: the lack of a filtering algorithm that can specifically weaken noise in the three-dimensional orthogonal directions while balancing data smoothness and preservation of structural surface details, resulting in insufficient accuracy in structural surface extraction. Therefore, developing an algorithm that can accurately reduce noise and significantly improve the accuracy of structural surface extraction has become a critical technical problem urgently needing to be solved in the field of rock slope engineering, and is of significant practical importance for ensuring the safety of mountainous engineering projects and promoting technological progress in the industry. Summary of the Invention

[0014] The purpose of this invention is to provide an algorithm for extracting structural surfaces of rock slopes to solve the problems mentioned in the background art.

[0015] To achieve the above objectives, the present invention provides the following technical solution:

[0016] An algorithm for extracting structural surfaces of rock slopes includes the following steps:

[0017] Step 1: 3D point cloud data acquisition

[0018] A laser scanner was used to acquire 3D point cloud data of a designated area of ​​a rock slope, obtaining the original point cloud dataset of the target rock slope. ,in The first point cloud The three-dimensional coordinates of the points This represents the total number of point clouds in the original point cloud dataset.

[0019] Furthermore, when the target rock slope is elongated (length to width ratio greater than 3:1), the limited single data acquisition window of the laser scanner (typically 30-100m) prevents it from covering the entire slope area at once. Therefore, image control points (ADCs) need to be set during data acquisition. ADCs are calibrated using GPSRTK positioning technology, with a planar positioning accuracy of no less than ±5mm and an elevation positioning accuracy of no less than ±10mm. The spacing between ADCs is determined based on the slope length, generally 30-50m, and they must be evenly distributed at both ends and key locations in the middle of the slope. This ensures that the original point cloud data collected from multiple segments can be accurately stitched together using point cloud stitching software (such as CloudCompare or Cyclone) to form a complete point cloud dataset of the target area.

[0020] Furthermore, the selection of a laser scanner should be determined based on the scale of the slope and the required measurement accuracy. Priority should be given to equipment with a point cloud density of not less than 50 points / cm², a ranging accuracy of not less than ±2mm, and an angular resolution of not less than 0.001° to ensure the quality of the original point cloud data.

[0021] Step 2: Neighborhood point set search

[0022] For each point cloud in the original point cloud dataset Traverse and search its neighborhood point set to construct a neighborhood point cloud dataset. Specifically, this is achieved through the following formula:

[0023]

[0024]

[0025] in, Point cloud The neighborhood point set, Point cloud With point clouds The Euclidean distance between them is calculated using the following formula: ; A function that sorts point cloud distances in descending order; The constant represents the point cloud. The number of neighboring point clouds ranges from 15 to 30 (too small a value will result in insufficient representativeness of the neighboring point clouds, while too large a value will introduce point clouds that are not on the same structural surface, affecting the accuracy of subsequent processing). This represents the total number of point clouds in the original point cloud dataset.

[0026] The above formula can be used for each target point cloud. Filter for the closest These point clouds form the neighborhood point set of the point cloud, providing basic data for subsequent filtering and smoothing processes.

[0027] Step 3: Three-dimensional orthogonal anisotropic Gaussian filtering smoothing process

[0028] A three-dimensional orthogonal anisotropic Gaussian filtering smoothing algorithm is used to smooth the neighborhood point cloud dataset. The processing is performed to specifically weaken the noise signals on the three orthogonal components of x, y, and z. The smoothed neighborhood point cloud dataset is obtained. The specific processing procedure is as follows:

[0029] 3.1 Calculate the Gaussian weights for the three orthogonal directions:

[0030] For neighborhood point sets Each point cloud Calculate its relationship with the target point cloud Coordinate differences in the three orthogonal directions x, y, and z , , Then, calculate the Gaussian weights for the three directions using the following formula:

[0031]

[0032]

[0033]

[0034] in, , , These are the smoothing control parameters for the three orthogonal directions (x, y, and z), with values ​​ranging from 0.5 to 1.5. Larger parameter values ​​result in a more significant smoothing effect during Gaussian filtering; conversely, smaller parameter values ​​preserve more detail in the point cloud data. In practical applications, the parameters can be adjusted based on noise distribution characteristics: if the rock surface roughness is large (significant noise in the z-direction), the parameters can be adjusted accordingly. The value should be set between 1.0 and 1.5; if the instrument measurement error results in significant noise in the x and y directions, the value can be increased appropriately. The value of .

[0035] 3.2 Gaussian weight normalization processing:

[0036] Since the absolute values ​​of the Gaussian weights in the three orthogonal directions may differ, directly using them for point cloud updates will lead to an imbalance in the smoothing intensity of each direction. Therefore, the calculated Gaussian weights need to be adjusted. , , Normalization is performed separately to ensure that the sum of the weights in each direction is 1. The normalization formula is:

[0037]

[0038]

[0039]

[0040] in, , , These are the normalized Gaussian weights for the x, y, and z directions, respectively. Through normalization, it can be ensured that the smooth adjustment in each direction is based on the same weight benchmark, avoiding filtering distortion caused by weight differences.

[0041] 3.3 Neighborhood point cloud dataset update:

[0042] Based on the normalized Gaussian weights, the neighborhood point cloud dataset Each point cloud coordinate in the dataset is updated with weights to obtain a smoothed neighborhood point cloud dataset. The updated formula is:

[0043]

[0044]

[0045]

[0046] in, For smoothed point clouds coordinates For neighborhood point set Midpoint cloud The original coordinates are obtained. Through the above weighted update, the noise signals in the three orthogonal directions can be effectively weakened while preserving the overall shape of the structural surface, making the point cloud data distribution more uniform and consistent.

[0047] Step 4: Covariance Matrix Construction and Eigenvalue Decomposition

[0048] For the smoothed neighborhood point cloud dataset The covariance matrix is ​​constructed, and the normal vector of the dataset is obtained through eigenvalue decomposition. The specific process is as follows:

[0049] 4.1 Calculate the coordinates of the center point of the smoothed neighborhood point cloud dataset:

[0050] First calculate The coordinates of the center point are obtained by averaging the coordinates of all point clouds in the x, y, and z directions. The calculation formula is:

[0051]

[0052]

[0053]

[0054] in, For smoothed neighborhood point clouds coordinates The number of neighboring point clouds (as in step 2) (The values ​​are consistent). The coordinates of the center point are the basis for constructing the covariance matrix, which can reflect the overall distribution center of the neighborhood point cloud.

[0055] 4.2 Constructing the covariance matrix:

[0056] Based on the center point coordinates Calculate the covariance of the neighborhood point cloud dataset in the x, y, and z directions, and construct a 3×3 covariance matrix. The formula is:

[0057]

[0058] Among them, covariance The calculation formula is:

[0059]

[0060] In the formula, , These are the coordinates of the smoothed neighborhood point cloud in the corresponding directions. These are the coordinates of the center point in the corresponding direction. The covariance matrix reflects the dispersion of the neighborhood point cloud in three-dimensional space and the correlation in each direction, and is the core input for subsequent eigenvalue decomposition.

[0061] 4.3 Eigenvalue decomposition of the covariance matrix:

[0062] The Jacobi iteration method is used to analyze the covariance matrix. Eigenvalue decomposition yields three eigenvalues. and the corresponding feature vectors The magnitude of the eigenvalue reflects the degree of dispersion of the neighborhood point cloud along the corresponding eigenvector direction: the largest eigenvalue... corresponding feature vector The principal direction of the point cloud distribution, the minimum eigenvalue corresponding feature vector The plane perpendicular to the fitted neighborhood point cloud, i.e., the normal vector of the smoothed neighborhood point cloud dataset.

[0063] The eigenvector corresponding to the smallest eigenvalue is chosen as the normal vector because this vector can most accurately reflect the normal direction of the structural surface, providing a key basis for solving the subsequent structural surface attitude information.

[0064] Step 5: Solving for structural surface attitude information

[0065] Based on the normal vectors of the smoothed neighborhood point cloud dataset Using a three-dimensional coordinate system (X-axis pointing east, Y-axis pointing north, Z-axis pointing sky), the dip direction of the structural plane of the rock slope is calculated using the following formula. and tilt angle Complete the structural surface extraction:

[0066] 5.1 Propensity Calculation:

[0067] The dip direction refers to the angle between the projection line of the structural surface onto the horizontal plane and the due north direction (clockwise is positive). The calculation formula is:

[0068]

[0069] If the calculation yields If the value is negative, it needs to be added to convert 360° to an angle value within the range of 0°-360° to ensure that the inclination representation conforms to engineering specifications.

[0070] 5.2 Inclination Calculation:

[0071] The tilt angle refers to the angle between the structural surface and the horizontal plane (ranging from 0° to 90°), and the calculation formula is:

[0072]

[0073] in, Normal vector The coordinate components, and the dip and inclination angles calculated using the above formulas, can accurately reflect the true geometric characteristics of the rock slope structure.

[0074] Compared with the prior art, the beneficial effects of the present invention are:

[0075] This invention innovatively proposes a three-dimensional orthogonal anisotropic Gaussian filtering smoothing algorithm, overcoming the technical limitation of traditional homogenization filtering algorithms that cannot distinguish noise in three-dimensional directions. By setting independent smoothing control parameters for the three orthogonal directions (x, y, and z), the smoothing intensity can be adjusted specifically according to the noise characteristics in different directions (such as rock roughness noise in the z-direction and instrument measurement noise in the x / y directions), achieving precise weakening of noise in each direction. Practical engineering verification shows that this filtering algorithm can reduce the "jump" amplitude of point cloud data by more than 80%, significantly improving the smoothness and consistency of the point cloud data, laying a high-quality data foundation for subsequent structural surface parameter calculation.

[0076] This invention employs a fully optimized design encompassing "directional filtering - covariance decomposition - normal vector solution - attitude calculation," ensuring precision control at every stage from data source to result output.

[0077] Directional filtering ensures that point cloud data is "de-noiseed without distortion" and preserves the true shape of structural surfaces;

[0078] The covariance matrix construction and eigenvalue decomposition can accurately capture the distribution pattern of the neighborhood point cloud, and the normal vector solution has high accuracy;

[0079] The attitude calculation is based on standardized coordinate system definitions and mathematical formulas, avoiding human error.

[0080] Comparative experiments have verified that the structural surface inclination error extracted by the algorithm of this invention is ≤0.5%, the tilt angle error is ≤1.5%, and the overall error is controlled within 2%, which is far better than the 5%-8% error level of existing intelligent extraction algorithms. It can provide engineers with accurate structural surface parameters, effectively avoid misjudgment of slope stability and unreasonable support scheme design caused by data errors, and significantly reduce engineering safety risks.

[0081] This invention fully considers the diverse morphologies of rock slopes in actual engineering projects. By setting image control points during the data acquisition phase, it effectively solves the technical challenge of stitching together point clouds for long, narrow slopes. In terms of algorithm parameter settings, smoothing control parameters are implemented. and the number of neighboring points All algorithms can be flexibly adjusted according to slope type (such as long, steep, or fractured) and noise characteristics to adapt to different terrain conditions and engineering needs. Furthermore, the algorithms do not rely on special hardware and can seamlessly integrate with existing mainstream 3D laser scanners, UAV photogrammetry systems, and point cloud processing software. This eliminates the need for modifications or upgrades to existing equipment, lowering the barrier to entry and cost for engineering applications, and thus possessing broad application value.

[0082] This invention achieves semi-automated processing of the entire process from 3D point cloud data acquisition to structural surface attitude information extraction:

[0083] During the data acquisition phase, laser scanners / drones can quickly acquire point cloud data of large-area slopes, improving acquisition efficiency by more than 50 times compared to traditional manual measurement.

[0084] In the data processing stage, neighborhood point set search, directional filtering, covariance decomposition and other steps can all be automated through programming without manual intervention. The extraction time for a single structural surface is only a few seconds, which greatly improves work efficiency compared to the traditional manual measurement of "point by point and surface by surface".

[0085] In addition, the standardized process of the algorithm avoids subjective errors caused by manual operation, and the data results are highly reliable and consistent, which can meet the needs of large-scale and high-intensity engineering surveys, shorten project timelines, and reduce labor costs.

[0086] The algorithm of this invention has strong adaptability to complex environments and can effectively resist the influence of various factors such as rock surface roughness, instrument measurement error, illumination changes, and environmental interference.

[0087] • Directional filtering algorithms can adapt to noise of different intensities and directions by adjusting weights;

[0088] • Neighborhood point set search uses Euclidean distance sorting and filtering to ensure the representativeness of the neighborhood point cloud and reduce interference from point clouds with different structural surfaces;

[0089] The covariance matrix eigenvalue decomposition uses the mature Jacobi iteration method, which has high computational stability and is not easily affected by the discreteness of the data.

[0090] Applications in multiple complex engineering scenarios have shown that even under conditions of fractured rock mass, intense sunlight, and harsh measurement environments, the algorithm of this invention can still stably output high-precision structural surface parameters, demonstrating significantly better robustness than existing algorithms. Attached Figure Description

[0091] Figure 1 This is a flowchart of the algorithm for extracting structural surfaces of rock slopes;

[0092] It demonstrates the five core steps from 3D point cloud data acquisition, neighborhood point set search, 3D orthogonal anisotropic Gaussian filtering smoothing, covariance matrix construction and eigenvalue decomposition to solving for structural surface attitude information, as well as the logical connections between each step.

[0093] Figure 2 This is a schematic diagram of smoothing the graph neighborhood point cloud dataset;

[0094] The left side shows the original 3D point cloud data of the rock slope (obvious "jumping" phenomenon is visible, and the noise signal is significant), while the right side shows the neighborhood point cloud dataset after being processed by the 3D orthogonal anisotropic Gaussian filtering smoothing algorithm (the point cloud distribution is uniform and consistent, and the noise signal is effectively weakened), which intuitively demonstrates the noise reduction effect of the filtering algorithm.

[0095] Figure 3 This is a schematic diagram of solving the attitude information of structural surfaces based on a smoothed neighborhood point cloud dataset.

[0096] The coordinate system is defined as follows: the X-axis points east, the Y-axis points north, and the Z-axis points to the sky. The normal vector of the smoothed neighborhood point cloud dataset. The slope of the structural plane (the angle between the horizontal projection line and the due north direction). The dip angle of the structural plane (the angle between the structural plane and the horizontal plane) clearly shows the correspondence between the normal vector and the attitude parameters of the structural plane. Detailed Implementation

[0097] The technical solution of the present invention will be further described in detail below with reference to specific embodiments.

[0098] An algorithm for extracting structural surfaces of rock slopes includes the following steps:

[0099] Step 1: 3D point cloud data acquisition

[0100] A laser scanner was used to acquire 3D point cloud data of a designated area of ​​a rock slope, obtaining the original point cloud dataset of the target rock slope. ,in The first point cloud The three-dimensional coordinates of the points This represents the total number of point clouds in the original point cloud dataset.

[0101] Furthermore, when the target rock slope is elongated (length to width ratio greater than 3:1), the limited single data acquisition window of the laser scanner (typically 30-100m) prevents it from covering the entire slope area at once. Therefore, image control points (ADCs) need to be set during data acquisition. ADCs are calibrated using GPSRTK positioning technology, with a planar positioning accuracy of no less than ±5mm and an elevation positioning accuracy of no less than ±10mm. The spacing between ADCs is determined based on the slope length, generally 30-50m, and they must be evenly distributed at both ends and key locations in the middle of the slope. This ensures that the original point cloud data collected from multiple segments can be accurately stitched together using point cloud stitching software (such as CloudCompare or Cyclone) to form a complete point cloud dataset of the target area.

[0102] Furthermore, the selection of a laser scanner should be determined based on the scale of the slope and the required measurement accuracy. Priority should be given to equipment with a point cloud density of not less than 50 points / cm², a ranging accuracy of not less than ±2mm, and an angular resolution of not less than 0.001° to ensure the quality of the original point cloud data.

[0103] Step 2: Neighborhood point set search

[0104] For each point cloud in the original point cloud dataset Traverse and search its neighborhood point set to construct a neighborhood point cloud dataset. Specifically, this is achieved through the following formula:

[0105]

[0106]

[0107] in, Point cloud The neighborhood point set, Point cloud With point clouds The Euclidean distance between them is calculated using the following formula: ; A function that sorts point cloud distances in descending order; The constant represents the point cloud. The number of neighboring point clouds ranges from 15 to 30 (too small a value will result in insufficient representativeness of the neighboring point clouds, while too large a value will introduce point clouds that are not on the same structural surface, affecting the accuracy of subsequent processing). This represents the total number of point clouds in the original point cloud dataset.

[0108] The above formula can be used for each target point cloud. Filter for the closest These point clouds form the neighborhood point set of the point cloud, providing basic data for subsequent filtering and smoothing processes.

[0109] Step 3: Three-dimensional orthogonal anisotropic Gaussian filtering smoothing process

[0110] A three-dimensional orthogonal anisotropic Gaussian filtering smoothing algorithm is used to smooth the neighborhood point cloud dataset. The processing is performed to specifically weaken the noise signals on the three orthogonal components of x, y, and z. The smoothed neighborhood point cloud dataset is obtained. The specific processing procedure is as follows:

[0111] 3.1 Calculate the Gaussian weights for the three orthogonal directions:

[0112] For neighborhood point sets Each point cloud Calculate its relationship with the target point cloud Coordinate differences in the three orthogonal directions x, y, and z , , Then, calculate the Gaussian weights for the three directions using the following formula:

[0113]

[0114]

[0115]

[0116] in, , , These are the smoothing control parameters for the three orthogonal directions (x, y, and z), with values ​​ranging from 0.5 to 1.5. Larger parameter values ​​result in a more significant smoothing effect during Gaussian filtering; conversely, smaller parameter values ​​preserve more detail in the point cloud data. In practical applications, the parameters can be adjusted based on noise distribution characteristics: if the rock surface roughness is large (significant noise in the z-direction), the parameters can be adjusted accordingly. The value should be set between 1.0 and 1.5; if the instrument measurement error results in significant noise in the x and y directions, the value can be increased appropriately. The value of .

[0117] 3.2 Gaussian weight normalization processing:

[0118] Since the absolute values ​​of the Gaussian weights in the three orthogonal directions may differ, directly using them for point cloud updates will lead to an imbalance in the smoothing intensity of each direction. Therefore, the calculated Gaussian weights need to be adjusted. , , Normalization is performed separately to ensure that the sum of the weights in each direction is 1. The normalization formula is:

[0119]

[0120]

[0121]

[0122] in, , , These are the normalized Gaussian weights for the x, y, and z directions, respectively. Through normalization, it can be ensured that the smooth adjustment in each direction is based on the same weight benchmark, avoiding filtering distortion caused by weight differences.

[0123] 3.3 Neighborhood point cloud dataset update:

[0124] Based on the normalized Gaussian weights, the neighborhood point cloud dataset Each point cloud coordinate in the dataset is updated with weights to obtain a smoothed neighborhood point cloud dataset. The updated formula is:

[0125]

[0126]

[0127]

[0128] in, For smoothed point clouds coordinates For neighborhood point set Midpoint cloud The original coordinates are obtained. Through the above weighted update, the noise signals in the three orthogonal directions can be effectively weakened while preserving the overall shape of the structural surface, making the point cloud data distribution more uniform and consistent.

[0129] Step 4: Covariance Matrix Construction and Eigenvalue Decomposition

[0130] For the smoothed neighborhood point cloud dataset The covariance matrix is ​​constructed, and the normal vector of the dataset is obtained through eigenvalue decomposition. The specific process is as follows:

[0131] 4.1 Calculate the coordinates of the center point of the smoothed neighborhood point cloud dataset:

[0132] First calculate The coordinates of the center point are obtained by averaging the coordinates of all point clouds in the x, y, and z directions. The calculation formula is:

[0133]

[0134]

[0135]

[0136] in, For smoothed neighborhood point clouds coordinates The number of neighboring point clouds (as in step 2) (The values ​​are consistent). The coordinates of the center point are the basis for constructing the covariance matrix, which can reflect the overall distribution center of the neighborhood point cloud.

[0137] 4.2 Constructing the covariance matrix:

[0138] Based on the center point coordinates Calculate the covariance of the neighborhood point cloud dataset in the x, y, and z directions, and construct a 3×3 covariance matrix. The formula is:

[0139]

[0140] Among them, covariance The calculation formula is:

[0141]

[0142] In the formula, , These are the coordinates of the smoothed neighborhood point cloud in the corresponding directions. These are the coordinates of the center point in the corresponding direction. The covariance matrix reflects the dispersion of the neighborhood point cloud in three-dimensional space and the correlation in each direction, and is the core input for subsequent eigenvalue decomposition.

[0143] 4.3 Eigenvalue decomposition of the covariance matrix:

[0144] The Jacobi iteration method is used to analyze the covariance matrix. Eigenvalue decomposition yields three eigenvalues. and the corresponding feature vectors The magnitude of the eigenvalue reflects the degree of dispersion of the neighborhood point cloud along the corresponding eigenvector direction: the largest eigenvalue... corresponding feature vector The principal direction of the point cloud distribution, the minimum eigenvalue corresponding feature vector The plane perpendicular to the fitted neighborhood point cloud, i.e., the normal vector of the smoothed neighborhood point cloud dataset.

[0145] The eigenvector corresponding to the smallest eigenvalue is chosen as the normal vector because this vector can most accurately reflect the normal direction of the structural surface, providing a key basis for solving the subsequent structural surface attitude information.

[0146] Step 5: Solving for structural surface attitude information

[0147] Based on the normal vectors of the smoothed neighborhood point cloud dataset Using a three-dimensional coordinate system (X-axis pointing east, Y-axis pointing north, Z-axis pointing sky), the dip direction of the structural plane of the rock slope is calculated using the following formula. and tilt angle Complete the structural surface extraction:

[0148] 5.1 Propensity Calculation:

[0149] The dip direction refers to the angle between the projection line of the structural surface onto the horizontal plane and the due north direction (clockwise is positive). The calculation formula is:

[0150]

[0151] If the calculation yields If the value is negative, it needs to be added to convert 360° to an angle value within the range of 0°-360° to ensure that the inclination representation conforms to engineering specifications.

[0152] 5.2 Inclination Calculation:

[0153] The tilt angle refers to the angle between the structural surface and the horizontal plane (ranging from 0° to 90°), and the calculation formula is:

[0154]

[0155] in, Normal vector The coordinate components, and the dip and inclination angles calculated using the above formulas, can accurately reflect the true geometric characteristics of the rock slope structure.

[0156] Example 1

[0157] A long, narrow rock slope, 600m long and 100m high, is part of a highway construction project in a mountainous area. The slope runs north-south and exhibits a typical elongated shape (length to width ratio approximately 6:1). The rock mass is granite with a relatively rough surface and localized fracture zones. It is necessary to extract the orientation information of the main structural planes of the slope to provide a basis for the design of the support scheme.

[0158] 1. 3D point cloud data acquisition:

[0159] A Leica ScanStation P50 3D laser scanner was selected, which boasts a point cloud density of 120 points / cm², a ranging accuracy of ±2mm, and an angular resolution of 0.001°, meeting the engineering accuracy requirements. Due to the elongated shape of the slope, a control point was established every 40m along the slope's direction, for a total of 15 control points. GPS RTK positioning technology was used for coordinate calibration (planar accuracy ±3mm, elevation accuracy ±8mm). Point cloud data was collected from five sections of the slope, each approximately 120m in length. After collection, the point clouds were stitched together using CloudCompare software to obtain a complete original point cloud dataset. Total number of point clouds .

[0160] 2. Neighborhood point set search:

[0161] Set the number of neighboring point clouds (Balancing representativeness and computational efficiency), for each point cloud in the original point cloud dataset The Euclidean distance between it and all other point clouds is calculated using formulas (1) and (2), and then... The function sorts the point cloud distances in descending order and selects the point cloud corresponding to the top 25 distances as the neighborhood point set. Construct a neighborhood point cloud dataset.

[0162] 3. Three-dimensional orthogonal anisotropic Gaussian filtering smoothing:

[0163] 3.1 Gaussian weight calculation: Considering the large surface roughness of the slope rock mass, a smoothing control parameter in the x-direction is set. y-direction smoothing control parameters z-direction smoothing control parameters (Noise is more significant in the z-direction, so increase the smoothing intensity). For each neighborhood point set... Calculate the coordinate difference between each point and the target point. Substituting into formulas (3)-(5), we obtain the Gaussian weights in the three directions. .

[0164] 3.2 Weight Normalization: The Gaussian weights are normalized using formulas (6)-(8) to ensure that the sum of the weights in each direction is 1, thus obtaining the normalized weights. .

[0165] 3.3 Point Cloud Update: According to formulas (9)-(11), the coordinates of the neighborhood point cloud are updated by weighting the normalized weights to obtain the smoothed neighborhood point cloud dataset. After processing, the "jumping" phenomenon in the point cloud data was basically eliminated, and the uniformity of point cloud distribution was improved by 85%.

[0166] 4. Covariance matrix construction and eigenvalue decomposition:

[0167] 4.1 Calculation of center point coordinates: The center point coordinates of each smoothed neighborhood point cloud dataset are calculated using formulas (12)-(14). .

[0168] 4.2 Covariance Matrix Construction: Based on formulas (15)-(16), calculate the covariance in each direction and construct the covariance matrix. For example, for a certain neighborhood point cloud dataset, the calculated covariance matrix is:

[0169]

[0170] 4.3 Eigenvalue Decomposition: The Jacobi iteration method is used to decompose the covariance matrix. Perform eigenvalue decomposition to obtain eigenvalues. , , The corresponding feature vector , , Select the smallest eigenvalue corresponding feature vector As a normal vector.

[0171] 5. Solving for the attitude information of structural surfaces:

[0172] normal vector Substituting into formulas (17) and (18), the dip direction of the structural surface is calculated. (Converted to 0°-360° range), tilt angle .

[0173] Verification results

[0174] The dip and inclination angle of the same structural surface of the slope were measured using a traditional manual compass, yielding a dip of 325.1° and an inclination angle of 47.9°. The dip error of the algorithm in this invention is... The tilt angle error is The overall error is controlled within 2%, meeting the requirements of high-precision engineering surveys. Meanwhile, the overall data processing cycle of this algorithm is 2.5 hours, significantly improving work efficiency compared to traditional manual measurement (which takes one day to measure a single structural surface).

[0175] Example 2

[0176] A steep rock slope in a certain water conservancy project has a height of 150m and a slope of 75°. The rock mass is sandstone, and there are many cracks on the slope surface. There is a risk of small-scale collapse in some areas. It is necessary to quickly and accurately extract the structural surface attitude information to provide data support for geological disaster early warning.

[0177] Implementation steps

[0178] 1. 3D point cloud data acquisition:

[0179] The RIEGL VZ-6000 3D laser scanner was selected. This device has long-distance measurement capabilities (maximum measurement distance 3000m), a point cloud density of 100 points / cm², and a ranging accuracy of ±1mm, making it suitable for non-contact measurement of steep slopes. Due to the steepness and height of the slope, segmented data acquisition was unnecessary; a single acquisition could cover the entire slope area, therefore no image control points were set. The resulting raw point cloud dataset was obtained. Total number of point clouds .

[0180] 2. Neighborhood point set search:

[0181] Considering the presence of numerous cracks on the slope, the point cloud data may exhibit local discreteness. Therefore, a set number of neighboring point clouds is determined. (To reduce interference from point clouds with different structural surfaces), the neighborhood point cloud dataset is constructed by searching the neighborhood point set of each point cloud using formulas (1) and (2). .

[0182] 3. Three-dimensional orthogonal anisotropic Gaussian filtering smoothing:

[0183] 3.1 Gaussian Weight Calculation: Slope cracks mainly cause point cloud noise in the x and y directions (planar direction), while the noise in the z direction is relatively small. Therefore, smoothing control parameters are set. , , Calculate the coordinate difference between each neighboring point and the target point, and substitute it into formulas (3)-(5) to obtain the Gaussian weights. .

[0184] 3.2 Weight normalization: The Gaussian weights are normalized using formulas (6)-(8).

[0185] 3.3 Point Cloud Update: Update the coordinates of the neighborhood point cloud according to formulas (9)-(11) to obtain the smoothed neighborhood point cloud dataset. The point cloud discrete noise caused by the cracks is effectively weakened.

[0186] 4. Covariance matrix construction and eigenvalue decomposition:

[0187] 4.1 Calculate the coordinates of the center point of the smoothed neighborhood point cloud .

[0188] 4.2 Constructing the covariance matrix The covariance matrix of a certain neighborhood point cloud is shown below:

[0189]

[0190] 4.3 Eigenvalue decomposition yields eigenvalues , , The corresponding feature vector , as the normal vector.

[0191] 5. Solving for the attitude information of structural surfaces:

[0192] Substituting into formulas (17) and (18), the dip direction of the structural surface is calculated. ,inclination .

[0193] Verification results

[0194] Verification was conducted using a combination of drone aerial photography and manual verification. The measurement results were a tilt of 327.9° and a tilt angle of 44.9°. The tilt error of the algorithm in this invention is... The tilt angle error is The error is significantly lower than that of existing algorithms, fully demonstrating the high-precision advantage of this invention in steep slope scenarios. Furthermore, the algorithm's data processing cycle is only 1.8 hours, rapidly extracting structural surfaces in high-risk areas and saving valuable time for geological disaster early warning.

[0195] Summary of Implementation Examples

[0196] The two embodiments described above, targeting two typical rock slope types—long and steep—respectively demonstrate the practicality and high accuracy of the algorithm of this invention. This is further enhanced by flexibly adjusting the number of neighboring points. and smooth control parameters This algorithm can adapt to different types of rock slopes, and the extraction accuracy is controlled within 2%, significantly improving work efficiency compared to traditional techniques. Furthermore, the algorithm does not rely on special hardware and can be seamlessly integrated with existing surveying technologies, demonstrating broad engineering application prospects.

[0197] The above embodiments are only used to illustrate the technical solutions of the present invention and are not intended to limit the present invention. Although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that modifications can still be made to the technical solutions described in the foregoing embodiments, or equivalent substitutions can be made to some of the technical features; and these modifications or substitutions can make various changes without departing from the spirit of the present invention.

Claims

1. An algorithm for extracting structural surfaces of rock slopes, characterized in that, Includes the following steps: Step 1: 3D point cloud data acquisition: Use a laser scanner to acquire 3D point cloud data of a designated area of ​​the rock slope to obtain the original point cloud dataset; Step 2: Neighborhood point set search: For each point cloud in the original point cloud dataset, traverse and search its neighborhood point set to construct a neighborhood point cloud dataset; Step 3: 3D Orthogonal Anisotropic Gaussian Filter Smoothing Processing: The neighborhood point cloud dataset is processed using a 3D orthogonal anisotropic Gaussian filter smoothing algorithm. The 3D orthogonal anisotropic Gaussian filter smoothing algorithm refers to: setting independent smoothing control parameters for the three orthogonal directions x, y, and z; calculating the Gaussian weights for the three orthogonal directions in sequence; normalizing the Gaussian weights; updating the neighborhood point cloud coordinates based on the normalized Gaussian weights; weakening the noise signals on the three orthogonal components x, y, and z; and obtaining the smoothed neighborhood point cloud dataset. Step 4: Covariance Matrix Construction and Eigenvalue Decomposition: Construct a covariance matrix for the smoothed neighborhood point cloud dataset, and perform eigenvalue decomposition on the covariance matrix to obtain the normal vector of the smoothed neighborhood point cloud dataset. Step 5: Solving for structural surface attitude information: Based on the normal vector, calculate the dip and dip angle of the rock slope structural surface to complete the structural surface extraction.

2. The algorithm for extracting structural surfaces of rock slopes according to claim 1, characterized in that, In step 1, when the target rock slope is distributed in a long strip shape, image control points are set during the data acquisition process. The image control points are calibrated using GPS RTK positioning technology, with a planar positioning accuracy of not less than ±5mm and an elevation positioning accuracy of not less than ±10mm.

3. The algorithm for extracting structural surfaces of rock slopes according to claim 1, characterized in that, In step 2, the search for the neighborhood point set is achieved using the following formula: in, For point clouds The neighborhood point set, For point clouds With point clouds Euclidean distance, A function to sort point cloud distances in descending order. The constant representing the number of neighboring point clouds, with a value ranging from 15 to 30. This represents the total number of points in the original point cloud dataset.

4. The algorithm for extracting structural surfaces of rock slopes according to claim 1, characterized in that, In step 3, the formula for calculating the Gaussian weights in the three orthogonal directions using the three-dimensional orthogonal anisotropic Gaussian filtering smoothing algorithm is as follows: in, , , These are the smoothing control parameters in the x, y, and z directions, with values ​​ranging from 0.5 to 1.

5. , , This represents the coordinate difference between the neighboring point and the target point in the corresponding direction.

5. The algorithm for extracting structural surfaces of rock slopes according to claim 3, characterized in that, In step 3, the Gaussian weights are calculated. , , Then, they are normalized separately. The normalization formula is as follows: in, , , These are the normalized Gaussian weights.

6. The algorithm for extracting structural surfaces of rock slopes according to claim 5, characterized in that, In step 3, the update formula for the neighborhood point cloud dataset is: in, The coordinates of the smoothed point cloud, The coordinates of the point cloud in the neighborhood point set.

7. The algorithm for extracting structural surfaces of rock slopes according to claim 1, characterized in that, In step 4, the formula for constructing the covariance matrix is: in, , , The coordinates of the center point of the smoothed neighborhood point cloud dataset. This is the smoothed neighborhood point cloud.

8. The algorithm for extracting structural surfaces of rock slopes according to claim 1, characterized in that, In step 4, the Jacobi iteration method is used to decompose the covariance matrix into eigenvalues, and the eigenvector corresponding to the smallest eigenvalue is selected as the normal vector of the smoothed neighborhood point cloud dataset.

9. The algorithm for extracting structural surfaces of rock slopes according to claim 1, characterized in that, In step 5, the formula for calculating the structural plane dip is: in, Let be the coordinate components of the normal vector, if If the value is negative, 360° needs to be added to convert it to an angle value within the range of 0°-360°.

10. The algorithm for extracting structural surfaces of rock slopes according to claim 1, characterized in that, In step 5, the formula for calculating the inclination angle of the structural surface is: in, These are the coordinate components of the normal vector. The value range is 0°-90°.

Images