A shield tunnel ring joint segmentation fitting method based on gradient value feature map density clustering

Abstract

A method for segmented fitting of shield tunnel circumferential joints based on gradient value feature map density clustering includes the following steps: unfolding the tunnel's three-dimensional laser point cloud data to generate a two-dimensional grayscale image; calculating the vertical gradient values ​​of all pixels in the grayscale image, and generating a positive gradient value feature map A, a negative gradient value feature map B, and a comprehensive gradient value feature map C based on the magnitude and sign of the gradient values; determining the initial upper and lower edge row coordinates of the circumferential joint based on Figures A and B, and taking the middle row as the initial circumferential joint position; extracting pixels from the circumferential joint region in Figure C based on the initial circumferential joint position, performing cluster analysis, and dividing the circumferential joint region into several regions along the horizontal direction; identifying the upper and lower edges of the circumferential joint in each region, and calculating the center point as the center point of the circumferential joint in that region; sequentially connecting each point and extending to the image edge to form a broken line segment, thereby achieving segmented fitting. This invention achieves high-precision, adaptive segmented fitting of non-linear circumferential joints in shield tunnels, improving the accuracy and robustness of circumferential joint detection.

Description

Technical Field

[0001] This invention relates to a method for fitting the segmented joints of shield tunnels based on gradient value feature map density clustering, belonging to the field of tunnel engineering detection technology. Background Technology

[0002] Currently, as a key component of urban rail transit systems, the structural safety of subway shield tunnels directly impacts operational safety and long-term reliability. Shield tunnels are typically assembled from several precast concrete segments, and the joints formed between adjacent segments are a typical structural feature and also prone to common defects such as leakage, misalignment, and deformation. Therefore, the accurate identification and location of these joints has become a fundamental aspect of tunnel structural health monitoring, defect diagnosis, and scientific maintenance management, playing a crucial role in ensuring the safe operation of tunnels throughout their entire lifespan.

[0003] Traditional circumferential joint detection mainly relies on manual inspection, which has limitations such as low efficiency, strong subjectivity, and difficulty in achieving quantitative and systematic evaluation. In recent years, with the popularization of 3D laser scanning technology, automated circumferential joint detection methods based on point cloud data have become mainstream. Existing methods mainly rely on 3D tunnel point cloud data, which is usually unfolded into a 2D planar point cloud or further converted into a grayscale image, and then the circumferential joint is identified in the following ways: locating the circumferential joint by statistically analyzing the local extrema of intensity values ​​or elevation data in the point cloud cross-section; indirectly inferring the location of the circumferential joint based on the arrangement characteristics of bolt holes in the tunnel inner wall; and using edge detection operators to extract the edges of the grayscale image and combining them with line detection methods such as Hough transform to identify the circumferential joint.

[0004] While existing methods have achieved automated processing, they generally suffer from the flaw of treating the unfolded circumferential joint as an ideal straight line. In actual engineering, due to factors such as segment assembly errors or local deformation, the unfolded circumferential joint is not a perfect straight line. Therefore, when dealing with the non-straight circumferential joints commonly found in actual tunnels, the problem of identification position offset often occurs, making it difficult to meet the requirements of high-precision and robust engineering inspection. Summary of the Invention

[0005] The purpose of this invention is to provide a segmented fitting method for shield tunnel circumferential joints based on gradient value feature map density clustering. This method can achieve high-precision, adaptive segmented fitting of non-linear circumferential joints in shield tunnels, improve the accuracy and robustness of circumferential joint detection in complex engineering environments, and ensure the safe service of the tunnel throughout its entire life cycle.

[0006] To achieve the above objectives, this invention provides a method for fitting the segmented joint of a shield tunnel annular joint based on gradient value feature map density clustering, comprising the following steps:

[0007] S1. Unfold the tunnel's three-dimensional laser point cloud data and generate a two-dimensional grayscale image;

[0008] S2. Calculate the vertical gradient value of all pixels in the grayscale image, and generate three gradient value feature maps based on the magnitude and sign of the gradient value, namely positive gradient value feature map A, negative gradient value feature map B, and comprehensive gradient value feature map C.

[0009] S3. Perform gray-level accumulation statistics, interval peak detection, and global consistency verification on the positive gradient value feature map A and the negative gradient value feature map B respectively to determine the row coordinates of the initial upper and lower edges of the ring seam, and take the middle row as the initial ring seam position.

[0010] S4. Based on the initial circumferential seam position, extract the pixels of the circumferential seam region in the comprehensive gradient value feature map C, and use the density clustering algorithm to perform clustering analysis on it, adaptively dividing the circumferential seam region into several regions along the horizontal direction.

[0011] S5. Within each region divided in S4, the upper and lower edges of the circumferential seam are identified based on the positive gradient value feature map A and the negative gradient value feature map B. The center point is calculated as the center point of the circumferential seam in that region. The points are connected sequentially according to the principle of similar vertical coordinates and extended to the edge of the image to form a polyline segment, thereby realizing the segmented fitting of the circumferential seam.

[0012] Furthermore, the specific process of S1 is as follows:

[0013] S1.1 Projecting the 3D point cloud data onto the tunnel cross-section From the plane, we obtain the corresponding two-dimensional point set. In the coordinate system of this section, the origin (0, 0) corresponds to the center of the scanner, and the center point of the tunnel cross section is set. The coordinates are (0, By calculating the mean of the vertical coordinates of each point in this two-dimensional point set, the vertical coordinates of the center point of the cross-section can be obtained. ;

[0014] S1.2 The 3D tunnel point cloud is acquired by a laser scanner mounted on a track trolley. During unfolding, the 3D tunnel point cloud is unfolded along the top axis of the tunnel according to the principle of cylindrical projection and under the constraint of the scanner's viewpoint. The specific steps are as follows:

[0015] S1.2-1. In the cross-sectional coordinate system, connect the origin. A ray is formed with any point, and its parametric equation is:

[0016] ;

[0017] Based on the center coordinates of the cross section and tunnel design radius Establish the equation of the tunnel cross section circle Substituting the ray equation into the cross-sectional circle equation yields the following about... The quadratic equation: ,in , , Solve the quadratic equation and choose The solution calculates the intersection point of the ray and the cross-sectional circle. coordinates ;

[0018] S1.2-2, Let line segment along The axis extends in the positive direction and intersects the cross-sectional circle at point . Calculation points With point The connecting line is the same as the straight line. The included angle ;

[0019] Calculate the arc length at any point And calculate the point. Intersection with the cross-section circle Distance between Modify the coordinates of any point in the 3D tunnel point cloud to This allows the 3D tunnel point cloud to unfold along the top of the tunnel into a 2.5D point cloud;

[0020] S1.3. Convert the 2.5D point cloud into a grayscale image using the raster method. The specific steps are as follows:

[0021] S1.3-1 Read the planar coordinates of each point in the unfolded point cloud and the corresponding strength value Constructing a two-dimensional point cloud dataset Based on the set physical size of the pixel unit The unit is meters per pixel, which determines the size of the grayscale image: , ,in , , , Point clouds coordinate, The maximum and minimum values ​​of the coordinates. Indicates the image width. Indicates the image height. Indicates rounding up;

[0022] S1.3-2, Establish the mapping relationship from point cloud to pixel: column index row index ,in, , , Indicates rounding down;

[0023] S1.3-3. Based on the mapping relationship from point cloud to pixel, assign the point cloud intensity value to the pixel. If multiple point clouds fall into the same pixel, calculate the average of the intensity values ​​of these multiple point clouds. Assign this pixel, where, Is falling into pixels The number of point clouds, if no point cloud falls into a pixel If the strength value is 0, then the strength value matrix is ​​formed;

[0024] S1.3-4. Normalize the intensity value matrix to the interval [0, 255]. The calculation formula is as follows: ,in, Represents pixels Normalized grayscale values This indicates rounding to the nearest integer. and These represent the maximum and minimum values ​​of the intensity value matrix, respectively, and the generated size is... An 8-bit grayscale image.

[0025] Furthermore, the specific process of generating the gradient value feature map in S2 is as follows:

[0026] S2.1, Using Gaussian kernels Perform a convolution operation on the generated grayscale image to obtain the filtered image:

[0027] ;

[0028] In the formula, For Gaussian kernel function in The weight value at that point, This is the offset within the kernel relative to the center pixel. To customize the standard deviation, Where is the kernel radius, and the kernel size is... ;

[0029] S2.2. Use a vertical edge detection operator to calculate the gradient of the filtered image to obtain the vertical gradient features that reflect the gray-level changes of the image. The calculation formula is as follows:

[0030] ;

[0031] in, For pixels The vertical gradient value, This is the offset within the kernel relative to the center pixel. For edge detection operators in the vertical direction, the convolution kernel is in The weight value at that point, The radius of the convolution kernel;

[0032] S2.3, Based on the mean of the absolute values ​​of the gradient values Set a gradient magnitude threshold Where scale is a user-defined parameter; three feature maps are generated based on the sign and magnitude of the gradient values, as follows:

[0033] (1) Only retain For each pixel, its gradient value is used as the gray value of that pixel, and the gray values ​​of other pixels are all set to 0. The gradient values ​​are then normalized and mapped to the [0,255] interval to generate a positive gradient value feature map A, which is used to identify the lower edge of the circumferential seam.

[0034] (2) Only retain For each pixel, the absolute value of its gradient is used as the gray value of that pixel, and the gray values ​​of other pixels are all set to 0. The gradient values ​​are normalized and mapped to the interval [0,255] to generate a negative gradient value feature map B, which is used to identify the upper edge of the circumferential seam.

[0035] (3) Retain For all pixels, the absolute value of their gradient values ​​is taken as the gray value of that pixel, and the gray values ​​of other pixels are set to 0. The gradient values ​​are then normalized and mapped to the interval [0, 255] to generate a comprehensive gradient value feature map C.

[0036] Furthermore, the specific process of S3 is as follows:

[0037] S3.1 Calculate the sum of gray values ​​of pixels in each row of the positive gradient feature map A and the negative gradient feature map B. , ,in Indicates the image at position The grayscale value at that location;

[0038] S3.2, Based on the design width of the tunnel segment ring Divide the image vertically into A series of continuous equal-height intervals, based on the set physical size of the pixel unit. Determine the pixel height of each interval as Total number of intervals ;

[0039] S3.3, Find the intervals respectively Largest row As candidate edge rows for this interval, Perform this operation on the positive gradient value feature map A and the negative gradient value feature map B to obtain the candidate lower edge row set and the candidate upper edge row set;

[0040] S3.4, Calculate the first... Candidate edge rows of each interval absolute difference from other candidate edge rows Count the number of intervals whose differences satisfy the following conditions. :

[0041] , making ;in, It is the set of positive integers. The pixel tolerance is the absolute difference between candidate edge rows. Defined allowable physical distance deviation;

[0042] like If the candidate behavior is verified by the edge row in that interval, then the candidate behavior is verified by the edge row. ;like Then within that interval, according to Starting from the largest to the smallest, try other rows as candidates in turn, repeat the above verification process until a row that meets the conditions is found;

[0043] S3.5 For each interval, check its valid upper edge row. With the bottom edge line Take the arithmetic mean to obtain the row coordinates of the circumferential seam center in this interval. As the initial circumferential seam position.

[0044] Furthermore, the specific process of S4 is as follows:

[0045] S4.1, Using the initial circumferential seam position as the coordinate system. Centered on the design width of the circumferential seam Set a physical distance threshold in the vertical direction. (usually) ), and its corresponding pixel distance threshold Only retain image C Within the range of pixels, the grayscale value of all other pixels is set to 0, generating a gradient value feature map D;

