A method for detecting the flatness of tunnel initial support based on three-dimensional point cloud

Abstract

The present invention relates to a method for detecting the flatness of a tunnel primary branch based on a three-dimensional point cloud, comprising the following steps: step S1: collecting and obtaining three-dimensional point cloud data in the tunnel primary branch area, and preprocessing the point cloud data; step S2: obtaining a fitting surface by fitting processing through an improved B cubic spline interpolation method; step S3: limiting the point cloud curvature threshold range on the fitting surface, limiting the curvature of the corresponding regional points on the fitting surface through an S-G filtering smoothing processing method, and obtaining a flatness calculation reference surface; step S4: establishing an overall and local flatness calculation method, and calculating the overall and local flatness of the tunnel primary branch surface; step S5: drawing a flatness distribution map, and comprehensively analyzing the flatness of the tunnel primary branch surface; step S7: establishing a tunnel engineering flatness detection discrimination standard based on three-dimensional laser scanning technology. The present invention improves the surface fitting accuracy of the tunnel primary branch surface, and can also ensure the accuracy of flatness detection.

Description

Technical Field

[0001] The present invention relates to the technical field of tunnel engineering and relates to a method for detecting the flatness of a tunnel primary support surface during a construction process based on three-dimensional point clouds. Background Art

[0002] As 3D technology becomes increasingly sophisticated, the emerging 3D laser scanning measurement is quietly entering the transportation field. 3D laser scanning emits laser signals at a certain frequency. When the signal encounters an obstacle, it can be reflected and fed back to the receiving device of the 3D scanner. The scanner will scan millions of point clouds per second, with 360° full coverage, and conduct a full-scale scan of the object. The data obtained after scanning is presented in the form of points, which is called a 3D point cloud. As the name suggests, the 3D point cloud obtained by scanning is calculated in the order of 100 million. After being processed by data technology (thinning and de-noising), it can fit the tunnel surface on the basis of meeting the requirements of engineering construction and provide guidance for construction.

[0003] The flatness of the initial support of a tunnel refers to the degree of unevenness of the concrete lining on the initial support surface during tunnel excavation, and is one of the important criteria for the quality acceptance of tunnel projects. Its result directly affects the laying of the later waterproofing membrane, the quality of the second-phase lining, and the safety of the tunnel project; excessive flatness of the initial support surface of the tunnel directly leads to a decrease in the construction quality of the laying of the waterproofing membrane and the second-phase lining, which indirectly affects the permeability of the tunnel surface, making it more likely to seep in and leak in the tunnel. Since tunnels are generally large in scale, have a long service life, and are in complex environmental conditions, it is necessary to perform surface fitting on the surface of the initial support section of the tunnel, and then analyze the overall flatness of the initial support surface of the tunnel.

[0004] Traditional methods for detecting the flatness of the initial support of a tunnel mainly include the 2m ruler combined with a wedge-shaped feeler gauge detection method and the total station detection method. The 2m ruler combined with a wedge-shaped feeler gauge detection method is to place the ruler close to the wall to be detected, then move the slider of the wedge-shaped feeler gauge to the top of the ruler, and then insert the top of the wedge-shaped feeler gauge into the gap between the ruler and the wall to be detected to read the slider scale at this time, which is the deviation value of the flatness. According to the error index data specified in the specification, if the detection value is greater than the specified error, the detection point will be marked as an unqualified point, and finally the flatness of the wall will be determined based on the percentage of unqualified points in the total detection points. However, due to the randomness of the selection of detection points, this detection method has low detection accuracy, slow detection speed and low efficiency; when using a total station to detect the flatness of the initial support surface of the tunnel, due to the influence of the prism-free mode on the measurement accuracy, the deviation between the fitted plane and the actual wall surface will be too large, and the accurate flatness of the initial support of the tunnel cannot be obtained. In addition, the selection of measurement points on the tunnel arch wall is subjective, and the situation of the initial support surface of the tunnel cannot be fully and truly reflected. Summary of the invention

[0005] In order to solve the above technical problems, the purpose of the present invention is to provide a tunnel primary support flatness detection method based on three-dimensional point cloud, which improves the surface fitting accuracy of the tunnel primary support surface, makes the final flatness calculation reference surface closer to the actual engineering situation, and also ensures the accuracy of flatness detection.

[0006] In order to achieve the above-mentioned object of the invention, the present invention adopts the following technical solutions:

[0007] A method for detecting the flatness of a tunnel primary support based on a three-dimensional point cloud comprises the following steps:

[0008] Step S1: using a 3D laser scanner to collect 3D point cloud data in the initial branch area of ​​the tunnel, and preprocessing the 3D point cloud data;

[0009] Step S2: The pre-processed three-dimensional point cloud data is fitted by the improved B cubic spline interpolation method to obtain a fitting surface;

[0010] Step S3: limiting the curvature threshold range of the point cloud on the fitting surface, and limiting the curvature of the points in the corresponding area on the fitting surface by the SG filtering smoothing method, so as to obtain the final flatness calculation reference surface suitable for flatness calculation;

