A three-dimensional laser scanning data processing method and system

Abstract

The application relates to the technical field of three-dimensional laser scanning data processing, and discloses a three-dimensional laser scanning data processing method and system, wherein the method comprises the following steps: acquiring multi-period scanning point cloud data and environmental temperature record values and denoising; performing cylindrical surface fitting and structure feature extraction on each station denoised point cloud, identifying a registration degeneration degree of freedom direction; performing coarse registration and ICP precise registration enhanced by degeneration direction regularization based on a structure feature set, generating a registration transformation matrix with a directional covariance; performing pose graph global optimization and fusion to generate a full-tunnel point cloud; performing grid-by-grid difference calculation and propagating registration uncertainty to generate a differential displacement field; accumulating the differential displacement field by period and propagating cumulative uncertainty; performing weighted least squares temperature-deformation regression with the reciprocal of the cumulative uncertainty variance as a weight, deducting a temperature effect deformation component, and extracting a deformation rate.

Description

Technical Field

[0001] This invention relates to the field of three-dimensional laser scanning data processing technology, and more specifically, to a three-dimensional laser scanning data processing method and system. Background Technology

[0002] In the completion acceptance and initial operation deformation monitoring of shield tunnels, three-dimensional laser scanning technology is used to continuously scan the inner wall of the tunnel in multiple phases to obtain three-dimensional point cloud data of the entire tunnel in each monitoring period. After registration, fusion, difference and regression processing, the structural deformation rate is extracted to track the deformation trend of the tunnel structure.

[0003] In existing technologies, the Iterative Closest Point (ICP) algorithm is typically used to perform fine registration of multi-site cloud data, and then the temperature effect deformation component is separated by linear regression to extract the true deformation rate of the structure.

[0004] However, the inner wall of a shield tunnel is assembled from precast concrete segments. The surface of these segments is highly uniform and lacks textural features, exhibiting an approximate smooth cylindrical surface. During ICP registration, the translational degrees of freedom along the cylindrical axis and the rotational degrees of freedom around the cylindrical axis constitute the registration degradation direction. The registration residuals in this direction are systematically higher than those in non-degradation directions such as the radial direction, leading to the superposition of registration pseudo-deformation components dependent on the direction in the differential displacement field. Simultaneously, the temperature inside the tunnel varies periodically due to seasonal climate influences. The thermal expansion effect of the segments results in the superposition of periodic temperature-effect deformation components in the cumulative deformation, with a magnitude comparable to the actual structural deformation. Existing linear regression methods assume uniform observation accuracy in all directions, treating observations with high registration noise in degradation directions and high-precision observations in non-degradation directions equally in temperature-deformation regression. This leads to systematic confusion between registration pseudo-deformation and temperature-effect deformation in degradation directions, incomplete separation of temperature effects, and deformation rate extraction results in degradation directions being doubly contaminated by registration noise and residual temperature effects, failing to reliably identify the true structural deformation trend. Summary of the Invention

[0005] This invention provides a three-dimensional laser scanning data processing method and system, which solves the technical problems in the long-term deformation monitoring of shield tunnels, such as point cloud registration degradation, accumulation of errors in multi-period data fusion, difficulty in separating the coupling between temperature effect and structural deformation, and difficulty in quantitatively assessing the reliability of deformation results.

[0006] This invention discloses a three-dimensional laser scanning data processing method, which includes the following steps: acquiring three-dimensional laser scanning point cloud data and environmental temperature records for each period of multiple consecutive monitoring cycles of the same shield tunnel; performing statistical filtering and denoising processing on the original point cloud of each station; and generating a denoised point cloud dataset for each station.

[0007] Cylindrical surface fitting is performed on the denoised point cloud of each station to obtain the unit vector in the axial direction. The translational degree of freedom along the axial direction and the rotational degree of freedom around the axial direction are identified as the registration degradation degree of freedom direction. The joint line segments of the pipe segments and the center coordinates of the bolt holes are extracted as the structural feature set.

[0008] Based on the set of structural features, coarse registration is performed to obtain the initial transformation matrix. Regularized constraint terms are constructed based on the correspondence of structural features and superimposed on the point-to-surface ICP method equations to perform fine registration. The fine registration transformation matrix and the residual covariance matrix reflecting the registration uncertainty in each direction are output.

[0009] Using the inverse of the residual covariance matrix as the information weight matrix, global optimization of the pose graph is performed to obtain the global optimal pose parameters of each station and fuse them to generate full tunnel point clouds for each phase.

[0010] The differential displacement field is obtained by performing grid-by-grid differential calculation on the fused point clouds of two adjacent periods, and the directional variance of the differential displacement of each grid cell is obtained through uncertainty propagation;

[0011] The cumulative deformation field is obtained by accumulating the differential displacement field period by period, and the cumulative uncertainty variance is obtained by accumulating the directional variance simultaneously.

[0012] Using the reciprocal of the cumulative uncertainty variance of each period as the observation weight, with the temperature time series as the independent variable and the cumulative deformation time series as the dependent variable, weighted least squares regression is performed to deduct the deformation component due to temperature effect, extract the deformation rate, and generate the deformation tracking results of the shield tunnel structure.

[0013] Furthermore, the statistical filtering denoising process includes: calculating the average distance from each point to its nearest neighbor set, identifying points whose average distance exceeds a preset multiple of the global mean and standard deviation as outliers and removing them; after denoising, performing Z-score standardization on the environmental temperature record values, performing mean normalization based on the range on the point cloud coordinate data, and performing joint regression operations involving temperature and deformation in subsequent steps based on the preprocessed values, and finally transforming the deformation rate result back to the original physical dimensions when outputting the deformation rate.

[0014] Further, the extraction of the joint line segments and bolt hole center coordinates as a structural feature set includes: estimating the local normal vector for each point in the denoised point cloud using principal component analysis with a set of nearest neighbors; calculating the angle between the normal vectors of adjacent points; marking points with normal vector angles exceeding a preset angle threshold as normal vector mutation points; aggregating spatially adjacent normal vector mutation points into a connected region as the joint line region; projecting the point cloud of the joint line region onto a local tangent plane and performing piecewise linear fitting to obtain the joint line segment set; and further processing the denoised point cloud... Each point in the model is estimated for local curvature based on its nearest neighbor set through principal component analysis. The ratio of the minimum eigenvalue to the sum of the eigenvalues ​​of the covariance matrix of the nearest neighbor set is used as the curvature measure. Points with local curvature exceeding the mean plus a preset multiple of the standard deviation are marked as high curvature points. Spatially adjacent high curvature points are aggregated into a connected region as the region around the bolt hole. The point cloud of the region around the bolt hole is projected onto the local tangent plane and then circle fitting is performed to obtain the set of coordinates of the bolt hole center. The set of joint line segments and the set of coordinates of the bolt hole center are merged into the set of structural features.

[0015] Furthermore, the step of constructing regularization constraint terms based on structural feature correspondences and superimposing them onto the point-to-surface ICP method equations to perform fine registration includes:

[0016] The source station cloud is transformed to the target station coordinate system based on the current pose estimation, a point-to-point correspondence is established, the point-to-surface residual is calculated based on the local normal vector of the target point, and the normal equation coefficient matrix and gradient vector of the standard point-to-surface ICP are constructed.

[0017] The spatial deviation between each corresponding feature pair in the structural feature correspondence is used as the constraint residual to construct the coefficient matrix and gradient vector of the structural feature constraint. The effect of the coefficient matrix of the structural feature constraint is concentrated on the component corresponding to the direction of the degenerate degree of freedom.

[0018] The ratio of the eigenvalues ​​of the coefficient matrix of the normal equation in the direction of degradation to the eigenvalues ​​of the coefficient matrix of the structural feature constraints in the direction of degradation is calculated, and this ratio is used as the regularization weight parameter.

[0019] The coefficient matrix of the normal equation is obtained by multiplying the coefficient matrix of the regularization weight parameter by the coefficient matrix of the structural feature constraint. The gradient vector is obtained by multiplying the gradient vector of the regularization weight parameter by the gradient vector of the structural feature constraint. The six-degree-of-freedom pose increment vector is obtained by solving the enhanced normal equation.

[0020] Iteratively update the pose estimation until the pose increment magnitude is lower than the convergence threshold.

