A tunnel point cloud adaptive elliptic cylinder filtering method based on affine constraint

Abstract

This invention provides an adaptive elliptical cylinder filtering method for tunnel point clouds based on affine constraints, belonging to the field of point cloud data processing technology. The method includes: extracting the tunnel centerline and setting cross-sections; initializing the elliptical cylinder model parameters using principal component analysis; setting an initial screening threshold and dynamically estimating the optimal threshold using the RANSAC algorithm; constructing a multi-objective optimization function containing data fitting error terms, parameter smoothing constraint terms, and affine transformation constraint terms to optimize the model parameters; and filtering the point cloud data based on the optimized model and the dynamic threshold. This invention innovatively introduces affine transformation constraints and adaptive parameter constraint mechanisms, dynamically adjusting the constraint threshold according to local curvature and deformation degree, and developing specialized identification and processing mechanisms for special regions. This significantly improves the geometric continuity of the model and its adaptability to complex environments, improves the separation effect between noise points and structural points, and provides a reliable data foundation for tunnel 3D reconstruction and deformation analysis.

Description

Technical Field

[0001] This invention relates to the field of point cloud data processing technology, and specifically to an adaptive elliptic cylindrical filtering method for tunnel point clouds based on affine constraints. Background Technology

[0002] In the field of tunnel point cloud data processing, accurate filtering is a crucial prerequisite for achieving 3D tunnel reconstruction and deformation monitoring. With the widespread application of 3D laser scanning technology in tunnel engineering, how to extract effective information from high-noise, high-density tunnel point cloud data and accurately identify tunnel structural surfaces has become a research hotspot. However, existing technologies still have many limitations when dealing with complex tunnel environments.

[0003] Existing methods for processing tunnel point clouds mainly include statistical analysis-based filtering methods, geometric feature extraction methods, and model fitting methods. Statistical analysis-based filtering methods, such as the Statistical Outlier Removal (SOR) model, typically use fixed parameters to identify and remove noise points. However, these methods are sensitive to parameter selection and struggle to adapt to dynamic changes in the shape and size of the tunnel cross-section, easily leading to under-filtering or over-filtering. Especially when dealing with localized tunnel deformations or special structures (such as refuges and cross passages), fixed-parameter methods often fail to accurately identify structural features, resulting in the loss of useful information or misjudgment.

[0004] Geometric feature extraction methods, including those based on normal vector analysis and local curvature estimation, can capture the geometric characteristics of the tunnel surface well, but their computational efficiency is low when processing high-noise, high-density point cloud data. Furthermore, these methods require high uniformity of point cloud density; in actual tunnel point cloud data, density non-uniformity caused by factors such as occlusion and scanning angle often affects the accuracy of feature extraction.

[0005] Model fitting methods, such as the Random Sample Consensus (RANSAC) algorithm and the Hough transform, identify tunnel structures by fitting predefined geometric models (such as cylinders or elliptical cylinders). These methods perform well under ideal conditions, but when faced with complex tunnel structures (such as twin-track tunnels or tunnels with varying cross-sections), a single model struggles to accurately describe the entire tunnel's geometry. Furthermore, these methods typically employ global optimization strategies, resulting in high computational complexity and making them unsuitable for real-time processing of large-scale point cloud data.

[0006] In existing technologies, such as Chinese patent application CN119090728A, an adaptive elliptic cylinder filtering method based on tunnel point cloud data is disclosed. This method combines the RANSAC algorithm with dynamic threshold estimation and uses PCA initialization and iterative least squares to optimize the elliptic cylinder model parameters. However, this method still has the following shortcomings: On the one hand, it lacks an effective affine transformation constraint mechanism, making it difficult to guarantee the rationality and continuity of parameter changes between adjacent sections. Traditional methods often use fixed thresholds or simple smoothing for constraints, which cannot adapt to the characteristics of curvature changes in different regions. Especially in the curved and deformed sections of tunnels, fixed threshold constraints are prone to over-constraint or under-constraint problems.

[0007] On the other hand, existing technologies for optimizing model parameters often employ a single error metric (such as the distance from a point to the model), lacking multi-objective optimization strategies that consider global geometric continuity. This results in a lack of reasonable variation patterns in model parameters along the longitudinal direction when dealing with long tunnels or tunnels under complex geological conditions, affecting the overall accuracy and reliability of the model.

[0008] Furthermore, existing technologies lack specialized identification and processing mechanisms for special regions such as portal sections, making it difficult to cope with changes in the geometric characteristics of these areas. At the same time, existing methods rely on limited approaches to smooth parameter transitions, failing to effectively handle abrupt changes in cross-sectional parameters and impacting the stability and reliability of the model.

[0009] In summary, existing technologies lack a comprehensive solution that fully considers model fitting accuracy, geometric continuity, local characteristic adaptability, and computational efficiency when processing tunnel point cloud data. There is an urgent need for a new method that can take all these aspects into account to improve the overall performance of tunnel point cloud data processing and provide a more reliable data foundation for tunnel 3D reconstruction and deformation monitoring. Summary of the Invention

