Rock surface erosion volume calculation method based on curved surface fitting optimization

Abstract

The invention provides a rock surface erosion volume calculation method based on curved surface fitting optimization. The rock surface erosion volume calculation method comprises the following steps: fully grinding and polishing the surface of a rock sample; performing three-dimensional laser scanning on the surface of the rock sample by using a laser confocal microscope to obtain three-dimensional point cloud data; carrying out model construction in the Matlab script; determining a reference surface of the surface of the rock sample before erosion treatment; grid division is carried out on the three-dimensional point cloud data according to the initial block size, and quadratic polynomial curved surface fitting is carried out on each block; calculating and accumulating the recess volumes of all the blocks to obtain the total recess volume before erosion; the total recess volume after erosion is obtained according to the steps, the difference value of the total recess volume before and after erosion is calculated, and the erosion volume is obtained. The method solves the problems that in the prior art, cross-scale data fusion errors are remarkable, and the volume calculation reliability is affected.

Description

Technical Field

[0001] The present invention relates to the technical field related to rock surface erosion volume calculation, and in particular to a rock surface erosion volume calculation method based on surface fitting optimization. Background Art

[0002] Driven by the global "dual carbon" strategy, CO2 geological storage technology has become a key means of achieving carbon neutrality. However, the complex water-rock-gas interactions triggered by CO2 injection into reservoirs can lead to dissolution and reprecipitation of rock minerals, pore structure reconstruction, and mechanical property degradation, seriously threatening the integrity of the storage system. For example, experiments have demonstrated that CO2-saturated water can reduce the elastic modulus of sandstone by up to 30% and induce the expansion of microfracture networks. Geochemical simulations show that caprock mineral dissolution can increase permeability by two orders of magnitude. These physicochemical processes urgently require the development of high-precision methods for quantifying rock surface erosion.

[0003] Although recent studies have revealed the advantages of surface fitting technology in materials science, the following technical gaps still exist in the field of rock erosion quantification: the existing surface fitting algorithms are not optimized for the multi-scale morphology of rock surfaces (millimeter to nanometer scale), resulting in significant errors in cross-scale data fusion.

[0004] Prior art discloses a Chinese patent with publication number CN 117969545B, which describes a method for calculating rock surface erosion volume based on 3D point cloud data. The method includes the following steps: polishing the rock sample surface to ensure consistent surface roughness before erosion; performing 3D laser scanning on the rock surface before and after erosion to obtain 3D point cloud data; using a Matlab script to create a visual model of the rock sample's true surface topography based on the 3D point cloud data; defining the concave volume and determining the reference planes before and after erosion in the constructed rock surface model; and calculating the concave volume and erosion volume of the rock surface. However, this method for calculating rock surface erosion volume (i.e., the line-surface method) relies on a fixed grid, making it less adaptable to steep surfaces and small pits. It also lacks a dynamic reference plane calibration mechanism, making it difficult to distinguish between eroded and non-eroded areas, impacting the reliability of the volume calculation. Summary of the Invention

[0005] The present invention aims to provide a rock surface erosion volume calculation method based on surface fitting optimization, so as to solve the above-mentioned problem of significant cross-scale data fusion errors and impact on the reliability of volume calculation.

[0006] To this end, the technical solution adopted by the present invention is: a method for calculating rock surface erosion volume based on surface fitting optimization, comprising the following steps:

[0007] S1. The surface of the rock sample is fully polished to ensure that the roughness of the rock sample surface is consistent before treatment;

[0008] S2. Performing three-dimensional laser scanning on the surface of the rock sample using a laser confocal microscope to obtain three-dimensional point cloud data of the rock sample surface having a preset resolution, and performing noise reduction processing on the three-dimensional point cloud data according to actual conditions, wherein the preset resolution is dynamically adjusted according to the surface morphology of the rock sample and the performance of the equipment;

[0009] S3, using the coordinates (X i , Y j , Z i,j ) constructing a model in a Matlab script to obtain a visualization model of the actual surface morphology of the rock sample;

[0010] S4. determining a reference plane of the surface of the rock sample before erosion treatment according to the constructed visualization model;

[0011] S5. Meshing the 3D point cloud data according to the initial block size, performing quadratic polynomial surface fitting on each block, and dynamically adjusting the block size based on the goodness of fit R² value; the blocks support multi-level recursive subdivision from the initial block size to a minimum block size of 2×2 until the goodness of fit within the block satisfies R² ≥ 0.99 or the minimum block size is reached;