[0021] Further, the residual covariance matrix is ​​generated as follows: the residual covariance matrix is ​​obtained by multiplying the variance estimate of the registration residual at convergence by the inverse matrix of the enhanced normal equation coefficient matrix, wherein the variance estimate of the registration residual is calculated by dividing the weighted sum of squares of all point-to-point residuals at convergence by the effective degrees of freedom, and the effective degrees of freedom is the difference between the total number of point pairs participating in registration and the dimension of the pose parameters; the degradation direction eigenvalue is extracted by performing eigenvalue decomposition on the normal equation coefficient matrix, arranging the eigenvalues ​​in ascending order, and taking the average of the two smallest eigenvalues ​​as the representative value of the degradation direction eigenvalue.

[0022] Furthermore, the objective function for global optimization of the pose graph is the weighted sum of squares of the relative pose constraint residuals of all adjacent station pairs. For each pair of adjacent stations, the relative pose constraint residual is defined as a six-dimensional deviation vector, expressed in Lie algebra, between the relative pose given by the current optimization variable and the fine registration transformation matrix. The inverse of the residual covariance matrix is ​​used as the information weight matrix of the residual. The objective function is the sum of the squares of the Mahalanobis distances of the residuals of all station pairs, and the Levenberg-Marquardt algorithm is used to solve iteratively. Based on the globally optimal pose parameters, the point clouds of each station are transformed to a unified coordinate system, and voxel downsampling and fusion are performed on the overlapping areas to generate the full tunnel fusion point cloud for each phase.

[0023] Further, the step of obtaining the directional variance of the differential displacement of each grid cell through uncertainty propagation includes: determining the residual covariance matrix corresponding to the scanning station to which the grid cell belongs and the registration link of the adjacent station participating in the differential calculation; calculating the Jacobian matrix of the grid cell displacement relative to the station pose parameters; multiplying the Jacobian matrix with the covariance matrix of the globally optimized station pose parameters and the transpose of the Jacobian matrix; and extracting the diagonal elements of the product matrix as the directional variance vector containing variance components in the axial, radial, and circumferential directions. The covariance matrix of the globally optimized station pose parameters is obtained from the inverse matrix of the normal equation coefficient matrix when the global optimization of the pose graph converges. For cases where the differential displacement involves different station poses in two data periods, the Jacobian matrix contains the partial derivative components of the pose parameters of the corresponding stations in the two periods, and the covariance matrix is ​​correspondingly expanded to a joint covariance matrix containing the covariance of the pose parameters of the corresponding stations in the two periods.

[0024] Furthermore, after obtaining the residual cumulative deformation time series by subtracting the temperature effect deformation component, the process further includes: performing a discrete Fourier transform on the residual cumulative deformation time series to obtain the spectral amplitude distribution, and checking whether there are frequency components with amplitudes exceeding a preset significance threshold at the annual cycle frequency; if the amplitude at the annual cycle frequency does not exceed the preset significance threshold, the deformation rate result is retained; if the amplitude at the annual cycle frequency exceeds the preset significance threshold, a second-order weighted regression correction is performed by introducing a quadratic temperature term, the temperature effect deformation component is recalculated and subtracted, and the deformation rate is updated to the corrected deformation rate slope; wherein, the preset significance threshold is determined based on a preset multiple of the mean amplitude of frequency components outside the annual cycle frequency neighborhood.

[0025] Furthermore, after obtaining the deformation rate, the method further includes: for each directional component of each grid cell, comparing the final cumulative uncertainty variance with a preset variance threshold; if the final cumulative uncertainty variance does not exceed the preset variance threshold, then marking the deformation rate reliability level of the grid cell in that directional component as high level; if the final cumulative uncertainty variance exceeds the preset variance threshold, then marking it as low level. The reliability level label is then attached to the deformation rate results of the corresponding grid cell and directional component, and the shield tunnel structure deformation tracking results with directional reliability labels are output.

[0026] This invention provides a three-dimensional laser scanning data processing system, comprising:

[0027] The data acquisition and denoising module is used to acquire three-dimensional laser scanning point cloud data and environmental temperature records for each period of multiple consecutive monitoring cycles of the same shield tunnel, and to perform statistical filtering and denoising processing on the original point cloud of each station to generate denoised point cloud datasets for each station.

[0028] The structural feature extraction and degradation direction identification module is used to perform cylindrical surface fitting on the denoised point cloud of each station, obtain the unit vector of the axis direction, identify the direction of the registered degradation degree of freedom, and extract the joint line segments of the pipe segments and the center coordinates of the bolt holes as a set of structural features.

[0029] The degradation direction regularization registration module is used to perform coarse registration based on the set of structural features to obtain the initial transformation matrix, construct regularization constraint terms based on the correspondence of structural features and superimpose them into the point-to-surface ICP method equation to perform fine registration, and output the fine registration transformation matrix and residual covariance matrix;

[0030] The pose graph global optimization and fusion module is used to perform global optimization of the pose graph with the inverse of the residual covariance matrix as the information weight matrix, obtain the global optimal pose parameters of each station, and fuse them to generate full tunnel point clouds for each phase.

[0031] The differential calculation and uncertainty propagation module is used to perform grid-by-grid differential calculations on the fused point clouds of two adjacent periods to obtain the differential displacement field, and to obtain the directional variance of the differential displacement of each grid cell through uncertainty propagation;

[0032] The cumulative deformation and cumulative uncertainty calculation module is used to accumulate the differential displacement field periodically to obtain the cumulative deformation field and simultaneously accumulate the directional variance to obtain the cumulative uncertainty variance.

[0033] The weighted regression and deformation rate extraction module is used to perform weighted least squares regression with the inverse of the cumulative uncertainty variance of each period as the observation weight, the temperature time series as the independent variable and the cumulative deformation time series as the dependent variable, to deduct the deformation component of temperature effect, extract the deformation rate, and generate the deformation tracking results of the shield tunnel structure.

[0034] This invention introduces degradation direction regularization constraints based on the structural features of segment joint lines and bolt holes in ICP registration, thereby generating a directional covariance matrix reflecting the differences in actual registration accuracy in each direction. The directional uncertainty information is propagated step by step to the observation weights of temperature-deformation regression. The reciprocal of the cumulative uncertainty variance is used to differentiate the weights of observations in each direction. This solves the technical problem of systematic confusion between degradation direction uncertainty and seasonal temperature effect deformation in ICP registration in low-texture shield tunnel scenarios in traditional equal-weight regression. It achieves the technical effects of suppressing the interference of pseudo-deformation in degradation direction registration on the estimation of temperature effect coefficients, effectively distinguishing the temperature effect deformation component from the actual structural deformation, and attaching reliability level labels to the deformation rate extraction results in each direction. Attached Figure Description

[0035] Figure 1 This is a flowchart of the three-dimensional laser scanning data processing method provided in the embodiments of the present invention;

[0036] Figure 2 This is a schematic diagram of the ambient temperature and Z-score standardized values ​​for each monitoring period provided in the embodiments of the present invention;

[0037] Figure 3 This is a schematic diagram illustrating the number of structural features extracted from station 3 provided in an embodiment of the present invention;

[0038] Figure 4 This is a schematic diagram comparing the diagonal components of the residual covariance matrix of the station pair (3, 4) in each direction according to an embodiment of the present invention.

[0039] Figure 5 This is a schematic diagram of the axial and radial cumulative deformation time series of the mesh element M-0473 provided in an embodiment of the present invention;

[0040] Figure 6This is a schematic diagram illustrating the evolution of the axial and radial cumulative uncertainty variance of the mesh element M-0473 provided in an embodiment of the present invention;

[0041] Figure 7 This is a schematic diagram of the axial and radial observation weight changes of the mesh cell M-0473 provided in an embodiment of the present invention;

[0042] Figure 8 This is a schematic diagram comparing the temperature linearity coefficients of each directional component with the deformation rate provided in the embodiments of the present invention. Detailed Implementation

[0043] In the completion acceptance and initial operation deformation monitoring of shield tunnels, 3D laser scanning technology is used to continuously scan the tunnel inner wall in multiple phases to obtain 3D point cloud data of the entire tunnel in each monitoring period. After registration, fusion, difference, and regression processing, the structural deformation rate is extracted to track the deformation trend of the tunnel structure. However, the inner wall of the shield tunnel is assembled from precast concrete segments, and the surface of the segments is highly uniform and lacks texture features, exhibiting an approximately smooth cylindrical surface characteristic. During ICP fine registration, the point-to-surface distance has extremely low sensitivity to changes in pose parameters in the translational degree of freedom along the cylindrical axis and the rotational degree of freedom around the cylindrical axis, constituting the registration degradation degree of freedom direction. The registration residual system in the registration degradation degree of freedom direction is systematically higher than that in the radial and other non-degradation directions, resulting in the superposition of direction-dependent registration pseudo-deformation components in the differential displacement field.

