Rock structure surface roughness generation method based on power spectrum density

Abstract

The invention belongs to the technical field of rock physical simulation, and particularly relates to a rock structure surface roughness generation method based on power spectral density. A rock structure surface roughness generation method based on power spectral density comprises the following steps: establishing an initialized Fourier coefficient storage array, generating a Fourier matrix through basic calculation parameters, and generating surface data X by using inverse Fourier transform; sampling the obtained surface data, and calculating a joint roughness value JRC through different roughness parameters; establishing a target joint roughness value optimization function, and iteratively obtaining an optimal parameter solution under various characteristic parameter standards through global optimization and local optimization; and for the surface data obtained according to different target joint roughness, 3D printing is used for achieving generation of the surface with the specific roughness. According to the method, repeatable corresponding joint roughness surface generation under various roughness parameter evaluation standards of the rock surface is realized, and the natural rock joint surface can be simulated more accurately.

Description

Technical Field

[0001] The present invention belongs to the technical field of rock physics simulation, and in particular relates to a method for generating rock structure surface roughness based on power spectrum density. Background Art

[0002] The surface roughness of rock joints refers to the degree of unevenness of the joint surface. It involves research in many fields such as geological engineering, earthquake engineering, and water resources engineering, and has a significant impact on the mechanical and hydraulic properties of rock joints.

[0003] Existing characterization methods for the roughness of rock joints are generally categorized as either qualitative or quantitative. For example, the Joint Roughness Coefficient (JRC) proposed by Barton and Choubey is a common qualitative characterization method that evaluates roughness by comparing the actual joint surface with a standard roughness curve. However, this method is susceptible to human subjectivity and has limited accuracy. Quantitative analysis methods primarily include laser scanning, image analysis, and profilometer measurement, but these methods are often highly dependent on equipment performance, are expensive, and require complex data processing. Some methods can only obtain local two-dimensional information and fail to reflect the overall three-dimensional morphology. Random samples obtained from three-dimensional laser scanning of joint surfaces in different studies have poor universality, making them unsuitable for controlled and repeated testing. Furthermore, further studies have struggled to accurately replicate and draw definitive conclusions, hindering their further application in engineering research.

[0004] In the existing technology, the generation of rock joint surface roughness relies on reasonable characteristics to characterize the joint roughness. The main characteristic parameters are: mean square value MS, root mean square RMS, centerline average value CLA, first-order derivative root mean square Z2, second-order inverse root mean square Z3, fluctuation amplitude mean parameter Z4, autocorrelation function ACF and structure function SF to characterize its roughness characteristics. However, some parameters are easily affected by the position of the baseline, and it is difficult to find the baseline. In addition, although the influence of roughness on engineering properties can be considered through rock joint contact surface model experiments and estimation of fluid mechanics behavior, repeated experiments are time-consuming and material-intensive, and lack good repeatability. In addition, due to the complex morphology of the rock structure surface, there are anisotropy and size effect characteristics. Combined with engineering practice, accurate characterization of rock joint roughness is still critical. Summary of the Invention

[0005] In response to the above problems, the purpose of the present invention is to provide a method for generating rock structure surface roughness based on power spectral density. The present invention uses inverse Fourier transform to generate initial fractal surface data, calculates joint roughness values ​​through different roughness parameters, and combines differential evolution algorithm and local optimization to find the optimal solution that meets the target roughness parameters. It realizes the repeatable generation of corresponding joint roughness surfaces under multiple roughness parameter evaluation criteria of rock surfaces, thereby more accurately simulating the natural rock joint surface.

[0006] The technical solution of the present invention is: a method for generating rock structure surface roughness based on power spectrum density, comprising the following steps:

[0007] S1: Establish and initialize the Fourier coefficient storage array, generate the Fourier matrix through basic calculation parameters, and generate surface data X using inverse Fourier transform;

[0008] S2: Sampling the obtained surface data, and calculating the joint roughness value JRC by using different roughness parameters;

[0009] S3: Establishing the target joint roughness value optimization function, and iteratively obtaining the optimal parameter solution under multiple characteristic parameter standards through global optimization and local optimization;