[0011] Step S4: Introduce the concept of overall flatness, establish an overall flatness calculation method, and calculate the overall flatness of the tunnel primary support surface;

[0012] Step S5: Introduce the concept of local flatness, establish a local flatness calculation method, and calculate the local flatness of the tunnel primary support surface;

[0013] Step S6: Draw a flatness distribution map based on the settlement results of the overall flatness and local flatness, and comprehensively analyze the flatness of the tunnel primary support surface;

[0014] Step S7: Establishing a tunnel engineering flatness detection and judgment standard based on three-dimensional laser scanning technology.

[0015] As a preferred solution: the steps of the improved B cubic spline interpolation method in step S2 are as follows:

[0016] First, the 3D point cloud data is sliced ​​according to the x direction, that is, the tunnel extension direction. That is, a small dx is taken as the threshold. The x coordinates of the 3D point cloud data within this range are considered to be points on a 2D slice.

[0017] Then, spline curve fitting is performed on this two-dimensional slice. First, the Cartesian coordinates are converted into polar coordinates. Since the cross section of the tunnel is a curve similar to an arc, after converting to polar coordinates, it is ensured that the depression angle and the polar radius can be in a one-to-one correspondence. Then, the polar coordinates are used to perform spline curve interpolation encryption, and then the polar coordinates are converted back to rectangular coordinates, thus completing the fitting of each cross section of the tunnel.

[0018] Follow the above steps again in the Y direction. After completion, curve interpolation is performed in the two orthogonal directions of the tunnel, x and y. After superimposing the slices in the two directions, a fitting surface of the entire tunnel primary branch surface is formed.

[0019] As a preferred solution: the specific process of the improved B cubic spline interpolation method in step S2 is as follows:

[0020] Step S21, B cubic spline interpolation method is a process of obtaining a curve function group by solving a set of three bending moment equations through a smooth curve of a series of shape value points; the definition of the function is: function F(x, y)∈C[a,b], and in each small interval [x j, y j+1 ] is a cubic polynomial, where a = x0 <x1<...<x n =b,a=y0 <y1<...<y n = b is a given node, then F(x) is called the node x0, x1, ...x n The cubic spline function on j Given function value Y j =f(X j ), (j=0,1,…,n), and S(x j )=y j , (j=0,1,…,n), and the function form on this three-dimensional function is F(x,y)=a+b*x+c*y+d*x^2+e*x*y+f*y^2;

[0021] Step S22: The function in the above step S21 is a piecewise defined formula. Given n+1 data points, there are n intervals in total, and the cubic spline equation satisfies the following conditions:

[0022] a. In each segment interval (i = 0, 1, ..., n-1, x increases), it is a cubic polynomial;

[0023] b. Satisfy (i=0,1,…,n);

[0024] c. Derivative, the second-order derivative is continuous in the interval [a, b], that is, the curve is smooth;

[0025] Therefore, n cubic polynomials can be written as: i = 0, 1, ..., n-1 where a i ,b i ,c i ,d i represents 4n unknown coefficients;

[0026] Step S23, boundary conditions: From the range of i, we know that there are n-1 formulas, but n+1 unknown quantities m; to solve the system of equations, we need two more formulas, so we need to solve the two endpoints x0 and x n and y0,y n The differential is restricted.

[0027] As a preferred solution: in step S3, the curvature threshold range of the point cloud on the fitting surface is limited, and the specific limiting steps are as follows:

[0028] Step S31, calculate the curvature of any point on the curve by using the curvature function. The curvature of a point is calculated based on the first-order and second-order derivatives of the two points before and after it. Since the surface in step S2 is fitted by the slicing method, the curvature in this step is also fitted in the same way. The two-dimensional curvature of each slice in step S2 is calculated, and the curvature of the entire surface is formed by superimposing them. The formula for calculating the curvature is as follows:

[0029]

[0030] In the formula, K is the curvature at a point on the curve, and y is the function of the corresponding curve;

[0031] Step S32: setting the upper limit of the curvature threshold of the tunnel point cloud to 0.395 and the lower limit to 0.104.

[0032] As a preferred solution: the SG filtering and smoothing processing method in step S3 is as follows:

[0033] Suppose we process a column of data x[n], select a set of 2M+1 consecutive data points centered at n=0, and use the following formula to fit the data:

[0034]

[0035] Therefore, the residual of the least squares fit is:

[0036]

[0037] When ε is the smallest, the partial derivatives of ε with respect to each parameter are all 0. At this time, an auxiliary matrix A with (2M+1) rows and (N+1) columns is introduced:

[0038] A={c n,i}, cn,i =n i , -M≤n≤M, 0≤i≤N

[0039] Introduce an auxiliary matrix B:

[0040] B=A T A

[0041] set up:

[0042] We can get: Ba=A T Aa=A T x

[0043] Then: a=(A T A) -1 A T x

[0044] Thus, we can get the coefficients of the polynomial a0, a1, ..., a N The value of , then, by translating with a window width of 2M+1, the filtered sequence y[n] can be obtained;