[0012] S6. For the blocks that satisfy the goodness of fit R²≥0.99, calculate the concave volume between the fitting surface and the reference surface by double integral; for the smallest block that does not satisfy the goodness of fit R²≥0.99, calculate the concave volume by linear interpolation; accumulate the concave volumes of all the blocks to obtain the total concave volume before erosion;

[0013] S7. Repeat steps S2-S3 for the rock sample after erosion treatment to obtain the visualization model, then determine the reference plane of the surface of the rock sample after erosion treatment based on the visualization model, and repeat steps S5-S6 to obtain the total depressed volume after erosion; the eroded volume of the rock sample surface can be calculated by subtracting the total depressed volume before erosion from the total depressed volume after erosion.

[0014] As a preferred embodiment of the above solution, in step S4, the average Z coordinate of each point in the three-dimensional point cloud data of the rock sample surface before erosion treatment is used as the reference plane height, so as to determine the reference plane of the rock sample surface before erosion treatment. The calculation formula is:

[0015]

[0016] Where zm0 is the base height of the rock sample surface before erosion treatment, n is the number of point clouds, z i is the z coordinate of the i-th point.

[0017] More preferably, in step S5, the initial block size is the original resolution size of the point cloud, and the quadratic polynomial surface fitting model adopts a standardized coordinate form, and its general formula is:

[0018]

[0019] Where a, b, c, d, e, and f are fitting coefficients, and are all solved by the least squares method. X, Y, and Z are the coordinate data of the point cloud.

[0020] More preferably, in step S5, the calculation formula of the goodness of fit of the block is as follows:

[0021]

[0022] Where R 2 is the goodness of fit of the block, Z obs is the observed value, Z fit is the fitted value.

[0023] Further preferably, in order to eliminate the dimension difference and ensure the stability of the fitting result, the region coordinates of each block are standardized, and the calculation formula is:

[0024]

[0025] Where u and v are the coordinates of the region after block standardization. 、 is the mean value of the coordinates within the block, 、 is the standard deviation, x i 、y j are the regional coordinates before block normalization processing;

[0026] The polynomial surface equation of each block is obtained by surface fitting, and its calculation formula is:

[0027]

[0028] Where f(u,v) is the independent variable of the binary function, c1, c2, c3, c4, c5, and c6 are the coefficients of the fitting equation, and u and v are the regional coordinates after block standardization.

[0029] More preferably, the concave volume V between the fitting surface and the reference surface in the block is calculated by the double integral block, the integration area D is the actual physical range of the current block, such as , the integral result is realized by Matlab's integratePolynomial function, and its calculation formula is:

[0030]

[0031] Where V block is the concave volume calculated by double integral, f(u,v) is the independent variable of the two-variable function, z m is the height of the reference plane, D is the integration area, and x and y are the actual physical ranges of the blocks respectively;

[0032] The calculation formula for calculating the concave volume using the linear interpolation method is as follows:

[0033]

[0034] Where V linear is the concave volume calculated by linear interpolation, is the average height of the point cloud in the block, and A is the area of ​​the block.

[0035] More preferably, the total concave volume is the sum of the concave volumes of each block, and its calculation formula is:

[0036]

[0037] Where V total is the total concave volume, V block is the concave volume calculated by double integral, V linear The volume of the depression calculated by linear interpolation.

[0038] Further preferably, in step S5, when the block is subdivided, if the goodness of fit R² of the block is less than 0.99, it is divided into 4 equal sub-blocks and the surface fitting is performed again.

[0039] Further preferably, in step S6, the calculation results of the concave volumes of all the blocks are globally optimized by the residual sum of squares to ensure that the overall goodness of fit R² is ≥ 0.99.

[0040] More preferably, in step S7, when determining the reference surface of the rock sample after erosion, it is necessary to randomly select specific areas with a grid size of p×q in the three-dimensional point cloud data after erosion, and calculate the root mean square roughness S of these areas. q , until a certain number t of root mean square roughness S are selected q Both are smaller than S q0 The area, S q0is the root mean square roughness of the rock sample surface before erosion treatment, and the arithmetic average of the Z coordinates of the points in these areas is used as the reference plane height z of the rock sample surface after erosion treatment. m1 , thereby determining the reference surface of the rock sample after erosion treatment, the root mean square roughness S q The calculation formula is as follows:

