*2.1. Comparison of Fitting Methods*

Common surface fitting methods include the meshing method [26], Poisson surface reconstruction method [27], Lagrangian interpolation method [28], and cubic spline interpolation method [29]. This experiment uses 3DLS The instrument obtains the surface point cloud data of the initial branch of the tunnel K109+870~K109+900 in the tunnel project, and takes the point cloud data of some areas as the analysis object, and compares and analyzes the fitting degree and fitting of the four methods when constructing the surface. Accuracy. In the schematic diagram of the degree of surface fitting by the four fitting methods, the *X* axis represents the direction along the central axis of the tunnel, the *Y* axis represents the horizontal direction of the tunnel cross-section, and the *Z* axis represents the vertical direction of the tunnel cross-section. In the schematic diagram, only the intercept Part of the area on the surface of the tunnel.

(1) The degree of fit of the four methods

After the tunnel surface is constructed by meshing (Figure 5), the surface is complete and has good continuity, but the smoothness of the surface is poor, and the details of the local area are not enough.

The tunnel surface obtained by the Poisson reconstruction method (Figure 6) is continuous and complete, which can reflect the unevenness of the tunnel surface, but the smoothness of the surface is not high, and the construction details of the point cloud cavity are insufficient.

The tunnel surface constructed by the Lagrangian interpolation method (Figure 7) can reflect the overall contour of the tunnel surface and can reflect the basic details of the tunnel surface in a local area. However, the smoothness of the surface constructed by this method is poor, and there are convex hulls. Phenomenon.

The tunnel surface constructed by the cubic B-spline interpolation method (Figure 8) is continuous and complete, with high smoothness, no local mutations, etc., and the details of the local area are rich, and the tunnel surface is better restored.

**Figure 5.** Partial area of tunnel surface constructed by meshing method.

**Figure 6.** Partial area of tunnel surface constructed by Poisson.

(2) Analysis of the fitting accuracy of the four methods

The data this time is a total of 40,000 point clouds. The data points are extracted 5 times from the point cloud data at 5 cm intervals and 100 data points are randomly selected each time, which is divided into 5 groups. After that, the *z* value corresponding to each point is stored in the order of arrangement, and then the surface is constructed using four methods for the 10 cm interval point cloud data. The *x* and *y* of each stored point on the surface correspond to the corresponding *z* value, 5 cm interval points are compressed to ge<sup>t</sup> 10 cm interval point clouds, so the *z* value on the surface obtained by the 10 cm interval point cloud fitting is different from the corresponding *z* values of the points stored in the 5 cm interval point cloud. According to the *x* and *y* of the stored point, the corresponding surface can be obtained and stored, and then the difference of the corresponding point is

calculated. This method is equivalent to the error calculation of the 10 cm interval point cloud, by calculating the sum of each point The fitting difference value corresponding to the fitting surface, the fitting difference value formula is

> *E*

$$=z\_p - z'\_p\tag{1}$$

where: *p* = 1, 2, 3, . . . , *n*.

**Figure 7.** Local area of tunnel surface constructed by Lagrangian interpolation method.

**Figure 8.** The obtained tunnel surface was constructed by cubic B-spline interpolation.

The fitting error generated by the constructed surface can also be called the root mean square error. By counting the root mean square error, you can ge<sup>t</sup> the error of the area where the surface other than the original point is located. The root mean square error formula is:

$$
\sigma = \sqrt{\frac{E\_1^2 + E\_2^2 + \dots + E\_N^2}{n}} = \sqrt{\frac{\sum E\_p^2}{n}} \tag{2}
$$

According to Formula (2), four interpolation methods are used to construct the surface, and the fitting error is calculated based on 5 sets of data, as shown in Table 2.


**Table 2.** Fitting errors of the four sets of data corresponding to the four methods.

Note: Method 1 means meshing method, method 2 means Poisson reconstruction method, method 3 means Lagrangian interpolation method, method 4 means bicubic spline interpolation method.

It can be seen from the tab that the point cloud fitting errors of the Poisson reconstruction method and cubic B-spline interpolation method are kept in a small range, and the fitting accuracy is high. The fitting accuracy of the Lagrangian interpolation method is not high. The grid division method has the largest fitting error and the lowest precision. The results show that the construction of tunnel surface by the Poisson reconstruction method and cubic B-spline interpolation method has higher fitting accuracy.

By comparing and analyzing the fitting degree and fitting accuracy of the four fitting methods when constructing the surface, the following conclusions are drawn: the smoothness of the surface constructed by the meshing method is poor, and the fitting accuracy is low; Poisson reconstruction The fitting accuracy of the method is high, but the surface lacks a certain degree of smoothness; the surface fitting effect constructed by the Lagrangian interpolation method is better, but the fitting accuracy is not high; the cubic B-spline interpolation method has high fitting accuracy, compared with other methods, the surface details are complete and the smoothness is higher. This method has more advantages in the construction of the tunnel surface and has the least influence on the calculation results of the surface flatness of the initial support of the tunnel.

