Method and system for rapid detection of subgrade compaction degree based on three-dimensional point cloud technology

Abstract

The application discloses a subgrade compaction degree rapid detection method and system based on three-dimensional point cloud technology and relates to the technical field of three-dimensional scanning. The method comprises the following steps: obtaining a normalized point cloud of an inner wall of a subgrade compaction degree test hole, combining macroscopic prior knowledge of a test hole drilling tool, adopting a process of fusing two algorithms to fit a base cylinder model parameter and calculate an inner point rate; calculating a normal distance of the point cloud to a model surface as a geometric residual, expanding a cylindrical surface to a two-dimensional plane and mapping the residual, and interpolating to generate an initial residual grid field; obtaining physical mechanics parameters of filling materials to deduce constraint conditions, thereby optimizing the initial residual grid field to obtain a physical compliance residual grid field; mapping the physical compliance residual grid field back to three dimensions and superimposing the physical compliance residual grid field to the base model to generate a triangular mesh model and calculate a volume through numerical integration, outputting a result after correction of the inner point rate, and combining filling material performance parameters to calculate the compaction degree. The application effectively improves the accuracy of subgrade compaction degree detection.

Description

Technical Field

[0001] This invention relates to the field of three-dimensional scanning technology, and more specifically to a method and system for rapid detection of roadbed compaction based on three-dimensional point cloud technology. Background Technology

[0002] With the rapid development of highway construction and the accelerating pace of construction, the requirements for real-time and efficiency in roadbed compaction testing are increasing. Traditional testing methods, such as the sand cone method, which are destructive and time-consuming, are no longer adequate for modern engineering needs. Against this backdrop, volume measurement methods based on three-dimensional point cloud technology, due to their non-contact and high-speed advantages, provide a new technical approach for rapid compaction testing.

[0003] However, existing technologies have not developed tailored solutions for the specific scenarios of highway subgrade compaction testing. Directly applying general 3D reconstruction algorithms to compaction test tunnels composed of loose subgrade fill faces numerous challenges. The point cloud on the inner wall of the subgrade compaction test tunnel suffers from sparseness, high noise, and data voids. The surfaces reconstructed by general algorithms are prone to geometric distortion and may produce shapes that violate soil mechanics principles, such as overhangs and sharp protrusions, leading to excessive calculation errors and insufficient accuracy. Summary of the Invention

[0004] This invention provides a rapid detection method and system for roadbed compaction based on three-dimensional point cloud technology, aiming to solve the technical problem of insufficient accuracy in roadbed compaction detection in existing technologies.

[0005] In view of the above problems, the present invention provides a method and system for rapid detection of roadbed compaction based on three-dimensional point cloud technology.

[0006] In a first aspect, the present invention provides a rapid method for detecting the compaction degree of roadbed based on three-dimensional point cloud technology, including:

[0007] Obtain the normalized point cloud of the inner wall of the roadbed compaction test tunnel; based on the macroscopic prior knowledge formed by the roadbed compaction test tunnel drilling tool, the normalized point cloud is fitted using a fitting process that combines the random sampling consensus algorithm and the minimum median flat method to solve the parameters of the matrix cylindrical model representing the macroscopic shape of the roadbed compaction test tunnel, and calculate the inner point ratio.

[0008] Calculate the normal distance from each point in the normalized point cloud to the surface of the base cylindrical model, and define it as the geometric residual; unfold the cylindrical surface of the base cylindrical model to a two-dimensional parameter plane, map each point and its corresponding geometric residual to the two-dimensional parameter plane to form a sparse residual distribution, and interpolate the sparse residual distribution to generate a continuous initial residual grid field.

[0009] Obtain the physical and mechanical parameters of the roadbed fill material, and derive the maximum surface slope constraint and the minimum surface undulation wavelength constraint based on the physical and mechanical parameters as boundary conditions. The physical and mechanical parameters include at least the internal friction angle and the maximum particle size.

[0010] Based on the initial residual mesh field, and with the boundary conditions as constraints, an optimization function is established with the smoothness of the residual field and the goodness of fit to the measurement data as objectives. The optimization function is solved to obtain the optimized physical compliance residual mesh field.

[0011] The physical compliance residual mesh field is mapped back to three-dimensional space and superimposed on the surface of the matrix cylindrical model to generate the final three-dimensional composite surface triangular mesh model. The calculated volume of the roadbed compaction test hole is obtained by numerical integration. The confidence level of the calculated volume is corrected by the interior point ratio, and the final volume value and the corresponding uncertainty range are output.

[0012] The performance parameters of the roadbed fill material are obtained, and a compaction degree calculation formula is built based on the performance parameters. The final volume value is substituted into the formula to calculate the roadbed compaction degree. The performance parameters include humidity, density, and moisture content.

[0013] Secondly, this invention provides a rapid roadbed compaction degree detection system based on three-dimensional point cloud technology, comprising:

[0014] The matrix cylindrical model acquisition module is used to acquire the normalized point cloud of the inner wall of the roadbed compaction test tunnel; based on the macroscopic prior knowledge formed by the roadbed compaction test tunnel drilling tool, the normalized point cloud is fitted using a fitting process that combines the random sampling consensus algorithm and the minimum median flat method to solve for the matrix cylindrical model parameters representing the macroscopic shape of the roadbed compaction test tunnel, and to calculate the inner point ratio.

[0015] The initial residual mesh acquisition module is used to calculate the normal distance from each point in the normalized point cloud to the surface of the base cylindrical model, which is defined as the geometric residual; unfold the cylindrical surface of the base cylindrical model to a two-dimensional parameter plane, map each point and its corresponding geometric residual to the two-dimensional parameter plane to form a sparse residual distribution, and interpolate the sparse residual distribution to generate a continuous initial residual mesh field;

[0016] The optimized boundary condition acquisition module is used to acquire the physical and mechanical parameters of the roadbed fill material, and derive the maximum surface slope constraint and the minimum surface undulation wavelength constraint based on the physical and mechanical parameters as boundary conditions. The physical and mechanical parameters include at least the internal friction angle and the maximum particle size.

[0017] The residual mesh optimization module is used to establish an optimization function based on the initial residual mesh field and the boundary conditions, with the smoothness of the residual field and the goodness of fit to the measurement data as objectives, solve the optimization function, and obtain the optimized physical compliance residual mesh field.

[0018] The final volume acquisition module is used to map the physical compliance residual mesh field back to three-dimensional space, superimpose it onto the surface of the matrix cylindrical model, generate the final three-dimensional composite surface triangular mesh model, and use the numerical integration method to obtain the calculated volume of the roadbed compaction test hole; use the interior point ratio to correct the confidence level of the calculated volume, and output the final volume value and the corresponding uncertainty range;

[0019] The subgrade compaction degree calculation module is used to obtain the performance parameters of the subgrade fill material, build a compaction degree calculation formula based on the performance parameters, and substitute the final volume value to calculate the subgrade compaction degree; wherein, the performance parameters include humidity, density, and moisture content.

[0020] One or more technical solutions provided in this invention have at least the following technical effects or advantages:

[0021] This invention provides a rapid detection method and system for subgrade compaction based on 3D point cloud technology. By combining prior macroscopic knowledge from test drilling tools with a fusion fitting algorithm, the accuracy and noise resistance of fitting the cylindrical model of the subgrade are improved. Through residual mapping, interpolation, and optimization using physical and mechanical constraints, distortions in surface reconstruction and morphologies that violate soil mechanics are avoided, ensuring the rationality of the residual grid field. Numerical integration and interior point rate confidence correction reduce errors in test hole volume calculation, clarify result uncertainty, and improve detection reliability. Finally, by accurately calculating compaction based on filler performance parameters, this method eliminates the cumbersome procedures of traditional methods, achieving rapid and accurate detection of subgrade compaction while balancing efficiency and accuracy. It meets the quality control requirements of modern highway subgrade construction and provides reliable technical support for subgrade bearing capacity and long-term stability. Attached Figure Description

[0022] Figure 1 This is a flowchart illustrating the rapid detection method for roadbed compaction based on three-dimensional point cloud technology provided in an embodiment of the present invention.

[0023] Figure 2 This is a schematic diagram of the structure of the rapid roadbed compaction detection system based on three-dimensional point cloud technology provided in an embodiment of the present invention;

[0024] The components represented by each number in the attached diagram are explained below:

[0025] The module includes: 11 for obtaining the cylindrical model of the base body; 12 for obtaining the initial residual mesh; 13 for obtaining the optimized boundary conditions; 14 for optimizing the residual mesh; 15 for obtaining the final volume; and 16 for calculating the compaction degree of the roadbed. Detailed Implementation

[0026] This invention provides a rapid detection method and system for roadbed compaction based on three-dimensional point cloud technology, which is used to address the technical problem of insufficient accuracy in roadbed compaction detection in existing technologies.

[0027] Example 1, as Figure 1 As shown, this invention provides a rapid method for detecting the compaction degree of roadbed based on three-dimensional point cloud technology, the method comprising:

[0028] S100: Obtain the normalized point cloud of the inner wall of the roadbed compaction test tunnel; based on the macroscopic prior knowledge formed by the roadbed compaction test tunnel drilling tool, the normalized point cloud is fitted using a fitting process that combines the random sampling consensus algorithm and the minimum median flat method to solve for the parameters of the matrix cylindrical model representing the macroscopic shape of the roadbed compaction test tunnel, and the inner point ratio is calculated.

[0029] In this embodiment of the invention, a normalized point cloud of the inner wall of a roadbed compaction test tunnel is obtained. Based on the macroscopic prior knowledge formed by the roadbed compaction test tunnel drilling tool, a fitting process combining the random sampling consensus algorithm and the minimum median square method is used to solve the parameters of the matrix cylindrical model representing the macroscopic shape of the roadbed compaction test tunnel, and the in-point ratio is calculated. The three-dimensional point cloud of the inner wall of the roadbed compaction test tunnel is affected by the loose roadbed fill and interference from the acquisition environment, resulting in a large number of noise points and anomalies. Direct geometric fitting is prone to model deviation and distortion. The test tunnel is formed by drilling with a special drilling tool and has fixed inherent geometric features. It is necessary to build macroscopic prior constraints based on these features and use a fusion fitting algorithm with strong anti-interference ability to accurately extract the macroscopic cylindrical shape of the test tunnel. Furthermore, the fitting quality needs to be evaluated through quantitative indicators to provide an accurate basic model for subsequent geometric residual calculation and surface reconstruction, avoiding deviations in subsequent compaction detection caused by matrix model errors.

[0030] Step S100 in the method provided in this embodiment of the invention includes:

[0031] First, a normalized point cloud of the inner wall of the roadbed compaction test tunnel is obtained. Normalized point cloud refers to the coordinate standardization processing of the original three-dimensional point cloud of the test tunnel inner wall, scaling the coordinate values ​​of all points to a uniform range, eliminating magnitude differences caused by acquisition equipment and distance, and providing a 3D point set that facilitates subsequent geometric calculations. Normalized point cloud means scaling the original point cloud coordinates to a uniform range of [0,1], eliminating magnitude differences in coordinate magnitude caused by acquisition equipment and scanning distance, facilitating subsequent geometric calculations, and eliminating additional errors. A handheld laser scanning device was used to scan the inner wall of the test cavity at a constant speed for one revolution to obtain the original three-dimensional point cloud data of the inner wall of the test cavity, ensuring that the point cloud covered the entire inner wall of the test cavity without obvious voids. A direct-pass filtering algorithm was used to remove outliers in the z-axis direction that exceeded the height range of the test cavity, retaining the valid original point cloud. Normalization calculation: The minimum and maximum values ​​of the x, y, and z coordinate axes of the original point cloud after preprocessing were statistically analyzed. The Min-Max normalization formula was used to normalize the x, y, and z coordinates of each point. After all points were normalized, the normalized point cloud was obtained.