[0010] The purpose of this invention is to provide an adaptive elliptic cylindrical filtering method for tunnel point clouds based on affine constraints, which solves the problem that existing technologies lack comprehensive consideration of model fitting accuracy, geometric continuity, and local characteristic adaptability when processing tunnel point cloud data, and provides a reliable data foundation for tunnel 3D reconstruction and deformation monitoring.

[0011] To achieve the above objectives, the present invention provides the following technical solution:

[0012] An adaptive elliptic cylindrical filtering method for tunnel point clouds based on affine constraints includes the following steps:

[0013] S1: Extract the tunnel centerline based on tunnel point cloud data, and set cross-sections at fixed intervals along the centerline axis;

[0014] S2: Initialize the parameters of the elliptical cylinder model using principal component analysis on each cross section to form the elliptical cylinder model;

[0015] S3: Set an initial screening threshold to initially screen outliers in the distance from all points within the cross section to the elliptical cylinder, and use the RANSAC algorithm to dynamically iteratively estimate the optimal screening threshold.

[0016] S4: Construct a multi-objective optimization function that includes data fitting error terms, parameter smoothing constraint terms, and affine transformation constraint terms, and use an optimization algorithm to optimize the parameters of the elliptical cylinder model until the parameters of the elliptical cylinder model converge.

[0017] S5: Based on the optimized elliptical cylinder model and dynamic threshold, the tunnel point cloud data is filtered, and the filtered point cloud data and elliptical cylinder model parameters are output.

[0018] Further: In step S1, the tunnel centerline is extracted from the preliminarily denoised tunnel point cloud data using the least squares method or the master curve algorithm, with a fixed interval of 0.5-3 meters.

[0019] Further: In step S2, the initialization of the elliptical cylinder model parameters using the principal component analysis method specifically includes: calculating the covariance matrix of the cross-sectional point cloud, performing eigenvalue decomposition on the covariance matrix; determining the directions of the major and minor axes of the elliptical cylinder model using the eigenvectors corresponding to the largest and second largest eigenvalues, and estimating the lengths of the major and minor axes based on the proportions of the eigenvalues.

[0020] Further: In step S3, the initial screening threshold is 20 cm, and the optimal screening threshold is the 95th percentile of the distance between interior points.

[0021] Further: In step S4, the multi-objective optimization function is ,in Let be the sum of the squares of the distances from the point cloud to the elliptical cylinder. It is the sum of squares of the differences in parameters between adjacent sections. λ1 and λ2 are the sum of squares of the differences between the affine transformation matrix and the desired transformation matrix between adjacent sections, and are weighting coefficients.

[0022] Furthermore, the feature is that the Levenberg-Marquardt algorithm is used to minimize the multi-objective optimization function, with an initial damping factor of 0.01 and a damping factor adjustment factor of 10.

[0023] Furthermore, the weighting coefficient λ1 ranges from 0.5 to 2.0, and the weighting coefficient λ2 ranges from 0.3 to 1.5.

[0024] Furthermore, step S4 also includes a step of dynamically adjusting the constraint threshold based on the local curvature and the degree of deformation, specifically including: dynamically calculating the orientation angle change threshold based on the local curvature, with a basic threshold of 2° and a maximum threshold of 8°; and dynamically calculating the axis length change threshold based on the degree of local deformation, with a basic threshold of 2% and a maximum threshold of 6%.

[0025] Furthermore, step S4 also includes a step of smoothing the affine transformation parameters, specifically including: calculating the affine transformation matrix and decomposing it into rotation, scaling and translation components; smoothing the transformation parameters using Gaussian weights; and applying cubic spline interpolation to achieve smooth parameter transition.

[0026] Furthermore: the curvature change rate threshold is set to 0.2, the optimization window for local re-optimization is set to 5 sections, the confidence level for abnormal change identification is set to 97%, and the weight decay rate for the local weighted average method is set to 0.7.

[0027] Further: In step S5, points within the optimal screening threshold range are retained as tunnel structure points, and points exceeding the optimal screening threshold are removed as unstructured points or noise.

[0028] Furthermore: the expected transformation matrix of the affine transformation constraint term is constructed according to the tunnel segment type, including: the expected transformation matrix of the straight line segment. ,in It is the identity matrix. This is the disturbance coefficient, with a value ranging from 0.01 to 0.05. The perturbation matrix; the desired transformation matrix of the curved segment contains rotational components. According to local curvature Construct; the expected transformation matrix of the deformed segment contains scaling components. Constructed based on the deformation rates of adjacent sections; the expected transformation matrix of the composite segment. .

[0029] Furthermore, it also includes a multi-level verification and correction mechanism, specifically including: physical feasibility verification to check whether the transformed ellipse meets the axis length variation limit. and shaft ratio limit The conflict constraint coordination mechanism determines the conflict when there is a conflict between data fitting and affine transformation constraints. The conflict criterion is as follows: Prioritize ensuring data fitting accuracy; the anomaly detection and handling mechanism marks and locally re-optimizes cases where abnormal changes occur in the model parameters after constraints.