[0044] Meanwhile, the temperature inside the tunnel varies periodically due to seasonal climate influences. The thermal expansion effect of the tunnel segments results in a cumulative deformation component with a periodic temperature effect deformation that is comparable in magnitude to the actual structural deformation. Existing linear regression methods assume uniform observation accuracy in all directions, treating observations with high registration noise in the degradation direction and high-precision observations in the non-degradation direction equally in temperature-deformation regression. This leads to systematic confusion between registration spurious deformation and temperature effect deformation in the degradation direction. Bias in the temperature-deformation coefficient estimation results in incomplete separation of the temperature effect. The deformation rate extraction results in the degradation direction are doubly contaminated by registration noise and residual temperature effects, making it impossible to reliably identify the true deformation trend of the structure.

[0045] Therefore, the technical problem to be solved by this embodiment is: in the low-texture scenario of shield tunnels, how to use the directional uncertainty information generated during the registration process to differentiate the weighting of temperature-deformation regression, thereby suppressing the contamination of temperature effect separation by pseudo-deformation of degraded direction registration, and accurately extracting the true deformation trend of the structure.

[0046] At least one embodiment of the present invention discloses a three-dimensional laser scanning data processing method, such as... Figure 1 As shown, it includes the following steps:

[0047] According to an embodiment of this invention, a three-dimensional laser scanning data processing method is provided. The execution entity of the three-dimensional laser scanning data processing method is a computer device equipped with point cloud processing software and numerical calculation programs. The computer device is communicatively connected to a three-dimensional laser scanner and a temperature sensor, and is capable of receiving point cloud data acquired by the three-dimensional laser scanner and ambient temperature values ​​recorded by the temperature sensor. The three-dimensional laser scanning data processing method includes the following steps:

[0048] Step 1: Obtain multi-period scan point cloud data and ambient temperature records, perform statistical filtering to denoise, and generate denoised point cloud datasets for each station;

[0049] Obtain the continuous shield tunnel The original point cloud data of the three-dimensional laser scan at each scanning station during each monitoring cycle, and the recorded values ​​of the ambient temperature inside the tunnel during each scanning period. ( ),in Total number of monitoring periods For the first Environmental temperature records for each monitoring period, This serves as an index for the monitoring period. Statistical filtering and denoising are performed on the original point cloud for each station: the average distance from each point to its nearest neighbor set is calculated, and points whose average distance exceeds a preset multiple of the global mean and standard deviation are identified as outliers and removed, generating a denoised point cloud dataset for each station.

[0050] It should be noted that the aforementioned nearest neighbor set refers to the set of surrounding points retrieved centered on the point to be processed, within a preset search radius, or according to a preset number of nearest neighbors. The aforementioned preset multiple range refers to the effective distance interval defined by a multiple of the global average distance plus or minus the standard deviation; for example, 1.0 times the standard deviation can be used as the rejection threshold.

[0051] To eliminate the impact of dimensional differences on subsequent multi-source parameter joint calculations, after denoising in step 1, data preprocessing is performed on the heterogeneous parameters participating in the joint calculations in subsequent steps: for environmental temperature records... Perform Z-score normalization to normalize the ambient temperature recording values. Deformation data with displacement dimensions are at a comparable scale in weighted regression calculations. Mean normalization based on the range is performed on the point cloud coordinate data to unify the coordinate values ​​from different scanning stations to the same range, eliminating the impact of absolute coordinate magnitude differences on the numerical stability of matrix operations. Subsequent steps involving joint regression calculations of temperature and deformation are performed based on the preprocessed values. Finally, the deformation rate is inversely transformed back to the original physical dimensions when outputting the deformation rate.

[0052] Step 2: Perform cylindrical surface fitting and structural feature extraction on the denoised point cloud of each station, identify the direction of the registration degradation degree of freedom, and generate the structural feature set of each station;

[0053] For the Each scanning station The denoised point cloud was fitted with cylindrical least squares to obtain the scanning station position. Unit vector of the axial direction of the local cylindrical surface and radius parameters Based on the unit vector of the axis direction The translational degrees of freedom along the axis and the rotational degrees of freedom about the axis are identified as the scanning station. The two registered degenerate degrees of freedom directions are denoted as the degenerate direction index set. ,in For the translational degenerate degrees of freedom along the axis, For the rotational degenerate degree of freedom about the axis.

[0054] In the denoised point cloud of each station, regions where the normal vector of the detected point cloud changes abruptly are taken as the segment joint line regions. After projecting the point cloud of the segment joint line regions onto a local tangent plane, piecewise linear fitting is performed to obtain the set of joint line segments. ,in , Equal to scanning stations Each joint line segment detected; regions with significantly higher local curvature than the surrounding cylindrical surface are identified as the bolt hole perimeter. The point cloud of the bolt hole perimeter is projected onto a local tangent plane, and circle fitting is performed to obtain the set of bolt hole center coordinates. ,in , Equal to scanning stations The coordinates of the center of each bolt hole detected in the image are used to collect the joint segments. Set of coordinates of the center of the bolt hole Merging, as scanning stations Structural feature set .

[0055] Furthermore, the above-mentioned least-squares fitting of the cylindrical surface is achieved in the following way: using a unit vector along the axial direction. Using the coordinates of a point on the axis as the cylindrical surface parameters, and taking the sum of the squared distances from all points in the point cloud to the cylindrical surface as the objective function, the Levenberg-Marquardt algorithm is used iteratively to find the cylindrical surface parameters that minimize the objective function. The distance from a point to the cylindrical surface is defined as the sum of the perpendicular distance from the point to the axis and the radius. The absolute value of the difference.

[0056] Furthermore, the detection of the above-mentioned point cloud normal vector abrupt change region is achieved in the following way: for each point in the denoised point cloud, the local normal vector is estimated by principal component analysis using the nearest neighbor set of each point; the angle between the normal vectors of adjacent points is calculated, and points whose normal vector angle exceeds a preset angle threshold are marked as normal vector abrupt change points; spatially adjacent normal vector abrupt change points are aggregated into connected regions as candidates for segment joint line regions.

[0057] Furthermore, the detection of the aforementioned regions with significantly higher local curvature than the surrounding cylindrical surface curvature is achieved through the following methods: For each point in the denoised point cloud, the local curvature is estimated based on the nearest neighbor set of each point through principal component analysis, specifically using the ratio of the minimum eigenvalue to the sum of eigenvalues ​​of the covariance matrix of the nearest neighbor set as the curvature measure; the mean and standard deviation of the curvature of the global cylindrical surface region are calculated, and points with local curvature exceeding the mean plus a preset multiple of the standard deviation are marked as high curvature points; spatially adjacent high curvature points are aggregated into connected regions as candidates for the region surrounding the bolt hole.

[0058] Step 3: Perform coarse registration based on the structural feature set, and perform ICP fine registration enhanced by degradation direction regularization to generate a registration transformation matrix with directional covariance;

[0059] For each pair of adjacent stations in each period of scanning data Based on their respective sets of structural features and Feature matching is performed on the spatial distribution patterns of the orientation angles, lengths, and bolt hole centers of the joint segments within the overlapping region to obtain the structural feature correspondences. Calculate the initial rigid body transformation matrix .

[0060] Furthermore, the above feature matching is achieved in the following ways: For seam segments within overlapping areas, a similarity metric is constructed using the difference in line angle and length, and line segment pairs with similarity exceeding a preset threshold are included in the candidate correspondence; for bolt hole centers, the Euclidean distance between centers is used as the similarity metric, and center pairs with distance below a preset threshold are included in the candidate correspondence; combining the above candidate correspondences, the Random Sampling Consensus Algorithm (RANSAC) is used to eliminate outliers, and the interior point correspondences that satisfy global geometric consistency are retained as the final structural feature correspondences. Based on the correspondence between interior points, the optimal rigid body transformation is solved through singular value decomposition (SVD) to obtain the initial rigid body transformation matrix. .

[0061] With the initial rigid body transformation matrix Using the initial values, perform point-to-surface ICP fine registration enhanced with degradation direction regularization, as follows:

[0062] (3-1) Transform the source station cloud to the target station coordinate system based on the current pose estimation, search for the nearest neighbor point in the target station cloud for each source point, establish the point-to-point correspondence, and calculate the point-to-surface residual based on the local normal vector of the target point;

[0063] Furthermore, the point-to-surface residual is defined as the projection of the difference between the source point's transformed coordinates and the target point's coordinates onto the direction of the target point's local normal vector, which is the inner product of the vector connecting the source point and the target point and the target point's normal vector.

