A method for identifying and calculating the volume of a rock joint shear failure zone

Abstract

The application discloses a rock joint surface shear failure zone identification and volume calculation method, point cloud data of rock joint surfaces before and after shearing is collected, and repeated points with the same coordinates are deleted; the center of the point cloud data is calculated and moved to the coordinate origin to obtain translated point cloud data; the point cloud data is rotated so that the normal direction is consistent with the z-axis direction according to the included angle of the normal vector and the z-axis of the three-dimensional coordinate system and the cross product result; after the point cloud data is projected on the xoy plane, the point cloud data is rotated according to the included angle of the boundary and the x-axis, and the boundary is adjusted to the coordinate axis; the two groups of point cloud data before and after shearing are accurately registered through a point cloud registration method; the shear failure zone is marked according to the height difference of the point cloud data, and the volume surrounded by the point cloud data before and after shearing is calculated, which is the shear failure volume. A series of optimization processes are performed on the point cloud data of the joint surfaces before and after shearing, and then accurate registration is performed, so that the accuracy and reliability of the shear failure zone identification and volume calculation are ensured, and the method is suitable for rock mass mechanical test result analysis.

Description

Technical Field

[0001] This application relates to the field of rock mechanics experimental technology, specifically to a method for identifying and calculating the volume of shear failure zones on rock joint surfaces. Background Technology

[0002] Due to long-term geological processes, numerous discontinuities such as joints and fissures inevitably exist in engineering rock masses. These discontinuities crisscross, significantly weakening the mechanical properties of intact rock. The surface morphology of joints is one of the important factors affecting their shear strength. Only by fully understanding and quantitatively characterizing the geometric features of joint surfaces can we accurately reveal their shear mechanical behavior. Therefore, accurate identification and volumetric calculation of shear failure zones at joint surfaces are helpful in better predicting the stability of engineering rock masses.

[0003] To address the aforementioned issues, several patents have proposed methods for calculating the shear failure area. For example, Chinese patent CN102749046B discloses a method for measuring the shear area in a direct shear test of a rock mass structural surface. Its main steps are: acquiring and preprocessing photographs of the rock mass structural surface; converting the structural surface photographs into black-and-white binary images; calculating the percentage of the area of ​​the white portion in the entire image; and calculating the actual total area of ​​the rock mass structural surface using a ruler used during photography as a reference. Another example is Chinese patent CN107144243A, which discloses a method and system for measuring the shear failure area of ​​a rock mass structural surface. This method first acquires point cloud data of the structural surface before and after shearing, then aligns the two and calculates the height deviation of the corresponding points. It then determines the number of point clouds with a height deviation greater than 0 and their percentage W of the total. Based on the percentage W and the area S of the rock mass structural surface before the direct shear test, the shear failure area of ​​the rock mass structural surface is calculated.

[0004] The two methods described above are simple to operate and have low cost, but neither can be directly used to calculate the shear failure volume of joint surfaces. Summary of the Invention

[0005] In order to solve the above-mentioned technical problems, this application proposes the following technical solution:

[0006] In a first aspect, embodiments of this application provide a method for identifying and calculating the volume of shear failure zones on rock joint surfaces, including:

[0007] Collect point cloud data of rock joint surfaces before and after shearing, and delete duplicate points with the same coordinates;

[0008] Calculate the point cloud data center and move it to the origin of the coordinate system to obtain the translated point cloud data;

[0009] Based on the angle between the normal vector and the z-axis of the three-dimensional coordinate system and the cross product, rotate the point cloud data so that its normal is aligned with the z-axis direction;

[0010] After projecting the point cloud data onto the xoy plane, rotate the point cloud data according to the angle between the boundary and the x-axis, and adjust the boundary to the coordinate axis;

[0011] The point cloud registration method is used to accurately register the two sets of point clouds before and after cutting.

[0012] The shear failure zone is marked based on the elevation difference of the point cloud data. The volume enclosed by the point cloud before and after shearing is calculated, which is the shear failure volume.

[0013] In one possible implementation, the acquisition of point cloud data of rock joint surfaces before and after shearing includes:

[0014] Clean the joint surfaces before and after shearing, and apply a material with extremely low reflectivity around the rock joint surfaces;

[0015] The joint surface was scanned multiple times from different locations using a 3D laser scanning device to obtain 3D point cloud data of the joint surface before and after shearing.

[0016] In one possible implementation, the calculation of the point cloud data center and its relocation to the origin to obtain the translated point cloud data includes:

[0017] Calculate the average value of the point cloud of rock joint surfaces on the x-axis. Average value on the y-axis Average value on the z-axis

[0018] Then all the coordinates (x) of the points i y i , z i )minus This involves moving the center of the 3D point cloud to the origin of the coordinate system to obtain the translated point cloud data.

[0019] In one possible implementation, rotating the point cloud data so that its normal vector aligns with the z-axis direction, based on the angle between the normal vector and the z-axis of the three-dimensional coordinate system and the cross product, includes:

[0020] Calculate the covariance matrix of the translated point cloud data. The matrix form is as follows:

[0021]

[0022] In the formula, Σ xx Σ represents the variance of column x; yy Σ represents the variance of column y; zz Σ represents the variance of column z; xy , Σ yx Σ is the covariance between columns x and y; xz , Σ zxIt is the covariance between columns x and z; Σ zy , Σ yz Let be the covariance between columns y and z. The formula for calculating the variance is:

[0023]

[0024] The formula for calculating covariance is:

[0025]

[0026] In the formula, n is the number of samples; These are the average values ​​of columns x and y, respectively.

[0027] The covariance matrix is ​​decomposed using QR decomposition. The elements on the diagonal after decomposition are eigenvalues, and the column vectors of the orthogonal matrix are eigenvectors. The eigenvector corresponding to the smallest eigenvalue is the normal vector of the point cloud plane.

[0028] Then, the angle between the normal vector and the vector (0, 0, 1) is calculated to obtain the rotation angle θ, and the cross product of the two vectors is used to obtain the rotation axis vector. The rotation axis vector is then normalized to obtain the unit vector of the rotation axis.

[0029] Then, the translated point cloud is rotated using the Rodriguez rotation formula to obtain a point cloud with its normal direction aligned with the z-axis direction, i.e., the rotated 3D point cloud. The Rodriguez rotation formula is as follows:

[0030] v′=vcosθ+(u·v)u(1-cosθ)+(u×v)sinθ

[0031] In the formula, v' is the rotated spatial vector; v is the spatial vector to be rotated; θ is the rotation angle; and u is the unit vector of the rotation axis.

[0032] In one possible implementation, the step of projecting the point cloud data onto the xoy plane, rotating the point cloud data according to the angle between the boundary and the x-axis, and adjusting the boundary to the coordinate axis includes:

[0033] Set the z value in the rotated 3D point cloud to 0 to obtain the 2D point cloud data projected onto the xoy plane, and then find the convex hull of the projected point cloud.

[0034] After obtaining the convex hull, sort the point cloud according to the angle of each point relative to the center of the convex hull;

