A fast method for extracting tunnel central axis based on point cloud data

Abstract

The present invention discloses a method for fast extraction of the central axis of a tunnel based on point cloud data, and belongs to the technical field of tunnel engineering. The present invention innovatively combines global principal component analysis, local segmentation processing, RANSAC robust fitting, and spline interpolation technologies to achieve a comprehensive analysis of the overall structure and local details of the tunnel, innovatively proposes an efficient local splicing strategy, effectively reduces the amount of data processing, and improves algorithm efficiency. It innovatively adopts segmentation processing and elliptical cylinder model fitting to improve adaptability to complex tunnel shapes. The present invention overcomes many limitations of the prior art, while paying attention to local and global features, and provides a more efficient, more accurate, and more reliable solution for central axis extraction in tunnel engineering.

Description

Technical Field

[0001] The invention belongs to the technical field of tunnel engineering, and in particular relates to a method for quickly extracting a tunnel central axis based on point cloud data. Background Art

[0002] In the field of tunnel engineering, central axis extraction is a key technology and is of great significance to tunnel design, construction and maintenance. However, the existing technology still has many shortcomings and it is difficult to meet the requirements of accuracy, efficiency and adaptability in engineering practice. Traditional tunnel central axis measurement mostly uses total stations, levels and other measuring instruments for manual measurement. Although this method is intuitive, it has the disadvantages of being time-consuming, inefficient, limited in accuracy and susceptible to human factors. Especially in complex terrain or long tunnels, it is difficult to ensure the continuity and accuracy of the measurement, and it is impossible to obtain high-precision continuous central axis data.

[0003] With the development of 3D laser scanning technology, the method of extracting the central axis of a tunnel based on point cloud data has gradually become a research hotspot and has been widely used. Traditional methods for extracting the central axis of a tunnel mainly include technologies such as the cross-section method, the random sample consensus (RANSAC) algorithm, and the principal component analysis (PCA). These methods have some common limitations in practical applications.

[0004] First, existing methods have limited ability to handle noise and outliers in point cloud data, which can easily lead to deviations in subsequent analysis results.

[0005] Secondly, when processing large-scale point cloud data, existing algorithms often require a lot of computing resources and are difficult to meet the needs of real-time or quasi-real-time processing.

[0006] In addition, most methods assume that the tunnel cross-section is a regular shape (such as circular or elliptical), which is difficult to adapt to complex or irregular tunnel structures.

[0007] Although global analysis methods (such as principal component analysis) can quickly obtain the general direction of the tunnel, they often ignore the changes in local details, resulting in low accuracy of the extracted central axis. Iteration-based methods are susceptible to cumulative errors, especially when dealing with long tunnels, where the error will gradually increase with the increase in the number of iterations. At the same time, existing methods often focus on the single analysis of local features or global features, and lack comprehensive consideration of the overall structure and local details of the tunnel. Summary of the invention

[0008] The present invention aims to solve the technical problems existing in the prior art and provides a method for fast extraction of the central axis of a tunnel based on point cloud data. It innovatively combines global principal component analysis, local segmentation processing, RANSAC robust fitting and spline interpolation technologies to achieve a comprehensive analysis of the overall structure and local details of the tunnel, innovatively proposes an efficient local splicing strategy, effectively reduces the amount of data processing, improves algorithm efficiency, and innovatively adopts segmentation processing and elliptical cylinder model fitting to improve adaptability to complex tunnel shapes. The present invention overcomes many limitations of the prior art, pays attention to local and global features at the same time, and provides a more efficient, more accurate and more reliable solution for central axis extraction in tunnel engineering.

[0009] The present invention provides a method for quickly extracting a tunnel central axis based on point cloud data, comprising the following steps:

[0010] Input tunnel point cloud dataset;

[0011] Perform global principal component analysis on the tunnel point cloud dataset to preliminarily determine the direction and scope of the tunnel;

[0012] According to the results of global principal component analysis, the tunnel point cloud dataset is segmented into the starting segment, the middle segment and the end segment;

[0013] Perform local principal component analysis on the starting segment and the ending segment, and obtain the starting point and the ending point of the center line according to the local principal component analysis results of each segment;

