Two-dimensional profile extraction method for maximum roughness of rock joint surface

Abstract

This application discloses a two-dimensional contour extraction method for the maximum roughness of rock joint surfaces. The method obtains a three-dimensional point cloud of a natural rock joint surface; denoises and rotates the point cloud so that the point cloud normal vector is collinear with the z-axis of the three-dimensional coordinate system; the geometric center of the rotated point cloud is moved to the origin of the three-dimensional coordinate system; the point cloud within an infinitely long cylinder with a base radius of R and the z-axis as its axis is selected and retained; given a profile equation and a set number of rotations n, the given profile is rotated counterclockwise by 180 / n degrees each time, using the z-axis as the rotation axis, to obtain the two-dimensional point cloud on the given profile; the two-dimensional point cloud is interpolated, and the roughness of the interpolated two-dimensional point cloud is calculated to find the maximum roughness and the corresponding two-dimensional point cloud; the interpolated two-dimensional point cloud with the maximum roughness is output in the DXF format commonly used by discrete element method (DEM) software. This method provides a powerful tool for two-dimensional joint roughness analysis and plays an important role in promoting the development of rock mechanics testing.

Description

Technical Field

[0001] This application relates to the field of rock mechanics experimental technology, specifically to a method for extracting the two-dimensional contour of the maximum roughness of rock joint surfaces. Background Technology

[0002] The geometric characteristics of rock joint surfaces play a decisive role in the mechanical properties and stability of rock masses. With the development of 3D scanning technology, point cloud data has become an important tool for characterizing the 3D roughness of rock joints. However, analyzing joint surface roughness using point cloud data is extremely complex and computationally inefficient due to the massive amount of data. Therefore, current research on rock joint surface roughness still mainly focuses on 2D morphology. However, how to scientifically and effectively convert 3D point clouds into 2D point clouds and extract the 2D contour of the joint surface under maximum roughness remains a key problem that needs to be solved.

[0003] Currently, existing patents provide methods for establishing three-dimensional discrete element models of jointed rock masses based on point cloud data. For example, Chinese patent CN112784403B discloses a numerical simulation method for establishing a discrete element model of jointed rock masses based on point cloud data. Its main steps are: acquiring a point cloud dataset of the jointed rock mass, performing denoising and filtering, calculating the attitude of the structural planes, and establishing a three-dimensional model of the structural planes using the point cloud data. The model is then processed into a data file readable by 3DEC software using Fracman and Rhino software, and then imported into 3DEC software for numerical analysis to generate the final three-dimensional model of the jointed rock mass.

[0004] The above method can generate a three-dimensional numerical model of jointed rock mass, but it cannot scientifically and reasonably generate a two-dimensional joint surface profile. 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 extracting the two-dimensional contour line of the maximum roughness of a rock joint surface, including:

[0007] A three-dimensional laser scanning device was used to scan the joint surfaces of natural rocks to obtain three-dimensional point clouds of the joint surfaces;

[0008] The three-dimensional point cloud of the joint surface is denoised to obtain a denoised three-dimensional point cloud;

[0009] Rotate the denoised 3D point cloud so that its normal vector is collinear with the unit vector of the z-axis of the 3D coordinate system;

[0010] The geometric center of the rotated 3D point cloud is moved to the origin of the 3D coordinate system to obtain the centered 3D point cloud.

[0011] The centered 3D point cloud is filtered and retained within the range of an infinitely long cylinder with the z-axis of the coordinate system as the axis and the base radius of R.

[0012] Given a profile equation, set the number of rotations n, rotate the given profile counterclockwise with the z-axis as the rotation axis, rotating 180 / n degrees each time, to obtain the two-dimensional point cloud of the filtered point cloud on the profile at each rotation.

[0013] For each of the n sets of two-dimensional point clouds, interpolate and calculate the roughness of the two-dimensional point cloud after each interpolation, and find the maximum roughness and the corresponding two-dimensional point cloud.

[0014] The two-dimensional point cloud after maximum roughness under-interpolation is output in the dxf format commonly used by discrete element software.