[0046] S4.2, Based on the physical width of the circumferential seam and pixel unit physical size Set the neighborhood radius and the minimum number of samples within the neighborhood radius for the DBSCAN algorithm. The DBSCAN algorithm is used to perform cluster analysis on the non-zero pixels in image D, and a minimum pixel count threshold for each cluster is set. The number of pixels deleted is less than Clusters that generate clusters containing Image F of the class result;

[0047] S4.3 For each cluster in image F Calculate the boundary of its enclosing rectangular region:

[0048] , ;

[0049] The circumferential seam area is adaptively divided into several rectangular segments along the transverse direction. .

[0050] Furthermore, the specific process of S5 is as follows:

[0051] S5.1, placing those located in the same Rectangular areas within the specified range are grouped into the same vertically similar group. ;

[0052] S5.2, in each longitudinally similar group Internally, all rectangular regions are arranged in ascending order according to their column coordinates, forming an ordered sequence. ,in , The number of rectangular regions within the group; for each rectangular region Define its adjacent segment set as: ;

[0053] S5.3, Based on the rectangular area Rectangular regions are extracted from the positive gradient feature map A and the negative gradient feature map B, respectively, and the cumulative grayscale value of each row of pixels in the rectangular region is calculated. Find each region Largest row As candidate edge rows for this region; perform this operation on all rectangular regions in the positive gradient value feature map A and the negative gradient value feature map B respectively to obtain the candidate lower edge row set and the candidate upper edge row set;

[0054] S5.4, According to the design width of the circumferential seam and pixel unit physical size Get the pixel width of the circumferential seam To ensure the spatial continuity and smoothness of the final fitted circumferential seam, edge position constraints are introduced between adjacent segments. For each longitudinally similar group... The candidate edge rows of the rectangular region within the area are optimized using the following steps:

[0055] (1) For segmentation The difference between its candidate edge row and the candidate edge rows of two adjacent segments satisfies the following condition: If not satisfied, then within that segment, according to... Starting from the largest to the smallest, try other rows as candidate rows one by one until a row that meets the criteria is found, and then update. If none of the rows in a segment meet the condition, then the segment is removed.

[0056] (2) After the above optimization, if the four adjacent segmented regions within the group Meet the conditions Then determine For the areas that need adjustment;

[0057] (3) For each detected The regions will be adjusted as follows:

[0058] exist Within the segmented area, according to Try other rows as candidate rows in descending order of size until a match is found.

[0059] satisfy The line, and update If no row meets the criteria, delete it. Segmented regions;

[0060] like If the segmented region is not deleted, then in Within the segmented area, according to Try other rows in descending order of size as candidate rows until a satisfactory row is found. The line, and update If no row meets the criteria, delete it. Segmented regions;

[0061] like If the segmented region has been deleted, then in Within the segmented area, according to Try other rows in descending order of size as candidate rows until a satisfactory row is found. The line, and update If no row meets the criteria, delete it. Segmented regions;

[0062] The above steps yield the final top edge row for each segmented region. and bottom edge row ;

[0063] S5.5, For each segment Use its final top edge row and bottom edge row Calculate the center point of the circumferential joint. Coordinates: , ;

[0064] S5.6, in each group Internally, adjacent suture points are connected sequentially to form a broken line segment that extends to the edge of the image, creating a continuous and smooth suture that runs through the entire image. , is represented as:

[0065] .