[0035] Then calculate the first point (x) in the sorted point cloud. j y j ) and the second point (x) j+1 y j+1 The slope of the line connecting (y) j+1 -y j) / (x j+1 -x j Continue to calculate arctan((y) j+1 -y j ) / (x j+1 -x j This value is the angle α between the line connecting the two axes and the x-axis.

[0036] Rotate all point cloud data by an angle -α around the z-axis, and subtract the minimum value on the y-axis from the y-value of all points in the rotated point cloud data, and subtract the minimum value on the x-axis from the x-value of all points. In other words, the left and lower boundaries of the point cloud are moved to the y-axis and x-axis respectively.

[0037] After performing the above processing on the point cloud data before and after cropping, two sets of point cloud data with coarse registration are obtained.

[0038] In one possible implementation, finding the convex hull of the projected point cloud includes:

[0039] Find the leftmost endpoint A and the rightmost endpoint B in the point cloud; these two points must lie on the convex hull.

[0040] A line segment is formed with A and B as endpoints, dividing all points into two groups: one above and one below the line segment AB.

[0041] For each group, find the point C that is farthest from line segment AB. Point C is also on the convex hull. This process requires using the formula for the distance from a point to a line segment:

[0042]

[0043] In the formula, (x1, y1) are the coordinates of point A, (x2, y2) are the coordinates of point B, and the coordinates of any point are (x, y).

[0044] Recursively apply the above steps to line segments AC and CB respectively, until no point is above the line segment.

[0045] In one possible implementation, the method of accurately registering the two sets of point clouds before and after cropping using point cloud registration includes:

[0046] The joint surface point cloud data after coarse registration and cut is used as the source point cloud, and the joint surface point cloud data before coarse registration and cut is used as the target point cloud.

[0047] The nearest point iterative algorithm is used to make the source point cloud continuously approach the target point cloud until convergence, resulting in two sets of finely registered point cloud data.

[0048] In one possible implementation, the step of calculating the volume enclosed by the point clouds before and after shearing, i.e., the shearing failure volume, based on the elevation difference marking of the shearing failure zone in the point cloud data, includes:

[0049] Interpolate the two sets of finely registered point cloud data respectively;

[0050] Then calculate the coordinates (x) of these two sets of point cloud data at the same coordinate point. i y i ) place z i The difference Δz;

[0051] A threshold is set; values ​​exceeding this threshold are considered to indicate shear failure and are marked accordingly.

[0052] The volume enclosed by the two sets of point clouds before and after shearing is then calculated by numerical integration, which is the shear failure volume.

[0053] In one possible implementation, interpolation is performed using cubic polynomial interpolation, as shown in the following formula:

[0054]

[0055] In the formula, a ij These are undetermined coefficients, which need to be determined using known data points and their derivatives at boundary points.

[0056] In this embodiment, a series of point cloud data optimization processes are performed on the joint surface point cloud data before and after shearing, followed by precise registration. This ensures the accuracy and reliability of shear failure zone identification and volume calculation, enabling precise identification of shear failure areas and calculation of their volumes. This method is more accurate than traditional calculation methods based on photographs or two-dimensional data. It is not only suitable for laboratory rock mechanics experiments but can also be extended to rock mass structural surface analysis in practical engineering. Attached Figure Description

[0057] Figure 1 A flowchart illustrating a method for identifying and calculating the volume of shear failure zones on rock joint surfaces, provided in an embodiment of this application;

[0058] Figure 2 This is an image of a joint surface in a natural rock.

[0059] Figure 3 Images of the 3D point cloud of the joint surfaces before and after cropping;

[0060] Figure 4 for Figure 3 Image of the translated joint surface;

[0061] Figure 5 for Figure 4 Image of the joint surface after rotation;

[0062] Figure 6 Images of the shear failure zone of the marked joint surfaces;

[0063] Figure 7Image showing the calculated shear failure volume of the joint surface. Detailed Implementation

[0064] The present solution will now be described in conjunction with the accompanying drawings and specific embodiments.

[0065] See Figure 1 The method for identifying and calculating the volume of shear failure zones at rock joint surfaces provided in this embodiment includes:

[0066] S101: Collect point cloud data of rock joint surfaces before and after shearing, and delete duplicate points with the same coordinates.

[0067] Clean the joint surface 1 before shearing, then place it in a bright place, and apply a layer of black material 2 around the rock joint surface, such as... Figure 2 As shown. Using a 3D laser scanning device, multiple partitioned scans of the joint surface were performed from different locations to obtain 3D point cloud data of the rock joint surface before shearing. The same method was used to obtain 3D point cloud data of the joint surface after shearing. See [link to documentation]. Figure 3 , Figure 3 No. 3 represents the point cloud of the joint surface before shearing, and No. 4 represents the point cloud of the joint surface after shearing. After obtaining the three-dimensional point cloud data of the joint surface after shearing, duplicate points with the same coordinates in the joint surface point cloud are deleted.

[0068] S102, calculate the point cloud data center and move it to the origin of the coordinate system to obtain the translated point cloud data.

[0069] Calculate the average value of all points on the joint surface before shearing along the x-axis. Average value on the y-axis Average value on the z-axis Then all the coordinates (x) of the points i y i , z i )minus This involves moving the center of the 3D point cloud to the origin of the coordinate system, resulting in translated point cloud data. The same method is used to process the point cloud data of the sheared joint surfaces, such as... Figure 4 As shown.