[0015] In one possible implementation, the denoising of the three-dimensional point cloud of the joint surface to obtain the denoised three-dimensional point cloud includes:

[0016] For each point p in the 3D point cloud i Calculate its K neighboring points (p i1 p i2 , ..., p iK mean and the standard deviation σ of the distance from the neighborhood points to the mean i The formulas are as follows:

[0017]

[0018] Set a threshold std, if Exceeded Then point p i It is considered noise and should be removed.

[0019] In one possible implementation, rotating the denoised 3D point cloud so that its normal vector is collinear with the unit vector of the 3D coordinate system z-axis includes:

[0020] Calculate the rotation axis *a* and rotation angle *α* between the point cloud normal vector *n* and the z-axis of the coordinate system. The rotation axis *a* is obtained by the cross product of the normal vector *n* and the z-axis unit vector *k*, with the formula: *a = n × k*. The rotation angle *α* is calculated by the dot product of the normal vector *n* and the z-axis unit vector *k*, with the formula:

[0021] The rotation matrix R is calculated using Rodrigues' rotation formula, which is:

[0022] R = I + sin(α)K + (1 - cos(α))K 2

[0023] In the formula, K is the antisymmetric matrix of the rotation axis, and α is the rotation angle;

[0024] Multiply the denoised 3D point cloud by the rotation matrix R to obtain the rotated point cloud.

[0025] In one possible implementation, the geometric center of the rotated 3D point cloud is moved to the origin of the 3D coordinate system to obtain a centered 3D point cloud, including:

[0026] Assume a point cloud consists of N points, each with coordinates p. i =(x i ,y i ,z i If the centroid C is 0, then the formula for calculating the centroid C is:

[0027]

[0028] Subtracting the centroid coordinates from the coordinates of each point in the point cloud yields the centered point cloud, using the following formula:

[0029] p i ′=p i -C=(x i ′,y i ′,z i ′)=(x i -C x ,y i -C y ,z i -C z ).

[0030] In one possible implementation, filtering and retaining the 3D point cloud within an infinitely long cylinder with a base radius of R and the z-axis as its axis includes: setting a radius R, and for each point (x... i y i , z i ), calculated value If d i If the value is greater than R, then delete the point.

[0031] In one possible implementation, the profile equation is:

[0032] y-tan(θ)x=0

[0033] In the formula, θ is the angle between the intersection of the given profile and the xoy plane and the positive x-axis.

[0034] In one possible implementation, obtaining the two-dimensional point cloud of the filtered point cloud on the cross-section at each rotation includes:

[0035] Set tolerance F;

[0036] Calculate the distance D from the filtered point cloud to the given profile. i The formula is:

[0037]

[0038] In the formula, x i y i Let x and y be the x-axis and y-axis coordinates of the i-th point in the point cloud;

[0039] If D i If ≤F, then the point is considered to be on the given profile.

[0040] In one possible implementation, the interpolation of each of the n sets of two-dimensional point clouds includes:

[0041] Define a new planar coordinate system on the given cross section: the origin of the new planar coordinate system is the same as that of the original three-dimensional coordinate system, the x-axis of the new planar coordinate system is denoted as the x' axis, which is the intersection of the given cross section and the xoy plane, with its positive direction being the positive direction of the original x-axis and its negative direction being the negative direction of the original x-axis, and the y-axis of the new planar coordinate system being the z-axis of the original coordinate system, denoted as the y' axis;

[0042] The formula for coordinate transformation between the new coordinate system and the original coordinate system is: y i ′=z i , where x i y i z i x represents the coordinates of the i-th point in the original coordinate system. i The sign of the plus or minus sign is determined by the x-axis of the original coordinate system. i The sign determines the sign;

[0043] Sort the x' axis values ​​from smallest to largest, then perform the following interpolation calculations:

[0044]

[0045] In the formula, x i y i These are the x and y coordinates of the i-th point in the original coordinate system.

[0046] In one possible implementation, the roughness of the 2D point cloud is preferentially expressed as the root mean square of the first derivative of the profile, Z², as shown in the following formula:

[0047]

[0048] In the formula, N is the number of interpolated points in the profile data; L is the projected length of the profile contour line; x i y i Let x and y be the x and y coordinates of the i-th point, respectively.

[0049] In this embodiment, based on three-dimensional point cloud data of joint surfaces, the two-dimensional joint contours and their coordinate values ​​on any cross-section are accurately obtained. From denoising the point cloud data and filtering the point cloud within the target range to finally outputting the two-dimensional point cloud after interpolation of the maximum roughness, a powerful tool is provided for two-dimensional joint roughness analysis, which plays an important role in promoting the development of the field of rock mechanics testing. Attached Figure Description

[0050] Figure 1 A flowchart illustrating a method for extracting the maximum roughness of a rock joint surface, provided in an embodiment of this application.

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

[0052] Figure 3 for Figure 2 Images of three-dimensional point clouds of joint surfaces in natural rocks;

[0053] Figure 4 To Figure 3 Rotated point cloud image;

[0054] Figure 5 An image of a point cloud after centralization;

[0055] Figure 6 A top view of the filtered point cloud;

[0056] Figure 7 The image is a two-dimensional point cloud image after interpolation with maximum roughness.

[0057] Figure 8 This is the two-dimensional contour line of the joint surface after interpolation of maximum roughness in AutoCAD. Detailed Implementation

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

[0059] See Figure 1 The method for extracting the two-dimensional contour of the maximum roughness of rock joint surfaces provided in this embodiment includes:

[0060] S101 uses a three-dimensional laser scanning device to scan the joint surfaces of natural rocks to obtain three-dimensional point clouds of the joint surfaces.

[0061] See Figure 2 To obtain an image of a natural rock joint surface, the surface of the natural joint is cleaned, and then a 3D laser scanning device is used to scan the key study area of ​​the natural rock joint surface to obtain 3D point cloud data of the rock joint surface, such as... Figure 3 The image shown is a 3D point cloud image of a natural rock joint surface.

[0062] S102, the three-dimensional point cloud of the joint surface is denoised to obtain the denoised three-dimensional point cloud.

[0063] Denoising of 3D point clouds on joint surfaces includes:

[0064] For each point p in the point cloud i Calculate its 30 neighboring points (p i1 p i2 , ..., p i30 mean and the standard deviation σ of the distance from the neighborhood points to the mean i The formulas are as follows:

[0065]

[0066] Set a threshold std = 1, if Exceeded Then point p i It is considered noise and should be removed.

[0067] S103, rotate the denoised 3D point cloud so that its normal vector is collinear with the unit vector of the z-axis of the 3D coordinate system.

[0068] 1) Calculate the rotation matrix

[0069] Calculate the rotation axis *a* and rotation angle *α* between the point cloud normal vector *n* and the z-axis of the coordinate system. The rotation axis *a* is obtained by the cross product of the normal vector *n* and the z-axis unit vector *k*(0, 0, 1), with the formula: *a = n × k*. The rotation angle *α* is calculated by the dot product of the normal vector *n* and the z-axis unit vector *k*(0, 0, 1), with the formula:

[0070]

[0071] 2) Calculate the rotation matrix R using Rodriguez's rotation formula, which is:

[0072] R = I + sin(α)K + (1 - cos(α))K 2

[0073] In the formula, K is the antisymmetric matrix of the rotation axis, and α is the rotation angle.

[0074] 3) Multiply the denoised 3D point cloud by the rotation matrix R to obtain the rotated point cloud. See [link to relevant documentation]. Figure 4 This is an image of the point cloud after rotation.