[0041]

[0042] Where S q is the root mean square roughness of the rock sample surface after erosion treatment, m and n are the number of point cloud rows and columns respectively, and Z m is the average height of the z coordinates of all 3D point clouds, Z i ,j is the z-coordinate of the point in the i-th row and j-th column.

[0043] Beneficial effects of the present invention:

[0044] 1. By integrating the block adaptive strategy, seamless connection of data at different resolution levels is achieved. Experiments show that this method reduces the error in cross-scale data volume calculation from 6% of traditional methods to 0.5%, thus reducing the error in cross-scale data fusion.

[0045] 2. Compared with the existing rock surface erosion volume calculation method based on three-dimensional point cloud data, the improved MLTS-LS algorithm significantly improves the accuracy and efficiency of rock surface erosion volume calculation, thereby improving the reliability of rock surface erosion volume calculation, more accurately quantifying the extent of rock surface erosion reaction during CO2 geological storage, and better restoring the actual situation of the rock surface.

[0046] 3. A mineral-morphology coupling model is used to establish a quantitative relationship between erosion volume and storage safety. This method breaks through the accuracy ceiling of existing technologies and provides key technical support for CO2 storage site screening and risk warning. BRIEF DESCRIPTION OF THE DRAWINGS

[0047] Figure 1 This is a visualization model diagram of the reconstruction of three-dimensional point cloud data through Matlab script in the present invention.

[0048] Figure 2 It is a two-dimensional schematic diagram of the rock surface morphology before and after treatment in the present invention.

[0049] Figure 3 This is a flow chart for determining the reference plane in the present invention.

[0050] Figure 4 This is a three-dimensional reconstruction model diagram of the rock surface before and after erosion after the reference plane is determined in the present invention.

[0051] Figure 5Schematic diagram of a local area of ​​a rough surface before and after flipping in the present invention.

[0052] Figure 6 This is a flow chart for calculating the volume of three-dimensional curved depressions on the rock surface in the present invention.

[0053] Figure 7 This is a distribution diagram of the block fitting size in the present invention.

[0054] Figure 8 Schematic diagram of local surface fitting in the present invention.

[0055] Figure 9 It is the distribution diagram of different fitting methods in the present invention.

[0056] Figure 10 Schematic diagram of the coal rock sample and the device for simulating erosion in the present invention.

[0057] Figure 11 This is a three-dimensional reconstruction model diagram of the coal rock sample before and after erosion in the present invention. DETAILED DESCRIPTION

[0058] The present invention will be further described below with reference to the accompanying drawings and embodiments.

[0059] like Figure 1-11 As shown, a method for calculating rock surface erosion volume based on surface fitting optimization includes the following steps:

[0060] S1. Thoroughly grind and polish the surface of the rock sample to ensure that the roughness of the rock sample surface remains consistent before processing.

[0061] This step lays the foundation for the subsequent acquisition of accurate three-dimensional point cloud data and the establishment of a reliable surface model, reducing the interference of large differences in surface roughness on the calculation results.

[0062] S2. Use a laser confocal microscope to perform three-dimensional laser scanning on the surface of the rock sample to obtain three-dimensional point cloud data with a preset resolution on the surface of the rock sample. Perform noise reduction on the three-dimensional point cloud data according to actual conditions. The preset resolution is dynamically adjusted according to the surface morphology characteristics of the rock sample and the performance of the equipment.

[0063] 3D laser scanning technology can quickly and accurately acquire spatial information about rock surfaces, providing rich data support for subsequent model construction and volume calculations. Noise reduction processing can remove noise interference from the data, improving data quality and making subsequent model construction and calculations more accurate.

[0064] S3, using the coordinates of each point in the three-dimensional point cloud data (X i , Y j , Z i,j) The model was constructed in Matlab script to obtain a visual model of the real surface morphology characteristics of the rock sample.

[0065] This visualization model can intuitively display the morphology of the rock surface, facilitating subsequent analysis and processing of rock surface characteristics.

[0066] S4. Determine the reference surface of the rock sample before erosion treatment based on the constructed visualization model.

[0067] After constructing a visual model of the rock surface's 3D point cloud data, the volume of the depression cannot be directly calculated. Instead, a reference plane is required to divide the constructed model into two parts, the upper and lower parts. The volume of the irregular shape formed by the part of the model below the reference plane and the part between the reference plane is called the rock's depression volume. To facilitate the calculation of the depression volume, the visualization model of the rock's actual surface topography is flipped upside down with the reference plane as the symmetry plane. After flipping, the volume of the space between the curved surface above the reference plane and the reference plane is the depression volume.