[0014] Calculate the angle that the center line needs to rotate to be parallel to the X-axis α and β ; Among them, the X-axis direction is along the preset tunnel central axis direction, the Y-axis direction is the direction perpendicular to the X-axis in the horizontal plane, and the Z-axis direction is the vertical upward direction; α and β Represents the angles of rotation around the Y-axis and Z-axis respectively;

[0015] The angle to be rotated is calculated α and β , construct a coordinate transformation matrix, and use the coordinate transformation matrix to transform the tunnel point cloud dataset;

[0016] The tunnel point cloud data is segmented along the X-axis direction, and the maximum and minimum values ​​of the Y coordinate and the maximum and minimum values ​​of the Z coordinate are searched in each segment of the tunnel point cloud data to determine the edge point set of the tunnel in the three-dimensional space;

[0017] The tunnel centerline is determined by the edge point set of each segment.

[0018] Preferably, a spline interpolation smoothing process is performed on the determined tunnel centerline.

[0019] Preferably, the tunnel point cloud data set is a tunnel point cloud segment data set obtained by locally splicing point cloud data of a plurality of adjacent stations through a sliding window.

[0020] Preferably, after inputting the tunnel point cloud data set, preliminary denoising is first performed, and then the tunnel point cloud data after preliminary denoising of multiple adjacent stations are locally spliced ​​in a sliding window manner using the ICP algorithm to form a tunnel point cloud segment data set.

[0021] Preferably, the number of sliding windows is one, two or four, and the optimal number and size of windows are selected according to the specific tunnel characteristics and the quality of the point cloud data; or a sliding window with an adjustable overlap rate is used.

[0022] Preferably, the first 1 / 5 of the tunnel point cloud data is the starting segment, the last 1 / 5 of the tunnel point cloud data is the ending segment, and the remaining part of the tunnel point cloud data is the middle segment.

[0023] Preferably, obtaining the starting point and end point of the center line according to the local principal component analysis results of each segment specifically includes: in the local principal component analysis results of each segment, finding a plane perpendicular to the X-axis direction, and projecting the point cloud of the segment onto this plane, and obtaining the starting point and end point of the center line by calculating the centroid of these projection points.

[0024] Preferably, the angle of rotation α and β The specific calculation method is as follows:

[0025] α =arctan((z_end-z_start) / (x_end-x_start));

[0026] β =arctan((y_end-y_start) / (x_end-x_start));

[0027] in,

[0028] z_end is the Z-axis coordinate of the end point of the center line;

[0029] z_start is the Z-axis coordinate of the starting point of the center line;

[0030] y_end is the Y-axis coordinate of the end point of the center line;

[0031] y_start is the Y-axis coordinate of the starting point of the center line;

[0032] x_end is the X-axis coordinate of the end point of the center line;

[0033] x_start is the X-axis coordinate of the starting point of the center line.

[0034] Preferably, segmenting the tunnel point cloud data along the X-axis direction specifically includes segmenting the tunnel point cloud data at preset intervals, or performing adaptive segmentation based on tunnel curvature.

[0035] Preferably, determining the tunnel central axis by using the edge point set of each segment specifically includes: fitting an elliptical cylinder model for each segment using the obtained edge point set, determining the ellipse center of each segment, and sorting the ellipse centers obtained for each segment according to their positions on the X-axis to form the tunnel central axis; wherein the elliptical cylinder model is expressed as follows:

[0036]

[0037] in,

[0038] (y0, z0) is the center of the ellipse;

[0039] a and b are the major and minor axes of the ellipse respectively.

[0040] Preferably, the RANSAC algorithm is used for iteration in the fitting process. In each iteration, a set of edge points is randomly selected to estimate the ellipse parameters, and then the degree of fit between the elliptical cylinder model and all edge points is evaluated, and the best fit result is retained.

[0041] Compared with the prior art, the present invention has the following beneficial effects:

[0042] 1. Improved accuracy of central axis extraction: The present invention significantly improves the accuracy of central axis extraction through the coordinate system transformation technology that combines global and local principal component analysis, as well as segmented elliptical cylinder fitting based on RANSAC. Especially when dealing with complex shapes or curved tunnels, the central axis deviation is significantly reduced compared with traditional methods.