[0010] S4: Based on the surface data obtained for different target joint roughness, 3D printing is used to generate a surface with a specific roughness.

[0011] In step S1, the process of generating surface data X by using inverse Fourier transform includes:

[0012] S11: According to the actual sampling interval and frame size, the number of grid points of the initialization Fourier coefficient storage array is defined, and the Fourier matrix A is generated through basic calculation parameters. The basic calculation parameters include N, H, IS, DF and ST, where N represents the array size. For the generated N×N two-dimensional array, it represents the number of grid points on the fractal surface; H is the fractal dimension, which controls the roughness of the surface and takes a value between 1 and 3; IS is the initial seed. The fractal surface generated by the same initial seed can be reproduced. The value of the initial seed will change the generation sequence of random numbers but will not affect the characteristics of the generated plane; DF is the directionality factor, which controls the phase change of the generated fracture surface in different directions. Considering increasing the change of the frequency component of the generated plane, it is usually within the range of 1±1; ST is the scale threshold, which controls the irregularity of the surface at a specific radius. It is applied to special processing of smaller radii to amplify the characteristics of the overall frequency component and is a real number greater than 0;

[0013] S12: Initialize the matrix A, traverse each point in the array, skip special points and assign values;

[0014] S13: Take the real parts of the values ​​corresponding to the special points (N / 2, 0), (0, N / 2), and (N / 2, N / 2) to meet the symmetry requirements of the inverse Fourier transform;

[0015] S14: Perform a two-dimensional inverse Fourier transform on the matrix A, convert the frequency domain data into the spatial domain and take the real part of the transformation result, and finally generate the surface data X of the fracture surface.

[0016] In step S2, the N×N matrix A is initialized to a zero matrix, and the calculation function for the surface data X is as follows:

[0017] A i，j =0

[0018]

[0019] Where A i,j The Fourier coefficients of each frequency point in the initialization Fourier coefficient storage array established in step S1, Rad i，j Represents the Euclidean distance from each matrix element (i, j) to the center point (N / 2, N / 2);

[0020]

[0021] Where, is the generated random phase, Temp is the generated random amplitude, and U(a, b) represents uniform distribution;

[0022] When Rad>Mismatch:A i，j =Temp

[0023] When Rad<Mismatch:

[0024] Where, is a random phase, used to increase the randomness of the surface; e i·Aniso·(i-j) The anisotropy parameter is introduced to make the surface different in different directions; as the radius increases, the amplitude decreases at the power of (3-H). As the fractal dimension H decreases, the details of the generated plane are smoother, that is, the high-frequency components of the components are suppressed;

[0025] A N / 2，0 =Re(A N / 2，0 ), A 0，N / 2 =Re(A 0，N / 2 ), A N / 2，N / 2 =Re(A N / 2，N / 2 )

[0026] Where Re represents the real part;

[0027]

[0028] Where, Represents a two-dimensional inverse Fourier transform. The returned X is a two-dimensional matrix representing the simulated fracture surface, where the element value represents the "height" of the surface.

[0029] In step S2, the surface data obtained are sampled and the joint roughness value JRC is calculated using different roughness parameters. The specific process is as follows:

[0030] S21: Randomly select 10 index rows and index columns from the obtained surface data to complete horizontal and vertical slicing of the data;

[0031] S22: Initialize a list to store the joint roughness value JRC calculated for each slice;

[0032] S23: Calculate the joint roughness values ​​JRC of each transverse and longitudinal slice based on different roughness parameters, and take the average joint roughness value JRC of each slice.

[0033] The roughness parameters include the root mean square first-order derivative value Z2, the structure function SF, the roughness profile index R p and 3D roughness parameters in:

[0034]

[0035]

[0036] Where Z i is the height value of the i-th point on the surface; Z i+1 is the height of the i+1th point on the surface; X i is the horizontal position of the i-th point; X i+1 is the horizontal position of the i+1th point; Δl is the horizontal distance between two adjacent points, i.e., the sampling interval; M x The number of points used to divide the region, l is the total length.

[0037] The calculation formula for the joint roughness value JRC is calculated by the root mean square first-order derivative value Z2:

[0038] JRC=51.16(Z2) 0.531 -11.44