[0075] S104, the geometric center of the rotated 3D point cloud is moved to the origin of the 3D coordinate system to obtain the centered 3D point cloud.

[0076] 1) Calculate the centroid of the point cloud

[0077] Assume a point cloud consists of N points, each with coordinates p. i =(x i ,y i ,z i If the centroid C is 0, then the formula for calculating the centroid C is:

[0078] 2) Centralize the point cloud

[0079] Subtract the centroid coordinates from the coordinates of each point in the point cloud to obtain the centered point cloud (see [link]). Figure 5 (For the image of the centered point cloud), the formula is:

[0080] p i ′=p i -C=(x i ′,y i ′,z i ′)=(x i -C x ,y i -C y ,z i -C z ).

[0081] S105, the centered 3D point cloud is filtered and retained within the range of an infinitely long cylinder with the z-axis of the coordinate system as the axis and the base radius of R.

[0082] Set the radius R to 0.05, for each point (x) in the point cloud i y i , z i ), calculated value If d i If the value is greater than R, then delete the point. See [link to relevant documentation]. Figure 6 This is a top view of the filtered point cloud.

[0083] S106, Given the profile equation, set the number of rotations n, rotate the given profile counterclockwise with the z-axis as the rotation axis, each rotation is 180 / n degrees, to obtain the two-dimensional point cloud of the filtered point cloud on the profile at each rotation.