[0068] In step S4, the average Z coordinate of each point in the three-dimensional point cloud data is used as the reference plane height for the rock sample surface before erosion treatment, thereby determining the reference plane of the rock sample surface before erosion treatment. The calculation formula is:

[0069]

[0070] Where z m0 is the base height of the rock sample surface before erosion treatment, n is the number of point clouds, z i is the z coordinate of the i-th point.

[0071] S5. Grid the 3D point cloud data according to the initial block size. Perform a quadratic polynomial surface fit on each block, and dynamically adjust the block size based on the goodness-of-fit R² value. Blocks are recursively subdivided from the initial block size to a minimum block size of 2×2 until the goodness-of-fit within the block meets R² ≥ 0.99 or the minimum block size is reached.

[0072] In step S5, the initial block size is the original resolution size of the point cloud, and the quadratic polynomial surface fitting model adopts the standardized coordinate form, and its general formula is:

[0073]

[0074] Where a, b, c, d, e, and f are fitting coefficients, and are all solved by the least squares method. X, Y, and Z are the coordinate data of the point cloud.

[0075] In step S5, the calculation formula for the goodness of fit of the blocks is as follows:

[0076]

[0077] Where R 2 is the goodness of fit of the block, Z obs is the observed value, Z fit is the fitted value.

[0078] In step S5, when the block is subdivided, if the goodness of fit R² of the block is less than 0.99, it is divided into 4 sub-blocks and the surface fitting is performed again.

[0079] The 3D point cloud data matrix is ​​calculated and processed using a Matlab script, and different blocks are fitted one by one. Its dynamic block strategy includes defining the physical space coordinate range (such as 0-256 um) and mapping the point cloud data to the actual size to ensure that the data is consistent with the actual physical size. The initial block size is set to the original resolution size, such as 1024×1024, and the minimum block size is 2×2. The entire point cloud data is divided into initial blocks, and each block is added to the processing queue. The block queue adopts a last-in-first-out (LIFO) management mechanism to ensure computational efficiency. For each block, its fitting accuracy R is calculated. 2 , if R 2 If the value is less than the threshold (e.g. 0.99), the current block is divided into 4 equal sub-blocks and the sub-blocks are added to the queue. The recursive division process continues until the block size reaches the minimum value or the fitting accuracy meets the requirements.

[0080] S6. For blocks that meet the goodness-of-fit R² ≥ 0.99, calculate the concave volume between the fitted surface and the reference surface using double integration. For the smallest block that does not meet the goodness-of-fit R² ≥ 0.99, calculate the concave volume using linear interpolation. Accumulate the concave volumes of all blocks to obtain the total concave volume before erosion.

[0081] In step S6, the calculation results of the concave volumes of all blocks are globally optimized by the residual sum of squares to ensure that the overall goodness of fit R² ≥ 0.99.

[0082] Based on Matlab's polynomial surface fitting technology, a quadratic polynomial is used to fit the surface of each block area. In order to eliminate dimensional differences and ensure the stability of the fitting results, the coordinates of each block area are standardized. The calculation formula is:

[0083]

[0084] Where u and v are the coordinates of the region after block standardization. 、 is the mean value of the coordinates within the block, 、 is the standard deviation, xi 、y j are the region coordinates before block normalization.

[0085] The polynomial surface equation of each block is obtained by surface fitting, and its calculation formula is:

[0086]

[0087] Where f(u,v) is the independent variable of the binary function, c1, c2, c3, c4, c5, and c6 are the coefficients of the fitting equation, and u and v are the regional coordinates after block standardization.

[0088] The coefficient matrix is ​​usually determined by the least squares method, and the fitting results are used to calculate the fitting residuals, which are calculated as follows:

[0089]

[0090] Where X is the normalized coordinate matrix, Z is the height value vector, T is the transpose of the coefficient matrix, and P is the vector of the coefficient matrix.

[0091] Calculate the concave volume V between the fitted surface and the reference surface in the block by double integral block , the integration area D is the actual physical range of the current block, such as , the integral result is realized by Matlab's integratePolynomial function, and its calculation formula is:

[0092]