[0039] The calculation formula for the joint roughness value JRC is calculated by the structure function SF:

[0040] JRC=73.95(SF) 0.266 -11.38

[0041] By roughness profile index Rp The calculation formula for the joint roughness value JRC is:

[0042] JRC=65.9(R p -1) 0.302 -9.65

[0043] Through 3D roughness parameters The process of calculating the joint roughness value JRC is:

[0044] S231: First, traverse the triangles formed by connecting three adjacent points in the generated surface, define the observation reference, and calculate the apparent inclination of the triangular surface using the azimuth and true inclination. The specific formula is:

[0045] tanθ * = -tanθ·cosα

[0046] Where θ * , θ, and α are the apparent inclination, true inclination, and azimuth of the surface connecting the three adjacent points on the rough surface, respectively;

[0047] S232: Then, by fitting the angle threshold-normalized area curve, the dimensionless parameter c is obtained. The specific formula is:

[0048]

[0049] In the formula is the normalized area, which is the surface area with an apparent inclination angle greater than the chosen threshold θ normalized to the total surface area; is the maximum apparent inclination of the surface in the selected analysis direction; C is a dimensionless fitting parameter calculated by nonlinear least squares regression analysis;

[0050] S233: Finally, substitute the joint roughness value function relationship Calculate the corresponding joint roughness value JRC value.

[0051] In step S3, the optimal parameter solution under various characteristic parameter standards is obtained through iteration through global optimization and local optimization. The specific process is as follows:

[0052] S31: define global parameter optimization boundaries;

[0053] S32: introducing the calculation formula for calculating the joint roughness value of different roughness parameters in step S2, defining the target joint roughness value, and calculating different roughness parameters and corresponding actual joint roughness values ​​for the surface data generated by simulation;

[0054] S33: Using the differential evolution algorithm for global optimization, iteratively searching for the optimal parameter solution with the minimum distance between the actual joint roughness value calculated by each roughness parameter and the target joint roughness value;

[0055] S34: The optimal parameters obtained from the global optimization are used as the initial guess values ​​for the local optimization, a new optimization boundary and number of iterations are set, and the search is refined to obtain the parameter combination closest to the target joint roughness value, thereby obtaining the optimal parameter solution.

[0056] In step S4, the surface data obtained for different target joint roughness are used to generate a surface with a specific roughness using 3D printing. The specific process is as follows:

[0057] S41: normalizing the surface data X to ensure that the actual size of the 3D printing meets the requirements;

[0058] S42: Convert the 2D height matrix into a 3D grid pattern and export it to STL file format;

[0059] S43: Use 3D printing slicing software to slice and set parameters according to the printer and material;

[0060] S44: Print and adjust the printing material and support structure as needed to ensure accurate surface details.

[0061] The technical effect of the present invention is as follows: the present invention introduces power spectrum density and fractal theory, uses Fourier transform to simulate the concept of power spectrum density for signal energy decomposition, synthesizes the overall plane space height through the point frequency components of the plane array, and generates surface data of random rock joint surfaces; by modifying the Fourier coefficients, a three-dimensional surface with different anisotropy and mismatch is generated, and the degree of suppression of high-frequency components can be changed by adjusting the value of the fractal dimension H, thereby increasing or decreasing the detailed roughness characteristics of the surface. Roughness characteristic parameters that have a good correlation with the joint roughness value JRC and are not affected by the baseline position are selected: the root mean square first-order derivative value (Z2), the structure function (SF), the roughness profile index (Rp), and the 3D roughness parameter The target joint roughness value is set, slices are sampled, and JRC values ​​are calculated and averaged using multiple roughness parameters. The sum of the absolute values ​​of the differences between the target JRC value and the calculated JRC values ​​of each parameter is taken as the final target. The differential evolution algorithm is introduced, and through a combination of global optimization and local optimization, the optimal solution of the fractal surface parameters that meets the target difference value minimization is obtained. The rock joint surface data that is consistent with the evaluation system of multiple main characteristic parameters of roughness is generated based on power spectral density. Then, the target surface is generated through 3D printing technology, and the repeatable generation of corresponding joint roughness surfaces under multiple roughness parameter evaluation standards of the rock surface is realized, thereby more accurately simulating the natural rock joint surface.