[0070] S103, based on the angle between the normal vector and the z-axis of the three-dimensional coordinate system and the cross product result, rotate the point cloud data so that its normal is consistent with the z-axis direction.

[0071] The covariance matrix of the translated point cloud data was calculated using principal component analysis, and the matrix form is as follows:

[0072]

[0073] In the formula, Σ xx Σ represents the variance of column x; yy Σ represents the variance of column y;zz Σ represents the variance of column z; xy , Σ yx Σ is the covariance between columns x and y; xz , Σ zx It is the covariance between columns x and z; Σ zy , Σ yz Let be the covariance between columns y and z. The formula for calculating the variance is:

[0074]

[0075] The formula for calculating covariance is:

[0076]

[0077] In the formula, n is the number of samples; These are the average values ​​of columns x and y, respectively.

[0078] The covariance matrix is ​​decomposed using QR decomposition. The elements on the diagonal after decomposition are eigenvalues, and the column vectors of the orthogonal matrix are eigenvectors. The eigenvector corresponding to the smallest eigenvalue is the normal vector of the point cloud plane. Then, the angle θ (rotation angle) between the normal vector and the vector (0, 0, 1) and the cross product of the two vectors (rotation axis) are calculated, and the rotation axis vector is normalized. Finally, the Rodriguez rotation formula is used to rotate the translated point cloud, resulting in a point cloud with the normal direction aligned with the z-axis direction, i.e., the rotated 3D point cloud. The Rodriguez rotation formula is as follows:

[0079] v′=vcosθ+(u·v)u(1-cosθ)+(u×v)sinθ

[0080] In the formula, v' is the rotated spatial vector; v is the spatial vector to be rotated; θ is the rotation angle; and u is the unit vector of the rotation axis. (See also...) Figure 5 for Figure 4 Image of the jointed surface after rotation.

[0081] S104: After projecting the point cloud data onto the xoy plane, rotate the point cloud data according to the angle between the boundary and the x-axis, and adjust the boundary to the coordinate axis.

[0082] Set the z-value of the rotated 3D point cloud to 0 to obtain the 2D point cloud data projected onto the xoy plane. Then, find the convex hull of the projected point cloud. The specific steps are as follows:

[0083] 1) Find the leftmost endpoint A and the rightmost endpoint B in the point cloud. These two points must lie on the convex hull.

[0084] 2) Using A and B as endpoints, form a line segment. Divide all points into two groups: one above and one below the line segment AB.

[0085] 3) For each group, find the point C that is farthest from line segment AB. Point C is also on the convex hull. This process requires using the formula for the distance from a point to a line segment:

[0086]

[0087] In the formula, (x1, y1) are the coordinates of point A, (x2, y2) are the coordinates of point B; the coordinates of any point are (x, y).

[0088] 4) Recursively apply steps 1) to 3) to line segments AC and CB respectively until no point is above the line segment.