### *2.2. Optimal Fitting of Cubic B-Spline Interpolation*

The main ideas for the optimal fitting of the original point cloud of the tunnel based on the cubic B-spline interpolation method are:

The point cloud data is divided into slices according to the *x*-direction, that is, the tunnel axis direction, which is equivalent to taking a tiny dx as the threshold. The *x* coordinate changes of the point cloud data within this range are considered to be on a 2-dimensional slice point. Then perform spline curve-fitting on this two-dimensional slice. First, convert Cartesian coordinates to polar coordinates. Since the cross-section of the tunnel is a curve similar to an arc, after changing to polar coordinates, it can be ensured that the depression angle and the polar radius can be in a one-to-one correspondence. Then use polar coordinates for interpolation and encryption based on the spline curve, and then convert the polar coordinates back to rectangular coordinates to complete the fitting of each section of the tunnel. For the same reason, perform the above processing again in the *Y* direction. After processing, the curve interpolation is carried out in the two orthogonal directions of the tunnel *x* and *y*, and the interpolation in the two directions is superimposed to form the fitting surface of the first branch surface of the tunnel [30].

The method of two-way slice complementary to the overall surface of the tunnel proposed in this study effectively eliminates the jagged layering effect of the one-way slice on the overall fitting surface of the tunnel, and the enlarged tunnel surface appears smoother. The cloud fitting operation speed is also faster than the overall point cloud fitting surface [31].

### *2.3. S-G Filter Smoothing Based on Curvature*

The curved surface after the fitting process by cubic B-spline interpolation has a high degree of smoothness, good continuity, and more complete and rich local details, which is more consistent with the actual engineering situation. However, the fitted surface still has many point clouds that deviate from the actual surface, and its accuracy cannot meet the requirements of flatness calculation. To make the final fitting surface that can meet the requirements of flatness calculation, the fitting surface should not only be close to the actual situation but also ensure that the fitting surface has sufficient smoothness. To further limit the smoothness and authenticity of the fitted surface, this study guarantees the reliability of the fitted surface by limiting the curvature of each point on the fitted surface. Based on the fitting processing of cubic B-spline interpolation, The fitted surface is again processed by SG filtering based on curvature limitation [32]. In all tunnel projects, the design parameters of the Leicaoshan Tunnel are universal, and among many tunnels, Leicaoshan Tunnel is the most typical. The flatness detection method in this study specifies the upper and lower limits of the tunnel point cloud curvature as the curvature parameter value of the Leicaoshan Tunnel. That is, the upper limit is specified as the maximum curvature of the vault in the design parameters of Leicaoshan Tunnel, 0.395, and the lower limit is specified as the minimum curvature of the arch bottom in the design parameters of Leicaoshan Tunnel, 0.104.

The specific plan for curvature limitation is as follows:

(1) First, the curvature of any point on the curve is calculated by the curve function. The curvature of a point is calculated based on the first and second derivatives of the two points before and after. Since the fitted surface obtained by the B-spline interpolation method is fitted by the slicing method, the curvature in this step is also used in the same way. The two-dimensional curvature of each tangent surface in the previous B-spline interpolation method is calculated. Superimposed to form the curvature of the entire surface, the formula for calculating the curvature is as follows:

$$K = \frac{|y^\*|}{\left(1 + y'^2\right)^{3/2}}\tag{3}$$

In the formula, *K* is the curvature at a point on the curve, and *y* is the function of the corresponding curve.

(2) There are points in the calculated surface point cloud curvature that are not within the limited range. At this time, it is necessary to perform smoothing and filtering again on the surface obtained by the B-spline interpolation method. For the cross-sectional direction and the longitudinal direction, the data in the two directions is smoothed and filtered again. Taking into account the symmetry of the tunnel and the requirements for the calculation of the surface flatness of the initial support of the tunnel, Savitzky–Golay filtering is used here. The processing of S-G filtering maintains the best shape of the original data, making the processed surface closer to the actual engineering situation.

(3) After removing the unqualified points through the above method, the remaining unqualified points are only the points at the bottom of the arch. Because the curvature calculation is calculated using the curvature of the front and back slices, the curvature calculated at the two ends of the point cloud data, namely the two arch bottoms on the left and right, is meaningless in itself, and it can be directly eliminated. At this point, the final flatness calculation datum (Figure 9) is obtained surface. In the figure, the *X* axis direction represents the direction along the central axis of the tunnel, the *Y* axis represents the horizontal direction of the tunnel cross-section, and the *Z* axis represents the vertical direction of the tunnel cross-section. Segment fitting surface.

**Figure 9.** Flatness calculation datum.