[0032] For example, a handheld laser scanning device with a resolution of 0.1 mm is used to scan a test hole with a depth of 500 mm to obtain the original point cloud. After preprocessing, 1000 valid points are retained. The coordinate range of the original point cloud is: x∈[1200,1300] mm, y∈[800,900] mm, z∈[0,500] mm. After performing Min-Max normalization on one of the original points (x=1250 mm, y=850 mm, z=250 mm), the coordinates of the point are (0.5,0.5,0.5). Similarly, the normalization of 1000 points is completed to obtain the normalized point cloud.

[0033] Secondly, based on the macroscopic prior knowledge formed by the roadbed compaction test hole drilling tool, the normalized point cloud is fitted using a fitting process that combines the random sampling consensus algorithm and the minimum median flat method to solve for the parameters of the matrix cylindrical model representing the macroscopic shape of the roadbed compaction test hole, and the in-point ratio is calculated.

[0034] Among them, based on the macroscopic prior knowledge formed by the roadbed compaction test tunnel drilling tool, the normalized point cloud is fitted using a process that integrates the random sampling consensus algorithm and the minimum median flat method to solve for the parameters of the matrix cylindrical model representing the macroscopic shape of the roadbed compaction test tunnel, and the in-point ratio is calculated, including:

[0035] Based on the parameters of the test hole drilling tool for roadbed compaction, the macroscopic prior knowledge is obtained, wherein the macroscopic prior knowledge includes at least the test hole diameter range, the approximate cylindrical shape, and the longitudinal axis direction of the test hole;

[0036] Based on the aforementioned macroscopic prior knowledge, a reasonable sampling range for the cylinder radius is set, and multiple initial cylinder hypotheses are generated by randomly selecting the minimum point set in the normalized point cloud.

[0037] For each initial cylindrical hypothesis, the distances from all points to the surface of the cylindrical hypothesis are calculated to obtain a set of distances. The median of the set of distances is selected as the goodness-of-fit index of the initial cylindrical hypothesis.

[0038] The initial cylinder assumption with the smallest goodness-of-fit index is selected as the optimal initial solution of the base cylinder model.

[0039] Based on the optimal initial solution, a distance threshold is set, and points in the normalized point cloud whose distance to the surface of the base cylindrical model is less than the distance threshold are determined as inner points, and points whose distance is greater than or equal to the distance threshold are determined as outer points.

[0040] Based on the point set consisting of all interior points, the least squares method is used for iterative optimization to solve for the parameters of the base cylinder model.

[0041] Calculate the inlier rate, which is the ratio of the number of inliers to the total number of points in the normalized point cloud.

[0042] First, based on the parameters of the test hole drilling tool for roadbed compaction, the macroscopic prior knowledge is obtained. This macroscopic prior knowledge includes at least the test hole diameter range, the approximate cylindrical shape, and the longitudinal axis direction of the test hole. Macroscopic prior knowledge refers to the inherent geometric characteristics of the test hole directly determined by the drilling tool specifications. After quantification, these characteristics serve as rigid constraints for subsequent fitting, requiring no additional derivation and can be directly applied using the drilling tool parameters.

[0043] Specifically, the specifications of the drilling tools used for this test borehole were retrieved from the construction log, clarifying two parameters: drill bit diameter D and borehole verticality requirements. Quantified macroscopic prior knowledge was established: Test borehole diameter range: based on drill bit diameter D, considering drilling error = [D×(1-3%), D×(1+3%)], with units standardized in mm; Approximate cylindrical shape: the inner wall of the test borehole was determined to be cylindrical, with no obvious irregular protrusions, and protrusion height ≤2mm, based on drilling tool accuracy; Test borehole longitudinal axis direction: based on borehole verticality requirements, the angle between the test borehole longitudinal axis and the plumb line was determined to be ≤0.5°, i.e., the longitudinal axis direction was approximately perpendicular to the ground. This macroscopic prior knowledge was compiled into written constraints for subsequent radius sampling range setting and model verification.

[0044] For example, this test hole uses a special drilling tool with a drill bit diameter D=150.0mm, and the drilling verticality requirement is ≤0.5°; quantified macroscopic prior knowledge: test hole diameter range: 150.0×(1-3%)=145.5mm, 150.0×(1+3%)=154.5mm, that is, diameter range [145.5mm, 154.5mm]; approximate cylindrical shape: the test hole is cylindrical as a whole, and the height of the hole wall protrusion is ≤2mm; longitudinal axis direction of the test hole: the angle between the longitudinal axis and the plumb line is ≤0.5°, approximately perpendicular to the ground.

[0045] Secondly, based on the aforementioned macroscopic prior knowledge, a reasonable sampling range for the cylinder radius is set. Multiple initial cylinder hypotheses are generated by randomly selecting the minimum point set from the normalized point cloud. The cylinder radius sampling range refers to the radius interval converted from the diameter range of the test hole, used to constrain the generation of initial hypotheses. The minimum point set consists of 3 non-collinear points, representing the minimum number of points required to fit the cylinder. The radius sampling range is converted from the diameter range and normalized. 100 sets of 3 non-collinear points are randomly selected from the normalized point cloud, and initial cylinder hypotheses are generated using the cylinder fitting formula, eliminating invalid hypotheses with radii exceeding the range. For example, the radius sampling range is 72.75mm~77.25mm, which after normalization is 0.1455~0.1545; 100 sets of non-collinear points are selected, and 100 valid initial cylinder hypotheses are generated through fitting.

[0046] Next, for each initial cylindrical hypothesis, the distances from all points to the surface of that hypothesis are calculated, resulting in a distance set. The median of this distance set is selected as the goodness-of-fit index for the initial cylindrical hypothesis. The goodness-of-fit index, the median distance from points to the cylindrical surface, measures the degree of fit between the cylindrical hypothesis and the point cloud; a smaller median indicates a better fit. Applying the formula for the distance from points to the cylindrical surface, the distances to all points under each hypothesis are calculated, generating a distance set. The distance sets are then sorted, and the median is taken as the goodness-of-fit index for that hypothesis. For example, calculating the distance values ​​of 1000 points for a certain initial cylindrical hypothesis, after sorting, the 500th point has a distance of 0.008, the 501st point has a distance of 0.009, and the goodness-of-fit index is 0.0085.

[0047] Furthermore, the initial cylindrical hypothesis with the smallest goodness-of-fit index is selected as the optimal initial solution for the base cylindrical model. The optimal initial solution refers to the cylindrical hypothesis with the smallest goodness-of-fit index and the highest fit to the point cloud among all initial cylindrical hypotheses, and it forms the basis for subsequent model optimization. By comparing the goodness-of-fit indices of all initial hypotheses, the hypothesis with the smallest index value is selected as the optimal initial solution for the base cylindrical model, and its parameters are recorded. For example, among 100 initial cylindrical hypotheses, the smallest goodness-of-fit index is 0.005, corresponding to the hypothesis parameters of radius 0.147 and axis (0.502, 0.498, 0.401), which is then designated as the optimal initial solution.

[0048] Furthermore, based on the optimal initial solution, a distance threshold is set. Points in the normalized point cloud whose distance to the surface of the base cylindrical model is less than the distance threshold are identified as interior points, and points whose distance is greater than or equal to the distance threshold are identified as exterior points. Interior points are valid points whose distance to the surface of the cylindrical model is less than the threshold; exterior points are abnormal noise points whose distance is greater than or equal to the threshold. The distance threshold is used to distinguish between valid points and abnormal points. Using four times the goodness-of-fit index of the optimal initial solution as the distance threshold, the distance from each point to the optimal solution is calculated one by one, and interior and exterior points are divided according to the threshold and their numbers are counted. For example, if the distance threshold is set to 0.02, 902 interior points and 98 exterior points are selected from 1000 points.

[0049] Subsequently, based on the point set composed of all interior points, the least squares method is used for iterative optimization to solve for the parameters of the base cylinder model. The least squares method is an optimization algorithm that aims to minimize the sum of squared errors from interior points to the model, iteratively correcting the cylinder parameters to improve model accuracy. Using the optimal initial solution as the initial value, the cylinder radius and axis coordinates are iteratively corrected using the least squares method. The process terminates after satisfying the error difference threshold or the required number of iterations, yielding the final base cylinder model parameters. For example, after 12 iterations, the error meets the termination condition, resulting in normalized parameters: radius 0.1473, axis (0.5023, 0.4978, 0.4012), which are then denormalized to obtain the actual engineering parameters.

[0050] Finally, the inlier rate is calculated, which is the ratio of the number of inliers to the total number of points in the normalized point cloud. The inlier rate, used to quantitatively evaluate the fitting quality of the cylindrical model, is calculated by dividing the number of inliers by the total number of points in the point cloud. The fitting quality is then judged according to industry standards. For example, if there are 902 inliers and 1000 total points, the inlier rate = 902 / 1000 × 100% = 90.2%, indicating a satisfactory fitting quality.

[0051] In this embodiment of the invention, relying on the macroscopic prior knowledge of the drilling tool, the random sampling consensus algorithm and the minimum median flat method are integrated to effectively filter out point cloud noise and outliers, and accurately solve for the base cylindrical model that fits the real shape of the test hole; by quantifying the fitting quality through the internal point rate, a precise benchmark is provided for subsequent geometric residual calculation and surface reconstruction, thereby reducing the compaction degree detection error from the source.

[0052] S200: Calculate the normal distance from each point in the normalized point cloud to the surface of the base cylindrical model, and define it as the geometric residual; unfold the cylindrical surface of the base cylindrical model to a two-dimensional parameter plane, map each point and its corresponding geometric residual to the two-dimensional parameter plane to form a sparse residual distribution, and interpolate the sparse residual distribution to generate a continuous initial residual grid field.

[0053] In this embodiment of the invention, the normal distance from each point in the normalized point cloud to the surface of the base cylindrical model is calculated and defined as the geometric residual. The cylindrical surface of the base cylindrical model is unfolded onto a two-dimensional parameter plane, and each point and its corresponding geometric residual are mapped onto this two-dimensional parameter plane to form a sparse residual distribution. The sparse residual distribution is then interpolated to generate a continuous initial residual mesh field. The base cylindrical model only reflects the macroscopic morphology of the test cavity, and there are slight deviations between the normalized point cloud and the model surface. These deviations are core data characterizing the true undulations of the inner wall of the test cavity. However, the point cloud is sparse and discretely distributed, and directly using it for subsequent optimization can easily lead to discontinuities and distortions in the surface reconstruction. It is necessary to accurately calculate the geometric residuals and distinguish the undulation directions first, and then convert the three-dimensional sparse deviations into a two-dimensional distribution through cylindrical surface unfolding and mapping. A continuous mesh field is generated through interpolation, providing a regular and continuous residual basis for subsequent physical constraint optimization, ensuring the continuity and accuracy of subsequent surface reconstruction.

[0054] Step S200 in the method provided in this embodiment of the invention includes:

[0055] The calculation of the normal distance from each point in the normalized point cloud to the surface of the base cylindrical model is defined as the geometric residual; including:

[0056] Based on the parameters of the base cylindrical model, the shortest distance vector from each normalized point to the surface of the base cylindrical model is calculated using analytical geometry methods.

[0057] Verify whether the shortest distance vector is consistent with the normal vector direction of the projection of the base cylindrical model surface at the corresponding point. If they are consistent, confirm that it is the true normal distance and use the absolute value of the true normal distance as the geometric residual. If they are inconsistent, mark the corresponding point as an invalid point and discard it.

[0058] The sign of the normal distance is also recorded, including both positive and negative signs.

[0059] First, based on the parameters of the base cylinder model, the shortest distance vector from each normalized point to the surface of the base cylinder model is calculated using analytical geometry. The shortest distance vector is a directed line segment pointing from the normalized point to the surface of the base cylinder model, with a length equal to the shortest distance from that point to the model surface, including both the magnitude and direction of the distance. Analytical geometry refers to the method of calculating the distance vector by applying fixed geometric formulas using known base cylinder model parameters. The normalized parameters of the base cylinder model obtained from S100 are retrieved: radius, axis coordinates, and vertical axis direction vector. For each normalized point (x, y, z), the shortest distance vector from that point to the model surface is calculated using analytical geometry formulas, with the direction pointing from the point to the model surface and the length being the shortest distance.