[0045] From the above derivation, it can be seen that the SG filtering performance is mainly determined by two filtering parameters: window width M and fitting order N. Therefore, if you want to process different data situations, you only need to make reasonable adjustments to these two filtering parameters.

[0046] As a preferred solution: the step S4 further includes:

[0047] Step S41: The normal line derived from the original point cloud intersects with the flatness calculation reference plane to form a normal vector distance;

[0048] Step S42: deriving the relationship between the standard deviation and the flatness; the overall flatness of the tunnel primary surface is determined by the degree of dispersion of the normal vector distance from the original point cloud to the fitted surface. A large dispersion degree indicates that the normal vector distance between the point cloud on the tunnel surface and the flatness calculation reference surface is greatly different, and the surface is rougher; a small dispersion degree indicates that the normal vector distance between the point cloud on the tunnel surface and the flatness calculation reference surface is less different, and the surface is flatter; the dispersion degree is calculated by the standard deviation formula;

[0049] Step S43: Substitute the normal vector distance and the number of point clouds as calculation parameters into the standard deviation formula to calculate and obtain the overall flatness.

[0050] As a preferred solution: the calculation formula for the overall flatness of the tunnel primary support surface is as follows:

[0051]

[0052] Where m0 is the overall flatness of the tunnel primary support surface, di is the normal vector distance, and n is the number of point clouds.

[0053] As a preferred solution: the step S5 further includes:

[0054] Step S51: drawing a flatness distribution map through normal vector distance;

[0055] Step S52: Calculate the gradient of each point through the three-dimensional coordinates of each point on the distribution map;

[0056] Step S53: Introduce the step distance x parameter in the gradient direction of each point, and obtain the point coordinates after the step distance x in the gradient direction of each point according to the Lagrange interpolation method, and calculate the height difference h between the front and rear points;

[0057] Step S54: Draw a schematic diagram for calculating local flatness. The step distance x and the height difference h between the front and rear points form a right triangle. The local flatness is defined as half of the step distance x multiplied by the sine of the elevation angle between the front and rear points.

[0058] As a preferred solution: the formula for defining the local flatness of the initial support surface of the tunnel engineering in step S5 is as follows:

[0059]

[0060] Where m1 is the local flatness of the tunnel primary support surface, x represents the stepping distance of the original point cloud along the gradient direction, h represents the height difference between the starting point and the end point in the stepping direction, and α represents the inclination angle between the two points.

[0061] As a preferred solution: the judgment standard in step S7 is: when the overall flatness calculation result of the tunnel primary support surface flatness is less than 30mm, the tunnel primary support surface flatness meets the requirements and meets the construction requirements. When the overall flatness calculation result of the tunnel primary support surface flatness is greater than 30mm, the tunnel primary support surface flatness does not meet the requirements, and it is necessary to calculate the local flatness of the tunnel primary support surface, and screen out the uneven area through the flatness distribution map, so that the construction team can perform later processing on the tunnel primary support section.

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

[0063] The present invention applies the emerging three-dimensional laser scanning technology to the detection of the flatness of the initial support of the tunnel. The three-dimensional laser scanner can quickly and accurately obtain a large amount of three-dimensional point cloud data on the initial support surface of the tunnel, and obtain the detailed situation of the surface of the initial support section of the tunnel. It is easy to operate and does not need to contact the surface of the initial support of the tunnel to be detected. The pre-processed three-dimensional point cloud data is processed by the improved B-spline interpolation method and the SG filter smoothing based on curvature limitation, so that the initial support surface of the tunnel can be accurately fitted, which meets the flatness calculation requirements. Finally, according to the normal vector distance between the original point cloud and the flatness calculation reference surface, a flatness calculation method is established. The calculation of flatness introduces two new concepts: overall flatness and local flatness. The calculation of overall flatness can directly reflect the overall flatness of the initial support surface of the tunnel, but it cannot reflect the uneven area of ​​the initial support surface of the tunnel and the flatness of the details on the initial support surface of the tunnel. The calculation of local flatness can be accurate to the flatness of every detail on the initial support surface of the tunnel. Overall flatness and local flatness complement each other and are indispensable. Both are very important for the flatness detection of the tunnel primary support surface. After completing the calculation of the overall flatness and local flatness, the local flatness size is colored into a map, which is presented in the form of a flatness distribution map, and a tunnel primary support surface flatness detection and discrimination standard based on three-dimensional laser scanning technology is established based on existing results. So far, the entire tunnel primary support surface flatness detection method system based on three-dimensional laser scanning technology has been fully established. This tunnel engineering primary support flatness detection system can not only determine the overall flatness of the tunnel primary support section, but also intuitively reflect the uneven areas of the tunnel engineering primary support section, providing theoretical methods and technical support for the detection of the flatness of the tunnel primary support surface during construction, and to a great extent improve the flatness detection system in tunnel engineering, so that the flatness detection system in tunnel engineering can meet the guidance expectations of future construction.