[0030] Furthermore, for straight tunnel segments exceeding 50 meters in length, principal directional stability constraints are also included, calculated through cumulative offset. and constraints ,in The length of the line segment is given, and a direction accumulation penalty term is added to the multi-objective optimization function. ,in This is the weighting coefficient, with a value ranging from 0.2 to 0.5.

[0031] Compared with the prior art, the present invention has the following advantages:

[0032] I. The method of this invention significantly improves the geometric continuity of the model by introducing affine transformation constraints and parameter smoothing transition processing. Affine transformation constraints ensure the physical rationality of parameter changes between adjacent sections; the adaptive parameter change threshold is dynamically adjusted according to the local curvature and deformation degree, ensuring deformation detection sensitivity while avoiding over-constraint; Gaussian weighted smoothing and cubic spline interpolation techniques ensure smooth parameter transition in the longitudinal direction, making the final model more realistically reflect the actual geometry of the tunnel.

[0033] Second, for areas with significant changes in geometric characteristics, such as tunnel portals, the specialized identification and processing mechanisms developed in this invention significantly improve the model's adaptability. Special regions are identified through curvature change rate analysis; abrupt changes are addressed through local re-optimization; abnormal changes are identified through a statistical model of parameter variations; and a local weighted averaging method ensures smooth transitions. These mechanisms work together to enable the model to accurately adapt to various special structures within the tunnel, improving the global accuracy of point cloud filtering.

[0034] Third, this invention constructs a comprehensive optimization objective function that includes data fitting error, parameter smoothness constraints, and affine transformation constraints, ensuring geometric continuity while maintaining fitting accuracy. Compared with traditional methods, this invention comprehensively considers data fitting accuracy, parameter smoothness, and geometric continuity through a multi-objective optimization function, avoiding the local optimum problem caused by focusing only on the distance from the point to the model, and ensuring the global optimum of the model parameters.

[0035] Fourth, this invention provides a method for dynamically adjusting constraint thresholds based on local curvature and deformation, enabling the change thresholds for orientation angle and axis length to adaptively adjust according to the local characteristics of the tunnel, thus solving the problem of poor adaptability of fixed threshold constraints in different regions. This dynamic parameter constraint mechanism can more accurately capture the local features of the tunnel, avoiding overfitting or underfitting while ensuring fitting accuracy.

[0036] V. This invention significantly improves the separation effect between noise points and structural points through dynamic thresholding and multi-stage optimization. The initial threshold setting provides a reasonable starting separation standard; the threshold is dynamically adjusted during RANSAC to adapt to the characteristics of different regions; multi-objective optimization improves the model fitting accuracy; the final filtering stage uses the optimized model and dynamic threshold for accurate separation, greatly improving the accuracy of structural point identification and providing a high-quality data foundation for subsequent tunnel 3D reconstruction and deformation analysis. Attached Figure Description

[0037] Figure 1 This is a flowchart illustrating an adaptive elliptical cylindrical filtering method for tunnel point clouds based on affine constraints, as described in one embodiment of the present invention. Detailed Implementation

[0038] The technical solution of the present invention will now be clearly and completely described with reference to the accompanying drawings. Obviously, the described embodiments are only some, not all, of the embodiments of the present invention. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0039] In the description of this invention, it should be noted that the terms "center," "upper," "lower," "left," "right," "vertical," "horizontal," "inner," and "outer," etc., indicate the orientation or positional relationship based on the orientation or positional relationship shown in the accompanying drawings. They are used only for the convenience of describing the invention and for simplifying the description, and do not indicate or imply that the device or element referred to must have a specific orientation, or be constructed and operated in a specific orientation. Therefore, they should not be construed as limitations on the invention. Furthermore, the terms "first," "second," and "third" are used for descriptive purposes only and should not be construed as indicating or implying relative importance.

[0040] This invention provides an adaptive elliptic cylindrical filtering method for tunnel point clouds based on affine constraints. This method achieves high-precision point cloud filtering in complex tunnel environments through dynamic parameter adjustment, multi-objective optimization, and special location processing. The implementation of this invention is described in detail below:

[0041] An adaptive elliptic cylindrical filtering method for tunnel point clouds based on affine constraints includes the following steps:

[0042] S1: Extract the tunnel centerline from the preliminary denoised point cloud using the least squares method or the master curve algorithm, and set cross sections at fixed intervals along the centerline axis.

[0043] In practical applications, the initial noise reduction threshold can be set to 2-5 times the standard deviation, preferably 3 times, to ensure the initial removal of obvious outliers. Using the extracted centerline as the reference axis, cross-sections are set along the axial direction at fixed intervals, ranging from 0.5-3m, preferably 1-2m. During implementation, the interval size can be adjusted according to the specific characteristics of the tunnel. For example, a larger interval (approximately 1.5-2m) can be used in straight sections of the tunnel; in curved sections or deformation areas, the interval should be appropriately reduced (approximately 1-1.5m) to improve accuracy. It is important to note that the cross-section normal should be kept as perpendicular to the tunnel direction as possible to ensure the accuracy and continuity of the model parameters.