[0062] The following is a further description with reference to the accompanying drawings. BRIEF DESCRIPTION OF THE DRAWINGS

[0063] Figure 1 This is a flow chart of a method for generating rock structure surface roughness based on power spectral density according to an embodiment of the present invention.

[0064] Figure 2 Schematic diagram of data slicing of generated surface data in an embodiment of the present invention.

[0065] Figure 3 Schematic diagram of rough joint data generated in an embodiment of the present invention. DETAILED DESCRIPTION

[0066] Example 1

[0067] like Figure 1 As shown, a method for generating rock structure surface roughness based on power spectrum density includes the following steps:

[0068] A method for generating rock structure surface roughness based on power spectral density comprises the following steps:

[0069] S1: Establish and initialize the Fourier coefficient storage array, generate the Fourier matrix through basic calculation parameters, and generate surface data X using inverse Fourier transform;

[0070] S2: Sampling the obtained surface data, and calculating the joint roughness value JRC by using different roughness parameters;

[0071] S3: Establishing the target joint roughness value optimization function, and iteratively obtaining the optimal parameter solution under multiple characteristic parameter standards through global optimization and local optimization;

[0072] S4: Based on the surface data obtained for different target joint roughness, 3D printing is used to generate a surface with a specific roughness.

[0073] In step S1, the process of generating surface data X by using inverse Fourier transform includes:

[0074] S11: According to the actual sampling interval and frame size, the number of grid points of the initialization Fourier coefficient storage array is defined, and the Fourier matrix A is generated through basic calculation parameters. The basic calculation parameters include N, H, IS, DF and ST, where N represents the array size. For the generated N×N two-dimensional array, it represents the number of grid points on the fractal surface; H is the fractal dimension, which controls the roughness of the surface and takes a value between 1 and 3; IS is the initial seed. The fractal surface generated by the same initial seed can be reproduced. The value of the initial seed will change the generation sequence of random numbers but will not affect the characteristics of the generated plane; DF is the directionality factor, which controls the phase change of the generated fracture surface in different directions. Considering increasing the change of the frequency component of the generated plane, it is usually within the range of 1±1; ST is the scale threshold, which controls the irregularity of the surface at a specific radius. It is applied to special processing of smaller radii to amplify the characteristics of the overall frequency component and is a real number greater than 0;

[0075] S12: Initialize the matrix A, traverse each point in the array, skip special points and assign values;

[0076] S13: Take the real parts of the values ​​corresponding to the special points (N / 2, 0), (0, N / 2), and (N / 2, N / 2) to meet the symmetry requirements of the inverse Fourier transform;

[0077] S14: Perform a two-dimensional inverse Fourier transform on the matrix A, convert the frequency domain data into the spatial domain and take the real part of the transformation result, and finally generate the surface data X of the fracture surface.

[0078] In step S2, the N×N matrix A is initialized to a zero matrix, and the calculation function for the surface data X is as follows:

[0079] A i，j =0

[0080]

[0081] Where A i,j The Fourier coefficients of each frequency point in the initialization Fourier coefficient storage array established in step S1, Rad i，j Represents the Euclidean distance from each matrix element (i, j) to the center point (N / 2, N / 2);

[0082]

[0083] Where, is the generated random phase, Temp is the generated random amplitude, and U(a, b) represents uniform distribution;

[0084] When Rad>Mismatch:A i，j =Temp

[0085] When Rad<Mismatch:

[0086] Where, is a random phase, used to increase the randomness of the surface; e i·Aniso·(i-j) The anisotropy parameter is introduced to make the surface different in different directions; as the radius increases, the amplitude decreases at the power of (3-H). As the fractal dimension H decreases, the details of the generated plane are smoother, that is, the high-frequency components of the components are suppressed;

[0087] A N / 2，0 =Re(A N / 2，0 ), A 0，N / 2 =Re(A 0，N / 2 ), A N / 2，N / 2 =Re(A N / 2，N / 2 )

[0088] Where Re represents the real part;

[0089]