[0064] The present invention introduces the improved B cubic spline interpolation method and the curvature-limited SG filter smoothing method into the surface fitting of the tunnel engineering point cloud, and establishes the concepts of overall flatness and local flatness of the tunnel based on the final fitting reference surface, obtains the flatness distribution map through the calculation of flatness, and uses the flatness distribution map to analyze the flatness of the initial support of the tunnel, thereby improving the detection system of the flatness of the initial support section of the tunnel engineering. The present invention solves the problems that the traditional detection method of tunnel flatness cannot cover the entire tunnel, the detection precision is low, the detection accuracy is low, and the tunnel flatness detection system is imperfect, improves the surface fitting accuracy of the initial support surface of the tunnel, makes the final flatness calculation reference surface closer to the actual engineering situation, and can also ensure the accuracy of flatness detection, and also provides theoretical guidance and technical support for the construction process in the actual tunnel engineering. At the same time, the application of the emerging technology of three-dimensional laser scanning technology in the detection of the flatness of the initial support of the tunnel can also meet the guidance expectations of the future construction of the tunnel engineering, and also conforms to the future development trend of the flatness detection of the tunnel engineering. The main innovations of the present invention include the following aspects:

[0065] (1) When performing surface fitting, the present invention selects different point cloud processing methods according to different regions to perform surface fitting, which has the advantages of strong adaptability, strong data integration, and high fault tolerance, and is suitable for a variety of tunnel surface fitting. It can comprehensively reflect the surface flatness and construction conditions of the initial support of the tunnel, and can also observe the conditions of the local area, which can improve the efficiency and effective range of the flatness detection work.

[0066] (2) The present invention is a fitting method based on the improved B cubic spline interpolation method combined with curvature-limited SG smoothing. The conventional B cubic spline interpolation method usually uses a one-way slicing method to perform curve interpolation encryption, and the final fitting surface will have a jagged layering effect. In order to eliminate the error effect of one-way slicing on the experiment, the present invention proposes a method of superimposing two-way slicing on each other to optimize the fitting of the initial branch section of the tunnel. This method greatly improves the fitting accuracy of the fitting surface and the degree of authenticity with the actual project, ensures the smoothness of the tunnel fitting surface and the continuity of the tunnel fitting surface, and ensures the accuracy of the flatness detection.

[0067] (3) The present invention organically combines the improved B cubic spline interpolation method and the SG filter smoothing method based on curvature limitation. According to the curvature of the three-dimensional point cloud of the tunnel primary branch surface, the curvature of the fitting surface is limited by the SG filter smoothing method to ensure the smoothness and neatness of the final fitting surface. Ensure that the final fitting surface obtained meets the calculation requirements of the tunnel surface flatness. At the same time, according to the curvature size of the tunnel surface point cloud data, the fitting area that does not meet the requirements can be effectively captured. The fitting method of the present invention has high optimization degree and high fitting accuracy, and is highly targeted to the optimization fitting processing of the point cloud. The obtained fitting surface is also more suitable for the tunnel flatness detection system.

[0068] (4) The present invention proposes the concept of the normal vector distance from the original point cloud to the flatness calculation reference surface, and judges the flatness of the initial support surface of the tunnel by the fluctuation of the normal vector distance. At the same time, a calculation method for overall flatness and local flatness is established based on the normal vector distance. The principle of the flatness calculation method based on three-dimensional laser scanning technology in the present invention is simple and easy to understand, well-reasoned, and can comprehensively reflect the flatness of the tunnel, ensuring the accuracy of the tunnel flatness detection, and also providing a theoretical basis and technical support for flatness detection in tunnel engineering.

[0069] (5) The present invention obtains the overall flatness and local flatness of the tunnel surface through a flatness calculation method based on three-dimensional laser scanning technology, and draws the flatness of the tunnel surface into a distribution map, and uses the flatness distribution map to observe the initial flatness of the tunnel project in more detail. The flatness distribution map of the present invention can not only reflect the overall flatness of the tunnel surface, but also reflect the details of the uneven areas of the tunnel surface, improve the entire flatness system, and has a good guiding role in the detection of flatness during the construction of the tunnel project.

[0070] (6) The present invention establishes a new standard for flatness detection and judgment of tunnel engineering based on three-dimensional laser scanning technology, improves the flatness detection system of tunnel engineering, and can also meet the guidance expectations for future construction of tunnel engineering. BRIEF DESCRIPTION OF THE DRAWINGS

[0071] The drawings in the specification, which constitute a part of the present application, are used to provide a further understanding of the present application. The illustrative embodiments of the present application and their descriptions are used to explain the present application and do not constitute a limitation on the present application.

[0072] Figure 1 It is a schematic diagram of the process of the present invention;

[0073] Figure 2 It is a three-dimensional distribution diagram of the flatness of the fitted surface after smoothing by the method of the present invention;

[0074] Figure 3 It is a gradient schematic diagram of the method of the present invention;