[0064] (3-2) Based on the point-to-surface residuals of all point pairs, construct the normal equation coefficient matrix of the standard point-to-surface ICP. and gradient vector ;

[0065] Furthermore, and The construction method is as follows: For each point pair, calculate the point-to-surface residual with respect to the six-degree-of-freedom pose increment vector. Jacobian row vectors The first three dimensions correspond to the partial derivatives of the translation components, and the last three dimensions correspond to the partial derivatives of the rotation components (rotation is represented by Lie algebra); the Jacobian row vectors of all point pairs are stacked row by row to form the Jacobian matrix. The corresponding residuals form a residual vector. ;but , ,in This indicates transpose.

[0066] (3-3) Correspondence of structural features The spatial deviation between corresponding feature pairs is used as the constraint residual. The deviation between joint segment pairs is the Euclidean distance between the endpoints of the transformed segment and the corresponding endpoints. The deviation between bolt hole center pairs is the Euclidean distance between the coordinates of the transformed center and the coordinates of the corresponding center. Based on these constraint residuals, a coefficient matrix for structural feature constraints is constructed. and gradient vector ,in The effect is concentrated in the direction of degenerate degrees of freedom. The corresponding components;

[0067] Furthermore, and The construction method and steps (3-2) and The construction methods are the same, the difference lies in that the residual source is the correspondence of structural features. The spatial deviation of each corresponding feature pair is not the point-to-surface residual; since the joint segments and bolt hole features are mainly distributed in discriminative positions along the degenerate degrees of freedom (axial and rotational directions), the partial derivatives of the spatial deviation of the corresponding feature pairs with respect to the pose components in the degenerate directions are significantly greater than those in the non-degenerate directions. Energy naturally concentrates in the direction of degenerate degrees of freedom. The corresponding components.

[0068] (3-4) Determine the regularization weight parameters : Calculation of the first iteration Degradation direction eigenvalues ​​and The ratio of the eigenvalues ​​in the degradation direction is used as the regularization weight parameter. The adaptive setting value makes the constraint strength of the regularization term on the degenerate direction and the constraint strength of the point-to-surface constraint on the non-degenerate direction of the same order of magnitude.

[0069] Furthermore, the extraction method for the aforementioned degradation direction feature values ​​is as follows: For Perform eigenvalue decomposition, sort the eigenvalues ​​in ascending order, and find the two smallest eigenvalues ​​corresponding to the degenerate degrees of freedom (axial translation and rotation about the axis, respectively). Take the average of these two eigenvalues ​​as the mean. The eigenvalues ​​representing the direction of degradation are denoted by the values; similarly, for Perform eigenvalue decomposition and take the mean of the eigenvalues ​​corresponding to the degradation direction as... The representative value of the degradation direction eigenvalue; regularization weight parameter Take the ratio of the two representative values ​​mentioned above.

[0070] (3-5) Solve the enhanced normal equations to obtain the six-degree-of-freedom pose increment vector. :

[0071]

[0072] in, This is the pose increment vector for six degrees of freedom. For regularization weight parameters, and These are the coefficient matrix and gradient vector of the point-to-surface ICP method equations, respectively. and These are the coefficient matrix and gradient vector, respectively, representing the structural feature constraints.

[0073] (3-6) Increment by pose Update the current pose estimate and determine whether the pose increment magnitude is lower than the convergence threshold; if it has not converged, return to step (3-1) to continue iterating; if it has converged, proceed to step (3-7).

[0074] (3-7) Calculate the residual covariance matrix using the coefficient matrix of the enhanced method equation at convergence. :

[0075]

[0076] in, The variance estimate of the registration residuals is calculated by dividing the weighted sum of squares of all point-to-point to surface residuals at convergence by the effective degrees of freedom; the residual covariance matrix is... The covariance components corresponding to the degradation direction truly reflect the actual registration uncertainty level of each direction after regularization enhancement; output the fine registration transformation matrix. and residual covariance matrix .

[0077] Furthermore, The effective degrees of freedom is the difference between the total number of point pairs participating in the registration and the dimension of the pose parameters (i.e., six).

[0078] Step 4: Perform global pose graph optimization to generate the global optimal pose parameters for each station, and fuse them to generate the full tunnel point cloud;

[0079] Using the pose parameters of all stations in this period as optimization variables, and the fine registration transformation matrix... For relative pose constraints, the residual covariance matrix is ​​used. The inverse of the matrix is ​​used as the information weight matrix of the corresponding constraint. Global optimization of the Levenberg-Marquardt pose graph is then performed to obtain the globally optimal pose parameters for all stations in this period. The Levenberg-Marquardt pose graph global optimization takes as input the relative pose constraints and their information weight matrices between stations, and outputs the optimal absolute pose parameters for each station under global consistency. This is because the residual covariance matrix... The Levenberg-Marquardt pose graph global optimization process automatically applies differentiated weights to each direction, reflecting the differences in actual registration accuracy in each direction. The relatively large covariance component in the degradation direction gives the degradation direction constraint a relatively low weight, thereby suppressing the propagation of residual error in the degradation direction to other stations.

[0080] Furthermore, the objective function for the global optimization of the Levenberg-Marquardt pose graph is the weighted sum of squared relative pose constraint residuals of all adjacent station pairs. Specifically, for each pair of adjacent stations... The relative pose constraint residual is defined as the transformation matrix between the relative pose and the fine registration given by the current optimization variable. The deviation between them (a six-dimensional vector expressed in Lie algebra), in Let be the information weight matrix of the residual, and let the objective function be the sum of the squared Mahalanobis distances of all stations to the residuals, i.e.:

[0081]

[0082] in, Adjacent stations The relative pose constraint residual vector, This is the corresponding information weight matrix. This indicates transpose, and summation iterates through all adjacent station pairs.

[0083] Based on global optimal pose parameters Transform the point clouds of each site to a unified coordinate system, perform voxel downsampling fusion on overlapping areas, and generate the full tunnel fused point cloud for this phase. .

[0084] Step 5: Perform grid-by-grid differential calculations on the fused point clouds of two adjacent periods, propagate the registration uncertainty to the grid-level differential displacement, and generate a differential displacement field with directional uncertainty;

[0085] Fusion of point clouds from two adjacent time periods and After spatial alignment, the tunnel space is divided into regular grid cells. Nearest neighbor search and differential calculation are performed on the point set in each grid cell to obtain the differential displacement components of the grid cell in the axial, radial, and circumferential directions in the local tunnel coordinate system, generating the first... Period and the Three-dimensional differential displacement field between periods ,in For grid cell indexing.

[0086] Furthermore, the specific implementation method of the above-mentioned grid-by-grid differential calculation is as follows: for grid cells... Inner Each point in the period set, at the th In the fused point cloud, the nearest neighbor point is searched, and the difference in coordinates between the two points is calculated as the displacement vector of that point pair. The average values ​​of the displacement vectors of all point pairs within the grid cell are taken in the axial, radial, and circumferential directions, respectively, as the differential displacement components of the grid cell, thus forming the differential displacement vector. .

[0087] For each grid cell Determine the grid cells Scanning station and the residual covariance matrix corresponding to the registration links of adjacent stations participating in the differential calculation. Calculate the Jacobian matrix of the mesh element displacement relative to the station pose parameters. The directional variance estimate of the differential displacement of the mesh elements is obtained through uncertainty propagation:

[0088]

[0089] in, This is the covariance matrix of the globally optimized station pose parameters. This indicates the extraction of diagonal elements from a matrix. Indicates transpose. It is a vector containing variance components in the axial, radial, and circumferential directions. Jacobian matrix. For grid cells Spatial coordinates about grid cells The partial derivatives of the pose parameters of the station location; specifically, let the mesh element be... Spatial coordinates are , grid cell The pose parameters of the station are (Including three-dimensional translation components and three-dimensional rotation components, a total of six dimensions), then The Listed as about No. The partial derivatives of each component, i.e. ,in This is the index of the pose parameter components; for cases where the differential displacement involves poses at different stations in two data periods, the Jacobian matrix... It includes the partial derivative components of the pose parameters of the corresponding stations in the two periods, that is, the partial derivatives of the pose parameters of the two stations are concatenated column by column to form a complete Jacobian matrix. , The corresponding extension is a joint covariance matrix containing the covariances of the pose parameters of the corresponding stations in the two periods. The directional variance estimation inherits information about the differences in directional accuracy during the registration process, and the differential displacement variance in the degenerate direction is systematically higher than that in the non-degenerate direction.