[0090] Where, Represents a two-dimensional inverse Fourier transform. The returned X is a two-dimensional matrix representing the simulated fracture surface, where the element value represents the "height" of the surface.

[0091] In the actual calculation process, the array size is set to a 256*256 point matrix, and the Euclidean distance Rad from the frequency point to the origin in the plane is calculated. i，j , and generate random phase and amplitude, and calculate the Fourier coefficient A corresponding to each frequency point ij The size of the Fourier coefficient is taken as the spatial height value of the point, and the overall elevation data of the plane is finally obtained.

[0092] In step S2, the surface data obtained are sampled and the joint roughness value JRC is calculated using different roughness parameters. The specific process is as follows:

[0093] S21: Randomly select 10 index rows and index columns for the obtained surface data, such as Figure 2 As shown, the horizontal and vertical slicing of data are completed;

[0094] S22: Initialize a list to store the joint roughness value JRC calculated for each slice;

[0095] S23: Calculate the joint roughness values ​​JRC of each transverse and longitudinal slice based on different roughness parameters, and take the average joint roughness value JRC of each slice.

[0096] The roughness parameters include the root mean square first-order derivative value Z2, the structure function SF, the roughness profile index R p and 3D roughness parameters in:

[0097]

[0098] Where Z i is the height value of the i-th point on the surface; Z i+1 is the height of the i+1th point on the surface; X i is the horizontal position of the i-th point; X i+1 is the horizontal position of the i+1th point; Δl is the horizontal distance between two adjacent points, i.e., the sampling interval; M x The number of points used to divide the region, l is the total length.

[0099] The calculation formula for the joint roughness value JRC is calculated by the root mean square first-order derivative value Z2:

[0100] JRC=51.16(Z2) 0.531 -11.44

[0101] The calculation formula for the joint roughness value JRC is calculated by the structure function SF:

[0102] JRC=73.95(SF) 0.266 -11.38

[0103] By roughness profile index R p The calculation formula for the joint roughness value JRC is:

[0104] JRC=65.9(R p -1) 0.302 -9.65

[0105] Through 3D roughness parameters The process of calculating the joint roughness value JRC is:

[0106] S231: First, traverse the triangles formed by connecting three adjacent points in the generated surface, define the observation reference, and calculate the apparent inclination of the triangular surface using the azimuth and true inclination. The specific formula is:

[0107] tanθ * = -tanθ·cosα

[0108] Where θ * , θ, and α are the apparent inclination, true inclination, and azimuth of the surface connecting the three adjacent points on the rough surface, respectively;

[0109] S232: Then, by fitting the angle threshold-normalized area curve, the dimensionless parameter c is obtained. The specific formula is:

[0110]

[0111] In the formula is the normalized area, which is the surface area with an apparent inclination angle greater than the chosen threshold θ normalized to the total surface area; is the maximum apparent inclination of the surface in the selected analysis direction; C is a dimensionless fitting parameter calculated by nonlinear least squares regression analysis;

[0112] S233: Finally, substitute the joint roughness value function relationship Calculate the corresponding joint roughness value JRC value.

[0113] In step S3, the optimal parameter solution under various characteristic parameter standards is obtained through iteration through global optimization and local optimization. The specific process is as follows:

[0114] S31: define global parameter optimization boundaries;

[0115] S32: introducing the calculation formula for calculating the joint roughness value of different roughness parameters in step S2, defining the target joint roughness value, and calculating different roughness parameters and corresponding actual joint roughness values ​​for the surface data generated by simulation;

[0116] S33: Using the differential evolution algorithm for global optimization, iteratively searching for the optimal parameter solution with the minimum distance between the actual joint roughness value calculated by each roughness parameter and the target joint roughness value;

[0117] S34: The optimal parameters obtained from global optimization are used as the initial guess values ​​for local optimization, new optimization boundaries and iteration times are set, and the search is refined to obtain the parameter combination closest to the target joint roughness value, and then the optimal parameter solution is obtained, such as Figure 3 shown.

[0118] The parameter combinations corresponding to different target joint roughness values ​​obtained by preliminary optimization calculation are shown in Table 1:

[0119] Table 1: Preliminary optimization parameter solutions for some target joint roughness values