[0093] Where V block is the concave volume calculated by double integral, f(u,v) is the independent variable of the two-variable function, z m is the reference height, D is the integration area, and x and y are the actual physical ranges of the blocks.

[0094] The formula for calculating the concave volume using linear interpolation is as follows:

[0095]

[0096] Where V linear is the concave volume calculated by linear interpolation, is the average height of the point cloud in the block, and A is the area of ​​the block.

[0097] The total depression volume is the sum of the depression volumes of each block, and its calculation formula is:

[0098]

[0099] Where Vtotal is the total concave volume, V block is the concave volume calculated by double integral, V linear The volume of the depression calculated by linear interpolation.

[0100] By calculating the depression volume on the rock surface before and after the reaction, the difference between the two can be used to obtain the erosion volume. The erosion distribution intuitively reflects the degree of reaction of different mineral areas on the rock surface and also provides a basis for the change in coal rock porosity.

[0101] S7. Repeat steps S2-S3 above for the eroded rock sample to obtain a visualization model. The reference plane of the eroded rock sample surface is then determined based on the visualization model. Steps S5-S6 above are then repeated to obtain the total depression volume after erosion. The eroded volume of the rock sample surface can be calculated by subtracting the total depression volume before erosion from the total depression volume after erosion.

[0102] After erosion, some minerals on the rock surface are affected by physical and chemical reactions, resulting in a certain degree of erosion and deformation, and changes in the surface morphology. Some minerals are unaffected, and the surface morphology remains unchanged. Therefore, the post-erosion reference level can be determined based on the portion of the rock surface that has not been affected by erosion.

[0103] In step S7, when determining the reference surface of the rock sample after erosion treatment, it is necessary to randomly select specific areas with a grid size of p×q in the three-dimensional point cloud data after erosion treatment and calculate the root mean square roughness S of these areas. q , until a certain number t of root mean square roughness S are selected q Both are smaller than S q0 The area, S q0 is the root mean square roughness of the rock sample surface before erosion treatment, and the arithmetic average of the Z coordinates of the points in these areas is used as the reference plane height z of the rock sample surface after erosion treatment. m1 , thereby determining the reference surface of the rock sample after erosion treatment, the root mean square roughness S q The calculation formula is as follows:

[0104]

[0105] Where S q is the root mean square roughness of the rock sample surface after erosion treatment, m and n are the number of point cloud rows and columns respectively, and Z m is the average height of the z coordinates of all 3D point clouds, Z i ,j is the z-coordinate of the point in the i-th row and j-th column.

[0106] Root mean square roughness S of the rock sample surface before erosion treatment q0 The root mean square roughness S of the rock sample surface after erosion treatmentq The calculation formula is the same as that of . The selection of P, q, and t needs to be determined based on the mineral size of the reaction area on the rock surface.

[0107] Coal rock samples were taken from the Zhengtong Coal Mine in Changwu County, Xianyang City, Shaanxi Province. After removing the weathered layer, square specimens with a size of 2 cm × 2 cm × 2 cm were obtained by wire cutting. The surface of the coal rock was polished with 600, 1200, 2000, and 3000 mesh sandpaper respectively. The test samples are shown in Figure 10 (a).

[0108] A random area of ​​the polished coal rock surface before erosion was scanned, and three random areas of the coal rock surface after erosion were scanned. After obtaining the three-dimensional point cloud data of the area, the three-dimensional reconstruction of the coal rock surface morphology before and after erosion was performed using the Matlab program. Figure 11 shown.

[0109] The 3D point cloud data is obtained by using the LEXTOLS4000 laser confocal microscope 3D laser scanning equipment produced by Olympus. Figure 10 (b)

[0110] In order to simulate the interaction process of CO2-water-coal-rock in real reservoirs, a high temperature and high pressure immersion test device was used, such as Figure 10 As shown in (c), based on the target reservoir depth, the test temperature was set to 40°C and the CO2 pressure was set to 8 MPa to simulate the supercritical CO2-water-coal-rock erosion effect in the real reservoir, and the reaction time was set to 10 days.

[0111] According to the above steps, the three-dimensional point cloud data of the coal rock surface before and after erosion treatment were processed to obtain the base surface height, root mean square roughness, depression volume and erosion volume of each area, as shown in Table 1.

[0112] Table 1