[0090] Furthermore, the covariance matrix of the above-mentioned globally optimized station pose parameters The Hessian matrix is ​​obtained from the inverse of the normal equation coefficient matrix at the convergence of the global optimization of the Levenberg-Marquardt pose graph. Specifically, the Hessian matrix is ​​obtained by taking the second-order partial derivative of the global optimization objective function with respect to the joint vector of pose parameters at all stations, and the inverse of the Hessian matrix is ​​taken as... The diagonal elements of the Hessian matrix correspond to the variance of each pose parameter component at each station, while the off-diagonal elements reflect the covariance between different stations or different pose parameter components.

[0091] Step 6: Accumulate the differential displacement field period by period, propagate the accumulated uncertainty synchronously, and generate the accumulated deformation time series and temperature time series;

[0092] Using the first phase of fused point cloud as a baseline, the differential displacement field is accumulated phase by phase to obtain the cumulative deformation field of each phase relative to the baseline phase:

[0093]

[0094] in, For the first Periodic grid cells The cumulative deformation vector, For the first Period and the Inter-period grid cells The difference displacement vector, For the summation index, For the current monitoring period index, This represents the total number of monitoring periods.

[0095] Since the registration errors of the differential displacements in each period are independent, the variance of the directional cumulative uncertainty of the cumulative deformation in each period is calculated through variance linear propagation:

[0096]

[0097] in, For the first Periodic grid cells The cumulative uncertainty variance vector in each direction, For the first The directional variance vector of the periodic difference displacement. For summation index. Extract each grid cell. In all directions ( Cumulative deformation time series on ) (in ) and the corresponding cumulative uncertainty variance time series and the temperature time series composed of temperature records from each period. Among them, temperature time series The Z-score-normalized values ​​from step 1 are used.

[0098] Step 7: Perform weighted least squares temperature-deformation regression with the reciprocal of the cumulative uncertainty variance as the weight, subtract the temperature effect deformation component, extract the deformation rate, and generate the shield tunnel structure deformation tracking results.

[0099] For each grid cell Each directional component , in terms of directional components Cumulative uncertainty variance The reciprocal of the first digit is used as the second digit. Weight of periodic observations temperature time series Independent variable, directional component Cumulative deformation time series Using the dependent variable, perform weighted least squares regression fitting:

[0100]

[0101] in, Directional components Temperature-deformation linearity coefficient, The slope of the deformation rate, For the first Monitoring period For constant offset terms, For the first Expected directional components The observation weights on the surface Total number of monitoring periods For summation index. In the above regression... The dimensionless temperature value after Z-score normalization. The normalized dimensionless deformation value is... The time values ​​are dimensionless after normalization. Since the three values ​​have the same dimensions, they can be directly used in weighted regression calculations.

[0102] Furthermore, the solution method for the above weighted least squares regression is as follows: for grid cells... Directional components Construct the design matrix Design matrix The Behavior ;

[0103] Construct the weight diagonal matrix Weight diagonal matrix The The diagonal elements are Construct observation vectors Observation vector The The elements are Then the regression coefficient vector From the weighted equation The solution is obtained, where This indicates transpose.

[0104] Due to the cumulative uncertainty variance in each period along the degradation direction The systematic aspect is greater than the non-degenerate aspect; the weight corresponding to the observation in the degenerate aspect is greater. The weight of observations in the degenerate direction is lower than that in the non-degenerate direction. Within the weighted least squares framework, observations in the degenerate direction containing higher registration spurious deformation components receive lower weights due to their high uncertainty. The temperature-deformation linearity coefficient of degenerate direction observations is affected. The estimated contribution is suppressed; high-precision observations in the non-degenerate direction are given higher weight and fully participate in the accurate estimation of the temperature effect coefficient, thus making the temperature-deformation linear coefficient more accurate. It is not affected by systemic contamination from registration pseudo-deformation in the degradation direction.

[0105] Based on temperature-deformation linear coefficient and temperature records for each period Calculate the temperature effect deformation components in each direction and period for each grid cell. The residual cumulative deformation time series is generated by subtracting it from the cumulative deformation in each direction and period.

[0106]

[0107] in, After deducting the temperature effect, the first Periodic grid cells In directional components The residual cumulative deformation value on.

[0108] Based on the slope of deformation rate The deformation rate and confidence interval of each grid cell in each direction are output from the residual cumulative deformation time series. After the deformation rate is inversely transformed to the original physical dimensions, the deformation tracking results of the shield tunnel structure are generated.

[0109] Furthermore, the aforementioned confidence interval is determined by the slope of the deformation rate in the weighted least squares regression. The standard error of the estimate is determined by the coefficient matrix of the weighted normal equation. The slope of the deformation rate in the inverse matrix The square root of the diagonal element is obtained, where Indicates transpose, confidence interval is expressed as slope of deformation rate. The estimated value is given in the form of a factor of several standard errors.

[0110] In this embodiment of the application, in order to verify whether the temperature effect component has been sufficiently removed, after obtaining the residual cumulative deformation time series in step 7, the following steps are also included:

[0111] Residual cumulative deformation time series Perform a discrete Fourier transform to obtain the spectral amplitude distribution and examine the frequency distribution within the annual cycle. Does the frequency component still exist at this location, with an amplitude exceeding the preset significance threshold? If the annual cycle frequency... If the amplitude at the location does not exceed the preset significance threshold, it is determined that the temperature effect has been sufficiently removed, and the deformation rate result obtained in step 7 is retained.

[0112] Furthermore, the aforementioned preset significance threshold is determined based on the spectral background noise level; specifically, the annual cycle frequency is... The mean amplitude of frequency components outside the neighborhood is used as the background noise benchmark, and a preset multiple of the background noise benchmark is used as the preset significance threshold. The annual cycle frequency... When the amplitude exceeds a preset significance threshold, a significant residual periodic component is determined to exist.

[0113] If the annual cycle frequency If the magnitude of the temperature effect exceeds a preset significance threshold, the linear temperature model is deemed to have failed to fully capture the nonlinear components of the temperature effect, and a second-order weighted regression correction is performed by introducing a quadratic temperature term.

[0114]

[0115] in, and These are the coefficients of the first and second terms of temperature, respectively. This is the corrected slope of the deformation rate. This is the corrected constant offset term. Total number of monitoring periods For summation indexing, the temperature effect deformation component is recalculated and subtracted to generate a refined residual cumulative deformation time series. The deformation rate is then updated to the corrected deformation rate slope. and the corrected deformation rate slope The output is obtained after inverse transformation to the original physical dimensions.

[0116] In this embodiment of the application, in order to quantify the reliability of deformation tracking results in each direction, after obtaining the deformation rate, the following steps are also included:

[0117] For each grid cell Each directional component Calculate the cumulative uncertainty variance at the end of the period. The cumulative uncertainty variance at the end of the period Compared with the preset variance threshold Compare. If Then label the grid cells. In directional components The deformation rate reliability level is high; if If the reliability level is low, it is labeled as low-level. The reliability level label is appended to the deformation rate results of the corresponding mesh element and directional component, outputting the shield tunnel structure deformation tracking results with directional reliability labels. A preset variance threshold is used. The required accuracy for monitoring tunnel structural deformation is determined based on the specific requirements. For example, a certain percentage of the allowable deformation in the tunnel design can be converted into a variance dimension and used as a preset variance threshold. .

[0118] This implementation addresses the technical problem of systematic confusion between ICP registration degradation direction uncertainty and seasonal temperature effect deformation in traditional equal-weight regression under low-texture shield tunnel scenarios. A targeted solution is achieved through the synergistic effect of the following technical means:

[0119] First, in step 2, the directions of registration degradation degrees of freedom caused by the cylindrical surface geometry were identified, and structural features of non-smooth surface regions such as segment joint lines and bolt holes were extracted. In step 3, the correspondences of structural features in the structural feature set were introduced as constraint residuals into the regularization penalty term of the ICP method equations to compensate for near-zero eigenvalues ​​in the degradation direction. This effectively constrains the pose parameters in the degradation direction, improving registration convergence accuracy, and simultaneously increasing the residual covariance matrix. It accurately reflects the actual registration accuracy differences in each direction after regularization enhancement.

[0120] Secondly, in step 5, the Jacobian matrix is ​​used. The station-level residual covariance matrix is ​​propagated to the directional variance estimation of the differential displacement at the grid cell level; in step 6, the directional cumulative uncertainty variance of the cumulative deformation in each period is obtained through linear variance propagation. Therefore, the observation reliability in each direction of each grid cell is quantitatively characterized, and the higher systematic uncertainty variance in the degradation direction and the lower uncertainty variance in the non-degradation direction are fully preserved and propagated step by step.