[0044] S2: On each cross section, the parameters of the elliptical cylinder model are initialized using the principal component analysis method to form the elliptical cylinder model.

[0045] Principal component analysis (PCA) can identify the main distribution directions of a point cloud, thereby estimating the major and minor axis ratios and orientations of the cross-section. The specific steps include calculating the covariance matrix of the point cloud cross-section, performing eigenvalue decomposition on the covariance matrix, using the eigenvectors corresponding to the largest and second-largest eigenvalues ​​to determine the major and minor axis orientations of the ellipse, and estimating the lengths of the major and minor axes based on the proportions of the eigenvalues. In practical applications, the eigenvalue screening threshold can be set to 5% of the largest eigenvalue, the principal axis orientation tolerance to be 2°, and the outlier removal threshold to be 3σ. For cases of non-uniform point cloud density, a weight compensation mechanism can be used to ensure that the initialization parameters are reasonable and reliable. The weight compensation mechanism can be calculated based on the local density of the point cloud, assigning higher weights to points in lower-density regions and lower weights to points in higher-density regions, thus balancing the contributions of different regions to the model fitting.

[0046] S3: Set an initial screening threshold to initially screen outliers in the distance from all points within the cross section to the elliptical cylinder, and use the RANSAC algorithm to dynamically iteratively estimate the optimal screening threshold.

[0047] Specifically, after initialization, the distance from all points within the cross-section to the elliptical cylinder is calculated. An initial screening threshold d = 20 cm is set for preliminary screening of outliers. The RANSAC algorithm is used to dynamically estimate the optimal threshold to adapt to the point cloud characteristics of different regions. RANSAC iterates 50-500 times, each time randomly selecting a minimum sample set of 5-15 points to fit the elliptical model, calculating the distance from all points to the model, and counting the number of inliers. The model with the most inliers is selected as the best fit result, and the threshold d is dynamically adjusted based on the best fit result, taking the 95th percentile of the inlier distance. This dynamic threshold strategy significantly improves the method's adaptability to changes in local characteristics. In practical applications, the number of RANSAC iterations can be appropriately adjusted according to the point cloud quality and tunnel characteristics; the number of iterations can be increased in areas with high noise and decreased in areas with low noise, to balance accuracy and efficiency.

[0048] S4: Construct a multi-objective optimization function that includes data fitting error terms, parameter smoothing constraint terms, and affine transformation constraint terms, and use an optimization algorithm to optimize the parameters of the elliptical cylinder model until the parameters of the elliptical cylinder model converge.

[0049] S401: Construction of Multi-Objective Optimization Function

[0050] The core of this invention lies in constructing a multi-objective optimization function to achieve comprehensive optimization of model parameters, the mathematical expression of which is:

[0051]

[0052] in, The sum of squared distances from the point cloud to the elliptical cylinder, to ensure data fitting accuracy, is calculated using the following formula:

[0053]

[0054] in, Indicates the first j On the first cross section i Each point corresponds to an elliptical cylinder model. The distance;

[0055] The sum of squares of the differences in parameters between adjacent sections, ensuring a smooth transition of parameters, is calculated using the following formula:

[0056]

[0057] Among them, They represent the first j The orientation angle, major axis length, and minor axis length of the cross-sectional ellipse;

[0058] The sum of squares of the differences between the affine transformation matrix and the desired transformation matrix between adjacent sections, ensuring geometric continuity, is calculated using the following formula:

[0059]

[0060] Among them, Indicates from the first j 1 The section to the first j The actual affine transformation matrix of each cross section Denotes the desired affine transformation matrix. This represents the Frobenius norm.

[0061] The weighting coefficients λ1 range from 0.5 to 2.0, and λ2 range from 0.3 to 1.5. In the embodiments of the present invention, λ1 = 1.0 and λ2 = 0.8 can be taken.

[0062] S402: The calculation and decomposition process of constructing the affine transformation matrix for a multi-objective optimization function.

[0063] Affine transformation matrix calculation:

[0064] When establishing affine transformation relationships between adjacent sections, the elliptical parameters of the two adjacent sections S1 and S2 are first extracted using PCA analysis. The major axis of the ellipse of section S1 is a unit vector v1, and the minor axis is a unit vector v2; the major axis of the ellipse of section S2 is a unit vector w1, and the minor axis is a unit vector w2.

[0065] Construct a system of linear equations:

[0066] A·v1=w1

[0067] A·v2=w2

[0068] Let V = [v1, v2] and W = [w1, w2]. Since v1 and v2 are orthogonal unit vectors (the principal and secondary axes of the ellipse), matrix V is invertible. Therefore, the affine transformation matrix A = W·V is calculated. -1 .