[0075] Figure 4 It is a schematic diagram of the local flatness definition based on three-dimensional laser scanning technology in the method of the present invention. DETAILED DESCRIPTION

[0076] It should be noted that the following detailed descriptions are illustrative and are intended to provide further explanation of the present application. Unless otherwise specified, all technical and scientific terms used herein have the same meanings as those commonly understood by those skilled in the art to which the present application belongs.

[0077] The present invention will be further described below in conjunction with the accompanying drawings and embodiments:

[0078] like Figure 1 As shown, the present invention is a method for detecting the flatness of the surface of the tunnel primary branch based on three-dimensional point cloud data, comprising the following steps:

[0079] (1) Data collection and processing

[0080] (1a) A 3D laser scanner is used to scan the tunnel. The spatial system of the 3D scanner is based on the existing spatial coordinate system. The existing coordinate system is introduced into the scanner system before laser scanning. Therefore, the tunnel 3D point cloud data collection selects appropriate spatial coordinate points for coordinate matching. Then the scanner instrument is positioned after the total station passes through the rear intersection control point, and the scanning is finally spliced ​​through software.

[0081] (1b) After obtaining the 3D point cloud data obtained by scanning, the corresponding point cloud thinning and denoising are performed according to the hardware computing power. Since two adjacent points may overlap, thinning and denoising are required. In this process, Amberg Tunnel software is mainly used to process these outliers. In the actual acquisition process, due to errors and various reasons, a small number of spatial points that are separated from the main point group will be generated. These outliers will affect the curvature of the subsequent local fitting of the point cloud.

[0082] (2) Division of the fitting calculation area: According to the design and actual situation of the tunnel, the three-dimensional point cloud data is collected, and the center point a of the tunnel bottom surface is taken as the coordinate origin. On the tunnel cross section, 180° of the tunnel is taken as a whole, and 2m is taken in the longitudinal direction according to the calculation capacity, and this area segment is used as the calculation basis;

[0083] (3) In this region, the B cubic spline interpolation method is used to perform surface fitting on the preprocessed three-dimensional point cloud data to obtain the fitted surface after the first processing, which can be compared with the three-dimensional point cloud data.

[0084] The principle of cubic spline interpolation method is introduced as follows:

[0085] 3a) B cubic spline interpolation is a process of obtaining a curve function group by solving the three bending moment equations through a smooth curve of a series of shape value points. The definition of this function is: function F(x, y)∈C[a,b], and it is a cubic polynomial on each small interval [xj, yj+1], where a=x0 <x1<...<x n =b,a=y0 <y1<...<y n = b is a given node, then F(x) is called the node x0, x1, ...x n The cubic spline function on . If at the node x j Given function value Y j =f(x j ), (j=0,1,…,n), and S(x j )=y j , (j=0,1,…,n), and the function form on this three-dimensional function is F(x,y)=a+b*x+c*y+d*x^2+e*x*y+f*y^2.

[0086] 3b) The function is a piecewise defined formula. Given n+1 data points, there are n intervals, and the cubic spline equation satisfies the following conditions:

[0087] a. In each segmented interval (i=0,1,…,n-1, x increases), it is a cubic polynomial.

[0088] b.Satisfy (i=0,1,…,n)

[0089] c. Derivative, the second-order derivative is continuous in the interval [a,b], that is, the curve is smooth.

[0090] So n cubic polynomials can be written piecewise as: i = 0, 1, ..., n-1 where a i ,b i ,c i ,d i Represents 4n unknown coefficients.

[0091] 3c) Boundary conditions

[0092] From the range of i, we can see that there are n-1 formulas, but n+1 unknowns m. To solve this system of equations, we need two more equations. So we need to calculate the two endpoints x0 and x n and y0,y n Add some restrictions to the differential.

[0093] The main idea of ​​fitting the tunnel point cloud based on the B cubic spline interpolation method is to slice and divide the three-dimensional point cloud data according to the x direction, that is, the extension direction of the tunnel, that is, take a small dx as the threshold, and the x coordinate changes of the three-dimensional point cloud data within this range are considered to be points on a two-dimensional slice. Then perform spline curve fitting on this two-dimensional slice. First, convert the Cartesian coordinates to polar coordinates. Since the cross-section of the tunnel is a curve similar to an arc, after converting to polar coordinates, it can be guaranteed that the depression angle and the polar radius can be in a one-to-one correspondence. Then use polar coordinates to perform spline curve-based interpolation encryption, and then convert the polar coordinates back to rectangular coordinates, and the fitting of each section of the tunnel is completed. The same reason is repeated in the Y direction, that is, the cutting direction this time is along the x direction instead of perpendicular to the x direction. After this is completed, curve interpolation is performed in the two orthogonal directions of the tunnel, x and y. After superimposing the slices in the two directions, a fitting surface of the entire tunnel primary support surface is formed.

[0094] (4) In order to ensure the smoothness of the initial branch fitting surface of the tunnel after fitting and its closeness to the actual engineering, it is also necessary to meet this requirement by limiting the curvature of all points on the fitting surface. The specific restriction steps are as follows:

[0095] 4a) First, calculate the curvature of any point on the curve using the curvature function. The curvature of a point is calculated based on the first-order and second-order derivatives of the two points before and after it. Since the surface in step (3) is fitted using the slicing method, the curvature in this step is also calculated in the same way. The two-dimensional curvature of each slice in step (3) is calculated, and the curvature of the entire surface is formed by superimposing them. The formula for calculating the curvature is as follows:

[0096]

[0097] Where K is the curvature of a point on the curve, and y is the function of the corresponding curve.

[0098] 4b) In order to make the final fitting surface usable for flatness calculation, the final fitting surface must not only be very close to the actual situation, but also ensure that it is smooth enough. In order to limit its smoothness and authenticity, the curvature of each point is limited to a certain range. Among all tunnel projects, the design parameters of the Leicaoshan Tunnel are universal, and among many tunnels, the Leicaoshan Tunnel is the most typical. Therefore, the curvature threshold of the tunnel point cloud is defined as the parameters of the Leicaoshan Tunnel. That is, the upper limit is defined as the maximum curvature of the arch in the design parameters of the Leicaoshan Tunnel, 0.395, and the lower limit is defined as the minimum curvature of the arch bottom in the design parameters of the Leicaoshan Tunnel, 0.104.

[0099] There are some points in the curvature of the tunnel point cloud calculated by the above method that are not within the specified range. At this time, it is necessary to perform another smoothing filter on the fitting surface obtained by the B cubic spline interpolation method in step (3). Here, Savitzky-Golay filtering is used.

[0100] The G filtering principle is introduced as follows:

[0101] Suppose we process a column of data x[n], select a set of 2M+1 consecutive data points centered at n=0, and use the following formula to fit the data:

[0102]

[0103] Therefore, the residual of the least squares fit is:

[0104]

[0105] When ε is the smallest, the partial derivatives of ε with respect to each parameter are all 0. At this time, an auxiliary matrix A with (2M+1) rows and (N+1) columns is introduced:

[0106] A={c n,i}, c n,i =n i , -M≤n≤M, 0≤i≤N

[0107] Introduce an auxiliary matrix B:

[0108] B=A T A

[0109] set up:

[0110] We can get: Ba=A T Aa=A T x

[0111] Then: a=(A T A) -1 A T x

[0112] Thus, we can get the coefficients of the polynomial a0, a1, ..., a N Then, by shifting with a window width of 2M+1, the filtered sequence y[n] can be obtained.

[0113] From the above derivation, it can be seen that the SG filtering performance is mainly determined by two filtering parameters: window width M and fitting order N. Therefore, if you want to process different data situations, you only need to make reasonable adjustments to these two filtering parameters.

[0114] This filtering method is based on the local polynomial least squares fitting. The biggest feature of this filter is that it can ensure that the shape and width of the signal remain unchanged while filtering out noise, that is, it keeps the shape of the original data as best as possible, making the processed surface closer to the actual engineering situation. After eliminating the points that do not meet the conditions through the above method, the remaining points that do not meet the conditions are only the points at the two arch bottom positions. Because the curvature calculation is to use the curvature of the front and back slices, the curvature calculated from the two sections of the three-dimensional point cloud data, i.e. the left and right arch bottoms, is meaningless in itself and can be directly eliminated. So far, the final flatness calculation reference surface is obtained.

[0115] (6) Through the above method, a flatness calculation reference surface suitable for calculating the flatness of the tunnel primary support surface has been obtained. The flatness calculation method in the present invention is based on the normal vector distance from the original point cloud of the tunnel primary support surface to the fitting surface. The calculation of the normal vector distance only requires drawing a normal line perpendicular to the fitting surface through the original point cloud. The normal line intersects the fitting surface at one point, and the distance between this point and the starting point of the ray is the normal vector distance di.

[0116] (7) The overall flatness of the tunnel primary surface is mainly determined by the degree of dispersion of the normal vector distance from the original point cloud to the fitted surface. If the degree of dispersion is large, it means that the normal vector distance between the point cloud on the tunnel surface and the flatness calculation reference surface is large, and the surface is rougher. If the degree of dispersion is small, it means that the normal vector distance between the point cloud on the tunnel surface and the flatness calculation reference surface is small, and the surface is flatter.

[0117] In order to intuitively reflect the overall flatness through the calculation formula and obtain the overall flatness result, the present invention introduces the concept of standard deviation. Simply put, standard deviation is a measurement concept of the degree to which a group of values ​​are dispersed from the average value. A larger standard deviation means that the difference between most values ​​and their average value is larger; a smaller standard deviation means that these values ​​are closer to the average value.

[0118] Substituting the normal vector distance di and the number of point clouds n into the standard deviation formula, the overall flatness calculation formula of the tunnel primary support surface can be obtained as follows:

[0119]

[0120] Where m0 is the overall flatness of the tunnel primary support surface, d i is the normal vector distance, and n is the number of point clouds.