[0060] For example, continuing with the base cylinder model of S100, the normalized parameters are: radius = 0.1473, axis (0.5023, 0.4978, 0.4012), and vertical axis vector (0, 0, 1.0008). Taking the normalized point (0.5, 0.5, 0.5), the shortest distance vector is calculated using the formula: the direction points to the model surface, and the length is 0.0027.

[0061] Next, verify whether the shortest distance vector is consistent with the normal vector direction of the projection of the base cylindrical model surface at the corresponding point. If consistent, confirm it as the true normal distance, and use the absolute value of the true normal distance as the geometric residual. If inconsistent, mark the corresponding point as invalid and discard it. The normal vector is a vector perpendicular to the tangent direction at a point on the base cylindrical model surface, with the positive direction pointing outwards from the model. The true normal distance is the distance when the shortest distance vector and the normal vector direction of the projection point are consistent, and its absolute value is the geometric residual, representing the magnitude of the deviation between the point and the model surface. Calculate the projection point of the normalized point on the base cylindrical model surface, solve for the normal vector at the projection point; verify whether the shortest distance vector and the normal vector direction are consistent, with a vector dot product ≥ 0 indicating consistency; if consistent, confirm the shortest distance as the true normal distance, and use its absolute value as the geometric residual of that point.

[0062] For example, the projection of point (0.5, 0.5, 0.5) onto the model surface is calculated. Its normal vector direction is consistent with the shortest distance vector direction, and the dot product = 0.98 ≥ 0. The true normal distance is 0.0027. Therefore, the geometric residual of this point = |0.0027| = 0.0027. If they are inconsistent, the corresponding point is marked as an invalid point and discarded.

[0063] In addition, the signs of the normal distances are recorded, including positive and negative signs. A positive sign (+) indicates that the normalized point is located on the theoretical outer side of the base cylindrical model; a negative sign (-) indicates that it is located on the theoretical inner side, corresponding to a concave inner wall of the test hole. Based on the true normal distance, its sign is recorded: the direction of the normal distance is consistent with the direction of the normal vector of the projected point, pointing outwards from the model, and is recorded as positive; the direction is opposite, and is recorded as negative. For example, the direction of the true normal distance of a point is consistent with the normal vector, and is recorded as positive (+); the geometric residual is 0.0027, and the sign is +, indicating that the inner wall of the test hole corresponding to this point is slightly convex.

[0064] Specifically, the cylindrical surface of the base cylindrical model is unfolded onto a two-dimensional parametric plane, and each point and its corresponding geometric residual are mapped onto this two-dimensional parametric plane to form a sparse residual distribution. The sparse residual distribution is then interpolated to generate a continuous initial residual mesh field, including:

[0065] Establish a two-dimensional parametric plane, where the U-axis of the two-dimensional parametric plane represents the circumferential angle of the cylindrical surface, and the V-axis represents the axial height of the cylindrical surface;

[0066] For each point in the normalized point cloud, calculate the axial height and circumferential angle of its projection point on the surface of the base cylindrical model to determine the coordinate position of the corresponding point on the two-dimensional parameter plane; at the same time, use the geometric residual of the corresponding point as the attribute value of the coordinate position.

[0067] All points in the normalized point cloud are sequentially mapped to the two-dimensional parameter plane to form a discrete set of points with attribute values ​​of geometric residuals, constituting the sparse residual distribution;

[0068] Generate a two-dimensional parametric planar regular mesh based on preset rules;

[0069] A preset interpolation algorithm is used to estimate the residual value at each grid node of the two-dimensional parametric plane regular grid, thereby generating the initial residual grid field that covers the entire parameter domain.

[0070] First, a two-dimensional parametric plane is established. The U-axis of this plane represents the circumferential angle of the cylindrical surface, and the V-axis represents the axial height of the cylindrical surface. The two-dimensional parametric plane refers to the two-dimensional coordinate system used to map the three-dimensional point cloud and its geometric residuals, where the surface of the three-dimensional base cylindrical model is unfolded into a plane. The U-axis and V-axis are its coordinate axes. Specifically, the U-axis is the circumferential angle around the longitudinal axis of the base cylindrical model, ranging from 0 to 360° or 0 to 2π radians; the V-axis is the height along the longitudinal axis of the base cylindrical model, consistent with the axial height of the model. A two-dimensional rectangular coordinate system is established as the parametric plane, with the coordinate axes defined as follows: the U-axis represents the circumferential angle of the base cylindrical surface, and the V-axis represents the axial height of the cylindrical surface. Normalized coordinates are used, consistent with the normalized system of the point cloud. The origin (0,0) corresponds to the circumferential starting point and axial lowest point of the center of the cylinder's base.

[0071] Secondly, for each point in the normalized point cloud, the axial height and circumferential angle of its projection point on the surface of the base cylindrical model are calculated to determine the coordinate position of the corresponding point on the two-dimensional parameter plane; simultaneously, the geometric residual of the corresponding point is used as the attribute value of that coordinate position. The axial height of the projection point refers to the height of the normalized point's projection point on the surface of the base cylindrical model along the longitudinal axis of the model, which corresponds to the V-axis after normalization. The circumferential angle of the projection point refers to the circumferential angle of the projection point around the longitudinal axis of the model, corresponding to the U-axis. The attribute value is the geometric residual corresponding to the two-dimensional coordinates. The axial height (V value) of the projection point of each normalized point on the surface of the base cylindrical model is calculated by directly taking the normalized z-coordinate of the projection point; the circumferential angle (U value) of the projection point is calculated by the horizontal angle between the projection point and the model axis; (U value, V value) is used as the coordinates of that point on the two-dimensional parameter plane, and the geometric residual is used as the attribute value of that coordinate.

[0072] For example, take the normalized point (0.5, 0.5, 0.5), and the normalized axial height V of its projection point is 0.5; calculate the circumferential angle U = 90°; therefore, the coordinates of this point in the two-dimensional parameter plane are (90°, 0.5), and the calculated geometric residual is 0.0027, that is, the attribute value is 0.0027.

[0073] Next, all points in the normalized point cloud are sequentially mapped to the two-dimensional parameter plane to form a discrete set of points with geometric residuals as attribute values, constituting the sparse residual distribution. A sparse residual distribution refers to the discrete coordinate point set formed after mapping all normalized points one by one to the two-dimensional parameter plane. Each point corresponds to a unique (U,V) coordinate and geometric residual attribute. Due to the discrete nature of the point cloud, the distribution is sparse. By mapping all normalized points in S100 one by one to the two-dimensional parameter plane and recording the (U,V) coordinates and corresponding geometric residuals of each point, all discrete points constitute a sparse residual distribution.

[0074] For example, 1000 normalized points are mapped sequentially to a two-dimensional parameter plane to obtain 1000 discrete points, such as (90°,0.5)-0.0027, (180°,0.3)-0.0015, (270°,0.7)-0.0032, etc. All points are discretely distributed on the plane, forming a sparse residual distribution.

[0075] Then, a two-dimensional parametric planar regular mesh is generated based on preset rules.

[0076] The generation of a two-dimensional parametric planar regular mesh based on preset rules includes:

[0077] Based on the height and circumference of the base cylindrical model, a rectangular computational domain covering the entire cylindrical surface is predefined on the two-dimensional parameter plane;

[0078] Within the rectangular computational domain, the rectangular computational domain is divided along the U-axis and V-axis at preset intervals to generate the two-dimensional parametric planar regular mesh.

[0079] First, based on the height and perimeter of the base cylindrical model, a rectangular computational domain covering the entire cylindrical surface is predefined on the two-dimensional parameter plane. The rectangular computational domain is a rectangular area covering the entire two-dimensional parameter plane, and its extent is determined by the perimeter (corresponding to the U-axis) and height (corresponding to the V-axis) of the base cylindrical model. The rectangular computational domain is defined as follows: U-axis range 0~360°, corresponding to the model perimeter; V-axis range 0~1, corresponding to the normalized height of the model, thus forming the rectangular computational domain.

[0080] Secondly, within the rectangular computational domain, the domain is divided along the U-axis and V-axis at preset intervals to generate the two-dimensional parametric planar regular mesh. The two-dimensional parametric planar regular mesh refers to a regular mesh formed by dividing the rectangular computational domain along the U and V axes at fixed intervals, with the mesh nodes serving as calculation points for subsequent interpolation. The division intervals are set as follows: U-axis interval 5°, V-axis interval 0.05. The rectangular computational domain is divided along the U and V axes at the set intervals to generate the regular mesh; the number of nodes is calculable, ensuring coverage of the entire plane. For example, if the rectangular computational domain U∈[0°,360°], V∈[0,1], the U-axis is divided at 5° intervals, resulting in 73 nodes, and the V-axis is divided at 0.05 intervals, resulting in 21 nodes, generating a regular mesh with 73×21=1533 mesh nodes.

[0081] Finally, a preset interpolation algorithm is applied to the sparse residual distribution to estimate the residual value at each grid node of the regular grid in the two-dimensional parameter plane, thereby generating a continuous initial residual grid field covering the entire parameter domain. The preset interpolation algorithm uses Kriging interpolation or radial basis functions to estimate the discrete values ​​of the sparse residual distribution as the residual values ​​of all nodes in the regular grid. The initial residual grid field covers the entire two-dimensional parameter plane, and each grid node has a continuous residual field corresponding to its residual value, used for subsequent optimization. Kriging interpolation is used; using the discrete points of the sparse residual distribution as samples, interpolation calculations are performed on each node of the two-dimensional regular grid to estimate the residual value of each node; after interpolation of all nodes is completed, a continuous initial residual grid field covering the entire parameter domain is formed.

[0082] For example, using the Kriging interpolation method, the geometric residuals of 1000 discrete points are used as samples to interpolate 1533 grid nodes; for example, the residual value of grid node (90°, 0.5) after interpolation is 0.0026. After all nodes are interpolated, a continuous initial residual grid field is generated.

[0083] In this embodiment of the invention, the geometric residuals between the point cloud and the base cylindrical model are accurately calculated, and the protrusions and depressions of the inner wall of the test hole are distinguished by symbols. By expanding and mapping the cylindrical surface, the three-dimensional sparse deviation is converted into a two-dimensional sparse residual distribution. Combined with the interpolation algorithm, a continuous initial residual mesh field is generated, which solves the problem that the sparse point cloud cannot be directly used for subsequent optimization. It provides regular and continuous basic data for subsequent physical constraint-based residual field optimization, ensuring the continuity and accuracy of subsequent surface reconstruction.

[0084] S300: Obtain the physical and mechanical parameters of the roadbed fill material, and derive the maximum surface slope constraint and the minimum surface undulation wavelength constraint based on the physical and mechanical parameters as boundary conditions. The physical and mechanical parameters include at least the internal friction angle and the maximum particle size.

[0085] In this embodiment of the invention, the physical and mechanical parameters of the subgrade fill material are obtained. Based on these parameters, the maximum surface slope constraint and the minimum surface undulation wavelength constraint are derived as boundary conditions. The physical and mechanical parameters include at least the internal friction angle and the maximum particle size. The initial residual mesh field is generated solely based on geometric data, without considering the physical and mechanical properties of the subgrade fill material. This can easily lead to non-physical high-frequency fluctuations caused by point cloud noise or individual abnormal particles, such as sharp protrusions or overhangs, violating soil mechanics laws and resulting in distortion of subsequent surface reconstruction and excessive volume calculation errors. The internal friction angle and maximum particle size of the subgrade fill material are core parameters determining the physical compliance morphology of the test tunnel inner wall. Quantified boundary constraints must be derived based on these parameters to filter out non-physical fluctuations, ensuring that the residual field conforms to soil mechanics laws and providing reasonable physical boundary conditions for subsequent residual field optimization.