[0121] In step 7, weighted least squares temperature-deformation regression is performed using the reciprocal of the cumulative uncertainty variance as the observation weight. Observations with higher registration spurious deformation components in the degradation direction receive lower weights due to their high uncertainty. The linear coefficient of the temperature-deformation regression is affected by the degradation direction observations. The contribution of the estimation is suppressed; high-precision observations in non-degenerate directions are given high weights and fully participate in the accurate estimation of the temperature effect coefficient. This differentiated weighting avoids the systematic contamination of the temperature effect coefficient estimation by pseudo-deformation of registration in the degradation direction, and enables accurate separation of the temperature effect deformation components.

[0122] Furthermore, spectral analysis and second-order correction steps further examine and eliminate residual components of nonlinear temperature effects that the linear temperature model failed to capture, and directional reliability labeling provides a confidence level reference for the deformation rate results in each direction for subsequent engineering decisions.

[0123] In summary, this implementation method overcomes the systematic confusion between pseudo-deformation due to degradation direction registration and temperature effect deformation in traditional equal-weight regression by introducing degradation direction registration regularization constraints based on structural features and generating a directional covariance matrix in the registration stage, and by propagating directional uncertainty information step by step to the observation weights of temperature-deformation regression. In the separation of temperature effects and extraction of deformation rate, it can accurately distinguish between temperature effect deformation and real structural deformation, and avoid misjudgment of deformation trend caused by degradation direction registration noise.

[0124] The shield tunnel of the AA section of a city's subway system was completed and accepted in 2023, entering the initial deformation monitoring phase during operation. The tunnel has an inner diameter of 6.2m and is assembled from precast concrete segments, each 1.2m wide, with each ring consisting of 6 segments. The tunnel is approximately 480m long, with 8 three-dimensional laser scanning stations (stations 1 to 8) deployed along the axial direction, with an overlap of approximately 15m between adjacent stations. The monitoring system consists of a three-dimensional laser scanner and tunnel wall temperature sensors. The temperature sensors record the ambient temperature inside the tunnel during each scan. Starting in January 2023, a full tunnel scan was performed approximately every two months, completing 8 monitoring cycles (K=8), covering a full year, to track the deformation trend of the tunnel structure during the initial operation phase.

[0125] like Figure 2-8 As shown, the method in this embodiment operates as follows:

[0126] Over eight monitoring cycles, a 3D laser scanner sequentially scanned stations 1 through 8, with each station's original point cloud containing approximately 2.8 million points. Simultaneously, temperature sensors recorded the ambient temperature inside the tunnel during each scan period. Statistical filtering denoising was performed on the original point clouds for each station: the average distance from each point to its 30 nearest neighbors was calculated, and points whose average distance exceeded the global mean plus or minus 1.0 standard deviation were identified as outliers and removed. On average, approximately 21,000 points were removed per station (approximately 0.75% of the original points), generating denoised point cloud datasets for each station. After denoising, Z-score standardization was performed on the ambient temperature records, and mean normalization based on the range was performed on the point cloud coordinate data.

[0127] Table 1. Recorded ambient temperature values ​​and Z-score standardized results for each monitoring period:

[0128]

[0129] The Z-score normalization formula is the average temperature. ℃, standard deviation ℃. Taking the first issue as an example:

[0130]

[0131] Taking the first phase of denoised point cloud at station 3 as an example, cylindrical least squares fitting is performed. After iterative convergence using the Levenberg-Marquardt algorithm, the unit vector along the local cylindrical axis of station 3 is obtained. (Approximately parallel to the tunnel axis), radius parameter m. Based on The translational degradation degree of freedom along the axis (axial translation) and the rotational degradation degree of freedom around the axis (rotation around the axis) are identified as the degradation direction index set.

[0132] In the normal vector mutation detection, points with a normal vector angle exceeding 25° were marked as mutation points. After aggregation, 14 segment joint lines were identified in the point cloud of station 3. In the high curvature point detection, points with local curvature exceeding the mean plus three times the standard deviation were aggregated to identify the center coordinates of 38 bolt holes. Finally, a set of structural features for station 3 was generated, including a set of 14 joint lines and a set of 38 bolt holes.

[0133] Table 2: Results of structural feature extraction at station 3:

[0134]

[0135] Taking the registration of adjacent station pairs (station 3, station 4) in the first phase of scanning as an example, based on the joint segments and bolt holes in the overlapping areas of the structural feature sets of station 3 and station 4, after removing outliers using RANSAC, 11 pairs of joint segment correspondences and 24 pairs of bolt hole center correspondences are retained. The initial rigid body transformation matrix is ​​obtained by solving using SVD.

[0136] Using this initial transformation matrix as the initial value, point-to-surface ICP fine registration is performed with degradation direction regularization enhancement. In the first iteration, eigenvalue decomposition is performed on the point-to-surface information matrix, and the mean of the two smallest eigenvalues ​​(corresponding to the axial translation and rotational degradation directions) is 0.0183; eigenvalue decomposition is performed on the structural feature information matrix, and the mean of the eigenvalues ​​corresponding to the degradation directions is 1.847. Based on this, the regularization weight parameters are adaptively set.

[0137]

[0138] After 18 iterations, the pose increment magnitude decreased below the convergence threshold, indicating registration convergence. Upon convergence, the estimated variance of the registration residuals was obtained by dividing the weighted sum of squared residuals from all point-to-point to surface relationships by the effective degrees of freedom (total number of point pairs minus 6). m² is then used to calculate the residual covariance matrix, where the diagonal component of the covariance corresponding to the degradation direction is significantly larger than the radial component.

[0139] Table 3. Diagonal components of the residual covariance matrix for station pair (3,4):

[0140]

[0141] Degeneration direction (axial translation) covariance diagonal components m² is approximately a radial translation The fact that the value is 61 times that of m² truly reflects the actual registration uncertainty of the degenerate direction, which, even after regularization enhancement, is still higher than that of the non-degenerate direction.

[0142] Taking the first phase of scanning as an example, the pose parameters of 8 stations are used as optimization variables, the fine registration transformation matrices of 7 pairs of adjacent stations are used as relative pose constraints, and the inverses of the corresponding residual covariance matrices are used as information weight matrices. Global optimization of the Levenberg-Marquardt pose graph is performed, with the objective function being:

[0143]

[0144] After global optimization and convergence, the globally optimal pose parameters for 8 stations were obtained. Due to the large covariance components in the degradation direction (e.g., the axial translation direction is approximately...), the optimal pose parameters were not obtained. The degradation direction component in the corresponding information weight matrix is ​​relatively small, and the degradation direction constraint automatically obtains a low weight in the global optimization, suppressing the cross-site propagation of the degradation direction residual error. Based on the globally optimal pose parameters, the point clouds of each station are transformed to a unified coordinate system, and voxel downsampling (voxel size 0.02m) fusion is performed on the overlapping area to generate the first phase of the full tunnel fused point cloud, with a total of approximately 13.4 million points. The above process is repeated for 8 monitoring cycles to generate 8 phases of full tunnel fused point clouds from phase 1 to phase 8.

[0145] Taking the differential calculation between the merged point clouds of Phase 2 and Phase 1 as an example, the tunnel space is divided into regular grid cells of 0.5m × 0.5m × 0.5m, generating approximately 14,600 effective grid cells. For grid cell M-0473 (located in the coverage area of ​​station 3, at approximately 143.5m along the tunnel axis), the nearest neighbor of each point in the Phase 2 point set within the grid is searched in the merged point cloud of Phase 1, forming 187 effective point pairs. The displacement averages are taken in the axial, radial, and circumferential directions to obtain the differential displacement vector:

[0146]

[0147] The three components correspond to the axial, radial, and circumferential directions, respectively. Using the Jacobian matrix (which is a 3×12 matrix obtained by taking the partial derivatives of the spatial coordinates of the grid cells with respect to the pose parameters of their respective station 3 and adjacent station 4) and the joint covariance matrix of the globally optimized pose parameters, the directional differential displacement variance estimate of M-0473 in the second phase is obtained through propagation:

[0148]

[0149] Axial variance mm² is approximately the radial variance It is 63 times the size of mm², inheriting the characteristic from step 3 that the covariance component of the degradation direction is significantly higher than that of the non-degradation direction.

[0150] For grid cell M-0473, the differential displacement fields of 8 monitoring periods are accumulated period by period to generate the cumulative deformation time series of each period relative to the base period (period 1). Simultaneously, the cumulative uncertainty variance time series is calculated through variance linear propagation, and the reciprocal of the cumulative variance is taken as the observation weight of each period.