[0084] Given the profile equation: y-tan(θ)x=0, set the number of rotations to 180, rotate the given profile counterclockwise with the z-axis as the rotation axis, rotating 1 degree each time, that is, from 0 to 180°, with an interval of 1°, to obtain the two-dimensional point cloud on the profile after each rotation, for a total of 180 sets of two-dimensional point clouds.

[0085] Set tolerance F = 0.0002;

[0086] Calculate the distance D from the filtered point cloud to the given profile. i The formula is:

[0087]

[0088] In the formula, x i y i Let x and y be the coordinates of the i-th point in the point cloud.

[0089] S107, interpolate n groups of two-dimensional point clouds respectively, calculate the roughness of the two-dimensional point cloud after each interpolation, and find the maximum roughness and the corresponding two-dimensional point cloud.

[0090] The steps for 2D point cloud interpolation are as follows:

[0091] 1) Define a new planar coordinate system on the given cross section: The origin of the new planar coordinate system is the same as that of the original three-dimensional coordinate system. The x-axis of the new planar coordinate system is denoted as the x' axis, which is the intersection of the given cross section and the xoy plane. Its positive direction is the positive direction of the original x-axis, and its negative direction is the negative direction of the original x-axis. The y-axis of the new planar coordinate system is the z-axis of the original coordinate system.

[0092] 2) The coordinate transformation formula between the new coordinate system and the original coordinate system is: y i ′=z i , where x i y i z i x represents the coordinates of the i-th point in the original coordinate system. i The sign of the plus or minus sign is determined by the x-axis of the original coordinate system. i The sign determines the sign.

[0093] 3) Sort the x' axis values ​​from smallest to largest, and then perform the following interpolation calculations:

[0094]

[0095] In the formula, x i y i Let x and y be the x and y coordinates of the i-th point in the original coordinate system.

[0096] Roughness is calculated using the root mean square of the first derivative of the profile, Z². The formula is as follows:

[0097]

[0098] In the formula, N is the number of interpolated points in the profile data, L is the projected length of the profile contour line, and x i y i These are the x and y coordinates of the i-th point, see [reference]. Figure 7 This is the two-dimensional point cloud after interpolation under maximum roughness.

[0099] S108, output the two-dimensional point cloud after the maximum roughness under-interpolation as the dxf format commonly used by discrete element software.

[0100] Using the Python third-party software library ezdxf, the 2D point cloud after maximum roughness interpolation is output in the dxf format commonly used by discrete element method (DEM) software. See [link to DEM software]. Figure 8 This is the two-dimensional contour line of the joint surface after interpolation of maximum roughness in AutoCAD.

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