[0086] Step S300 in the method provided in this embodiment of the invention includes:

[0087] The internal friction angle and maximum particle size of the roadbed fill material are obtained as physical and mechanical parameters.

[0088] The tangent of the internal friction angle is set as the maximum allowable surface slope ratio. Based on the radius of the base cylindrical model and the maximum surface slope ratio, the maximum allowable change in geometric residual between any adjacent mesh nodes is calculated in the two-dimensional parametric plane unfolded on the cylindrical surface. The maximum surface slope constraint is constructed based on the maximum change in geometric residual.

[0089] Based on the maximum particle size, a minimum surface feature wavelength is set, and a spatial low-pass filter corresponding to the minimum surface feature wavelength is established in the two-dimensional parameter plane. By constraining the change amplitude of the initial residual grid field after passing through the spatial low-pass filter, the minimum surface undulation wavelength constraint is constructed.

[0090] First, obtain the internal friction angle and maximum particle size of the subgrade fill material as physical and mechanical parameters. The internal friction angle refers to the frictional characteristic parameter between subgrade fill material particles, which determines the maximum stable slope of the fill material accumulation. The maximum particle size refers to the maximum size of the particles in the subgrade fill material, which determines the minimum physical wavelength of the surface undulation of the test tunnel inner wall. Prioritize retrieving the test report of the fill material for the test tunnel section from the subgrade construction test log to extract the internal friction angle and maximum particle size; if the log is unavailable, use a direct shear test to determine the internal friction angle and a sieve analysis test to determine the maximum particle size to ensure parameter accuracy. For example, retrieve the subgrade fill material test data for the test tunnel section: internal friction angle φ = 30°, maximum particle size D_max = 20mm.

[0091] Secondly, the tangent of the internal friction angle is set as the maximum allowable surface slope ratio. Based on the radius of the matrix cylindrical model and the maximum surface slope ratio, the maximum allowable change in geometric residual between any adjacent grid nodes is calculated in the two-dimensional parametric plane unfolded from the cylindrical surface. The maximum surface slope constraint is constructed based on the maximum change in geometric residual. The maximum surface slope ratio is the tangent of the internal friction angle, tanφ, representing the maximum stable slope of the roadbed fill material. Surface morphology exceeding this value violates the laws of soil mechanics. The maximum change in geometric residual refers to the maximum allowable difference in geometric residual between adjacent grid nodes in the two-dimensional parametric plane, derived from the maximum surface slope ratio and the radius of the matrix cylindrical model. The maximum surface slope constraint is a boundary condition constructed based on the maximum change in geometric residual, limiting the residual variation amplitude of adjacent nodes in the residual field and avoiding the morphology of overly stable slopes.

[0092] Specifically, the maximum surface slope ratio is calculated by converting the internal friction angle to radians and calculating tanφ to obtain the maximum surface slope ratio. The actual spacing of the 2D parametric plane mesh nodes is converted by converting the U-axis angular spacing to arc length based on the actual radius of the base cylindrical model, and the V-axis normalized spacing to actual height. The maximum change in geometric residual = maximum surface slope ratio × actual spacing of mesh nodes. A constraint is established: the absolute value of the difference in geometric residual between any two adjacent mesh nodes in the 2D parametric plane is ≤ the maximum change in geometric residual in the corresponding direction.

[0093] For example, the actual radius of the base cylindrical model is R=73.65mm, the spacing between the U-axis of the two-dimensional parametric plane is 5°, and the normalized spacing between the V-axis is 0.05; the internal friction angle is 30°, tan30°≈0.577, that is, the maximum surface slope ratio = 0.577; the maximum variation of the residual is calculated as follows: U-axis ≈ 0.577 × 6.4 ≈ 3.7mm, which is 0.0074 after normalization, V-axis ≈ 0.577 × 25 ≈ 14.4mm, which is 0.0288 after normalization; the constraint is: the difference in residual between adjacent nodes of the U-axis ≤ 0.0074, and the difference in residual between adjacent nodes of the V-axis ≤ 0.0288.

[0094] Furthermore, a minimum surface characteristic wavelength is set based on the maximum particle size. A spatial low-pass filter corresponding to the minimum surface characteristic wavelength is established in the two-dimensional parameter plane. By constraining the variation amplitude of the initial residual grid field after passing through this spatial low-pass filter, the minimum surface undulation wavelength constraint is constructed. The minimum surface characteristic wavelength refers to the minimum physical undulation wavelength of the test cavity inner wall determined by the maximum particle size; undulations smaller than this value are non-physical. The spatial low-pass filter is a filtering tool that allows low-frequency signals to pass through while suppressing high-frequency signals, used to filter out non-physical high-frequency undulations in the residual field. The minimum surface undulation wavelength constraint means that by constraining the variation amplitude of the initial residual grid field after passing through the low-pass filter, non-physical high-frequency undulations are suppressed, ensuring that the residual field conforms to the filler characteristics.

[0095] Specifically, based on industry experience, a minimum surface feature wavelength is set: λ_min = 2 × D_max, to avoid non-physical fluctuations caused by single particles; λ_min is converted to the grid cutoff frequency of the two-dimensional parametric plane, and a Gaussian low-pass filter σ = λ_min / 6 is constructed; a constraint is constructed: after the initial residual grid field passes through this filter, the change amplitude of high-frequency components with wavelength < λ_min is ≤ 5%.

[0096] For example, the maximum particle size is 20 mm, and the minimum surface characteristic wavelength λ_min = 40 mm; the angle of λ_min to the U-axis is approximately 31.3°, and the cutoff frequency of the filter is determined; a Gaussian low-pass filter σ≈6.7 mm is constructed; the constraint is: after filtering, the amplitude change of high-frequency fluctuations with wavelength <40 mm in the initial residual grid field is ≤5%, to suppress noise or non-physical fluctuations caused by single particles.

[0097] In this embodiment of the invention, by obtaining the physical and mechanical parameters of the subgrade fill material, quantitative boundary constraints are derived and constructed: the maximum surface slope constraint limits the maximum stable slope of the test tunnel inner wall, avoiding steep slope morphology that violates the soil friction characteristics; the minimum surface undulation wavelength constraint filters out non-physical high-frequency undulations smaller than the minimum physical wavelength, eliminating point cloud noise or distortion caused by single particles. These two constraints serve as boundary conditions for subsequent residual field optimization, ensuring that the optimized residual field conforms to the laws of soil mechanics, and avoiding morphological distortion in subsequent surface reconstruction from a physical perspective.

[0098] S400: Based on the initial residual grid field, and with the boundary conditions as constraints, establish an optimization function with the smoothness of the residual field and the goodness of fit to the measurement data as objectives, solve the optimization function, and obtain the optimized physical compliance residual grid field.

[0099] In this embodiment of the invention, based on the initial residual mesh field and constrained by the boundary conditions, an optimization function is established with the smoothness of the residual field and the goodness of fit to the measured data as objectives. Solving the optimization function yields an optimized physically compliant residual mesh field. The initial residual mesh field achieves continuity only through interpolation, without considering the smoothness of the residual field and the goodness of fit to the original measured data. Furthermore, it is not included in the physical boundary constraints of S300, which can easily lead to problems such as severe local fluctuations, non-physical high-frequency fluctuations, or excessive deviations from the measured data. A dual-objective optimization function needs to be constructed to transform physical constraints into quantitative mathematical constraints. Through standardized mathematical solution methods, a physically compliant residual mesh field that conforms to the laws of soil mechanics and also considers smoothness and good data fit is obtained, providing a reliable residual foundation for subsequent three-dimensional surface reconstruction.

[0100] Step S400 in the method provided in this embodiment of the invention includes:

[0101] The optimization function includes two weighted sums, where the first term is a smoothing term, which consists of the second-order difference modulus of the grid node values, and the second term is a data fitting term, which consists of the mean square error of the interpolation residuals at the grid nodes and the geometric residuals of the corresponding grid nodes.

[0102] The maximum surface slope constraint and the minimum surface undulation wavelength constraint are respectively transformed into upper limit constraints on the gradient magnitude of the grid nodes and upper limit constraints on the second-order difference magnitude of the grid nodes.

[0103] The constrained optimization problem consisting of the optimization function and the boundary conditions is transformed into a Lagrange dual problem for solution, ultimately yielding the optimized physically compliant residual mesh field and the overall smoothness evaluation result of the physically compliant residual mesh field.

[0104] First, the optimization function comprises two weighted sums. The first term is a smoothing term, consisting of the second-order difference magnitude of the grid node values. The second term is a data fitting term, consisting of the mean square error of the interpolation residuals at the grid nodes and the corresponding geometric residuals. The optimization function is a weighted sum function aimed at achieving both residual field smoothness and data fit, balancing the two requirements of minimizing residual field fluctuations and ensuring close alignment with the original measurement data. The smoothing term, consisting of the second-order difference magnitude of the grid node values, penalizes severe fluctuations in the residual field; a smaller value results in a smoother residual field. The data fitting term, consisting of the mean square error of the interpolation residuals at the grid nodes and the corresponding geometric residuals, ensures that the optimized residual field matches the original measurement data; a smaller value indicates a higher degree of fit.

[0105] Specifically, the optimization function is defined as: F = α × S + β × D, where S is the smoothing term, D is the data fitting term, α is the weight of the smoothing term, and β is the weight of the data fitting term. Smoothing term S: ,in, Let be the residual value of grid node (i,j). For second-order difference operators; data fitting term D: Where N is the original number of point clouds, The interpolation residual for the grid node corresponding to the k-th point. Let be the geometric residual at the k-th point, and let the mean squared error represent the goodness of fit. Weighting coefficients are set: α=0.3, β=0.7, prioritizing good data fit. Based on the initial residual grid field of S200, the initial residual values ​​of all grid nodes are extracted and substituted into the formula to complete the construction of the optimization function.

[0106] For example, the smoothing term is calculated as follows: Take a group of adjacent nodes (90°, 0.5), (95°, 0.5), and (90°, 0.55), calculate the second-order difference modulus ≈ 0.0002, and sum the second-order difference moduli of all nodes to obtain an initial S ≈ 0.35; the data fitting term is calculated as follows: the geometric residual of the k-th point coordinates (90°, 0.5). =0.0027, corresponding to the interpolation residual at the grid nodes. =0.0026, single-group error (0.0026−0.0027) 2 =1×10 −8 Accumulate the errors of 1000 points and take the average, resulting in an initial D≈8×10. −9 Initial value of the optimization function: F = 0.3 × 0.35 + 0.7 × 8 × 10 −9 ≈0.105.

[0107] Secondly, the maximum surface slope constraint and the minimum surface undulation wavelength constraint are respectively transformed into upper limit constraints on the gradient magnitude of the grid nodes and upper limit constraints on the second-order difference magnitude of the grid nodes. The gradient magnitude of the grid nodes characterizes the variation range of the residual values ​​of adjacent grid nodes, corresponding to the surface slope. The upper limit constraint on the gradient magnitude is derived from the maximum surface slope constraint, limiting the slope of the residual field to not exceed the stable slope of the filler. The upper limit constraint on the second-order difference magnitude is derived from the minimum surface undulation wavelength constraint, suppressing high-frequency non-physical fluctuations in the residual field. Transformation of the maximum surface slope constraint: grid node gradient magnitude , The maximum change in geometric residuals derived for S300; the minimum surface undulation wavelength constraint: the second-order difference modulus of the mesh nodes. , Derived from the minimum surface characteristic wavelength, take =0.0005; The two upper limit constraints mentioned above are rearranged into mathematical inequalities and used as boundary conditions for the optimization function.