[0151] Table 4. Cumulative deformation and cumulative uncertainty variance (axial and radial components) of mesh element M-0473 for each period:

[0152]

[0153] The comparison between axial and radial directions in Table 4 intuitively reflects the difference in directional uncertainty: the radial observation weight (34364.0) in the second period is about 63 times that of the axial observation weight (546.4). This ratio is completely consistent with the variance propagation result in step 5, reflecting the systematic difference in observation reliability between the degenerate direction (axial) and the non-degenerate direction (radial).

[0154] For the axial component of grid cell M-0473, using the reciprocal of the cumulative variance of the axial direction in each period in Table 4 as the weight, and standardized temperature as the independent variable and cumulative axial deformation as the dependent variable (periods 2 to 8, a total of 7 observations), a design matrix (7×3, each row containing standardized temperature, standardized time, and constant terms), a weighted diagonal matrix (7×7, with diagonal elements representing the axial observation weights) and an observation vector (7-dimensional) were constructed. The weighted method equation was solved, and the linear coefficient of axial temperature-deformation was obtained as 0.412 mm / (standardized temperature unit), and the slope of the axial deformation rate was -0.071 mm / (standardized time unit).

[0155] For the radial component, a weighted least squares regression was performed similarly with the radial weights in Table 4, yielding a radial temperature-deformation linear coefficient of 0.198 mm / (standardized temperature units) and a radial deformation rate slope of -0.134 mm / (standardized time units).

[0156] Based on the axial temperature-deformation linear coefficient (0.412) and the standardized temperature for each period, the axial temperature effect deformation component is calculated period by period and subtracted from the axial cumulative deformation to generate a time series of axial residual cumulative deformation. A discrete Fourier transform is performed on this residual series. The amplitude at the annual period frequency is 0.031 mm, the mean background noise is 0.028 mm, and the ratio is 1.11, which does not exceed the preset significance threshold (2.0 times the mean background noise, i.e., 0.056 mm). Therefore, it is determined that the temperature effect has been sufficiently removed, and the current deformation rate result is retained.

[0157] After inversely transforming the deformation rate slope to the original physical dimensions, the axial deformation rate of mesh element M-0473 is -0.42 mm / year, and the radial deformation rate is -0.81 mm / year. A preset variance threshold is used. Based on mm², the axial final variance mm² exceeds the threshold, indicating a low reliability level for axial deformation rate; radial terminal variance... When the diameter of the radial deformation rate is below the threshold, the reliability of the radial deformation rate is rated as high.

[0158] Table 5 Summary of final deformation tracking results for mesh element M-0473:

[0159]