[0043] 2. Significantly improved processing efficiency: The present invention adopts a sliding window local splicing strategy, which significantly reduces the amount of data processed each time. For a tunnel point cloud data containing n points, the processing time of the present invention is reduced from O(n^2) or O(nlogn) required by the traditional global processing method to O(km^2) or O(kmlogm), where m is the number of points in each window and k is the number of windows.

[0044] 3. Improved robustness to noise and outliers: The application of the RANSAC algorithm enables the present invention to perform well when processing data containing noise and outliers.

[0045] 4. A continuous and smooth central axis is achieved: Through spline interpolation and smoothing, the central axis obtained by the present invention has better continuity and smoothness. Compared with the traditional piecewise straight line fitting method, it can more accurately reflect the actual curvature change of the tunnel, especially when dealing with tunnel sections with large curvature changes, the advantages of the present invention are more obvious.

[0046] 5. Reduced cumulative error: The sliding window local splicing strategy effectively avoids cumulative error. Compared with the traditional iterative method, the present invention shows significant advantages in processing long-distance tunnels, especially in error control near the end of the tunnel. Traditional methods often produce large deviations at the end of the tunnel due to error accumulation, while the present invention effectively controls the propagation and accumulation of errors by combining local processing with global optimization, thereby maintaining a high accuracy over the entire tunnel length. BRIEF DESCRIPTION OF THE DRAWINGS

[0047] Figure 1 The present invention is a flowchart of a method for quickly extracting a tunnel central axis based on point cloud data according to an embodiment of the present invention. DETAILED DESCRIPTION

[0048] The following will be combined with the drawings in the embodiments of the present invention to clearly and completely describe the technical solutions in the embodiments of the present invention. The described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

[0049] like Figure 1 The present invention provides a method for quickly extracting the central axis of a tunnel based on point cloud data, comprising the following steps:

[0050] Input tunnel point cloud dataset;

[0051] Perform global principal component analysis on the tunnel point cloud dataset to preliminarily determine the direction and scope of the tunnel;

[0052] According to the results of global principal component analysis, the tunnel point cloud dataset is segmented into the starting segment, the middle segment and the end segment;

[0053] Perform local principal component analysis on the starting segment and the ending segment, and obtain the starting point and the ending point of the center line according to the local principal component analysis results of each segment;

[0054] Calculate the angle that the center line needs to rotate to be parallel to the X-axis α and β ; Among them, the X-axis direction is along the preset tunnel central axis direction, the Y-axis direction is the direction perpendicular to the X-axis in the horizontal plane, and the Z-axis direction is the vertical upward direction;α and β Represents the angles of rotation around the Y-axis and Z-axis respectively;

[0055] The angle to be rotated is calculated α and β , construct a coordinate transformation matrix, and use the coordinate transformation matrix to transform the tunnel point cloud dataset;

[0056] The tunnel point cloud data is segmented along the X-axis direction, and the maximum and minimum values ​​of the Y coordinate and the maximum and minimum values ​​of the Z coordinate are searched in each segment of the tunnel point cloud data to determine the edge point set of the tunnel in the three-dimensional space;

[0057] The tunnel centerline is determined by the edge point set of each segment.

[0058] In certain specific embodiments of the present invention, a spline interpolation smoothing process is performed to determine the tunnel centerline.

[0059] In certain specific embodiments of the present invention, the tunnel point cloud dataset is a tunnel point cloud segment dataset obtained by locally splicing point cloud data of a plurality of adjacent stations through a sliding window.

[0060] In certain specific embodiments of the present invention, after the tunnel point cloud data set is input, preliminary denoising is first performed, and then the preliminary denoised tunnel point cloud data of multiple adjacent stations are locally spliced ​​in a sliding window manner using the ICP algorithm to form a tunnel point cloud segment data set.

[0061] In certain specific embodiments of the present invention, the number of sliding windows is one, two or four, and the optimal number and size of windows are selected according to specific tunnel characteristics and point cloud data quality; or a sliding window with adjustable overlap rate is used.

[0062] In certain specific embodiments of the present invention, the first 1 / 5 of the tunnel point cloud data is the starting segment, the last 1 / 5 of the tunnel point cloud data is the ending segment, and the remaining part of the tunnel point cloud data is the middle segment.