[0108] For example, gradient modulus constraint: the residual difference between adjacent nodes on the U-axis, such as (90°, 0.5) and (95°, 0.5), is ≤0.0074; the residual difference between adjacent nodes on the V-axis, such as (90°, 0.5) and (90°, 0.55), is ≤0.0288; second-order difference modulus constraint: the second-order difference modulus of any node is ≤0.0005, suppressing high-frequency fluctuations with wavelength <40mm; final constraint condition: ≤max(0.0074,0.0288) ≤0.0005.

[0109] Furthermore, the constrained optimization problem formed by the optimization function and the boundary conditions is transformed into a Lagrange dual problem for solution, ultimately yielding the optimized physically compliant residual grid field and the overall smoothness evaluation result of the physically compliant residual grid field. The Lagrange dual problem refers to transforming a constrained optimization problem into an unconstrained dual problem, facilitating numerical solution. The physically compliant residual grid field refers to the final residual grid field that satisfies all physical constraints and balances smoothness and fit. The overall smoothness evaluation result is based on the overall curvature energy of the residual field, representing the degree of smoothness; the smaller the value, the higher the smoothness.

[0110] Specifically, the two physical constraints mentioned above are incorporated as penalty terms into the optimization function. When the residual field parameters exceed the constraint range, the penalty term increases, thereby forcing the constraint conditions to be satisfied. The original problem of minimizing the optimization function under physical constraints is transformed into a Lagrangian dual problem, and the influence of the constraints is reflected through the dual variables, transforming the constrained solution problem into a more computationally efficient unconstrained solution problem. Mature numerical optimization algorithms used in engineering, such as the scipy library optimization module in Python and the constrained optimization function in Matlab, are selected to solve the transformed dual problem. An iterative convergence criterion is set: when the difference between the optimization function values ​​calculated in two adjacent iterations is less than 1 × 10⁻⁶, the problem is solved. -6 When the iteration stops, the iteration is stopped; after the iteration converges, the optimized residual values ​​of all grid nodes are extracted and integrated to generate a physically compliant residual grid field; at the same time, the sum of the second-order difference correlation values ​​of all grid nodes in the residual field is calculated and used as the overall smoothness evaluation result.

[0111] For example, penalty weights are set for the two physical constraints, and a Lagrangian function is constructed after incorporating them into the optimization function. This ensures that the function value increases significantly when the residual field parameters exceed the slope or wavelength constraints. After transforming the constrained optimization problem into a Lagrangian dual problem, the Python function `scipy.optimize.minimize` is called for numerical solution. For instance, after 15 iterations, the difference in the optimization function value between two adjacent iterations is approximately 8 × 10⁻⁶. -7 Once the convergence threshold requirement is met, the iteration stops; the optimized physically compliant residual grid field is output: the residual values ​​of all grid nodes not only meet the physical constraints that the variation of residuals between adjacent nodes does not exceed the maximum slope and there are no high-frequency fluctuations less than 40mm wavelength, but also maintain a high degree of fit with the original measurement data; at the same time, the sum of the second-order difference correlation values ​​of all nodes is approximately 0.12, which is the overall smoothness evaluation result, indicating that the residual field fluctuates gently and without drastic fluctuations.

[0112] In this embodiment of the invention, by constructing a dual-objective optimization function that takes into account both smoothness and data fit, the physical constraints of S300 are transformed into quantitative mathematical constraints. By solving the Lagrange dual problem, the resulting physically compliant residual grid field satisfies the mechanical laws of the roadbed filler and also takes into account the fit of the original measurement data, while possessing good smoothness. The overall smoothness evaluation result can quantitatively characterize the quality of the residual field, providing a physically compliant, data reliable, and shape smooth foundation for subsequent three-dimensional surface reconstruction.

[0113] S500: Map the physical compliance residual mesh field back to three-dimensional space, superimpose it onto the surface of the matrix cylindrical model to generate the final three-dimensional composite surface triangular mesh model, and use the numerical integration method to obtain the calculated volume of the roadbed compaction test hole; use the interior point ratio to correct the confidence level of the calculated volume, and output the final volume value and the corresponding uncertainty range.

[0114] In this embodiment of the invention, the physically compliant residual mesh field is mapped back to three-dimensional space and superimposed on the surface of the base cylindrical model to generate the final three-dimensional composite curved surface triangular mesh model. The calculated volume of the roadbed compaction test tunnel is obtained using a numerical integration method. The interior point ratio is used to correct the confidence level of the calculated volume, outputting the final volume value and the corresponding uncertainty range. The physically compliant residual mesh field is two-dimensional planar data and cannot directly represent the true three-dimensional morphology of the tunnel's inner wall. It needs to be mapped back to three-dimensional space and superimposed on the base cylindrical model to restore the actual undulations of the tunnel's inner wall. The tunnel volume needs to be accurately calculated through numerical integration. The interior point ratio reflects the validity of the original point cloud data, and the smoothness of the physically compliant residual mesh field reflects the quality of the residual field. The confidence level of the volume result needs to be corrected based on these two indicators, and the uncertainty range needs to be quantified to ensure the reliability and traceability of the volume result, providing an accurate volume basis for the final calculation of the roadbed compaction degree.

[0115] Step S500 in the method provided in this embodiment of the invention includes:

[0116] Specifically, the physically compliant residual mesh field is mapped back to three-dimensional space and superimposed onto the surface of the matrix cylindrical model to generate the final three-dimensional composite curved surface triangular mesh model. The calculated volume of the roadbed compaction test tunnel is then obtained using a numerical integration method, including:

[0117] Each grid node of the physically compliant residual grid field on the two-dimensional parameter plane is back-calculated to its corresponding three-dimensional spatial coordinates on the surface of the base cylindrical model, based on its coordinates on the two-dimensional parameter plane.

[0118] By translating the residual value corresponding to each mesh node along the direction of its normal vector at the corresponding point on the surface of the base cylindrical model, a three-dimensional point set of the composite surface is obtained.

[0119] The three-dimensional point set of the composite surface is subjected to Delaunay triangulation to generate a three-dimensional composite surface triangular mesh model representing the inner wall surface of the test hole.

[0120] Along the axial direction of the base cylindrical model, the three-dimensional composite surface triangular mesh model is sliced ​​at equal intervals with a preset spacing;

[0121] Calculate the area of ​​the polygonal cross-section formed by the intersection of each slice plane and the triangular mesh model;

[0122] The Simpson integral method is used to numerically integrate the areas of all polygonal cross sections along the axial direction to obtain the calculated volume of the roadbed compaction test hole.

[0123] First, for each grid node of the physically compliant residual grid field on the two-dimensional parameter plane, its corresponding three-dimensional spatial coordinates on the surface of the base cylindrical model are calculated based on its coordinates on the two-dimensional parameter plane. Calculating the three-dimensional spatial coordinates means inversely deriving the corresponding three-dimensional spatial position on the surface of the base cylindrical model from the (U,V) coordinates of the two-dimensional parameter plane, thus restoring the two-dimensional data to a three-dimensional form. The normal vector direction is the normal vector of the corresponding point on the surface of the base cylindrical model; the radial outward direction is the positive direction and serves as the reference direction for subsequent residual translation. The (U,V) coordinates and corresponding residual values ​​of all nodes in the physically compliant residual grid field of S400 are retrieved, along with the actual parameters of the base cylindrical model of S100. For each grid node, its radial orientation is calculated based on its circumferential angle along the U-axis, and combined with its axial height along the V-axis, it is inversely normalized to the actual height. The original three-dimensional spatial coordinates of the node on the surface of the base cylindrical model are then calculated inversely. The normal vector direction on the surface of the base cylindrical model at this three-dimensional coordinate point is determined and recorded for subsequent translation operations.

[0124] For example, the two-dimensional coordinates of a certain mesh node are U=90° and V=0.5. The original three-dimensional coordinates of the node on the surface of the base cylindrical model are (1250.23mm, 849.78mm, 250mm). The normal vector of this point is radially outward.

[0125] Secondly, each mesh node is translated along its normal vector direction at its corresponding point on the surface of the base cylindrical model by shifting the residual value corresponding to the mesh node, thus obtaining the three-dimensional point set of the composite surface. The residual value translation involves moving the three-dimensional coordinates along the normal vector direction; a positive residual value translates outwards from the model, while a negative residual value translates inwards, restoring the true undulations of the test tunnel's inner wall. The three-dimensional point set of the composite surface refers to the set of all three-dimensional coordinate points that, after superimposing the residuals, can completely represent the three-dimensional morphology of the test tunnel's inner wall. For each original three-dimensional coordinate point obtained through back-calculation, it is translated along its normal vector direction according to the physical compliance residual value corresponding to that node: a positive residual value translates outwards from the model by a corresponding distance; a negative residual value translates inwards from the model by a corresponding distance. All translated three-dimensional coordinates are recorded and integrated to form a three-dimensional point set of the composite surface covering the entire height and circumference of the test tunnel.

[0126] For example, the physical compliance residual value of the above node is 0.00265, which is 1.325mm after inverse normalization. The sign is positive. Therefore, it is translated outward along the normal vector by 1.325mm. After translation, the three-dimensional coordinates are (1251.555mm, 851.105mm, 250mm). The translation of 1533 mesh nodes is completed to obtain the complete three-dimensional point set of the composite surface.

[0127] Next, the three-dimensional point set of the composite surface is subjected to Delaunay triangulation to generate a three-dimensional composite surface triangular mesh model representing the inner wall surface of the test cavity. Delaunay triangulation is a mature method for generating three-dimensional surface meshes in the engineering field, ensuring that the triangular facets are uniform, without narrow triangles, and accurately fit the surface shape. The three-dimensional composite surface triangular mesh model is composed of triangular facets, completely representing the three-dimensional true shape of the inner wall of the test cavity, and is the core carrier for volume calculation. Commonly used 3D modeling software in engineering, such as CloudCompare or MeshLab, is used to import the three-dimensional point set of the composite surface; the built-in Delaunay triangulation function of the software is called to triangulate the three-dimensional point set to generate a triangular mesh; the mesh integrity is checked: invalid triangular facets that are too large / too small or overlapping are removed to ensure that the mesh completely covers the inner wall of the test cavity.

[0128] For example, 1533 translated 3D points are imported into CloudCompare, and Delaunay triangulation is performed to generate a 3D composite surface triangular mesh model with approximately 3000 valid triangular faces, with no invalid faces, which perfectly fits the undulating shape of the test hole's inner wall.

[0129] Furthermore, the three-dimensional composite curved surface triangular mesh model is sliced ​​at equal intervals along the axial direction of the base cylindrical model at a preset spacing. The axial direction refers to the longitudinal axis of the base cylindrical model, the direction of the test hole depth, and is perpendicular to the ground. Equal-interval slicing means cutting the triangular mesh model at fixed intervals along the axial direction to obtain parallel slicing planes. The industry standard preset spacing is 5mm, balancing computational accuracy and efficiency.

[0130] For example, the slicing range is set to cover the entire height along the test hole axis (z-axis), i.e., from z=0mm (bottom of the test hole) to z=500mm (top of the test hole); slicing planes are set at equal intervals of 5mm, generating a total of 100 slicing planes; the intersection area between each slicing plane and the three-dimensional composite surface triangular mesh model is calculated, and the closed intersection boundary is determined.

[0131] Subsequently, the area of ​​the polygonal cross-section formed by the intersection of each slice plane and the triangular mesh model is calculated. A polygonal cross-section refers to a closed polygon formed by the intersection of the slice plane and the triangular mesh model, and its area reflects the actual cross-sectional area of ​​the test hole at that height. The cross-sectional area calculation relies on the software's built-in algorithm, accurately calculating the area based on the polygon vertex coordinates, eliminating the need for manual calculation. For each slice plane, the coordinates of all vertices of the polygon formed by its intersection with the triangular mesh model are extracted; the software's polygon area calculation function is called to calculate the area of ​​each polygon, and all cross-sectional area values ​​are recorded in order of slice height.