[0160] Throughout the implementation process, the data underwent a complete step-by-step transformation from its initial state: In Step 1, the original point cloud from the 8 scans was statistically filtered and denoised, and the ambient temperature was Z-score standardized (e.g., the temperature in Phase 1 was standardized from 12.3℃ to -1.234), forming a unified dimensional basis for subsequent calculations; The set of structural features extracted in Step 2 (14 joint segments and 38 bolt holes at Station 3) was introduced into the ICP method equation as a regularization constraint in Step 3, making the diagonal component of the residual covariance in the degradation direction (axial direction) ( m²) is approximately radial ( The directional accuracy difference information is completely preserved, being 61 times that of m²; the pose parameter covariance matrix generated by global optimization in step 4 is propagated to mesh element M-0473 in step 5 through the Jacobian matrix, and the axial differential displacement variance (m²) is 61 times that of m². mm²) is approximately radial ( The axial cumulative variance was 63 times that of mm²; after step 6, the cumulative variance continued to increase until the end of the period. mm², corresponding to an axial observation weight (546.4) and a radial observation weight (34364.0) that differ by about 63 times; in step 7, this weight difference directly affects the weighted least squares regression, which systematically suppresses the observation contribution of the degradation direction (axial) containing high-registration pseudo-deformation, and the temperature effect coefficient estimation is not polluted by axial registration noise. Finally, after spectral verification and reliability labeling, the deformation rate result with directional reliability level under physical dimensions is output, realizing the full-link data flow from the original point cloud to the reliable deformation trend.

[0161] The embodiments of the present invention have been described above. However, the embodiments are not limited to the specific implementation methods described above. The specific implementation methods described above are merely illustrative and not restrictive. Those skilled in the art can make more equivalent embodiments under the guidance of the present embodiments, and all of them are within the protection scope of the present embodiments.

Claims

1. A method for processing three-dimensional laser scanning data, characterized in that, Includes the following steps: The three-dimensional laser scanning point cloud data and environmental temperature records for each period of the same shield tunnel are obtained for multiple consecutive monitoring cycles. Statistical filtering and denoising processing are performed on the original point cloud of each station to generate a denoised point cloud dataset for each station. Cylindrical surface fitting is performed on the denoised point cloud of each station to obtain the unit vector in the axial direction. The translational degree of freedom along the axial direction and the rotational degree of freedom around the axial direction are identified as the registration degradation degree of freedom direction. The coordinates of the pipe segment joint line segment and the center of the bolt hole are extracted as the structural feature set. Based on the set of structural features, coarse registration is performed to obtain the initial transformation matrix. Regularized constraint terms are constructed based on the correspondence of structural features and superimposed on the point-to-surface ICP method equations to perform fine registration. The fine registration transformation matrix and the residual covariance matrix reflecting the registration uncertainty in each direction are output. The ratio of the eigenvalues ​​of the coefficient matrix of the normal equation in the direction of degradation to the eigenvalues ​​of the coefficient matrix of the structural feature constraints in the direction of degradation is calculated, and this ratio is used as the regularization weight parameter. The residual covariance matrix is ​​generated by multiplying the variance estimate of the registration residuals at convergence by the inverse matrix of the coefficient matrix of the enhanced normal equation. Using the inverse of the residual covariance matrix as the information weight matrix, global optimization of the pose graph is performed to obtain the global optimal pose parameters of each station and fuse them to generate full tunnel point clouds for each phase. The differential displacement field is obtained by performing grid-by-grid differential calculation on the fused point clouds of two adjacent phases. The directional variance of the differential displacement of each grid cell is obtained through uncertainty propagation, including: determining the scanning station to which the grid cell belongs and the residual covariance matrix corresponding to the registration link of the adjacent station participating in the differential calculation; calculating the Jacobian matrix of the grid cell displacement relative to the station pose parameters; multiplying the Jacobian matrix with the covariance matrix of the globally optimized station pose parameters and the transpose of the Jacobian matrix; and extracting the diagonal elements of the product matrix as the directional variance vector containing variance components in the axial, radial and circumferential directions. The cumulative deformation field is obtained by accumulating the differential displacement field period by period, and the cumulative uncertainty variance is obtained by accumulating the directional variance simultaneously. The registration errors of the differential displacements in each period are independent, and the variance of the directional cumulative uncertainty of the cumulative deformation in each period is calculated by variance linear propagation. Using the reciprocal of the cumulative uncertainty variance of each period as the observation weight, with the temperature time series as the independent variable and the cumulative deformation time series as the dependent variable, weighted least squares regression is performed to deduct the deformation component of temperature effect, extract the deformation rate, and generate the deformation tracking results of the shield tunnel structure. Based on the temperature-deformation linearity coefficient and the temperature records for each period, the temperature effect deformation components for each grid cell in each direction and period are calculated and deducted from the cumulative deformation in each direction and period.

2. The three-dimensional laser scanning data processing method according to claim 1, characterized in that, The statistical filtering denoising process includes: calculating the average distance from each point to its nearest neighbor set, identifying points whose average distance exceeds the global mean and standard deviation by a preset multiple as outliers and removing them; after denoising, performing Z-score standardization on the ambient temperature record values, performing mean normalization based on the range on the point cloud coordinate data, and performing joint regression operations involving temperature and deformation in subsequent steps based on the preprocessed values, and finally transforming the deformation rate result back to the original physical dimensions when outputting the deformation rate.

3. The three-dimensional laser scanning data processing method according to claim 1, characterized in that, The extraction of the joint line segments and bolt hole center coordinates as a structural feature set includes: estimating the local normal vector for each point in the denoised point cloud using principal component analysis with a nearest neighbor set; calculating the angle between the normal vectors of adjacent points; marking points with normal vector angles exceeding a preset angle threshold as normal vector mutation points; aggregating spatially adjacent normal vector mutation points into a connected region as the joint line region; projecting the joint line region point cloud onto a local tangent plane and performing piecewise linear fitting to obtain the joint line segment set; and further processing each point in the denoised point cloud... The local curvature of each point is estimated using principal component analysis based on the nearest neighbor point set. The ratio of the smallest eigenvalue to the sum of eigenvalues ​​of the covariance matrix of the nearest neighbor point set is used as the curvature measure. Points with local curvature exceeding the mean plus a preset multiple of the standard deviation are marked as high curvature points. Spatially adjacent high curvature points are aggregated into a connected region as the region around the bolt hole. The point cloud of the region around the bolt hole is projected onto the local tangent plane and circle fitting is performed to obtain the set of coordinates of the bolt hole center. The set of joint line segments and the set of coordinates of the bolt hole center are merged into the set of structural features.

4. The three-dimensional laser scanning data processing method according to claim 1, characterized in that, The step of constructing regularization constraint terms based on structural feature correspondences and superimposing them into the point-to-surface ICP method equations to perform fine registration includes: The source station cloud is transformed to the target station coordinate system based on the current pose estimation, a point-to-point correspondence is established, the point-to-surface residual is calculated based on the local normal vector of the target point, and the normal equation coefficient matrix and gradient vector of the standard point-to-surface ICP are constructed. The spatial deviation between each corresponding feature pair in the structural feature correspondence is used as the constraint residual to construct the coefficient matrix and gradient vector of the structural feature constraint. The effect of the coefficient matrix of the structural feature constraint is concentrated on the component corresponding to the direction of the degenerate degree of freedom. The coefficient matrix of the normal equation is obtained by multiplying the coefficient matrix of the regularization weight parameter by the coefficient matrix of the structural feature constraint. The gradient vector is obtained by multiplying the gradient vector of the regularization weight parameter by the gradient vector of the structural feature constraint. The six-degree-of-freedom pose increment vector is obtained by solving the enhanced normal equation. Iteratively update the pose estimation until the pose increment magnitude is lower than the convergence threshold.

5. The three-dimensional laser scanning data processing method according to claim 1, characterized in that, The variance estimate of the registration residual is calculated by dividing the weighted sum of squares of all point-to-surface residuals at convergence by the effective degrees of freedom. The effective degrees of freedom is the difference between the total number of point pairs participating in the registration and the dimension of the pose parameters. The degradation direction eigenvalue is extracted by performing eigenvalue decomposition on the normal equation coefficient matrix, arranging the eigenvalues ​​in ascending order, and taking the average of the two eigenvalues ​​with the smallest values ​​as the representative value of the degradation direction eigenvalue.

6. The three-dimensional laser scanning data processing method according to claim 1, characterized in that, The objective function for global optimization of the pose map is the weighted sum of squared relative pose constraint residuals of all adjacent station pairs. For each pair of adjacent stations, the relative pose constraint residual is defined as a six-dimensional deviation vector, expressed in Lie algebra, between the relative pose given by the current optimization variable and the fine registration transformation matrix. The inverse of the residual covariance matrix is ​​used as the information weight matrix of the residual. The objective function is the sum of squared Mahalanobis distances of the residuals of all station pairs. The Levenberg-Marquardt algorithm is used to solve iteratively. Based on the globally optimal pose parameters, the point clouds of each station are transformed to a unified coordinate system. Voxel downsampling and fusion are performed on the overlapping areas to generate full tunnel fusion point clouds for each phase.

7. The three-dimensional laser scanning data processing method according to claim 1, characterized in that, The covariance matrix of the station pose parameters after global optimization is obtained by the inverse matrix of the normal equation coefficient matrix when the global optimization of the pose graph converges. For the case where the differential displacement involves different station poses in two periods of data, the Jacobian matrix contains the partial derivative components of the pose parameters of the corresponding stations in the two periods, and the covariance matrix is ​​correspondingly expanded to a joint covariance matrix containing the covariance of the pose parameters of the corresponding stations in the two periods.

8. The three-dimensional laser scanning data processing method according to claim 1, characterized in that, After obtaining the residual cumulative deformation time series by subtracting the temperature effect deformation component, the process further includes: performing a discrete Fourier transform on the residual cumulative deformation time series to obtain the spectral amplitude distribution, and checking whether there are frequency components with amplitudes exceeding a preset significance threshold at the annual cycle frequency; if the amplitude at the annual cycle frequency does not exceed the preset significance threshold, the deformation rate result is retained; if the amplitude at the annual cycle frequency exceeds the preset significance threshold, a second-order weighted regression correction is performed by introducing a quadratic temperature term, the temperature effect deformation component is recalculated and subtracted, and the deformation rate is updated to the corrected deformation rate slope; wherein, the preset significance threshold is determined based on a preset multiple of the mean amplitude of frequency components outside the annual cycle frequency neighborhood.

9. The three-dimensional laser scanning data processing method according to claim 1, characterized in that, After obtaining the deformation rate, the process also includes: for each directional component of each grid cell, comparing the final cumulative uncertainty variance with a preset variance threshold; if the final cumulative uncertainty variance does not exceed the preset variance threshold, then marking the deformation rate reliability level of the grid cell in that directional component as high level; if the final cumulative uncertainty variance exceeds the preset variance threshold, then marking it as low level. The reliability level label is then attached to the deformation rate results of the corresponding grid cell and directional component, and the shield tunnel structure deformation tracking results with directional reliability labels are output.

10. A three-dimensional laser scanning data processing system, used to execute the three-dimensional laser scanning data processing method according to any one of claims 1 to 9, characterized in that, include: The data acquisition and denoising module is used to acquire three-dimensional laser scanning point cloud data and environmental temperature records for each period of multiple consecutive monitoring cycles of the same shield tunnel, and to perform statistical filtering and denoising processing on the original point cloud of each station to generate denoised point cloud datasets for each station. The structural feature extraction and degradation direction identification module is used to perform cylindrical surface fitting on the denoised point cloud of each station, obtain the unit vector of the axis direction, identify the direction of the registered degradation degree of freedom, and extract the joint line segments of the pipe segments and the center coordinates of the bolt holes as a set of structural features. The degradation direction regularization registration module is used to perform coarse registration based on the set of structural features to obtain the initial transformation matrix, construct regularization constraint terms based on the correspondence of structural features and superimpose them into the point-to-surface ICP method equation to perform fine registration, and output the fine registration transformation matrix and residual covariance matrix; The ratio of the eigenvalues ​​of the coefficient matrix of the normal equation in the direction of degradation to the eigenvalues ​​of the coefficient matrix of the structural feature constraints in the direction of degradation is calculated, and this ratio is used as the regularization weight parameter. The residual covariance matrix is ​​generated by multiplying the variance estimate of the registration residuals at convergence by the inverse matrix of the coefficient matrix of the enhanced normal equation. The pose graph global optimization and fusion module is used to perform global optimization of the pose graph with the inverse of the residual covariance matrix as the information weight matrix, obtain the global optimal pose parameters of each station, and fuse them to generate full tunnel point clouds for each phase. The differential calculation and uncertainty propagation module is used to perform grid-by-grid differential calculations on the fused point clouds of two adjacent periods to obtain the differential displacement field. It obtains the directional variance of the differential displacement of each grid cell through uncertainty propagation, including: determining the scanning station to which the grid cell belongs and the residual covariance matrix corresponding to the registration link of the adjacent station participating in the differential calculation; calculating the Jacobian matrix of the grid cell displacement relative to the station pose parameters; multiplying the Jacobian matrix with the covariance matrix of the globally optimized station pose parameters and the transpose of the Jacobian matrix; and extracting the diagonal elements of the product matrix as the directional variance vector containing variance components in the axial, radial and circumferential directions. The cumulative deformation and cumulative uncertainty calculation module is used to accumulate the differential displacement field periodically to obtain the cumulative deformation field and simultaneously accumulate the directional variance to obtain the cumulative uncertainty variance. The registration errors of the differential displacements in each period are independent, and the variance of the directional cumulative uncertainty of the cumulative deformation in each period is calculated by variance linear propagation. The weighted regression and deformation rate extraction module is used to perform weighted least squares regression with the inverse of the cumulative uncertainty variance of each period as the observation weight, the temperature time series as the independent variable and the cumulative deformation time series as the dependent variable, to deduct the deformation component of temperature effect, extract the deformation rate, and generate the deformation tracking results of the shield tunnel structure. Based on the temperature-deformation linearity coefficient and the temperature records for each period, the temperature effect deformation components for each grid cell in each direction and period are calculated and deducted from the cumulative deformation in each direction and period.

Images