[0066] This invention unfolds the three-dimensional laser point cloud data of the tunnel and generates a two-dimensional grayscale image; calculates the vertical gradient value of all pixels in the grayscale image, and generates three gradient value feature maps based on the magnitude and sign of the gradient values: positive gradient value feature map A, negative gradient value feature map B, and comprehensive gradient value feature map C; performs grayscale accumulation statistics, interval peak detection, and global consistency verification on positive gradient value feature map A and negative gradient value feature map B respectively to determine the row coordinates of the initial upper and lower edges of the ring joint, and takes the middle row as the initial ring joint position; based on the initial ring joint position, extracts the pixels of the ring joint region in the comprehensive gradient value feature map C, and performs cluster analysis on them using a density clustering algorithm to adaptively divide the ring joint region into several regions along the horizontal direction; within each divided region, identifies the upper and lower edges of the ring joint in that region based on the positive gradient value feature map A and negative gradient value feature map B, and calculates the center point as the ring joint center point of that region; connects each point sequentially according to the principle of similar vertical coordinates and extends to the image edge to form a broken line segment, realizing the segmented fitting of the ring joint. This invention does not rely on prior assumptions about the morphology of the circumferential joint, and achieves effective identification of non-linear circumferential joints. This enables high-precision and robust automated detection in complex tunnel environments, providing a solid technical foundation for structural health monitoring and disease diagnosis, and ensuring safe service throughout the entire life cycle of the tunnel. Attached Figure Description

[0067] Figure 1 This is a flowchart of the method of the present invention;

[0068] Figure 2 This is a schematic diagram of the three-dimensional point cloud unfolding and two-dimensional grayscale image generation of the present invention;

[0069] Figure 3 This is the gradient value feature map generated in this invention;

[0070] Figure 4 This is a schematic diagram of the initial circumferential seam positioning result of the present invention;

[0071] Figure 5 This is a schematic diagram of the lateral division of the annular seam region based on density clustering in this invention;

[0072] Figure 6 This is a schematic diagram of the edge recognition results within the segment in this invention;

[0073] Figure 7 This is a schematic diagram of the fitting result of the circumferential seam segmentation in this invention;

[0074] Figure 8 This is a comparison of the results of this method and the image grayscale thresholding method in identifying non-linear circumferential seams. Detailed Implementation

[0075] The invention will now be further described with reference to the accompanying drawings.

[0076] like Figure 1 As shown, a method for fitting the segmented joint of a shield tunnel based on gradient value feature map density clustering includes the following steps:

[0077] S1. Unfold the tunnel's three-dimensional laser point cloud data and generate a two-dimensional grayscale image;

[0078] S2. Calculate the vertical gradient value of all pixels in the grayscale image, and generate three gradient value feature maps based on the magnitude and sign of the gradient value, namely positive gradient value feature map A, negative gradient value feature map B, and comprehensive gradient value feature map C.

[0079] S3. Perform grayscale accumulation statistics, interval peak detection, and global consistency verification on the positive gradient value feature map A and the negative gradient value feature map B respectively to determine the initial upper and lower edge row coordinates of the annular fracture, and take the middle row as the initial annular fracture position; in this embodiment, the finite element method is used to construct the geometric model of the rock, such as Figure 2 As shown, a cylindrical model was established, which was obtained by rotating a rectangular domain with a height of 100 mm and a diameter of 50 mm around the axis of symmetry, and has two coordinate directions, r and z.

[0080] S4. Based on the initial circumferential seam position, extract the pixels of the circumferential seam region in the comprehensive gradient value feature map C, and use the density clustering algorithm to perform clustering analysis on it, adaptively dividing the circumferential seam region into several regions along the horizontal direction.

[0081] S5. Within each region divided in S4, the upper and lower edges of the circumferential seam are identified based on the positive gradient value feature map A and the negative gradient value feature map B. The center point is calculated as the center point of the circumferential seam in that region. The points are connected sequentially according to the principle of similar vertical coordinates and extended to the edge of the image to form a polyline segment, thereby realizing the segmented fitting of the circumferential seam.

[0082] like Figure 2 As shown, the specific process of S1 is as follows:

[0083] S1.1 Projecting the 3D point cloud data onto the tunnel cross-section From the plane, we obtain the corresponding two-dimensional point set. In the coordinate system of this section, the origin (0, 0) corresponds to the center of the scanner, and the center point of the tunnel cross section is set. The coordinates are (0, By calculating the mean of the vertical coordinates of each point in this two-dimensional point set, the vertical coordinates of the center point of the cross-section can be obtained. ;

[0084] S1.2 The 3D tunnel point cloud is acquired by a laser scanner mounted on a track trolley. During unfolding, the 3D tunnel point cloud is unfolded along the top axis of the tunnel according to the principle of cylindrical projection and under the constraint of the scanner's viewpoint. The specific steps are as follows:

[0085] S1.2-1. In the cross-sectional coordinate system, connect the origin. A ray is formed with any point, and its parametric equation is:

[0086] ;

[0087] Based on the center coordinates of the cross section and tunnel design radius Establish the equation of the tunnel cross section circle Substituting the ray equation into the cross-sectional circle equation yields the following about... The quadratic equation: ,in , , Solve the quadratic equation and choose The solution calculates the intersection point of the ray and the cross-sectional circle. coordinates ;

[0088] S1.2-2, Let line segment along The axis extends in the positive direction and intersects the cross-sectional circle at point . Calculation points With point The connecting line is the same as the straight line. The included angle ;

[0089] Calculate the arc length at any point And calculate the point. Intersection with the cross-section circle Distance between Modify the coordinates of any point in the 3D tunnel point cloud to This allows the 3D tunnel point cloud to unfold along the top of the tunnel into a 2.5D point cloud;

[0090] S1.3. Convert the 2.5D point cloud into a grayscale image using the raster method. The specific steps are as follows:

[0091] S1.3-1 Read the planar coordinates of each point in the unfolded point cloud and the corresponding strength value Constructing a two-dimensional point cloud dataset Based on the set physical size of the pixel unit The unit is meters per pixel, which determines the size of the grayscale image: , ,in , , , Point clouds coordinate, The maximum and minimum values ​​of the coordinates. Indicates the image width. Indicates the image height. Indicates rounding up;

[0092] S1.3-2, Establish the mapping relationship from point cloud to pixel: column index row index ,in, , , Indicates rounding down;

[0093] S1.3-3. Based on the mapping relationship from point cloud to pixel, assign the point cloud intensity value to the pixel. If multiple point clouds fall into the same pixel, calculate the average of the intensity values ​​of these multiple point clouds. Assign this pixel, where, Is falling into pixels The number of point clouds, if no point cloud falls into a pixel If the strength value is 0, then the strength value matrix is ​​formed;

[0094] S1.3-4. Normalize the intensity value matrix to the interval [0, 255]. The calculation formula is as follows: ,in, Represents pixels Normalized grayscale values This indicates rounding to the nearest integer. and These represent the maximum and minimum values ​​of the intensity value matrix, respectively, and the generated size is... An 8-bit grayscale image.

[0095] Furthermore, the specific process of generating the gradient value feature map in S2 is as follows:

[0096] S2.1, Using Gaussian kernels Perform a convolution operation on the generated grayscale image to obtain the filtered image:

[0097] ;

[0098] In the formula, For Gaussian kernel function in The weight value at that point, This is the offset within the kernel relative to the center pixel. To customize the standard deviation, Where is the kernel radius, and the kernel size is... ;

[0099] S2.2. Use a vertical edge detection operator to calculate the gradient of the filtered image to obtain the vertical gradient features that reflect the gray-level changes of the image. The calculation formula is as follows:

[0100] ;

[0101] in, For pixels The vertical gradient value, This is the offset within the kernel relative to the center pixel. For edge detection operators in the vertical direction, the convolution kernel is in The weight value at that point, The radius of the convolution kernel;

[0102] S2.3, Based on the mean of the absolute values ​​of the gradient values Set a gradient magnitude threshold Where scale is a user-defined parameter; three feature maps are generated based on the sign and magnitude of the gradient values, as follows:

[0103] (1) Only retain For each pixel, its gradient value is used as its grayscale value, and the grayscale values ​​of other pixels are all set to 0. The gradient values ​​are then normalized and mapped to the [0, 255] interval to generate a positive gradient feature map A. Rotating it 90 degrees counterclockwise... Figure 3 (a) shows the lower edge of the circumferential seam;

[0104] (2) Only retain For each pixel, the absolute value of its gradient is used as its grayscale value, and the grayscale values ​​of other pixels are all set to 0. The gradient values ​​are then normalized and mapped to the [0, 255] interval to generate a negative gradient feature map B. Rotating it 90 degrees counterclockwise, as shown... Figure 3 (b) shows the upper edge used to identify the circumferential seam;

[0105] (3) Retain For all pixels, the absolute value of their gradient is used as the grayscale value of that pixel, and the grayscale values ​​of other pixels are set to 0. The gradient values ​​are then normalized and mapped to the interval [0, 255] to generate a comprehensive gradient feature map C. Rotating it 90 degrees counterclockwise, as shown... Figure 3 As shown in (c).

[0106] The specific process of S3 is as follows:

[0107] S3.1 Calculate the sum of gray values ​​of pixels in each row of the positive gradient feature map A and the negative gradient feature map B. , ,in Indicates the image at position The grayscale value at that location;

[0108] S3.2, Based on the design width of the tunnel segment ring Divide the image vertically into A series of continuous equal-height intervals, based on the set physical size of the pixel unit. Determine the pixel height of each interval as Total number of intervals ;

[0109] S3.3, Find the intervals respectively Largest row As candidate edge rows for this interval, Perform this operation on the positive gradient value feature map A and the negative gradient value feature map B to obtain the candidate lower edge row set and the candidate upper edge row set;

[0110] S3.4, Calculate the first... Candidate edge rows of each interval absolute difference from other candidate edge rows Count the number of intervals whose differences satisfy the following conditions. :

[0111] , making ;in, It is the set of positive integers. The pixel tolerance is the absolute difference between candidate edge rows. Defined allowable physical distance deviation;

[0112] like If the candidate behavior is verified by the edge row in that interval, then the candidate behavior is verified by the edge row. ;like Then within that interval, according to Starting from the largest to the smallest, try other rows as candidates in turn, repeat the above verification process until a row that meets the conditions is found;

[0113] S3.5 For each interval, check its valid upper edge row. With the bottom edge line Take the arithmetic mean to obtain the row coordinates of the circumferential seam center in this interval. As the initial circumferential seam position, such as Figure 4 The diagram shown is a schematic of the initial circumferential seam positioning result (rotated 90 degrees counterclockwise).

[0114] The specific process of S4 is as follows:

[0115] S4.1, Using the initial circumferential seam position as the coordinate system. Centered on the design width of the circumferential seam Set a physical distance threshold in the vertical direction. (usually) ), and its corresponding pixel distance threshold Only retain image C Within the range of pixels, the grayscale value of all other pixels is set to 0, generating a gradient value feature map D;