[0069] Matrix decomposition process:

[0070] Calculate matrix A^T·A and perform eigenvalue decomposition, solve the characteristic equation det(A^T·A-λI)=0, and obtain the two eigenvalues. λ 1 and λ 2 (all are positive real numbers, because A^T·A is a positive definite symmetric matrix), construct a diagonal matrix Simultaneously, the corresponding eigenvectors q1 and q2 are obtained, and orthogonalized to form an orthogonal matrix. .

[0071] Calculate the scaling matrix , represents the scaling of the ellipse along its major and minor axes, where the scaling factor is... and These are the scaling amounts along the major and minor axes, respectively.

[0072] Calculate the rotation matrix:

[0073]

[0074] Calculate the rotation angle .

[0075] S403: Construction of the expected affine transformation matrix for multi-objective optimization functions

[0076] Expected affine transformation matrix It is an ideal transformation matrix constructed based on the local geometric and physical properties of the tunnel, and is used as a reference standard for the actual affine transformation of adjacent sections.

[0077] Expected transformation of straight segments: For straight segments of a tunnel, the expected affine transformation is close to the identity matrix I, indicating that the change in cross-sectional shape is very small.

[0078]

[0079] in I It is the identity matrix. ε For small perturbation coefficients (values ​​ranging from 0.01 to 0.05), D This is the perturbation matrix.D The introduction of this allows for minute shape variations, avoiding over-constraint.

[0080] Desired Transformation of Curved Segments: For curved segments, the desired transformation mainly manifests as a rotational transformation, constructed based on the local curvature κ.

[0081]

[0082] in α The empirical proportionality coefficient between local curvature and rotation angle is determined experimentally to be in the range of 1.2-1.8, and is dynamically adjusted for different radii of curvature.

[0083] For large radii of curvature (>500m), take α =1.2-1.4

[0084] For medium radius of curvature (200-500m), take... α =1.4-1.6

[0085] For small radius of curvature (<200m), take α =1.6-1.8

[0086] Deformation Segment Expected Transformation: For the deformed segment, the expected transformation mainly manifests as a scaling transformation, constructing the expected scaling matrix based on the deformation rates of adjacent sections.

[0087]

[0088] a and b These represent the relative rates of change of the major and minor axes, respectively. β The smoothing factor (values ​​0.8-0.95) is used to calculate the deformation rate.

[0089]

[0090]

[0091] here and This is the estimated axis length of the next section based on the current deformation trend.

[0092] Expected Transformation of Composite Segment: For a composite segment that simultaneously exhibits bending and deformation, the complete expected affine transformation matrix is ​​calculated as follows:

[0093]

[0094] The weights of each component are dynamically adjusted based on the dominant transformation type of the current paragraph.

[0095] S404: Adaptive setting of constraints for constructing multi-objective optimization functions

[0096] Local geometric property analysis: First, calculate the local curvature of the centerline. Where r(s) is the equation of the centerline parameters, and the rate of change of the elliptic parameters of adjacent sections is also calculated. and Set curvature determination threshold =0.005 and deformation judgment threshold =0.03.

[0097] Paragraph type determination:

[0098] when κ < and a < and b < When it is determined to be a straight line segment

[0099] when κ ≥ and a < and b < It was determined to be a curved section.