[0063] In certain specific embodiments of the present invention, obtaining the starting point and end point of the center line based on the local principal component analysis results of each segment specifically includes: in the local principal component analysis results of each segment, finding a plane perpendicular to the X-axis direction, and projecting the point cloud of the segment onto this plane, and obtaining the starting point and end point of the center line by calculating the centroid of these projection points.

[0064] In certain embodiments of the present invention, the angle of rotation α and β The specific calculation method is as follows:

[0065] α=arctan((z_end-z_start) / (x_end-x_start));

[0066] β =arctan((y_end-y_start) / (x_end-x_start));

[0067] in,

[0068] z_end is the Z-axis coordinate of the end point of the center line;

[0069] z_start is the Z-axis coordinate of the starting point of the center line;

[0070] y_end is the Y-axis coordinate of the end point of the center line;

[0071] y_start is the Y-axis coordinate of the starting point of the center line;

[0072] x_end is the X-axis coordinate of the end point of the center line;

[0073] x_start is the X-axis coordinate of the starting point of the center line.

[0074] In some specific embodiments of the present invention, segmenting the tunnel point cloud data along the X-axis direction specifically includes segmenting the tunnel point cloud data at preset intervals, or performing adaptive segmentation based on the curvature of the tunnel.

[0075] In certain specific embodiments of the present invention, determining the tunnel central axis by using the edge point set of each segment specifically includes: fitting an elliptical cylinder model for each segment using the obtained edge point set, determining the ellipse center of each segment, and sorting the ellipse centers obtained for each segment according to their positions on the X-axis to form the tunnel central axis; wherein the elliptical cylinder model is represented as follows:

[0076]

[0077] in,

[0078] (y0, z0) is the center of the ellipse;

[0079] a and b are the major and minor axes of the ellipse respectively.

[0080] In certain specific embodiments of the present invention, the RANSAC algorithm is used to iterate during the fitting process. In each iteration, a set of edge points is randomly selected to estimate the ellipse parameters, and then the degree of fit between the elliptical cylinder model and all edge points is evaluated, and the best fit result is retained.

[0081] In certain specific embodiments of the present invention, the process of constructing a coordinate transformation matrix involves a combination of rotation matrices, and a common three-dimensional rotation formula, such as the Rodriguez rotation formula, can be used for specific implementation.

[0082] Example 1

[0083] A specific embodiment of the present invention is described in detail below.

[0084] Step 1: For the preliminary denoised data, call CloudCompare's built-in ICP algorithm to perform sliding window local splicing of the point cloud data of three adjacent stations to form a smaller point cloud segment, that is, first process stations 1-2-3, then 2-3-4, then 3-4-5, and so on. This method can reduce the amount of data processed each time and ensure the continuity and smooth transition of adjacent areas. At the same time, it also avoids cumulative errors and has better consistency. When calling the ICP algorithm, based on the experience of processing tunnel point cloud data, the maximum number of iterations is set to 100 times and the distance threshold is set to 0.5 meters.

[0085] Step 2: Perform global principal component analysis on the locally spliced ​​tunnel point cloud data to determine the approximate direction and range of the tunnel.

[0086] Step 3, according to the results of global PCA, the point cloud data are divided into three main parts: the starting segment (the first 1 / 5 of the tunnel point cloud data), the middle segment (the point cloud data has poor symmetry), and the end segment (the last 1 / 5 of the tunnel point cloud data).

[0087] Step 4: Perform local PCA analysis on the starting and ending segments. In the PCA results of each segment, find a plane perpendicular to the main axis and project the points of the segment onto this plane. By calculating the centroid of these projected points, we can get the points close to the center line. These points represent the starting and ending points of the tunnel center line.

[0088] Step 5, using the obtained centerline starting point and end point, calculate the rotation angles α and β required to make the tunnel axis roughly parallel to the X-axis. α and β represent the rotation around the Y-axis and Z-axis respectively. Among them, the X-axis direction is along the preset tunnel centerline direction, the Y-axis direction is perpendicular to the X-axis direction in the horizontal plane, and the Z-axis direction is vertically upward. The calculation of the rotation angles α and β is based on the coordinate difference between the starting point and the end point. The specific calculation method is as follows:

[0089] α = arctan((z_end - z_start) / (x_end - x_start));

[0090] β= arctan((y_end - y_start) / (x_end - x_start));

[0091] in,

[0092] z_end is the Z-axis coordinate of the end point of the center line;

[0093] z_start is the Z-axis coordinate of the starting point of the center line;

[0094] y_end is the Y-axis coordinate of the end point of the center line;

[0095] y_start is the Y-axis coordinate of the starting point of the center line;

[0096] x_end is the X-axis coordinate of the end point of the center line;

[0097] x_start is the X-axis coordinate of the starting point of the center line.

[0098] Step 6: Based on the calculated rotation angle α and β , construct a coordinate transformation matrix. This matrix is ​​used to rotate the entire point cloud so that the tunnel axis is roughly parallel to the X-axis. The construction process of the coordinate transformation matrix involves the combination of rotation matrices, and the specific implementation can use common three-dimensional rotation formulas. Finally, use this coordinate transformation matrix to transform the entire point cloud. This step will retain the original three-dimensional coordinate information of the point cloud, while adjusting the tunnel axis to be roughly parallel to the X-axis.

[0099] Step 7: Segment the point cloud at a certain interval of L=1.5m along the X-axis direction, and search not only the maximum and minimum values ​​of the Y coordinate in each segmented area, but also the maximum and minimum values ​​of the Z coordinate. In this way, we obtain the edge point set of the tunnel in three-dimensional space, including edge points in the four directions of up, down, left, and right, and determine the central axis of the tunnel through these edge points, including:

[0100] Step 71, use the obtained edge point set to fit an elliptical cylinder model for each segment. This model can be expressed as

[0101]

[0102] in,

[0103] (y0, z0) is the center of the ellipse;

[0104] a and b are the major and minor axes of the ellipse respectively.

[0105] The RANSAC algorithm is used iteratively during the fitting process to improve the robustness to noise and outliers.

[0106] In each iteration, a subset of edge points is randomly selected to estimate the ellipse parameters, and then the fit of the model to all edge points is evaluated, and the best fit result is retained, effectively handling the noise and outliers that may exist in the point cloud. For each segment, the ellipse center (y0, z0) of the segment is determined by fitting the ellipse equation, which represents the center of the tunnel cross section of the segment.

[0107] Step 72, sorting the ellipse center points obtained from each segment according to their positions on the X-axis, and the sorted center point sequences form the three-dimensional central axis of the tunnel.

[0108] In order to further improve the quality of the central axis, the obtained central axis is smoothed by spline interpolation to reduce the possible noise influence. The final central axis is represented as a series of three-dimensional points (x, y, z), where x is the coordinate along the length of the tunnel, and y and z are the center coordinates obtained by ellipse fitting.

[0109] This invention innovatively combines global PCA with local PCA, which not only captures the overall direction of the tunnel, but also adapts to local subtle changes. This method obtains a more accurate tunnel axis direction by performing local PCA analysis on the starting and ending segments, thereby achieving accurate coordinate system conversion of the entire point cloud data.

[0110] The present invention organically combines the steps of local splicing, coordinate system conversion, elliptical cylinder fitting and spline interpolation smoothing to form a complete and efficient processing flow. This multi-step fusion method can effectively process large-scale and complex structure tunnel point cloud data, significantly improving processing efficiency while ensuring high accuracy.

[0111] The above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. For those skilled in the art, the present invention may have various modifications and variations. Any modification, equivalent replacement, improvement, etc. made within the spirit and principle of the present invention are included in the protection scope of the present invention.

Claims