[0132] For example, the slice plane at z=250mm intersects with the triangular mesh model to form an approximately circular polygon. After extracting the vertex coordinates, the area of ​​this cross-section is calculated to be approximately 8500mm². 2 ; Calculate the cross-sectional area of ​​100 slices in sequence to form an area sequence from z=0 to z=500mm.

[0133] Finally, Simpson's integral method was used to numerically integrate the areas of all polygonal cross sections along the axial direction to obtain the calculated volume of the roadbed compaction test hole. Simpson's integral method is a commonly used numerical integration method in engineering. By fitting the variation law of the cross-sectional area along the axial direction, it accurately calculates the overall volume of the test hole. The calculated volume refers to the original value of the test hole volume without confidence correction. The cross-sectional area values ​​of all slices were arranged into a sequence according to the axial height; numerical calculation tools such as Excel or Python's scipy library were used to call the Simpson's integral function, using the axial height as the integration variable and the cross-sectional area as the integrand, to perform integration calculations along the entire height of the test hole; the output integration result is the calculated volume of the roadbed compaction test hole.

[0134] For example, by importing the cross-sectional area sequence of 100 slices into Python, the calculated volume of the test hole using Simpson's integral method is approximately 4,250,000 mm². 3.

[0135] The confidence level of the calculated volume is corrected using the interior point ratio, and the final volume value and corresponding uncertainty range are output, including:

[0136] Based on the inlier rate, a preset confidence correction coefficient table is queried to obtain the volume correction coefficient corresponding to the current inlier rate; wherein, the confidence correction coefficient table defines a monotonically increasing mapping relationship between the inlier rate and correction coefficients greater than 0 and less than or equal to 1.

[0137] Multiply the calculated volume by the volume correction factor to obtain the preliminary corrected volume value;

[0138] Based on the current in-point rate and the overall smoothness evaluation results of the physical compliance residual grid field, a comprehensive uncertainty factor is calculated; the comprehensive uncertainty factor is negatively correlated with the in-point rate and negatively correlated with the overall smoothness of the physical compliance residual grid field.

[0139] Based on the comprehensive uncertainty factor and the preliminary corrected volume value, the uncertainty range corresponding to the final volume value is calculated and output, wherein the width of the uncertainty range is positively correlated with the comprehensive uncertainty factor.

[0140] First, based on the inner point rate, a preset confidence correction coefficient table is consulted to obtain the volume correction coefficient corresponding to the current inner point rate. The confidence correction coefficient table defines a monotonically increasing mapping relationship between the inner point rate and correction coefficients greater than 0 and less than or equal to 1. The confidence correction coefficient table is a standardized table pre-defined in the engineering field, establishing a correspondence between the inner point rate and the volume correction coefficient. The correction coefficient ranges from 0 < correction coefficient ≤ 1, and the two exhibit a monotonically increasing mapping relationship. The higher the inner point rate, the closer the correction coefficient is to 1, representing a higher confidence level in the volume result. The volume correction coefficient is a coefficient matched from the correction coefficient table based on the current inner point rate, used to correct the confidence level of the calculated volume of the test tunnel, offsetting the result deviation caused by insufficient inner point rate.

[0141] Specifically, the in-point rate of the test hole point cloud calculated in S100 is retrieved; a preset confidence correction coefficient table is retrieved, for example: the correction coefficient for an in-point rate of 80% is 0.85, the correction coefficient for an in-point rate of 85% is 0.9, ..., and the correction coefficient for an in-point rate of 100% is 1; if the in-point rate is not an integer in the table, the corresponding correction coefficient is calculated using linear interpolation; the final volume correction coefficient is determined and recorded. For example, if the in-point rate is 90.2%, the correction coefficient calculated by linear interpolation is 0.951, that is, the volume correction coefficient corresponding to the current in-point rate is 0.951.

[0142] Next, the calculated volume is multiplied by the volume correction factor to obtain the preliminary corrected volume value. The preliminary corrected volume value is the product of the calculated test tunnel volume and the volume correction factor; it is the initial volume result after correction for the interior point rate confidence level, balancing calculation accuracy and data validity. The calculated test tunnel volume and volume correction factor obtained through Simpson's integral method are retrieved, and the preliminary corrected volume value = calculated volume × volume correction factor; the preliminary corrected volume value is recorded, maintaining the same accuracy as the calculated volume. For example, the preliminary corrected volume value = 4250000 mm. 3 ×0.951=4041750mm 3 .

[0143] Next, based on the aforementioned interior point rate and the overall smoothness evaluation result of the physically compliant residual grid field, a comprehensive uncertainty factor is calculated. This comprehensive uncertainty factor is negatively correlated with both the interior point rate and the overall smoothness of the physically compliant residual grid field. The comprehensive uncertainty factor is used to quantify the error magnitude of the volumetric result. The comprehensive uncertainty factor is negatively correlated with the interior point rate; the higher the interior point rate, the smaller the comprehensive uncertainty factor and the smaller the error. It is also negatively correlated with the overall smoothness of the physically compliant residual grid field; the smaller the smoothness evaluation value, the smoother the residual field, the smaller the comprehensive uncertainty factor, and the smaller the error. The overall smoothness evaluation result is the sum of the second-order difference correlation values ​​of all nodes in the residual field statistically analyzed in S400; the smaller the value, the better the smoothness.

[0144] Specifically, retrieve the current in-point rate and the overall smoothness assessment result output by S400; calculate the factor using industry-preset quantification rules: Basic rule: Comprehensive uncertainty factor = Basic factor - In-point rate correction term - Smoothness correction term; where, Basic factor = 0.09, In-point rate correction term = In-point rate × 0.05, Smoothness correction term = (0.5 / Smoothness assessment value) × 0.01. The Basic factor is an empirical engineering value and can be adjusted within the range of 0.08 to 0.1 according to project accuracy requirements.

[0145] For example, the internal point ratio is 90.2%, and the smoothness evaluation value is 0.12. The internal point ratio correction term = 0.902 × 0.05 = 0.0451; the smoothness correction term = (0.5 / 0.12) × 0.01 ≈ 0.0417; the comprehensive uncertainty factor = 0.09 - 0.0451 - 0.0417 = 0.032.

[0146] Furthermore, based on the comprehensive uncertainty factor and the preliminary corrected volume value, the uncertainty range corresponding to the final volume value is calculated and output. The width of the uncertainty range is positively correlated with the comprehensive uncertainty factor. The uncertainty range refers to the error interval centered on the final volume value. The interval width is positively correlated with the comprehensive uncertainty factor; the larger the factor, the wider the interval and the larger the error range. The final volume value refers to the final volume value using the preliminary corrected volume value, which is the core result after confidence level correction. The preliminary corrected volume value is directly used as the final volume value; uncertainty magnitude = final volume value × comprehensive uncertainty factor; lower bound = final volume value - uncertainty magnitude; upper bound = final volume value + uncertainty magnitude; the final volume value and the uncertainty range presented in the form of lower bound to upper bound are clearly indicated.

[0147] For example, the final volume value is 4041750 mm. 3 The amplitude is calculated as follows: 4041750 × 0.0032 = 12933.6 mm. 3 The lower bound is 4041750 − 12933.6 = 4028816.4 mm. 3 The upper bound = 4041750 + 12933.6 = 4054683.6 mm 3 Output: Final volume value 4041750mm 3 The uncertainty range is 4028816.4 mm³ to 4054683.6 mm³. 3 .

[0148] In this embodiment of the invention, a two-dimensional physical compliance residual mesh field is accurately mapped back to three-dimensional space and superimposed onto a base cylindrical model to generate a three-dimensional composite curved surface triangular mesh model that fits the actual shape of the test tunnel. The volume of the test tunnel is accurately calculated using the Simpson integral method. The confidence correction based on the interior point rate improves the reliability of the volume results. The uncertainty range of the quantification of the interior point rate and the residual field smoothness is combined to clarify the error boundary of the volume results. All steps have clear parameters and mature methods, which can be directly implemented, providing accurate and verifiable volume data for the final calculation of the roadbed compaction degree.

[0149] S600: Obtain the performance parameters of the subgrade fill material, build a compaction degree calculation formula based on the performance parameters, substitute the final volume value, and calculate the subgrade compaction degree; wherein, the performance parameters include humidity, density, and moisture content.

[0150] In this embodiment of the invention, performance parameters of the subgrade filler are obtained, and a compaction degree calculation formula is constructed based on these performance parameters. The final volume value is then substituted into the formula to calculate the subgrade compaction degree. The performance parameters include humidity, density, and moisture content. The preceding steps have already completed the calculation of the final volume and uncertainty of the test tunnel. The calculation of the subgrade compaction degree requires the combined use of the test tunnel volume and the performance parameters of the subgrade filler itself: volume reflects the actual size of the test tunnel, while humidity, density, and moisture content directly determine the compaction characteristics of the subgrade filler. If volume alone cannot quantify the compaction degree of the filler, it is necessary to obtain the filler performance parameters, construct a standardized compaction degree calculation formula, and substitute the final volume into the calculation to accurately obtain the subgrade compaction degree result. This meets the requirements for subgrade construction quality inspection and acceptance, completing the closed loop of the entire subgrade compaction degree testing process.

[0151] The performance parameters include humidity, density, and moisture content. Moisture content (ω) is the ratio of the mass of water in the subgrade fill to the mass of dry soil, reflecting the moisture content of the fill and directly affecting the compaction effect. Density (ρ_d) is the dry density of the subgrade fill, referring to the ratio of the mass of the fill after drying to its corresponding volume. Humidity (H) is the moisture content state in the fill, quantified as a humidity coefficient, used to correct the influence of moisture content on the compaction degree calculation. Compaction degree (K) is a quantitative index characterizing the actual compacted quality of the subgrade, which is the ratio of the actual dry density of the subgrade to the maximum dry density of the fill.

[0152] Specifically, three key data points are extracted from construction testing records or field tests: Test cavity fill dry mass (m): the dry soil mass of the test cavity fill measured by the drying method; final test cavity volume (V): the test cavity volume after confidence correction in S500; maximum dry density of the fill (ρ_max): the maximum dry density of the fill obtained through standard compaction tests. The compaction degree is calculated using the formula: K = (ρ_field / ρ_max) × 100%, where the measured dry density ρ_field = m / V. Substituting parameters: Substituting the final test cavity volume (V), fill dry mass (m), and maximum dry density (ρ_max) into the formula, the subgrade compaction degree is calculated.

[0153] For example, the dried mass of the test cavity filler is m = 7328g, and the final volume of the test cavity is V = 4041.75cm³. 3 The maximum dry density of the filler is ρ_max = 1.95 g / cm³. 3 Actual measured dry density: ρ_field = 7328 / 4041.75 ≈ 1.813 g / cm³ 3 The compaction degree is calculated as follows: K = (1.813 / 1.95) × 100% ≈ 93.0%. Industry personnel can perform subsequent operations based on the calculation results.

[0154] In this embodiment of the invention, by acquiring the performance parameters of the subgrade fill material, a standardized and practical compaction calculation formula is established. This organically combines quantified volume data with fill material characteristic parameters to accurately calculate the subgrade compaction degree. The final output compaction degree result can be directly used for subgrade construction quality assessment and acceptance, completing a closed-loop process from test tunnel scanning, data processing, volume calculation to compaction degree assessment, meeting the practical needs and accuracy requirements of subgrade compaction degree testing in geotechnical engineering.

[0155] Through the specific implementation methods described above, the embodiments of the present invention achieve the following technical effects:

[0156] This invention provides a rapid method and system for detecting subgrade compaction based on 3D point cloud technology. It effectively overcomes the shortcomings of traditional detection methods, such as low accuracy, poor anti-interference, and inability to accurately reflect the physical properties of the subgrade filler. By relying on macroscopic prior constraints and noise-resistant fitting algorithms, the accuracy of the basic model is ensured. Physical and mechanical constraints filter out unrealistic undulations, ensuring that the reconstructed surface strictly conforms to the soil mechanics laws of the subgrade filler. Numerical integration and confidence level correction enable high-precision calculation of test tunnel volume and quantification of the error range. Finally, by combining the filler performance parameters, accurate calculation of compaction degree is completed. This method combines geometric accuracy with physical rationality. The detection process is reproducible, and the results are traceable, improving the efficiency, accuracy, and reliability of subgrade compaction detection. It can meet the actual engineering needs of highway subgrade construction quality inspection and acceptance.