[0113] Scan area Reference surface height z / μm <![CDATA[Root mean square roughness S q > <![CDATA[Sunken volume V / μm 3 > <![CDATA[Erosion volume V e / V / μm 3 > <![CDATA[A0]]> 1.7060 0.1254 <![CDATA[3.2223×10 3 ]]> / <![CDATA[A1]]> 11.1643 0.5023 <![CDATA[1.1173×10 4 ]]> <![CDATA[7.9507×10 3 ]]> <![CDATA[A2]]> 11.4238 0.8496 <![CDATA[2.1853×10 4 ]]> <![CDATA[1.8631×10 4 ]]> <![CDATA[A3]]> 16.3953 1.7172 <![CDATA[4.2743×10 4 ]]> <![CDATA[3.9521×10 4 ]]>

[0114] The comparison between the traditional three-dimensional point cloud line-surface method and the surface fitting optimization method of the present invention is shown in Table 2.

[0115] Table 2

[0116] Comparison Dimension Traditional 3D point cloud line and surface method The present invention (surface fitting optimization method) Volume calculation method Fixed grid line-surface segmentation method Adaptive block surface fitting method (dynamic grid + polynomial fitting) Precision control Depends on grid uniformity, complex topography is prone to cumulative errors R²≥0.99 guarantees local fitting accuracy Adaptability to complex shapes Insufficient sensitivity to steep surfaces / small pits Block recursive refinement + residual optimization, supporting sub-millimeter accuracy Computational efficiency Uniform block calculation amount is fixed Dynamic block balance calculation accuracy and efficiency Determination of datum Single RMS roughness screening Multi-area screening + dynamic block reference surface calibration

[0117] Compared with the SEM rock surface morphology observation and quantitative characterization method that uses surface roughness and fractal number to measure the degree of erosion, the erosion volume proposed in the present invention has stronger reliability. The surface fitting-based depression volume calculation method proposed in the present invention is not only applicable to rock samples, but can also be applied to the calculation of erosion volume on the surfaces of various materials after erosion reactions according to actual needs. The quantitative relationship between erosion volume and storage safety is: by laser scanning to identify different mineral areas, and then obtain point cloud data of different minerals, the surface of different minerals can be imaged before and after the reaction, and the erosion volume of different mineral-filled areas can be obtained. The erosion volume of different mineral-filled areas can reflect the degree of erosion they have suffered, and different degrees of erosion can correspond to the safety of storage.

[0118] While embodiments of the present invention have been shown and described, it will be appreciated by those skilled in the art that various changes, modifications, substitutions, and variations may be made to the embodiments without departing from the principles and spirit of the invention, and that the scope of the invention is defined by the claims and their equivalents.

Claims

1. A method for calculating rock surface erosion volume based on surface fitting optimization, characterized in that: The following steps are involved: S1. The surface of the rock sample is fully polished to ensure that the roughness of the rock sample surface is consistent before treatment; S2. Performing three-dimensional laser scanning on the surface of the rock sample using a laser confocal microscope to obtain three-dimensional point cloud data of the rock sample surface having a preset resolution, and performing noise reduction processing on the three-dimensional point cloud data according to actual conditions, wherein the preset resolution is dynamically adjusted according to the surface morphology of the rock sample and the performance of the equipment; S3, using the coordinates (X i , Y j , Z i,j ) constructing a model in a Matlab script to obtain a visualization model of the actual surface morphology of the rock sample; S4. determining a reference plane of the surface of the rock sample before erosion treatment according to the constructed visualization model; S5. Divide the three-dimensional point cloud data into grids according to the initial block size, perform quadratic polynomial surface fitting on each block, and dynamically adjust the block granularity according to the goodness of fit R² value; The block supports multi-level recursive subdivision from the initial block size to the minimum block size of 2×2, until the goodness of fit within the block meets R²≥0.99 or the minimum block size is reached; S6. For the blocks that satisfy the goodness of fit R²≥0.99, calculate the concave volume between the fitting surface and the reference surface by double integral; for the smallest block that does not satisfy the goodness of fit R²≥0.99, calculate the concave volume by linear interpolation; accumulate the concave volumes of all the blocks to obtain the total concave volume before erosion; S7. Repeat steps S2-S3 for the rock sample after erosion treatment to obtain the visualization model, then determine the reference plane of the surface of the rock sample after erosion treatment based on the visualization model, and repeat steps S5-S6 to obtain the total depressed volume after erosion; the eroded volume of the rock sample surface can be calculated by subtracting the total depressed volume before erosion from the total depressed volume after erosion.

2. The method for calculating rock surface erosion volume based on surface fitting optimization according to claim 1, characterized in that: In step S4, the average Z coordinate of each point in the three-dimensional point cloud data of the rock sample surface before erosion treatment is used as the reference plane height to determine the reference plane of the rock sample surface before erosion treatment. The calculation formula is: ， Where z m0 is the base height of the rock sample surface before erosion treatment, n is the number of point clouds, z i is the z coordinate of the i-th point.

3. The method for calculating rock surface erosion volume based on surface fitting optimization according to claim 1, characterized in that: In step S5, the initial block size is the original resolution size of the point cloud, and the quadratic polynomial surface fitting model adopts a standardized coordinate form, and its general formula is: ， Where a, b, c, d, e, and f are fitting coefficients, and are all solved by the least squares method. X, Y, and Z are the coordinate data of the point cloud.

4. The method for calculating rock surface erosion volume based on surface fitting optimization according to claim 1, characterized in that: In step S5, the calculation formula for the goodness of fit of the block is as follows: ， Where R 2 is the goodness of fit of the block, Z obs is the observed value, Z fit is the fitted value.

5. The method for calculating rock surface erosion volume based on surface fitting optimization according to claim 4, characterized in that: In order to eliminate the dimension difference and ensure the stability of the fitting results, the regional coordinates of each block are standardized, and the calculation formula is: ， Where u and v are the coordinates of the region after block standardization. 、 is the mean value of the coordinates within the block, 、 is the standard deviation, x i 、y j are the regional coordinates before block normalization processing; The polynomial surface equation of each block is obtained by surface fitting, and its calculation formula is: ， Where f(u,v) is the independent variable of the binary function, c1, c2, c3, c4, c5, and c6 are the coefficients of the fitting equation, and u and v are the regional coordinates after block standardization.

6. The method for calculating rock surface erosion volume based on surface fitting optimization according to claim 5, characterized in that: The concave volume V between the fitting surface and the reference surface in the block is calculated by the double integral block , the integration area D is the actual physical range of the current block, such as , the integral result is realized by Matlab's integratePolynomial function, and its calculation formula is: ， Where V block is the concave volume calculated by double integral, f(u,v) is the independent variable of the two-variable function, z m is the height of the reference plane, D is the integration area, and x and y are the actual physical ranges of the blocks respectively; The calculation formula for calculating the concave volume using the linear interpolation method is as follows: ， Where V linear is the concave volume calculated by linear interpolation, is the average height of the point cloud in the block, and A is the area of ​​the block.

7. The method for calculating rock surface erosion volume based on surface fitting optimization according to claim 6, characterized in that: The total concave volume is the sum of the concave volumes of each block, and its calculation formula is: ， Where V total is the total concave volume, V block is the concave volume calculated by double integral, V linear The volume of the depression calculated by linear interpolation.

8. The method for calculating rock surface erosion volume based on surface fitting optimization according to claim 1, characterized in that: In step S5, when the block is subdivided, if the goodness of fit R² of the block is less than 0.99, it is divided into 4 equal sub-blocks and the surface fitting is performed again.

9. The method for calculating rock surface erosion volume based on surface fitting optimization according to claim 1, characterized in that: In step S6, the calculation results of the concave volumes of all the blocks are globally optimized by the residual sum of squares to ensure that the overall goodness of fit R² is ≥ 0.

99.

10. The method for calculating rock surface erosion volume based on surface fitting optimization according to claim 1, characterized in that: In step S7, when determining the reference surface of the rock sample after erosion, it is necessary to randomly select specific areas with a grid size of p×q in the three-dimensional point cloud data after erosion, and calculate the root mean square roughness S of these areas. q , until a certain number t of root mean square roughness S are selected q Both are smaller than S q0 The area, S q0 is the root mean square roughness of the rock sample surface before erosion treatment, and the arithmetic average of the Z coordinates of the points in these areas is used as the reference plane height z of the rock sample surface after erosion treatment. m1 , thereby determining the reference surface of the rock sample after erosion treatment, the root mean square roughness S q The calculation formula is as follows: ， Where S q is the root mean square roughness of the rock sample surface after erosion treatment, m and n are the number of point cloud rows and columns respectively, and Z m is the average height of the z coordinates of all 3D point clouds, Z i ,j is the z-coordinate of the point in the i-th row and j-th column.

Images