[0089] After obtaining the convex hull, sort the point cloud according to the angle of each point relative to the center of the convex hull. Then calculate the first point (x) in the sorted point cloud. j y j ) and the second point (x) j+1 y j+1 The slope of the line connecting (y) j+1 -y j ) / (x j+1 -x j Continue to calculate arctan((y) j+1 -y j ) / (x j+1 -x j The value is the angle α between the connecting line and the x-axis. Rotate all point cloud data around the z-axis by an angle -α. Then, subtract the minimum y-value from the minimum x-value for all points in the rotated point cloud data, and subtract the minimum x-value from the minimum x-value for all points. This essentially moves the left and lower boundaries of the point cloud onto the y-axis and x-axis, respectively. After performing this process on the point cloud data before and after cropping, two sets of coarsely registered point cloud data are obtained.

[0090] S105 uses a point cloud registration method to accurately register the two sets of point clouds before and after cutting.

[0091] Using the coarsely registered joint surface point cloud data after cutting as the source point cloud and the coarsely registered joint surface point cloud data before cutting as the target point cloud, the nearest point iteration algorithm is used to make the source point cloud continuously approach the target point cloud until convergence, thus obtaining two sets of finely registered point cloud data.

[0092] S106. Based on the elevation difference of the point cloud data, mark the shear failure zone and calculate the volume enclosed by the point cloud before and after shearing, which is the shear failure volume.

[0093] Interpolate the two sets of point cloud data that have been precisely registered, and then calculate the coordinates (x, y) of the two sets of point cloud data at the same point. i y i ) place z iThe difference Δz. A threshold of 0.3 mm is set; Δz exceeding this threshold is considered shear failure and is marked accordingly. Figure 6 The central failure zone is shown as number 5. The volume enclosed by the two sets of point clouds before and after shearing is then calculated; this is the shear failure volume. See [reference needed]. Figure 7 , Figure 7 The region indicated by the number 6 is the shear failure volume region.

[0094] In this embodiment, cubic polynomial interpolation is used, and the formula is as follows:

[0095]

[0096] In the formula, a ij These are undetermined coefficients, which need to be determined using known data points and their derivatives at boundary points.

[0097] In this application embodiment, "at least one" refers to one or more, and "more than one" refers to two or more. "And / or" describes the relationship between related objects, indicating that three relationships can exist. For example, A and / or B can represent the existence of A alone, the simultaneous existence of A and B, or the existence of B alone. A and B can be singular or plural. The character " / " generally indicates that the preceding and following related objects are in an "or" relationship. "At least one of the following" and similar expressions refer to any combination of these items, including any combination of single or plural items. For example, at least one of a, b, and c can represent: a, b, c, ab, ac, bc, or abc, where a, b, and c can be single or multiple.

[0098] The above description is merely a specific embodiment of this application. Any variations or substitutions that can be easily conceived by those skilled in the art within the scope of the technology disclosed in this application should be included within the protection scope of this application. The protection scope of this application should be determined by the protection scope of the claims.

Claims

1. A method for identifying and calculating the volume of shear failure zones at rock joint surfaces, characterized in that, include: Collect point cloud data of rock joint surfaces before and after shearing, and delete duplicate points with the same coordinates; The point cloud data center is calculated and moved to the origin to obtain the translated point cloud data, including: Calculate the point cloud of rock joint surfaces x Average value on axis ,exist y Average value on axis ,exist z Average value on axis ; Then take the coordinates of all points ( x i , y i , z i )minus( , , That is, moving the center of the 3D point cloud to the origin of the coordinate system to obtain the translated point cloud data; Based on the normal vector and the three-dimensional coordinate system z The angle between the axes and the cross product result, rotate the point cloud data so that its normal is aligned with the axis. z The axial direction is consistent, including: Calculate the covariance matrix of the translated point cloud data. The matrix form is as follows: In the formula, for x The variance of the column; for y The variance of the column; for z The variance of the column; , for x and y Covariance between columns; , yes x and z Covariance between columns; , for y and z The covariance between columns, where the variance is calculated using the following formula: The formula for calculating covariance is: In the formula, n The number of samples; , Samples x column sum y The average of the column; The covariance matrix is ​​decomposed using QR decomposition. The elements on the diagonal after decomposition are eigenvalues, and the column vectors of the orthogonal matrix are eigenvectors. The eigenvector corresponding to the smallest eigenvalue is the normal vector of the point cloud plane. Then calculate the angle between the normal vector and the vector (0, 0, 1) to obtain the rotation angle. θ The rotation axis vector is obtained by taking the cross product of the two vectors and normalizing the rotation axis vector to obtain the unit vector of the rotation axis. Then, using the Rodriguez rotation formula, the translated point cloud is rotated to obtain the normal direction and... z For point clouds with consistent axial directions, i.e., rotated 3D point clouds, the Rodriguez rotation formula is as follows: In the formula, v' The space vector after rotation; v Let be the space vector to be rotated; θ The rotation angle; u The unit vector is the axis of rotation; Point cloud data in xoy After planar projection, based on the boundary and x Rotate the point cloud data by the included angle of the axes and adjust the boundaries to the coordinate axes; The point cloud registration method is used to accurately register the two sets of point clouds before and after cutting. The shear failure zone is marked based on the elevation difference of the point cloud data. The volume enclosed by the point cloud before and after shearing is calculated, which is the shear failure volume.

2. The method for identifying and calculating the volume of shear failure zones on rock joint surfaces according to claim 1, characterized in that, The collected point cloud data of rock joint surfaces before and after shearing includes: Clean the joint surfaces before and after shearing, and apply a material with extremely low reflectivity around the rock joint surfaces; The joint surface was scanned multiple times from different locations using a 3D laser scanning device to obtain 3D point cloud data of the joint surface before and after shearing.

3. The method for identifying and calculating the volume of shear failure zones on rock joint surfaces according to claim 1, characterized in that, The point cloud data is in xoy After planar projection, based on the boundary and x Rotate the point cloud data by the included angle of the axes and adjust the boundaries to the coordinate axes, including: The rotated 3D point cloud z Setting the value to 0 will result in projection onto... xoy Two-dimensional point cloud data in a plane, and then find the convex hull of the projected point cloud; After obtaining the convex hull, sort the point cloud according to the angle of each point relative to the center of the convex hull; Then calculate the first point in the sorted point cloud. x j , y j ) and the second point ( x j+1 , y j+1 The slope of the connecting line Continue to find This value is the connection between x Angle between axes α ; by z The axis is the rotation axis, rotating all point cloud data. α Angle, and rotate all points of the point cloud data. y value minus y Minimum value on the axis, all points x value minus x The minimum value on the axis refers to moving the left and lower boundaries of the point cloud to the minimum values ​​on the axis. y shaft and x On the axis; After performing the above processing on the point cloud data before and after cropping, two sets of point cloud data with coarse registration are obtained.

4. The method for identifying and calculating the volume of shear failure zones at rock joint surfaces according to claim 2, characterized in that, The process of finding the convex hull of the projected point cloud includes: Find the leftmost endpoint A and the rightmost endpoint B in the point cloud. These two points must lie on the convex hull. A line segment is formed with A and B as endpoints, dividing all points into two groups: one above and one below the line segment AB. For each group, find the point C that is farthest from line segment AB. Point C is also on the convex hull. This process requires using the formula for the distance from a point to a line segment: In the formula, ( x 1, y 1) are the coordinates of point A, ( x 2, y 2) Let B be the coordinates of point B, and let the coordinates of any point be ( ). x , y ); Recursively apply the above steps to line segments AC and CB respectively, until no point is above the line segment.

5. The method for identifying and calculating the volume of shear failure zones on rock joint surfaces according to claim 3 or 4, characterized in that, The method of accurately registering two sets of point clouds before and after cropping using point cloud registration includes: The joint surface point cloud data after coarse registration and cut is used as the source point cloud, and the joint surface point cloud data before coarse registration and cut is used as the target point cloud. The nearest point iterative algorithm is used to make the source point cloud continuously approach the target point cloud until convergence, resulting in two sets of finely registered point cloud data.

6. The method for identifying and calculating the volume of shear failure zones on rock joint surfaces according to claim 5, characterized in that, The step of marking the shearing failure zone based on the elevation difference of the point cloud data and calculating the volume enclosed by the point cloud before and after shearing, which is the shearing failure volume, includes: Interpolate the two sets of finely registered point cloud data respectively; Then calculate the coordinates of these two sets of point cloud data at the same coordinate point ( x i , y i ) z i The difference Δ z ; Set a threshold, Δ z Anything exceeding this threshold is considered to have experienced shear failure and is marked accordingly; The volume enclosed by the two sets of point clouds before and after shearing is then calculated by numerical integration, which is the shear failure volume.

7. The method for identifying and calculating the volume of shear failure zones on rock joint surfaces according to claim 6, characterized in that, The interpolation uses cubic polynomial interpolation, and the formula is as follows: In the formula, These are undetermined coefficients, which need to be determined using known data points and their derivatives at boundary points.

Images