1. A method for quickly extracting the central axis of a tunnel based on point cloud data, characterized in that: The steps include: Input tunnel point cloud dataset; Perform global principal component analysis on the tunnel point cloud dataset to preliminarily determine the direction and scope of the tunnel; According to the results of global principal component analysis, the tunnel point cloud dataset is segmented into the starting segment, the middle segment and the end segment; Perform local principal component analysis on the starting segment and the ending segment, and obtain the starting point and the ending point of the center line according to the local principal component analysis results of each segment; Calculate the angle that needs to be rotated to make the center line parallel to the X-axis α and β ; Among them, the X-axis direction is along the preset tunnel central axis direction, the Y-axis direction is the direction perpendicular to the X-axis in the horizontal plane, and the Z-axis direction is the vertical upward direction; α and β Represents the angles of rotation around the Y-axis and Z-axis respectively; The angle to be rotated is calculated α and β , construct a coordinate transformation matrix, and use the coordinate transformation matrix to transform the tunnel point cloud dataset; The tunnel point cloud data is segmented along the X-axis direction, and the maximum and minimum values ​​of the Y coordinate and the maximum and minimum values ​​of the Z coordinate are searched in each segment of the tunnel point cloud data to determine the edge point set of the tunnel in the three-dimensional space; The tunnel centerline is determined by the edge point set of each segment.

2. The method for quickly extracting the central axis of a tunnel based on point cloud data according to claim 1 is characterized in that: Spline interpolation smoothing is performed to determine the tunnel centerline.

3. The method for quickly extracting the central axis of a tunnel based on point cloud data according to claim 1 is characterized in that: After the tunnel point cloud dataset is input, preliminary denoising is performed first, and then the ICP algorithm is used to perform sliding window local splicing of the preliminary denoised tunnel point cloud data of multiple adjacent stations to form a tunnel point cloud segment dataset.

4. The method for quickly extracting the central axis of a tunnel based on point cloud data according to claim 3 is characterized in that: The sliding windows are one, two or four, and the optimal number and size of windows are selected according to the specific tunnel characteristics and point cloud data quality; or a sliding window with adjustable overlap rate is used.

5. The method for quickly extracting the central axis of a tunnel based on point cloud data according to claim 1 is characterized in that: The first 1 / 5 of the tunnel point cloud data is the starting segment, the last 1 / 5 of the tunnel point cloud data is the ending segment, and the rest of the tunnel point cloud data is the middle segment.

6. The method for quickly extracting the central axis of a tunnel based on point cloud data according to claim 5 is characterized in that: The starting point and the end point of the center line are obtained according to the local principal component analysis results of each segment. Specifically, in the local principal component analysis results of each segment, a plane perpendicular to the X-axis direction is found, and the point cloud of the segment is projected onto the plane. The starting point and the end point of the center line are obtained by calculating the centroid of these projection points.

7. The method for quickly extracting the central axis of a tunnel based on point cloud data according to claim 1 is characterized in that: Angle of rotation α and β The specific calculation method is as follows: α =arctan((z_end-z_start) / (x_end-x_start)); β =arctan((y_end-y_start) / (x_end-x_start)); in, z_end is the Z-axis coordinate of the end point of the center line; z_start is the Z-axis coordinate of the starting point of the center line; y_end is the Y-axis coordinate of the end point of the center line; y_start is the Y-axis coordinate of the starting point of the center line; x_end is the X-axis coordinate of the end point of the center line; x_start is the X-axis coordinate of the starting point of the center line.

8. The method for quickly extracting the central axis of a tunnel based on point cloud data according to claim 1, characterized in that: Segmenting the tunnel point cloud data along the X-axis direction specifically includes segmenting the tunnel point cloud data at preset intervals, or performing adaptive segmentation based on the curvature of the tunnel.

9. The method for quickly extracting the central axis of a tunnel based on point cloud data according to claim 1, characterized in that: Determining the tunnel center axis by using the edge point set of each segment specifically includes: fitting an elliptical cylinder model for each segment using the obtained edge point set, determining the ellipse center of each segment, and sorting the ellipse centers obtained for each segment according to their positions on the X-axis to form the tunnel center axis; wherein the elliptical cylinder model is expressed as follows: in, (y0, z0) is the center of the ellipse; a and b are the major and minor axes of the ellipse respectively.

10. The method for quickly extracting the central axis of a tunnel based on point cloud data according to claim 9, characterized in that: The RANSAC algorithm is used to iterate during the fitting process. In each iteration, a set of edge points is randomly selected to estimate the ellipse parameters, and then the degree of fit between the elliptical cylinder model and all edge points is evaluated, and the best fitting result is retained.

Images