[0157] Example 2, as Figure 2 As shown, this invention provides a rapid detection system for roadbed compaction based on three-dimensional point cloud technology, the system comprising:

[0158] The matrix cylindrical model acquisition module 11 is used to acquire the normalized point cloud of the inner wall of the roadbed compaction test tunnel; based on the macroscopic prior knowledge formed by the roadbed compaction test tunnel drilling tool, the normalized point cloud is fitted by a fitting process that combines the random sampling consistency algorithm and the minimum median flat method to solve the matrix cylindrical model parameters representing the macroscopic shape of the roadbed compaction test tunnel, and to calculate the inner point ratio.

[0159] The initial residual mesh acquisition module 12 is used to calculate the normal distance from each point in the normalized point cloud to the surface of the base cylindrical model, which is defined as the geometric residual; unfold the cylindrical surface of the base cylindrical model to a two-dimensional parameter plane, map each point and its corresponding geometric residual to the two-dimensional parameter plane to form a sparse residual distribution, and interpolate the sparse residual distribution to generate a continuous initial residual mesh field;

[0160] The boundary condition acquisition module 13 is used to acquire the physical and mechanical parameters of the roadbed fill material, and derive the maximum surface slope constraint and the minimum surface undulation wavelength constraint based on the physical and mechanical parameters as boundary conditions. The physical and mechanical parameters include at least the internal friction angle and the maximum particle size.

[0161] The residual mesh optimization module 14 is used to establish an optimization function based on the initial residual mesh field and the boundary conditions as constraints, with the smoothness of the residual field and the goodness of fit to the measurement data as objectives, solve the optimization function, and obtain the optimized physical compliance residual mesh field.

[0162] The final volume acquisition module 15 is used to map the physical compliance residual mesh field back to three-dimensional space, superimpose it onto the surface of the matrix cylindrical model, generate the final three-dimensional composite surface triangular mesh model, and use the numerical integration method to obtain the calculated volume of the roadbed compaction test hole; use the interior point ratio to correct the confidence of the calculated volume, and output the final volume value and the corresponding uncertainty range.

[0163] The roadbed compaction degree calculation module 16 is used to obtain the performance parameters of the roadbed fill material, build a compaction degree calculation formula based on the performance parameters, substitute the final volume value, and calculate the roadbed compaction degree; wherein, the performance parameters include humidity, density, and moisture content.

[0164] In one embodiment, the base cylinder model acquisition module 11 is further configured to:

[0165] Among them, based on the macroscopic prior knowledge formed by the roadbed compaction test tunnel drilling tool, the normalized point cloud is fitted using a process that integrates the random sampling consensus algorithm and the minimum median flat method to solve for the parameters of the matrix cylindrical model representing the macroscopic shape of the roadbed compaction test tunnel, and the in-point ratio is calculated, including:

[0166] Based on the parameters of the test hole drilling tool for roadbed compaction, the macroscopic prior knowledge is obtained, wherein the macroscopic prior knowledge includes at least the test hole diameter range, the approximate cylindrical shape, and the longitudinal axis direction of the test hole;

[0167] Based on the aforementioned macroscopic prior knowledge, a reasonable sampling range for the cylinder radius is set, and multiple initial cylinder hypotheses are generated by randomly selecting the minimum point set in the normalized point cloud.

[0168] For each initial cylindrical hypothesis, the distances from all points to the surface of the cylindrical hypothesis are calculated to obtain a set of distances. The median of the set of distances is selected as the goodness-of-fit index of the initial cylindrical hypothesis.

[0169] The initial cylinder assumption with the smallest goodness-of-fit index is selected as the optimal initial solution of the base cylinder model.

[0170] Based on the optimal initial solution, a distance threshold is set, and points in the normalized point cloud whose distance to the surface of the base cylindrical model is less than the distance threshold are determined as inner points, and points whose distance is greater than or equal to the distance threshold are determined as outer points.

[0171] Based on the point set consisting of all interior points, the least squares method is used for iterative optimization to solve for the parameters of the base cylinder model.

[0172] Calculate the inlier rate, which is the ratio of the number of inliers to the total number of points in the normalized point cloud.

[0173] In one embodiment, the initial residual mesh acquisition module 12 is further configured to:

[0174] The calculation of the normal distance from each point in the normalized point cloud to the surface of the base cylindrical model is defined as the geometric residual; including:

[0175] Based on the parameters of the base cylindrical model, the shortest distance vector from each normalized point to the surface of the base cylindrical model is calculated using analytical geometry methods.

[0176] Verify whether the shortest distance vector is consistent with the normal vector direction of the projection of the base cylindrical model surface at the corresponding point. If they are consistent, confirm that it is the true normal distance and use the absolute value of the true normal distance as the geometric residual. If they are inconsistent, mark the corresponding point as an invalid point and discard it.

[0177] The sign of the normal distance is also recorded, including both positive and negative signs.

[0178] Specifically, the cylindrical surface of the base cylindrical model is unfolded onto a two-dimensional parametric plane, and each point and its corresponding geometric residual are mapped onto this two-dimensional parametric plane to form a sparse residual distribution. The sparse residual distribution is then interpolated to generate a continuous initial residual mesh field, including:

[0179] Establish a two-dimensional parametric plane, where the U-axis of the two-dimensional parametric plane represents the circumferential angle of the cylindrical surface, and the V-axis represents the axial height of the cylindrical surface;

[0180] For each point in the normalized point cloud, calculate the axial height and circumferential angle of its projection point on the surface of the base cylindrical model to determine the coordinate position of the corresponding point on the two-dimensional parameter plane; at the same time, use the geometric residual of the corresponding point as the attribute value of the coordinate position.

[0181] All points in the normalized point cloud are sequentially mapped to the two-dimensional parameter plane to form a discrete set of points with attribute values ​​of geometric residuals, constituting the sparse residual distribution;

[0182] Generate a two-dimensional parametric planar regular mesh based on preset rules;

[0183] A preset interpolation algorithm is used to estimate the residual value at each grid node of the two-dimensional parametric plane regular grid, thereby generating the initial residual grid field that covers the entire parameter domain.

[0184] The generation of a two-dimensional parametric planar regular mesh based on preset rules includes:

[0185] Based on the height and circumference of the base cylindrical model, a rectangular computational domain covering the entire cylindrical surface is predefined on the two-dimensional parameter plane;

[0186] Within the rectangular computational domain, the rectangular computational domain is divided along the U-axis and V-axis at preset intervals to generate the two-dimensional parametric planar regular mesh.

[0187] In one embodiment, the boundary condition acquisition module 13 is further configured to:

[0188] The internal friction angle and maximum particle size of the roadbed fill material are obtained as physical and mechanical parameters.

[0189] The tangent of the internal friction angle is set as the maximum allowable surface slope ratio. Based on the radius of the base cylindrical model and the maximum surface slope ratio, the maximum allowable change in geometric residual between any adjacent mesh nodes is calculated in the two-dimensional parametric plane unfolded on the cylindrical surface. The maximum surface slope constraint is constructed based on the maximum change in geometric residual.

[0190] Based on the maximum particle size, a minimum surface feature wavelength is set, and a spatial low-pass filter corresponding to the minimum surface feature wavelength is established in the two-dimensional parameter plane. By constraining the change amplitude of the initial residual grid field after passing through the spatial low-pass filter, the minimum surface undulation wavelength constraint is constructed.

[0191] In one embodiment, the residual mesh optimization module 14 is further configured to:

[0192] The optimization function includes two weighted sums, where the first term is a smoothing term, which consists of the second-order difference modulus of the grid node values, and the second term is a data fitting term, which consists of the mean square error of the interpolation residuals at the grid nodes and the geometric residuals of the corresponding grid nodes.

[0193] The maximum surface slope constraint and the minimum surface undulation wavelength constraint are respectively transformed into upper limit constraints on the gradient magnitude of the grid nodes and upper limit constraints on the second-order difference magnitude of the grid nodes.

[0194] The constrained optimization problem consisting of the optimization function and the boundary conditions is transformed into a Lagrange dual problem for solution, ultimately yielding the optimized physically compliant residual mesh field and the overall smoothness evaluation result of the physically compliant residual mesh field.

[0195] In one embodiment, the final volume acquisition module 15 is further configured to:

[0196] Specifically, the physically compliant residual mesh field is mapped back to three-dimensional space and superimposed onto the surface of the matrix cylindrical model to generate the final three-dimensional composite curved surface triangular mesh model. The calculated volume of the roadbed compaction test tunnel is then obtained using a numerical integration method, including:

[0197] Each grid node of the physically compliant residual grid field on the two-dimensional parameter plane is back-calculated to its corresponding three-dimensional spatial coordinates on the surface of the base cylindrical model, based on its coordinates on the two-dimensional parameter plane.

[0198] By translating the residual value corresponding to each mesh node along the direction of its normal vector at the corresponding point on the surface of the base cylindrical model, a three-dimensional point set of the composite surface is obtained.

[0199] The three-dimensional point set of the composite surface is subjected to Delaunay triangulation to generate a three-dimensional composite surface triangular mesh model representing the inner wall surface of the test hole.

[0200] Along the axial direction of the base cylindrical model, the three-dimensional composite surface triangular mesh model is sliced ​​at equal intervals with a preset spacing;

[0201] Calculate the area of ​​the polygonal cross-section formed by the intersection of each slice plane and the triangular mesh model;

[0202] The Simpson integral method is used to numerically integrate the areas of all polygonal cross sections along the axial direction to obtain the calculated volume of the roadbed compaction test hole.

[0203] The confidence level of the calculated volume is corrected using the interior point ratio, and the final volume value and corresponding uncertainty range are output, including:

[0204] Based on the inlier rate, a preset confidence correction coefficient table is queried to obtain the volume correction coefficient corresponding to the current inlier rate; wherein, the confidence correction coefficient table defines a monotonically increasing mapping relationship between the inlier rate and correction coefficients greater than 0 and less than or equal to 1.

[0205] Multiply the calculated volume by the volume correction factor to obtain the preliminary corrected volume value;

[0206] Based on the current in-point rate and the overall smoothness evaluation results of the physical compliance residual grid field, a comprehensive uncertainty factor is calculated; the comprehensive uncertainty factor is negatively correlated with the in-point rate and negatively correlated with the overall smoothness of the physical compliance residual grid field.

[0207] Based on the comprehensive uncertainty factor and the preliminary corrected volume value, the uncertainty range corresponding to the final volume value is calculated and output, wherein the width of the uncertainty range is positively correlated with the comprehensive uncertainty factor.

[0208] It should be noted that the above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the protection scope of the present invention.

Claims