[0121] (8) Through MATLAB programming, we first obtain the flatness distribution map, which is the normal vector distance distribution map. Figure 2 shown.

[0122] 8a) The flatness data in the figure constitutes a three-dimensional scalar field. First, the gradient value of this scalar data point at each point is calculated. However, the size of this gradient has no practical significance for the calculation of local flatness here, because only the direction of the gradient needs to be known. Theoretically, every point on a three-dimensional surface has countless directions that can change. The gradient of each point is calculated only to obtain the direction with the largest change in flatness.

[0123] 8b) Since it is the gradient of a binary function, the values ​​obtained correspond to the components in the x and y directions, which is the following formula:

[0124]

[0125] In this way, the direction of the gradient vector can be determined based on the x component and the y component, and the direction angle of the gradient can also be calculated.

[0126] 8c) The flatness distribution map obtained based on the normal vector distance is similar to a topographic map. Here, the principle of slope in the topographic map is borrowed. Slope = height difference / horizontal distance. In the gradient direction of each point, the slope i and the inclination angle α are calculated. The horizontal distance X can be adjusted according to the accuracy requirements.

[0127] 8d) Once the horizontal distance is determined, the height difference between the two points can be known. It can be imagined as a circle with a certain point as the center and the horizontal distance as the radius. The gradient direction is unique. In the gradient direction, the height changes fastest, and the height difference is the largest when the horizontal distance is constant.

[0128] Note: What needs to be noted here is how to find the corresponding horizontal distance according to the direction of the gradient. Figure 3 As shown in the figure: Assume that the red arrow is the calculated gradient vector direction angle, and the black dot is the normal vector distance data of a certain point, that is, the normal vector distance data is in the direction indicated by the red arrow, which is the direction of the fastest change (gradient direction). The original coordinates of the point can be marked in the figure, and the length X along the gradient direction can be represented as the green point in the figure.

[0129] Therefore, the coordinates of the green point, which is the data point that meets the requirements, can be calculated through the X value and the direction angle. Because the coordinate values ​​of all the original point cloud data on the plane are known, the normal vector distance corresponding to the green point coordinates can be calculated by interpolation. At the same time, the height difference h between the two points (black and green points) can also be calculated.

[0130] 8e) In order to make the definition of local flatness meet the measurement requirements and comply with the technical specifications for tunnel engineering construction, the definition formula for the flatness of the initial support surface of the tunnel engineering based on 3D laser scanning technology is as follows:

[0131]

[0132] In the formula, m1 is the local flatness of the tunnel primary surface, x represents the step distance of the original point cloud along the gradient direction, h represents the height difference between the starting point and the end point in the step direction, and α represents the inclination angle between the two points. The schematic diagram of the definition of local flatness is shown in Figure 4 shown.

[0133] (9) The overall flatness calculation results and the local flatness calculation results are plotted into flatness distribution maps and analyzed comprehensively with the respective flatness calculation results. Based on the existing results, a standard for judging the flatness of the tunnel primary support surface based on three-dimensional laser scanning technology is established. The judgment standard is: when the overall flatness calculation result of the tunnel primary support surface flatness is less than 30mm, the flatness of the tunnel primary support surface meets the requirements and meets the construction requirements. When the overall flatness calculation result of the tunnel primary support surface flatness is greater than 30mm, the flatness of the tunnel primary support surface does not meet the requirements, and it is necessary to calculate the local flatness of the tunnel primary support surface, and screen out the uneven areas through the flatness distribution map, so that the construction team can perform later processing on the tunnel primary support section.

[0134] The detection of the flatness of the initial support of the tunnel is a key step in the entire construction process of the tunnel project. The method of the present invention solves the problem that the current traditional detection method of the flatness of the initial support of the tunnel cannot meet the expectations of future construction guidance and the situation of the existing tunnel project flatness detection system is not optimistic.

[0135] In the description of this specification, the description with reference to the terms "one embodiment", "some embodiments", "examples", "specific examples", or "some examples" means that the specific features, structures, materials or characteristics described in conjunction with the embodiment or example are included in at least one embodiment or example of the present invention. In this specification, the schematic representation of the above terms does not necessarily refer to the same embodiment or example. Moreover, the specific features, structures, materials or characteristics described may be combined in any one or more embodiments or examples in a suitable manner.

[0136] Although the embodiments of the present invention have been shown and described above, it can be understood that the above embodiments are exemplary and cannot be understood as limiting the present invention. Ordinary technicians in the field can change, modify, replace and modify the above embodiments within the scope of the present invention without departing from the principles and purpose of the present invention. Any simple modifications, equivalent changes and modifications made to the above embodiments based on the technical essence of the present invention are still within the scope of the technical solution of the present invention.

Claims