[0116] S4.2, Based on the physical width of the circumferential seam and pixel unit physical size Set the neighborhood radius of the DBSCAN algorithm to Minimum number of samples within the neighborhood radius The DBSCAN algorithm is used to perform cluster analysis on the non-zero pixels in image D, and a minimum pixel count threshold for each cluster is set. The number of pixels deleted is less than Clusters that generate clusters containing Image F of the class result;

[0117] S4.3 For each cluster in image F Calculate the boundary of its enclosing rectangular region:

[0118] , ;

[0119] The circumferential seam area is adaptively divided into several rectangular segments along the transverse direction. , to obtain Figure 5 The diagram shown illustrates the lateral division of the circumferential seam region based on density clustering (rotated 90 degrees counterclockwise).

[0120] The specific process of S5 is as follows:

[0121] S5.1, placing those located in the same Rectangular areas within the specified range are grouped into the same vertically similar group. ;

[0122] S5.2, in each longitudinally similar group Internally, all rectangular regions are arranged in ascending order according to their column coordinates, forming an ordered sequence. ,in , The number of rectangular regions within the group; for each rectangular region

[0123] Define its adjacent segment set as: ;

[0124] S5.3, Based on the rectangular area Rectangular regions are extracted from the positive gradient feature map A and the negative gradient feature map B, respectively, and the cumulative grayscale value of each row of pixels in the rectangular region is calculated. Find each region Largest row As candidate edge rows for this region; perform this operation on all rectangular regions in the positive gradient value feature map A and the negative gradient value feature map B respectively to obtain the candidate lower edge row set and the candidate upper edge row set;

[0125] S5.4, According to the design width of the circumferential seam and pixel unit physical size Get the pixel width of the circumferential seam To ensure the spatial continuity and smoothness of the final fitted circumferential seam, edge position constraints are introduced between adjacent segments. For each longitudinally similar group... The candidate edge rows of the rectangular region within the area are optimized using the following steps:

[0126] (1) For segmentation The difference between its candidate edge row and the candidate edge rows of two adjacent segments satisfies the following condition: If not satisfied, then within that segment, according to... Starting from the largest to the smallest, try other rows as candidate rows one by one until a row that meets the criteria is found, and then update. If none of the rows in a segment meet the condition, then the segment is removed.

[0127] (2) After the above optimization, if the four adjacent segmented regions within the group Meet the conditions Then determine For the areas that need adjustment;

[0128] (3) For each detected The regions will be adjusted as follows:

[0129] exist Within the segmented area, according to Try other rows in descending order of size as candidate rows until a satisfactory row is found. The line, and update If no row meets the criteria, delete it. Segmented regions;

[0130] like If the segmented region is not deleted, then in Within the segmented area, according to Try other rows in descending order of size as candidate rows until a satisfactory row is found. The line, and update If no row meets the criteria, delete it. Segmented regions;

[0131] like If the segmented region has been deleted, then in Within the segmented area, according to Try other rows in descending order of size as candidate rows until a satisfactory row is found. The line, and update If no row meets the criteria, delete it. Segmented regions;

[0132] The above steps yield the final top edge row for each segmented region. and bottom edge row ,like Figure 6 As shown, where, Figure 6 In the image, (a) is the result of the lower edge recognition (rotated 90 degrees counterclockwise), and (b) is the result of the upper edge recognition (rotated 90 degrees counterclockwise).

[0133] S5.5, For each segment Use its final top edge row and bottom edge row Calculate the center point of the circumferential joint. Coordinates: , ;

[0134] S5.6, such as Figure 7 As shown, in each group Internally, adjacent suture points are connected sequentially to form a broken line segment that extends to the edge of the image, creating a continuous and smooth suture that runs through the entire image. , is represented as:

[0135] .

[0136] To verify the effectiveness and superiority of the method of this invention, it was compared with the traditional image grayscale thresholding method on the same tunnel point cloud grayscale image for ring seam identification, such as... Figure 8 As shown, where, Figure 8 Image (a) shows the circumferential seam identification result obtained using the method of the present invention (rotated 90 degrees counterclockwise). Figure 8 (b) shows the recognition result of the image grayscale thresholding method (rotated 90 degrees counterclockwise). As can be seen from the comparison, the method of the present invention can more accurately locate the position of the non-linear circumferential seam after unfolding, thus improving the accuracy and robustness of circumferential seam detection in complex engineering environments.

Claims

1. A method for segmented fitting of shield tunnel annular joints based on gradient value feature map density clustering, characterized in that, Includes the following steps: S1. Unfold the tunnel's three-dimensional laser point cloud data and generate a two-dimensional grayscale image; S2. Calculate the vertical gradient value of all pixels in the grayscale image, and generate three gradient value feature maps based on the magnitude and sign of the gradient value, namely positive gradient value feature map A, negative gradient value feature map B, and comprehensive gradient value feature map C. S3. Perform gray-level accumulation statistics, interval peak detection, and global consistency verification on the positive gradient value feature map A and the negative gradient value feature map B respectively to determine the row coordinates of the initial upper and lower edges of the ring seam, and take the middle row as the initial ring seam position. S4. Based on the initial circumferential seam position, extract the pixels of the circumferential seam region in the comprehensive gradient value feature map C, and use the density clustering algorithm to perform clustering analysis on it, adaptively dividing the circumferential seam region into several regions along the horizontal direction. S5. Within each region divided in S4, the upper and lower edges of the circumferential seam are identified based on the positive gradient value feature map A and the negative gradient value feature map B. The center point is calculated as the center point of the circumferential seam in that region. The points are connected sequentially according to the principle of similar vertical coordinates and extended to the edge of the image to form a polyline segment, thereby realizing the segmented fitting of the circumferential seam.