1. A rapid detection method for roadbed compaction degree based on three-dimensional point cloud technology, characterized in that, The method includes: Obtain the normalized point cloud of the inner wall of the roadbed compaction test tunnel; based on the macroscopic prior knowledge formed by the roadbed compaction test tunnel drilling tool, the normalized point cloud is fitted using a fitting process that combines the random sampling consensus algorithm and the minimum median flat method to solve the parameters of the matrix cylindrical model representing the macroscopic shape of the roadbed compaction test tunnel, and calculate the inner point ratio. Calculate the normal distance from each point in the normalized point cloud to the surface of the base cylindrical model, and define it as the geometric residual; unfold the cylindrical surface of the base cylindrical model to a two-dimensional parameter plane, map each point and its corresponding geometric residual to the two-dimensional parameter plane to form a sparse residual distribution, and interpolate the sparse residual distribution to generate a continuous initial residual grid field. Obtain the physical and mechanical parameters of the roadbed fill material, and derive the maximum surface slope constraint and the minimum surface undulation wavelength constraint based on the physical and mechanical parameters as boundary conditions. The physical and mechanical parameters include at least the internal friction angle and the maximum particle size. Based on the initial residual mesh field, and with the boundary conditions as constraints, an optimization function is established with the smoothness of the residual field and the goodness of fit to the measurement data as objectives. The optimization function is solved to obtain the optimized physical compliance residual mesh field. The physical compliance residual mesh field is mapped back to three-dimensional space and superimposed on the surface of the matrix cylindrical model to generate the final three-dimensional composite surface triangular mesh model. The calculated volume of the roadbed compaction test hole is obtained by numerical integration. The confidence level of the calculated volume is corrected by the interior point ratio, and the final volume value and the corresponding uncertainty range are output. The performance parameters of the roadbed fill material are obtained, and a compaction degree calculation formula is built based on the performance parameters. The final volume value is substituted into the formula to calculate the roadbed compaction degree. The performance parameters include humidity, density, and moisture content.

2. The method for rapid detection of roadbed compaction based on three-dimensional point cloud technology according to claim 1, characterized in that, Based on the macroscopic prior knowledge formed by the roadbed compaction test tunnel drilling tool, the normalized point cloud is fitted using a process that integrates the random sampling consensus algorithm and the minimum median flat method to solve for the parameters of the matrix cylindrical model representing the macroscopic shape of the roadbed compaction test tunnel, and the in-point ratio is calculated, including: Based on the parameters of the test hole drilling tool for roadbed compaction, the macroscopic prior knowledge is obtained, wherein the macroscopic prior knowledge includes at least the test hole diameter range, the approximate cylindrical shape, and the longitudinal axis direction of the test hole; Based on the aforementioned macroscopic prior knowledge, a reasonable sampling range for the cylinder radius is set, and multiple initial cylinder hypotheses are generated by randomly selecting the minimum point set in the normalized point cloud. For each initial cylindrical hypothesis, the distances from all points to the surface of the cylindrical hypothesis are calculated to obtain a set of distances. The median of the set of distances is selected as the goodness-of-fit index of the initial cylindrical hypothesis. The initial cylinder assumption with the smallest goodness-of-fit index is selected as the optimal initial solution of the base cylinder model. Based on the optimal initial solution, a distance threshold is set, and points in the normalized point cloud whose distance to the surface of the base cylindrical model is less than the distance threshold are determined as inner points, and points whose distance is greater than or equal to the distance threshold are determined as outer points. Based on the point set consisting of all interior points, the least squares method is used for iterative optimization to solve for the parameters of the base cylinder model. Calculate the inlier rate, which is the ratio of the number of inliers to the total number of points in the normalized point cloud.

3. The method for rapid detection of roadbed compaction based on three-dimensional point cloud technology according to claim 1, characterized in that, Calculate the normal distance from each point in the normalized point cloud to the surface of the base cylindrical model, defined as the geometric residual; including: Based on the parameters of the base cylindrical model, the shortest distance vector from each normalized point to the surface of the base cylindrical model is calculated using analytical geometry methods. Verify whether the shortest distance vector is consistent with the normal vector direction of the projection of the base cylindrical model surface at the corresponding point. If they are consistent, confirm that it is the true normal distance and use the absolute value of the true normal distance as the geometric residual. If they are inconsistent, mark the corresponding point as an invalid point and discard it. The sign of the normal distance is also recorded, including both positive and negative signs.

4. The method for rapid detection of roadbed compaction based on three-dimensional point cloud technology according to claim 1, characterized in that, The cylindrical surface of the base cylindrical model is unfolded onto a two-dimensional parametric plane. Each point and its corresponding geometric residual are mapped onto this two-dimensional parametric plane to form a sparse residual distribution. The sparse residual distribution is then interpolated to generate a continuous initial residual mesh field, including: Establish a two-dimensional parametric plane, where the U-axis of the two-dimensional parametric plane represents the circumferential angle of the cylindrical surface, and the V-axis represents the axial height of the cylindrical surface; For each point in the normalized point cloud, calculate the axial height and circumferential angle of its projection point on the surface of the base cylindrical model to determine the coordinate position of the corresponding point on the two-dimensional parameter plane; at the same time, use the geometric residual of the corresponding point as the attribute value of the coordinate position. All points in the normalized point cloud are sequentially mapped to the two-dimensional parameter plane to form a discrete set of points with attribute values ​​of geometric residuals, constituting the sparse residual distribution; Generate a two-dimensional parametric planar regular mesh based on preset rules; A preset interpolation algorithm is used to estimate the residual value at each grid node of the two-dimensional parametric plane regular grid, thereby generating the initial residual grid field that covers the entire parameter domain.

5. The rapid detection method for roadbed compaction based on three-dimensional point cloud technology according to claim 4, characterized in that, Generate a two-dimensional parametric planar regular mesh based on preset rules, including: Based on the height and circumference of the base cylindrical model, a rectangular computational domain covering the entire cylindrical surface is predefined on the two-dimensional parameter plane; Within the rectangular computational domain, the rectangular computational domain is divided along the U-axis and V-axis at preset intervals to generate the two-dimensional parametric planar regular mesh.

6. The method for rapid detection of roadbed compaction based on three-dimensional point cloud technology according to claim 1, characterized in that, Obtain the physical and mechanical parameters of the roadbed fill material, and derive the maximum surface slope constraint and minimum surface undulation wavelength constraint based on these parameters as boundary conditions, including: The internal friction angle and maximum particle size of the roadbed fill material are obtained as physical and mechanical parameters. The tangent of the internal friction angle is set as the maximum allowable surface slope ratio. Based on the radius of the base cylindrical model and the maximum surface slope ratio, the maximum allowable change in geometric residual between any adjacent mesh nodes is calculated in the two-dimensional parametric plane unfolded on the cylindrical surface. The maximum surface slope constraint is constructed based on the maximum change in geometric residual. Based on the maximum particle size, a minimum surface feature wavelength is set, and a spatial low-pass filter corresponding to the minimum surface feature wavelength is established in the two-dimensional parameter plane. By constraining the change amplitude of the initial residual grid field after passing through the spatial low-pass filter, the minimum surface undulation wavelength constraint is constructed.

7. The rapid detection method for roadbed compaction based on three-dimensional point cloud technology according to claim 6, characterized in that, Based on the initial residual mesh field, and constrained by the boundary conditions, an optimization function is established with the smoothness of the residual field and the goodness of fit to the measurement data as objectives. Solving the optimization function yields the optimized physically compliant residual mesh field, including: The optimization function includes two weighted sums, where the first term is a smoothing term, which consists of the second-order difference modulus of the grid node values, and the second term is a data fitting term, which consists of the mean square error of the interpolation residuals at the grid nodes and the geometric residuals of the corresponding grid nodes. The maximum surface slope constraint and the minimum surface undulation wavelength constraint are respectively transformed into upper limit constraints on the gradient magnitude of the grid nodes and upper limit constraints on the second-order difference magnitude of the grid nodes. The constrained optimization problem consisting of the optimization function and the boundary conditions is transformed into a Lagrange dual problem for solution, ultimately yielding the optimized physically compliant residual mesh field and the overall smoothness evaluation result of the physically compliant residual mesh field.

8. The method for rapid detection of roadbed compaction based on three-dimensional point cloud technology according to claim 1, characterized in that, The physically compliant residual mesh field is mapped back to three-dimensional space and superimposed onto the surface of the matrix cylindrical model to generate the final three-dimensional composite curved surface triangular mesh model. The calculated volume of the roadbed compaction test tunnel is then obtained using a numerical integration method, including: Each grid node of the physically compliant residual grid field on the two-dimensional parameter plane is back-calculated to its corresponding three-dimensional spatial coordinates on the surface of the base cylindrical model, based on its coordinates on the two-dimensional parameter plane. By translating the residual value corresponding to each mesh node along the direction of its normal vector at the corresponding point on the surface of the base cylindrical model, a three-dimensional point set of the composite surface is obtained. The three-dimensional point set of the composite surface is subjected to Delaunay triangulation to generate a three-dimensional composite surface triangular mesh model representing the inner wall surface of the test hole. Along the axial direction of the base cylindrical model, the three-dimensional composite surface triangular mesh model is sliced ​​at equal intervals with a preset spacing; Calculate the area of ​​the polygonal cross-section formed by the intersection of each slice plane and the triangular mesh model; The Simpson integral method is used to numerically integrate the areas of all polygonal cross sections along the axial direction to obtain the calculated volume of the roadbed compaction test hole.

9. The method for rapid detection of roadbed compaction based on three-dimensional point cloud technology according to claim 7, characterized in that, The calculated volume is corrected for confidence level using the interior point ratio, and the final volume value and corresponding uncertainty range are output, including: Based on the inlier rate, a preset confidence correction coefficient table is queried to obtain the volume correction coefficient corresponding to the current inlier rate; wherein, the confidence correction coefficient table defines a monotonically increasing mapping relationship between the inlier rate and correction coefficients greater than 0 and less than or equal to 1. Multiply the calculated volume by the volume correction factor to obtain the preliminary corrected volume value; Based on the current in-point rate and the overall smoothness evaluation results of the physical compliance residual grid field, a comprehensive uncertainty factor is calculated; the comprehensive uncertainty factor is negatively correlated with the in-point rate and negatively correlated with the overall smoothness of the physical compliance residual grid field. Based on the comprehensive uncertainty factor and the preliminary corrected volume value, the uncertainty range corresponding to the final volume value is calculated and output, wherein the width of the uncertainty range is positively correlated with the comprehensive uncertainty factor.

10. A rapid detection system for roadbed compaction based on three-dimensional point cloud technology, characterized in that, The system is used to implement the rapid detection method for roadbed compaction based on three-dimensional point cloud technology as described in any one of claims 1-9. The system includes: The matrix cylindrical model acquisition module is used to acquire the normalized point cloud of the inner wall of the roadbed compaction test tunnel; based on the macroscopic prior knowledge formed by the roadbed compaction test tunnel drilling tool, the normalized point cloud is fitted using a fitting process that combines the random sampling consensus algorithm and the minimum median flat method to solve for the matrix cylindrical model parameters representing the macroscopic shape of the roadbed compaction test tunnel, and to calculate the inner point ratio. The initial residual mesh acquisition module is used to calculate the normal distance from each point in the normalized point cloud to the surface of the base cylindrical model, which is defined as the geometric residual; unfold the cylindrical surface of the base cylindrical model to a two-dimensional parameter plane, map each point and its corresponding geometric residual to the two-dimensional parameter plane to form a sparse residual distribution, and interpolate the sparse residual distribution to generate a continuous initial residual mesh field; The optimized boundary condition acquisition module is used to acquire the physical and mechanical parameters of the roadbed fill material, and derive the maximum surface slope constraint and the minimum surface undulation wavelength constraint based on the physical and mechanical parameters as boundary conditions. The physical and mechanical parameters include at least the internal friction angle and the maximum particle size. The residual mesh optimization module is used to establish an optimization function based on the initial residual mesh field and the boundary conditions, with the smoothness of the residual field and the goodness of fit to the measurement data as objectives, solve the optimization function, and obtain the optimized physical compliance residual mesh field. The final volume acquisition module is used to map the physical compliance residual mesh field back to three-dimensional space, superimpose it onto the surface of the matrix cylindrical model, generate the final three-dimensional composite surface triangular mesh model, and use the numerical integration method to obtain the calculated volume of the roadbed compaction test hole; use the interior point ratio to correct the confidence level of the calculated volume, and output the final volume value and the corresponding uncertainty range; The subgrade compaction degree calculation module is used to obtain the performance parameters of the subgrade fill material, build a compaction degree calculation formula based on the performance parameters, and substitute the final volume value to calculate the subgrade compaction degree; wherein, the performance parameters include humidity, density, and moisture content.

Images