[0120]

[0121] In step S4, the surface data obtained for different target joint roughness are used to generate a surface with a specific roughness using 3D printing. The specific process is as follows:

[0122] S41: normalizing the surface data X to ensure that the actual size of the 3D printing meets the requirements;

[0123] S42: Convert the 2D height matrix into a 3D grid pattern and export it to STL file format;

[0124] S43: Use 3D printing slicing software to slice and set parameters according to the printer and material;

[0125] S44: Print and adjust the printing material and support structure as needed to ensure the accuracy of surface details. The above description is only a preferred embodiment of the present invention, but the scope of protection of the present invention is not limited to this. Any changes or substitutions that can be easily conceived by a person skilled in the art within the technical scope disclosed in the present invention should be covered by the scope of protection of the present invention.

Claims

1. A method for generating rock structure surface roughness based on power spectral density, characterized by: The following steps are involved: S1: Establish and initialize the Fourier coefficient storage array, generate the Fourier matrix through basic calculation parameters, and generate surface data X using inverse Fourier transform; S2: Sampling the obtained surface data, and calculating the joint roughness value JRC by using different roughness parameters; S3: Establishing the target joint roughness value optimization function, and iteratively obtaining the optimal parameter solution under multiple characteristic parameter standards through global optimization and local optimization; S4: Based on the surface data obtained for different target joint roughness, 3D printing is used to generate a surface with a specific roughness.

2. The method for generating rock structure surface roughness based on power spectrum density according to claim 1, characterized in that: In step S1, the process of generating surface data X by using inverse Fourier transform includes: S11: According to the actual sampling interval and frame size, the number of grid points of the initialization Fourier coefficient storage array is defined, and the Fourier matrix A is generated through basic calculation parameters. The basic calculation parameters include N, H, IS, DF and ST, where N represents the array size. For the generated N×N two-dimensional array, it represents the number of grid points on the fractal surface; H is the fractal dimension, which controls the roughness of the surface and takes a value between 1 and 3; IS is the initial seed. The fractal surface generated by the same initial seed can be reproduced. The value of the initial seed will change the generation sequence of random numbers but will not affect the characteristics of the generated plane; DF is the directionality factor, which controls the phase change of the generated fracture surface in different directions. Considering increasing the change of the frequency component of the generated plane, it is usually within the range of 1±1; ST is the scale threshold, which controls the irregularity of the surface at a specific radius. It is applied to special processing of smaller radii to amplify the characteristics of the overall frequency component and is a real number greater than 0; S12: Initialize the matrix A, traverse each point in the array, skip special points and assign values; S13: Take the real parts of the values ​​corresponding to the special points (N / 2, 0), (0, N / 2), and (N / 2, N / 2) to meet the symmetry requirements of the inverse Fourier transform; S14: Perform a two-dimensional inverse Fourier transform on the matrix A, convert the frequency domain data into the spatial domain and take the real part of the transformation result, and finally generate the surface data X of the fracture surface.

3. The method for generating rock structure surface roughness based on power spectrum density according to claim 2, characterized in that: In step S2, the N×N matrix A is initialized to a zero matrix, and the calculation function for the surface data X is as follows: A i，j ＝0 Where A i,j The Fourier coefficients of each frequency point in the initialization Fourier coefficient storage array established in step S1, Rad i，j Represents the Euclidean distance from each matrix element (i, j) to the center point (N / 2, N / 2); Where, is the generated random phase, Temp is the generated random amplitude, and U(a, b) represents uniform distribution; When Rad > Mismatch: A i，j = Temp when Where, is a random phase, used to increase the randomness of the surface; e i·Aniso·(i-j) The anisotropy parameter is introduced to make the surface different in different directions; as the radius increases, the amplitude decreases at the power of (3-H). As the fractal dimension H decreases, the details of the generated plane are smoother, that is, the high-frequency components of the components are suppressed; A N / 2,0 =Re(A N / 2,0 ), A 0，N / 2 =Re(A 0，N / 2 ), A N / 2，N / 2 =Re(A N / 2，N / 2 ) Where Re represents the real part; Where, Represents a two-dimensional inverse Fourier transform. The returned X is a two-dimensional matrix representing the simulated fracture surface, where the element value represents the "height" of the surface.