2. The method for segmented fitting of shield tunnel annular joints based on gradient value feature map density clustering according to claim 1, characterized in that, The specific process of S1 is as follows: S1.1 Projecting the 3D point cloud data onto the tunnel cross-section From the plane, we obtain the corresponding two-dimensional point set. In the coordinate system of this section, the origin (0, 0) corresponds to the center of the scanner, and the center point of the tunnel cross section is set. The coordinates are (0, By calculating the mean of the vertical coordinates of each point in this two-dimensional point set, the vertical coordinates of the center point of the cross-section can be obtained. ; S1.2 The 3D tunnel point cloud is acquired by a laser scanner mounted on a track trolley. During unfolding, the 3D tunnel point cloud is unfolded along the top axis of the tunnel according to the principle of cylindrical projection and under the constraint of the scanner's viewpoint. The specific steps are as follows: S1.2-1. In the cross-sectional coordinate system, connect the origin. A ray is formed with any point, and its parametric equation is: ； Based on the center coordinates of the cross section and tunnel design radius Establish the equation of the tunnel cross section circle Substituting the ray equation into the cross-sectional circle equation yields the following about... The quadratic equation: ,in , , Solve the quadratic equation and choose The solution calculates the intersection point of the ray and the cross-sectional circle. coordinates ; S1.2-2, Let line segment along The axis extends in the positive direction and intersects the cross-sectional circle at point . Calculation points With point The connecting line is the same as the straight line. The included angle ; Calculate the arc length at any point And calculate the point. Intersection with the cross-section circle Distance between Modify the coordinates of any point in the 3D tunnel point cloud to This allows the 3D tunnel point cloud to unfold along the top of the tunnel into a 2.5D point cloud; S1.

3. Convert the 2.5D point cloud into a grayscale image using the raster method. The specific steps are as follows: S1.3-1 Read the planar coordinates of each point in the unfolded point cloud and the corresponding strength value Constructing a two-dimensional point cloud dataset Based on the set physical size of the pixel unit The unit is meters per pixel, which determines the size of the grayscale image: , ,in , , , Point clouds coordinate, The maximum and minimum values ​​of the coordinates. Indicates the image width. Indicates the image height. Indicates rounding up; S1.3-2, Establish the mapping relationship from point cloud to pixel: column index row index ,in, , , Indicates rounding down; S1.3-3. Based on the mapping relationship from point cloud to pixel, assign the point cloud intensity value to the pixel. If multiple point clouds fall into the same pixel, calculate the average of the intensity values ​​of these multiple point clouds. Assign this pixel, where, Is falling into pixels The number of point clouds, if no point cloud falls into a pixel If the strength value is 0, then the strength value matrix is ​​formed; S1.3-4. Normalize the intensity value matrix to the interval [0, 255]. The calculation formula is as follows: ,in, Represents pixels Normalized grayscale values This indicates rounding to the nearest integer. and These represent the maximum and minimum values ​​of the intensity value matrix, respectively, and the generated size is... An 8-bit grayscale image.

3. The method for segmented fitting of shield tunnel annular joints based on gradient value feature map density clustering according to claim 2, characterized in that, The specific process of generating the gradient value feature map in S2 is as follows: S2.1, Using Gaussian kernels Perform a convolution operation on the generated grayscale image to obtain the filtered image: ； In the formula, For Gaussian kernel function in The weight value at that point, This is the offset within the kernel relative to the center pixel. To customize the standard deviation, Where is the kernel radius, and the kernel size is... ; S2.

2. Use a vertical edge detection operator to calculate the gradient of the filtered image to obtain the vertical gradient features that reflect the gray-level changes of the image. The calculation formula is as follows: ； in, For pixels The vertical gradient value, This is the offset within the kernel relative to the center pixel. For edge detection operators in the vertical direction, the convolution kernel is in The weight value at that point, The radius of the convolution kernel; S2.3, Based on the mean of the absolute values ​​of the gradient values Set a gradient magnitude threshold Where scale is a user-defined parameter; three feature maps are generated based on the sign and magnitude of the gradient values, as follows: (1) Only retain For each pixel, its gradient value is used as the gray value of that pixel, and the gray values ​​of other pixels are all set to 0. The gradient values ​​are then normalized and mapped to the [0,255] interval to generate a positive gradient value feature map A, which is used to identify the lower edge of the circumferential seam. (2) Only retain For each pixel, the absolute value of its gradient is used as the gray value of that pixel, and the gray values ​​of other pixels are all set to 0. The gradient values ​​are normalized and mapped to the interval [0,255] to generate a negative gradient value feature map B, which is used to identify the upper edge of the circumferential seam. (3) Retain For all pixels, the absolute value of their gradient values ​​is taken as the gray value of that pixel, and the gray values ​​of other pixels are set to 0. The gradient values ​​are then normalized and mapped to the interval [0, 255] to generate a comprehensive gradient value feature map C.

4. The method for segmented fitting of shield tunnel annular joints based on gradient value feature map density clustering according to claim 3, characterized in that, The specific process of S3 is as follows: S3.1 Calculate the sum of gray values ​​of pixels in each row of the positive gradient feature map A and the negative gradient feature map B. , ,in Indicates the image at position The grayscale value at that location; S3.2, Based on the design width of the tunnel segment ring Divide the image vertically into A series of continuous equal-height intervals, based on the set physical size of the pixel unit. Determine the pixel height of each interval as Total number of intervals ; S3.3, Find the intervals respectively Largest row As candidate edge rows for this interval, Perform this operation on the positive gradient value feature map A and the negative gradient value feature map B to obtain the candidate lower edge row set and the candidate upper edge row set; S3.4, Calculate the first... Candidate edge rows of each interval absolute difference from other candidate edge rows Count the number of intervals whose differences satisfy the following conditions. : , making ;in, It is the set of positive integers. The pixel tolerance is the absolute difference between candidate edge rows. Defined allowable physical distance deviation; like If the candidate behavior is verified by the edge row in that interval, then the candidate behavior is verified by the edge row. ;like Then within that interval, according to Starting from the largest to the smallest, try other rows as candidates in turn, repeat the above verification process until a row that meets the conditions is found; S3.5 For each interval, check its valid upper edge row. With the bottom edge line Take the arithmetic mean to obtain the row coordinates of the circumferential seam center in this interval. As the initial circumferential seam position.

5. The method for segmented fitting of shield tunnel annular joints based on gradient value feature map density clustering according to claim 4, characterized in that, The specific process of S4 is as follows: S4.1, Using the initial circumferential seam position as the coordinate system. Centered on the design width of the circumferential seam Set a physical distance threshold in the vertical direction. The corresponding pixel distance threshold Only retain image C Within the range of pixels, the grayscale value of all other pixels is set to 0, generating a gradient value feature map D; S4.2, Based on the physical width of the circumferential seam and pixel unit physical size Set the neighborhood radius and the minimum number of samples within the neighborhood radius for the DBSCAN algorithm. The DBSCAN algorithm is used to perform cluster analysis on the non-zero pixels in image D, and a minimum pixel count threshold for each cluster is set. The number of pixels deleted is less than Clusters that generate clusters containing Image F of the class result; S4.3 For each cluster in image F Calculate the boundary of its enclosing rectangular region: ， ； The circumferential seam area is adaptively divided into several rectangular segments along the transverse direction. .

6. The method for segmented fitting of shield tunnel annular joints based on gradient value feature map density clustering according to claim 5, characterized in that, The specific process of S5 is as follows: S5.1, placing those located in the same Rectangular areas within the specified range are grouped into the same vertically similar group. ; S5.2, in each longitudinally similar group Internally, all rectangular regions are arranged in ascending order according to their column coordinates, forming an ordered sequence. ,in , The number of rectangular regions within the group; for each rectangular region Define its adjacent segment set as: ; S5.3, Based on the rectangular area Rectangular regions are extracted from the positive gradient feature map A and the negative gradient feature map B, respectively, and the cumulative grayscale value of each row of pixels in the rectangular region is calculated. Find each region Largest row As candidate edge rows for this region; perform this operation on all rectangular regions in the positive gradient value feature map A and the negative gradient value feature map B respectively to obtain the candidate lower edge row set and the candidate upper edge row set; S5.4, According to the design width of the circumferential seam and pixel unit physical size Get the pixel width of the circumferential seam ; To ensure the spatial continuity and smoothness of the final fitted circumferential seam, edge position constraints are introduced between adjacent segments. For each longitudinally similar group... The candidate edge rows of the rectangular region within the area are optimized using the following steps: (1) For segmentation The difference between its candidate edge row and the candidate edge rows of two adjacent segments satisfies the following condition: If not satisfied, then within that segment, according to... Starting from the largest to the smallest, try other rows as candidate rows one by one until a row that meets the criteria is found, and then update. If none of the rows in a segment meet the condition, then the segment is removed. (2) After the above optimization, if the four adjacent segmented regions within the group Meet the conditions Then determine For the areas that need adjustment; (3) For each detected The regions will be adjusted as follows: exist Within the segmented area, according to Try other rows in descending order of size as candidate rows until a satisfactory row is found. The line, and update If no row meets the criteria, delete it. Segmented regions; like If the segmented region is not deleted, then in Within the segmented area, according to Try other rows in descending order of size as candidate rows until a satisfactory row is found. The line, and update If no row meets the criteria, delete it. Segmented regions; like If the segmented region has been deleted, then in Within the segmented area, according to Try other rows in descending order of size as candidate rows until a satisfactory row is found. The line, and update If no row meets the criteria, delete it. Segmented regions; The above steps yield the final top edge row for each segmented region. and bottom edge row ; S5.5, For each segment Use its final top edge row and bottom edge row Calculate the center point of the circumferential joint. Coordinates: , ; S5.6, in each group Internally, adjacent suture points are connected sequentially to form a broken line segment that extends to the edge of the image, creating a continuous and smooth suture that runs through the entire image. , is represented as: 。

Images