[0102] 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 extracting the two-dimensional contour of the maximum roughness of a rock joint surface, characterized in that, include: A three-dimensional laser scanning device was used to scan the joint surfaces of natural rocks to obtain three-dimensional point clouds of the joint surfaces; Denoising the 3D point cloud of the joint surface to obtain a denoised 3D point cloud includes: For each point in the 3D point cloud pi Calculate its K neighborhood points ( pi 1, pi 2, ..., piK mean and the standard deviation of the distance from the neighborhood points to the mean The formulas are as follows: ， Set a threshold std, if Exceeded Then point p i This is considered noise and should be removed. Rotate the denoised 3D point cloud so that its normal vector is aligned with the 3D coordinate system. z The unit vectors of the axes are collinear, including: Calculate point cloud normal vectors With coordinate system z Rotation axis between axes and rotation angle α Among them, the rotation axis Through the normal vector and z Axial unit vector The cross product is obtained, and the formula is: Rotation angle α Through the normal vector and z Axial unit vector The dot product is calculated using the following formula: ; The rotation matrix R is calculated using Rodrigues' rotation formula, which is: In the formula, Let be the antisymmetric matrix of the rotation axis. α The rotation angle; Multiply the denoised 3D point cloud by the rotation matrix R to obtain the rotated point cloud; Moving the geometric center of the rotated 3D point cloud to the origin of the 3D coordinate system to obtain a centered 3D point cloud includes: Assume a point cloud consists of N points, and the coordinates of each point are... Then the formula for calculating the centroid C is: Subtracting the centroid coordinates from the coordinates of each point in the point cloud yields the centered 3D point cloud, as shown in the formula: ； The centered 3D point cloud is filtered and retained in coordinate system. z The axis is the axis, and the radius of the base is... R A three-dimensional point cloud within the range of an infinitely long cylinder; Given the profile equation, set the number of rotations. n ,by z The axis is a rotation axis, rotating the given profile counterclockwise by 180 degrees each time. n The degree is used to obtain the two-dimensional point cloud on the cross-section after each rotation of the filtered point cloud; To each n Interpolate the two-dimensional point clouds in groups, calculate the roughness of the two-dimensional point cloud after each group of interpolation, and find the maximum roughness and the corresponding two-dimensional point cloud. The two-dimensional point cloud after maximum roughness under-interpolation is output in the dxf format commonly used by discrete element software.

2. The method for extracting the two-dimensional contour line of the maximum roughness of rock joint surfaces according to claim 1, characterized in that, The centered 3D point cloud is filtered and retained in coordinate system. z The axis is the axis, and the radius of the base is... R A 3D point cloud within the range of an infinitely long cylinder, including: a set radius R For each point in the point cloud ( x i , y i , z i ), calculated value ,like d i > R If so, then delete that point.

3. The method for extracting the two-dimensional contour of the maximum roughness of rock joint surfaces according to claim 1, characterized in that, The cross-sectional equation is: In the formula, θ For a given profile and xoy The intersection of planes and x The angle between the positive axes.

4. The method for extracting the two-dimensional contour of the maximum roughness of rock joint surfaces according to claim 3, characterized in that, Obtaining the two-dimensional point cloud on the cross-section after each rotation of the filtered point cloud includes: Set tolerance F ; Calculate the distance from the filtered point cloud to the given profile. D i The formula is: In the formula, x i , y i For the first point cloud i points x shaft and y Axis coordinates; like D i ≤ F If the point is on the given cross section, then it is considered to be on the given cross section.

5. The method for extracting the two-dimensional contour of the maximum roughness of rock joint surfaces according to claim 1, characterized in that, The respective to n Group 2D point cloud interpolation, including: Define a new planar coordinate system on the given cross-section: the origin of the new planar coordinate system is the same as that of the original three-dimensional coordinate system, and the new planar coordinate system... x Axis is denoted as x 'Axis, which is a given profile and xoy The line of intersection of planes, its positive direction is the original. x The positive direction of the axis and the negative direction are the original directions. x The negative direction of the axis, the new plane coordinate system y The axes are the original coordinate system z Axis, denoted as y 'axis; The formula for coordinate transformation between the new coordinate system and the original coordinate system is: , ,in, x i , y i , z i The first in the original coordinate system i The coordinates of each point; The sign changes from the original coordinate system. x shaft x i The sign determines the sign; according to x The values ​​of the axis are sorted from smallest to largest, and then the following interpolation calculations are performed: In the formula, x i , y i The first in the original coordinate system i The x and y coordinates of each point.

6. The method for extracting the two-dimensional contour of the maximum roughness of rock joint surfaces according to claim 5, characterized in that, For roughness of 2D point clouds, the root mean square of the first derivative of the profile is preferred. Z 2 The formula is as follows: In the formula, N This represents the number of points after interpolation of the profile data; L The projected length of the profile outline; x i , y i The first i The x and y coordinates of each point.

Images