1. A method for detecting the flatness of a tunnel primary support based on three-dimensional point cloud, characterized in that: The steps include: Step S1: using a 3D laser scanner to collect 3D point cloud data in the initial branch area of ​​the tunnel, and preprocessing the 3D point cloud data; Step S2: The pre-processed three-dimensional point cloud data is fitted by the improved B cubic spline interpolation method to obtain a fitting surface; Step S3: limiting the curvature threshold range of the point cloud on the fitting surface, and limiting the curvature of the points in the corresponding area on the fitting surface by the SG filtering smoothing method, so as to obtain the final flatness calculation reference surface suitable for flatness calculation; Step S4: Introduce the concept of overall flatness, establish an overall flatness calculation method, and calculate the overall flatness of the tunnel primary support surface; Step S5: Introduce the concept of local flatness, establish a local flatness calculation method, and calculate the local flatness of the tunnel primary support surface; Step S6: Draw a flatness distribution map based on the settlement results of the overall flatness and local flatness, and comprehensively analyze the flatness of the tunnel primary support surface; Step S7: Establishing a tunnel engineering flatness detection and discrimination standard based on three-dimensional laser scanning technology; The steps of the improved B cubic spline interpolation method in step S2 are as follows: First, the 3D point cloud data is sliced ​​according to the x direction, that is, the tunnel extension direction. That is, a threshold is taken. The x coordinate changes of the 3D point cloud data within this range are considered to be points on a 2D slice. Then, spline curve fitting is performed on this two-dimensional slice. First, the Cartesian coordinates are converted into polar coordinates. Since the cross section of the tunnel is an arc-shaped curve, after converting to polar coordinates, it is ensured that the depression angle and the polar radius are in a one-to-one correspondence. Then, the polar coordinates are used to perform spline curve interpolation encryption, and then the polar coordinates are converted back to rectangular coordinates to complete the fitting of each cross section of the tunnel. Follow the above steps again in the Y direction. After completion, curve interpolation is performed in the two orthogonal directions of the tunnel, x and y. After superimposing the slices in the two directions, a fitting surface of the entire tunnel primary branch surface is formed.

2. The method for detecting the flatness of a tunnel primary support based on three-dimensional point cloud according to claim 1, characterized in that: In step S3, the curvature threshold range of the point cloud on the fitting surface is limited, and the specific limiting steps are as follows: Step S31, calculate the curvature of any point on the curve by using the curvature function. The curvature of a point is calculated based on the first-order and second-order derivatives of the two points before and after it. Since the surface in step S2 is fitted by the slicing method, the curvature in this step is also fitted in the same way. The two-dimensional curvature of each slice in step S2 is calculated, and the curvature of the entire surface is formed by superimposing them. The formula for calculating the curvature is as follows: In the formula, K is the curvature at a point on the curve, and y is the function of the corresponding curve; Step S32: setting the upper limit of the curvature threshold of the tunnel point cloud to 0.395 and the lower limit to 0.

104.

3. The method for detecting the flatness of a tunnel primary support based on three-dimensional point cloud according to claim 1, characterized in that: The step S4 further comprises: Step S41: The normal line derived from the original point cloud intersects with the flatness calculation reference plane to form a normal vector distance; Step S42: deriving the relationship between the standard deviation and the flatness; the overall flatness of the tunnel primary surface is determined by the degree of dispersion of the normal vector distance from the original point cloud to the fitted surface. A large dispersion degree indicates that the normal vector distance between the point cloud on the tunnel surface and the flatness calculation reference surface is greatly different, and the surface is rougher; a small dispersion degree indicates that the normal vector distance between the point cloud on the tunnel surface and the flatness calculation reference surface is less different, and the surface is flatter; the dispersion degree is calculated by the standard deviation formula; Step S43: Substitute the normal vector distance and the number of point clouds as calculation parameters into the standard deviation formula to calculate and obtain the overall flatness.

4. The method for detecting the flatness of a tunnel primary support based on three-dimensional point cloud according to claim 3, characterized in that: The calculation formula for the overall flatness of the tunnel primary support surface is as follows: Where m0 is the overall flatness of the tunnel primary support surface, di is the normal vector distance, and n is the number of point clouds.

5. The method for detecting the flatness of a tunnel primary support based on three-dimensional point cloud according to claim 1, characterized in that: The step S5 further comprises: Step S51: drawing a flatness distribution map through normal vector distance; Step S52: Calculate the gradient of each point through the three-dimensional coordinates of each point on the flatness distribution map; Step S53: Introduce the parameter of step distance X in the gradient direction of each point, obtain the coordinates of the point after stepping distance X in the gradient direction of each point according to the Lagrange interpolation method, and calculate the height difference h between the front and rear points; Step S54: Draw a local flatness calculation diagram. The step distance X and the height difference h between the front and rear points form a right triangle. The local flatness is defined as half of the step distance X multiplied by the sine of the elevation angle between the front and rear points.

6. The method for detecting the flatness of the initial support of a tunnel based on three-dimensional point cloud according to claim 5, characterized in that: The calculation formula for the local flatness of the initial support surface of the tunnel engineering in step S5 is as follows: Where m1 is the local flatness of the tunnel primary support surface, x represents the stepping distance of the original point cloud along the gradient direction, h represents the height difference between the starting point and the end point in the stepping direction, and α represents the inclination angle between the two points.

Images