[0100] when κ < and( a ≥ or b ≥ When it is determined to be a deformed segment

[0101] when κ ≥ and( a ≥ or b ≥ When it is determined to be a comprehensive section

[0102] Constraint parameters are adaptively adjusted:

[0103] α Parameter adjustment formula: α =1.2+0.6·min(1, κ / 0.01)

[0104] β Parameter adjustment formula: β=0.8+0.15·min(1,max( a , b ) / 0.05)

[0105] Disturbance coefficient for straight segment: ε =0.01+0.04·max( a, b ) /

[0106] Dynamic control of constraint strength: introduction of regional adaptive weights Controlling constraint strength:

[0107]

[0108] in This indicates the degree of geometric change at the current cross-section, where k is the sensitivity coefficient (value 5-10). The constraint term is represented as .

[0109] S405: Multi-objective optimization function construction with multi-level verification and correction mechanism

[0110] Physical feasibility verification: Check whether the transformed ellipse meets the physical constraints, including:

[0111] Shaft length variation limit:

[0112] Shaft ratio limit:

[0113] If a constraint results in a violation of physical feasibility, the constraint strength is automatically reduced to [a lower value]. .

[0114] Conflict constraint coordination mechanism: In and When there are obvious conflicts in constraints, priority should be given to ensuring the accuracy of data fitting. The conflict determination criteria are as follows:

[0115]

[0116] The conflict coordination method is to temporarily reduce It is reduced to 50% of its original value and then linearly restored in subsequent iterations.

[0117] Anomaly detection and handling mechanism: When abnormal changes occur in the model parameters after constraints (relative change > 20%), they are marked. For outlier points, a local re-optimization strategy is adopted to reduce the constraint weight and increase the data fitting weight. After re-optimization, cubic spline smoothing is used to achieve a smooth transition with the front and rear sections.

[0118] S406: Multi-objective optimization function construction and optimization algorithm implementation

[0119] In the actual optimization process, the Levenberg-Marquardt algorithm is used to minimize the optimization function. The initial damping factor is 0.01, and the damping factor adjustment factor is 10. The optimization process includes calculating the Jacobian matrix and error vector, using QR decomposition to improve the stability of the solution, multiple random initializations to avoid local optima, and iteratively updating the model parameters until convergence (relative rate of change less than 0.1%) or reaching the preset number of iterations (50).

[0120] The Jacobian matrix can be calculated using numerical differentiation methods. For complex cross-sectional shapes, analytical expressions can also be considered to improve computational efficiency. The QR decomposition can utilize an improved Gram-Schmidt orthogonalization method to enhance numerical stability.

[0121] S407: Adaptive Parameter Constraint Mechanism for Constructing Multi-Objective Optimization Functions

[0122] Dynamic threshold calculation: The orientation angle change threshold is dynamically calculated based on the local curvature.

[0123]

[0124] in The baseline threshold is 2°. It is 8° (maximum threshold). The maximum curvature is taken as a reference (value 0.015).

[0125] The threshold for axis length variation is dynamically calculated based on the degree of local deformation.

[0126]

[0127] in The base threshold is 2%. The maximum threshold is 6%. As an indicator for evaluating the degree of deformation, It is 0.08.

[0128] Main directional stability constraints for long straight tunnel segments: For straight tunnel segments exceeding 50m in length, main directional stability constraints are set and calculated through cumulative offset. and constraints (Where L is the length of the line segment), add a direction accumulation penalty term to the multi-objective optimization function. ,in λ 3 is a weighting coefficient (valued between 0.2 and 0.5) to ensure the directional stability of long straight line segments.

[0129] S408: Smooth Transition Processing Mechanism for Multi-Objective Optimization Function Construction

[0130] Parameter sequence preprocessing: After optimization, to ensure a smooth transition of parameters between adjacent sections, parameter sequence preprocessing is performed first. Initialize the parameter sequence. ,in Indicates the first j Ellipse parameters for each cross section.

[0131] Outlier parameters are identified using the Median Absolute Deviation (MAD) method, by calculating... And mark as satisfied conditional These are outliers, and local median filtering is used to perform preliminary corrections on them.

[0132] Gaussian weighted smoothing: Set the window size w to 5%-10% of the total number of cross-sections (minimum 5 cross-sections), and apply Gaussian weighted smoothing to each parameter sequence. The smoothing formula is as follows:

[0133]

[0134] The weight calculation adopts , When the window exceeds the boundary, a mirror expansion method is used to supplement the boundary data. For special areas (such as portal sections), the smoothing coefficient is adaptively adjusted to 60% of the original value.

[0135] Cubic spline interpolation implementation: Constructing a cubic spline function based on smoothed parameter control points S(t) In each interval superior, S(t) Represented as a cubic polynomial:

[0136]

[0137] The following conditions must be met:

[0138] Interpolation conditions:

[0139] The first derivative is continuous:

[0140] The second derivative is continuous:

[0141] Natural boundary conditions:

[0142] By solving the tridiagonal linear equations, the polynomial coefficients for all intervals are obtained. Finally, an interpolation point is set every 10-20 cm between the original cross-sectional positions to generate a smooth and continuous parametric curve.

[0143] Change rate constraint verification: Verify whether the change rate between adjacent points in the interpolated parameter sequence meets the constraint conditions. If there is an interval that violates the constraint, adjust the control point of the interval and re-interpolate. Iterate this process until all intervals meet the constraint conditions or the maximum number of iterations (5 times) is reached, ensuring that the final generated parameter sequence meets both the smoothness requirement and the physical constraint conditions.

[0144] S409: Special Position Handling Mechanism for Constructing Multi-Objective Optimization Functions

[0145] Portal Segment Identification and Quantification Standards: A specialized identification and processing mechanism has been developed for areas with significant geometrical changes, such as portal segments. Portal segment identification employs multiple quantitative indicators:

[0146] Rate of change of curvature:

[0147] Rate of change of elliptic axis ratio:

[0148] Rate of change of cross-sectional area:

[0149] Regions that simultaneously meet these conditions are identified as areas with significant changes in geometric properties that require special handling.

[0150] Local re-optimization process: Local re-optimization is performed on the marked mutation points, with an optimization window of 5 cross-sections. The specific process is as follows:

[0151] Determine the optimization window: Using the marked point as the center, take two cross-sections before and after it to form a local region with a total of five cross-sections.

[0152] Adjusting constraint weights: The smoothing constraint weight λ1 is adjusted to 0.7 times its original value, and the affine constraint weight λ2 is adjusted to 0.5 times its original value, increasing the model's flexibility.

[0153] Construct a local optimization objective function:

[0154] Optimization Solution: The objective function is optimized using the Levenberg-Marquardt algorithm.

[0155] Parameter update: Only update the parameters of the marker point (i.e., the middle section of the window), keeping the parameters of other sections within the window unchanged.

[0156] Statistical Model Anomaly Identification: A statistical model for parameter changes is established to identify anomalies, with a confidence level of 97%. The specific implementation is as follows:

[0157] Calculate the rate of change of adjacent sections for each parameter.

[0158] Calculate the mean and standard deviation of the rate of change series.

[0159] Thresholds are set based on a 97% confidence level:

[0160] Points where the rate of change exceeds a threshold are marked as abnormal change points.

[0161] Local weighted averaging: A local weighted averaging method is used to process parameters in outlier regions, with a weight decay rate of 0.7. For the parameter sequence... The parameters after Gaussian smoothing are calculated as follows:

[0162]

[0163] Among them, weight σ is the smoothing window parameter, ranging from 2 to 5. This method achieves a smooth transition of parameters through weighted averaging, while preserving the main trends in parameter variation and avoiding the loss of details caused by over-smoothing.

[0164] S5: Final Filtering and Result Output

[0165] Based on the optimized elliptical cylinder model and dynamic threshold, the tunnel point cloud data is filtered, and the filtered point cloud data and elliptical cylinder model parameters are output.

[0166] Specifically, points within a threshold range are retained as tunnel structure points (with a confidence level of 98%), while points exceeding the threshold are discarded as unstructured points or noise. The model parameters are then smoothed (with a smoothing coefficient of 0.3) to ensure the overall continuity and consistency of the model. Finally, the filtered point cloud data and elliptical cylinder model parameters are output, achieving millimeter-level processing accuracy to meet the precision requirements for tunnel 3D reconstruction and deformation analysis.

[0167] S5: Based on the optimized elliptical cylinder model and dynamic threshold, perform final filtering on the point cloud and output the filtered point cloud data and elliptical cylinder model parameters.

[0168] Specifically, points within a threshold range are retained as tunnel structure points (with a confidence level of 98%), while points exceeding the threshold are discarded as unstructured points or noise. The model parameters are then smoothed (with a smoothing coefficient of 0.3) to ensure the overall continuity and consistency of the model. Finally, the filtered point cloud data and elliptical cylinder model parameters are output, achieving millimeter-level processing accuracy to meet the precision requirements for tunnel 3D reconstruction and deformation analysis.

[0169] In practical applications, the parameter values ​​can be adjusted according to the specific characteristics of different tunnels. For example, for ordinary highway tunnels, a larger cross-sectional interval (approximately 2m) and a more lenient constraint threshold can be used; for high-speed railway tunnels, a smaller cross-sectional interval (approximately 1m) and a more stringent constraint threshold can be used to meet higher accuracy requirements. For old tunnels with severe deformation, the parameter variation threshold can be appropriately relaxed to improve the model's adaptability to abnormal deformation; for newly built tunnels, geometric continuity constraints can be strengthened to obtain smoother model results.

[0170] The method of this invention has been verified in actual tunnel cases. Taking a highway tunnel as an example, the tunnel is approximately 1.2 kilometers long with a horseshoe-shaped cross-section. Point cloud data was acquired using mobile 3D laser scanning, with a total of approximately 250 million points. After processing using the method of this invention, the accuracy rate of noise point identification reached 96.5%, the retention rate of structural points reached 98.7%, and the processing efficiency was improved by approximately 35% compared with traditional methods. At the same time, the geometric continuity of the model was significantly improved, and the changes in parameters of adjacent sections were smoother and more reasonable. Especially in special areas such as the portal section, the model can accurately adapt to changes in geometric shape, providing a reliable data foundation for subsequent tunnel deformation analysis.

[0171] In non-real-time processing environments with sufficient computational resources, the method of this invention can be extended and replaced. For example, more complex geometric models such as hyperelliptical cylinders or B-spline surfaces can be used to replace the elliptical cylinder model; M-estimated sample consistency (MSAC) or progressive consistent sampling (PROSAC) can be used to replace the RANSAC algorithm; genetic algorithms, particle swarm optimization, or simulated annealing algorithms can be used to replace the iterative least squares method; and Bayesian optimization and other methods can be introduced to automatically adjust various parameters in the algorithm. These alternatives may provide better performance under certain conditions, but they usually increase computational complexity. In practical applications, a balance between accuracy and efficiency must be considered comprehensively.

[0172] In summary, the adaptive elliptic cylinder filtering method for tunnel point clouds based on affine constraints provided by this invention significantly improves the accuracy, stability, and adaptability to complex environments of tunnel point cloud filtering by introducing affine transformation constraints, multi-objective optimization, and adaptive parameter constraints, providing a reliable data foundation for tunnel 3D reconstruction and deformation monitoring.

[0173] The above embodiments are only for illustrating the technical concept and features of the present invention, and are intended to enable those skilled in the art to understand the content of the present invention and implement it accordingly. They should not be construed as limiting the scope of protection of the present invention. All equivalent transformations or modifications made in accordance with the spirit and essence of the present invention should be covered within the scope of protection of the present invention.

Claims

1. An affine constraint based tunnel point cloud adaptive elliptic cylinder filtering method, characterized in that, Includes the following steps: S1: Extract the tunnel centerline based on tunnel point cloud data, and set cross-sections at fixed intervals along the centerline axis; S2: Initialize the parameters of the elliptical cylinder model using principal component analysis on each cross section to form the elliptical cylinder model; S3: Set an initial screening threshold to initially screen outliers in the distance from all points within the cross section to the elliptical cylinder, and use the RANSAC algorithm to dynamically iteratively estimate the optimal screening threshold. S4: Construct a multi-objective optimization function that includes data fitting error terms, parameter smoothing constraint terms, and affine transformation constraint terms, and use an optimization algorithm to optimize the parameters of the elliptical cylinder model until the parameters of the elliptical cylinder model converge. S5: Based on the optimized elliptical cylinder model and the optimal filtering threshold, filter the tunnel point cloud data and output the filtered point cloud data and elliptical cylinder model parameters; Step S4 also includes a step of dynamically adjusting the constraint threshold based on the local curvature and the degree of deformation. Specifically, this includes: dynamically calculating the orientation angle change threshold based on the local curvature using the formula: Base threshold 2°, maximum threshold It is 8°. The local curvature of the centerline, the maximum curvature It is 0.015; The threshold for axis length change is dynamically calculated based on the degree of local deformation, using the following formula: Base threshold 2%, maximum threshold It is 6%. For evaluating the degree of deformation, the maximum degree of deformation is evaluated. It is 0.

08.

2. The method of claim 1, wherein the method is based on affine constraints. In step S1, the tunnel centerline is extracted from the preliminarily denoised tunnel point cloud data using the least squares method or the master curve algorithm, with a fixed interval of 0.5-3 meters.

3. The method of claim 1, wherein the method further comprises: In step S2, the initialization of the elliptical cylinder model parameters using principal component analysis specifically includes: calculating the covariance matrix of the cross-sectional point cloud, performing eigenvalue decomposition on the covariance matrix; determining the directions of the major and minor axes of the elliptical cylinder model using the eigenvectors corresponding to the largest and second largest eigenvalues, and estimating the lengths of the major and minor axes based on the proportions of the eigenvalues.

4. The affine constraint based tunnel point cloud adaptive elliptic cylinder filtering method according to claim 1, characterized in that, In step S3, the initial screening threshold is 20 cm, and the optimal screening threshold is the 95th percentile of the distance between interior points.

5. The affine constraint based tunnel point cloud adaptive elliptic cylinder filtering method according to claim 1, characterized in that, In step S4, the multi-objective optimization function is: ,in Let be the sum of the squares of the distances from the point cloud to the elliptical cylinder. It is the sum of the squares of the differences in orientation angle, major axis length, and minor axis length between adjacent sections. λ1 and λ2 are the sum of squares of the differences between the affine transformation matrix and the desired transformation matrix between adjacent sections, and are weighting coefficients. The Levenberg-Marquardt algorithm is used to minimize the multi-objective optimization function. The initial damping factor is 0.01, and the damping factor adjustment factor is 10. Weighting coefficient λ The value of 1 ranges from 0.5 to 2.0, and the weighting coefficient is... λ The value of 2 ranges from 0.3 to 1.

5.

6. The adaptive elliptic cylindrical filtering method for tunnel point clouds based on affine constraints according to claim 1, characterized in that, Step S4 also includes a step of smoothing the affine transformation parameters, specifically including: calculating the affine transformation matrix and decomposing it into rotation, scaling and translation components; smoothing the transformation parameters using Gaussian weights; and applying cubic spline interpolation to achieve smooth parameter transition.

7. The adaptive elliptic cylindrical filtering method for tunnel point clouds based on affine constraints according to claim 1, characterized in that, Step S4 also includes a step of special processing for regions with significant changes in geometric properties, specifically including: calculating the local curvature change rate and marking curvature abrupt change points; performing local re-optimization on the marked points; establishing a statistical model of parameter changes to identify abnormal changes; and using a local weighted average method to process parameters in abnormal regions.

8. The method of claim 5, wherein the method further comprises: The affine transformation constraint term also includes a dynamic constraint strength control mechanism, through regional adaptive weights. The formula for controlling constraint strength is as follows: Among them, the weighting coefficient The sensitivity coefficient ranges from 0.3 to 1.

5. k The value range is 5-10. Indicates the degree of geometric change of the current cross section; Affine transformation constraint terms Defined as , Indicates from the first j 1 The section to the first j The actual affine transformation matrix of each cross section Denotes the desired affine transformation matrix. This represents the Frobenius norm.

Images