4. The method for generating rock structure surface roughness based on power spectrum density according to claim 1, characterized in that: In step S2, the surface data obtained are sampled and the joint roughness value JRC is calculated using different roughness parameters. The specific process is as follows: S21: Randomly select 10 index rows and index columns from the obtained surface data to complete horizontal and vertical slicing of the data; S22: Initialize a list to store the joint roughness value JRC calculated for each slice; S23: Calculate the joint roughness values ​​JRC of each transverse and longitudinal slice based on different roughness parameters, and take the average joint roughness value JRC of each slice.

5. The method for generating rock structure surface roughness based on power spectrum density according to claim 4, characterized in that: The roughness parameters include the root mean square first-order derivative value Z2, the structure function SF, the roughness profile index R p and 3D roughness parameters in: Where Z i is the height of the i-th point on the surface; Z i+1 is the height of the i+1th point on the surface; X i is the horizontal position of the i-th point; X i+1 is the horizontal position of the i+1th point; Δl is the horizontal distance between two adjacent points, i.e., the sampling interval; M x The number of points used to divide the region, l is the total length.

6. The method for generating rock structure surface roughness based on power spectrum density according to claim 5, characterized in that: The calculation formula for the joint roughness value JRC is calculated by the root mean square first-order derivative value Z2: JRC=51.16(Z2) 0.531 -11.44 The calculation formula for the joint roughness value JRC is calculated by the structure function SF: JRC=73.95(SF) 0.266 -11.38 By roughness profile index R p The calculation formula for the joint roughness value JRC is: JRC=65.9(R p -1) 0.302 -9.65 Through 3D roughness parameters The process of calculating the joint roughness value JRC is: S231: First, traverse the triangles formed by connecting three adjacent points in the generated surface, define the observation reference, and calculate the apparent inclination of the triangular surface using the azimuth and true inclination. The specific formula is: tanθ * =-tanθ·cosα Where θ * , θ, and α are the apparent inclination, true inclination, and azimuth of the surface connecting the three adjacent points on the rough surface, respectively; S232: Then, by fitting the angle threshold-normalized area curve, the dimensionless parameter c is obtained. The specific formula is: In the formula is the normalized area, which is the surface area with an apparent inclination angle greater than the chosen threshold θ normalized to the total surface area; is the maximum apparent inclination of the surface in the selected analysis direction; C is a dimensionless fitting parameter calculated by nonlinear least squares regression analysis; S233: Finally, substitute the joint roughness value function relationship Calculate the corresponding joint roughness value JRC value.

7. The method for generating rock structure surface roughness based on power spectrum density according to claim 6, characterized in that: In step S3, the optimal parameter solution under various characteristic parameter standards is obtained through iteration through global optimization and local optimization. The specific process is as follows: S31: define global parameter optimization boundaries; S32: introducing the calculation formula for calculating the joint roughness value of different roughness parameters in step S2, defining the target joint roughness value, and calculating different roughness parameters and corresponding actual joint roughness values ​​for the surface data generated by simulation; S33: Using the differential evolution algorithm for global optimization, iteratively searching for the optimal parameter solution with the minimum distance between the actual joint roughness value calculated by each roughness parameter and the target joint roughness value; S34: The optimal parameters obtained from the global optimization are used as the initial guess values ​​for the local optimization, a new optimization boundary and number of iterations are set, and the search is refined to obtain the parameter combination closest to the target joint roughness value, thereby obtaining the optimal parameter solution.

8. The method for generating rock structure surface roughness based on power spectrum density according to claim 7, characterized in that: In step S4, the surface data obtained for different target joint roughness are used to generate a surface with a specific roughness using 3D printing. The specific process is as follows: S41: normalizing the surface data X to ensure that the actual size of the 3D print meets the requirements; S42: Convert the 2D height matrix into a 3D grid pattern and export it to STL file format; S43: Use 3D printing slicing software to slice and set parameters according to the printer and material; S44: Print and adjust the printing material and support structure as needed to ensure accurate surface details.

Images