# Research on Fitting Algorithm of Scattered Point Cloud Based on 3D Laser Scanning

Authors: Tianen Jiang; Xiaojun Lin; Liang Fang; Tianyu Yao; Yuzhen Zhao; Zheng Wei
Journal: Mobile Information Systems
 (2022)

---

## Abstract

3D laser scanning equipment is an important equipment for reverse engineering. It is used as an information sensing device for the Industrial IoT. After analyzing and processing the measured data, it can issue construction instructions to equipment such as TBMs and shotcrete manipulator. It can be used to solve the problems of initial support contour detection, flatness detection, monitoring measurement, and so on in tunnel engineering. This paper considers that 3D laser scanning point cloud is a scattered point cloud and cannot be used directly. In order to extract the point cloud information as much as possible, it is necessary to fit the point cloud obtained by scanning first. This paper uses the triangulation method, the Poisson reconstruction method, and the Lagrange method, etc. to fit the scattered point cloud, analyzes the advantages and disadvantages of several fitting methods, and analyzes the accuracy of the surfaces generated by several fitting methods. This creates conditions for the subsequent utilization of the tunnel 3D laser scanning point cloud.


---

## Body

## 1. Introduction

At present, equipment such as TBMs and shotcrete manipulator has been widely used in tunnel construction. In tunnel construction, there is a structure first, and then, loads and the structure need to be inspected and monitored. In the construction of the New Austrian Tunnel, the initial support and the surrounding rock share most of the load. The correct tunnel profile can reduce stress concentrations. Better flatness can prevent the waterproof membrane from being damaged and the secondary lining from being emptied. The deformation monitoring of the initial support can guide the on-site construction [1]. The tunnel engineer issues construction decisions based on the inspection and monitoring results. In order to make accurate decisions, it is necessary to conduct comprehensive and rapid inspection and monitoring of tunnels, and to analyze the results. 3D laser scanning equipment is an important Industrial IoT information sensing equipment [2–4]. Comprehensive and rapid detection and monitoring can be achieved. In the construction of the ENN tunnel, the detected point cloud contains the geometric information of the tunnel, and this technology is gradually adopted in tunnel contour detection, flatness detection, and monitoring measurement, etc. 3D laser scanning generated in the construction of the New Austrian Tunnel is a scatter cloud, which cannot be used directly. In order to solve the point cloud utilization problem and improve the calculation accuracy, it is necessary to extract as much point cloud information as possible, and to fit the point cloud [5].The author considers fitting the scattered point cloud formed by 3D laser scanning first, extracting more point cloud information, and eliminating or diluting the influence of noise on subsequent work. Therefore, point cloud surface fitting is a means to solve such problems. At the same time, using a point cloud to fit the surface can also extract surface information from the point cloud as much as possible. Realize the comparison of surface and surface and overcome the problem that it is difficult to directly compare scattered point clouds.The fitting method in this paper considers the fitting methods of the parameter reconstruction method, the triangular mesh method, and the implicit function method, respectively, and selects the commonly used Lagrange interpolation method, triangulation method, and the Poisson construction method. According to the fitting situation of various methods, the applicable scope of tunnel detection and monitoring is analyzed.


## 2. Related Works

The Amberg Tunnel software developed by Amberg Technologies uses the direct comparison method between point cloud and design surface to realize over-under-digging detection [6, 7], but this method is difficult to obtain high detection accuracy in the case of point cloud noise [8–10]. Ho et al. [11], Yao et al. [12] and others projected on the cross section through point cloud data [13–17], and compared the projected contour with the design contour, and obtained the over-undercut amount and over-undercut area [18]. These methods usually have a larger error than the point cloud spacing because 3D laser scanning point cloud is a scattered point cloud [19].Tang et al. used the heat map of the Amberg Tunnel software to analyze the smoothness of the tunnel [20]. Based on 3D laser scanning point cloud surface, China Railway Group Limited uses the pseudo-scale measurement method to realize the flatness detection in the circumferential and longitudinal directions of the tunnel, but this method is not accurate enough when the grid division is large, and when the grid division is small, the amount of calculation is too large [21, 22].The outline of shield tunnels is relatively regular, and 3D laser scanning technology has been widely used in construction and operation [23]. Liu et al. [24] of Tongji University generated point cloud slices by the projection method, found the center of the structure through the point cloud slice, and monitored the change in the radius of the tunnel structure to judge the tunnel deformation [11, 25]; Van Gosliga and Pfeifer [25]. The least squares method is used to fit the actual line shape of the tunnel; Li [26]. performed ellipse fitting on point cloud slices, hoping to judge the deformation of the tunnel through the ellipse change after fitting; Guangzhou Metro Design and Research Institute Co., Ltd. Zhang of the company [27] used the 3D laser scanning results to detect the ovality of shield tunnels [28, 29].


## 3. Mathematical Theory of Fitting

Surface reconstruction methods are mainly divided into three categories: parameter reconstruction, triangular mesh method, and implicit function reconstruction [30, 31].The interpolation function is one of the methods of parameter reconstruction. Given a known point, construct a curve or surface, requiring the curve or surface to pass through all known points [32]. It is an important method of discrete function approximation [33].The triangular mesh method usually refers to the Delaunay triangulation method [34], which generally constructs a triangular surface from the two nearest points at a certain starting point and constructs the entire area into a triangular surface network through continuous growth [35].Implicit function methods include radial basis [36], Poisson construction method, etc. the implicit surface is to first separate the object point cloud, and then perform the module on the processed data by default, and each module is fitted together and then joined together. Using the Poisson equation to solve this problem can reduce the above calculations [37]. Because of the difficulty of 3D laser scanning monitoring and measurement in the construction of the New Austrian Tunnel, few people are involved [38].


## 4. Delaunay Meshing Method

### 4.1. Theory of Delaunay Meshing Method

The meshing method is to construct a mesh close to the original surface after the known data control points are connected through a certain relationship method. It is a common method for processing 3D models. Most of the existing meshing methods are as follows. It is improved on the theoretical basis of the Delaunay triangulation method. The Delaunay triangulation method usually converts a point cloud into a point set first, and then maps it from space to a plane, and then makes the points on the plane form a closed line segment, and then, according to the formed line segment, extends outwards, so it is suitable for point cloud data with a small amount of point cloud data and a tighter topological relationship, such as the generation of finite element meshes [5]. Therefore, for different types of point cloud data forms, the meshing method needs to be changed accordingly.The Delaunay triangulation method needs to grow the grid unit to form a complete surface. There are usually three widely used growth algorithms, the triangulation growth method, the segmented area method, and the point-by-point insertion method [39].(1)
The principle of the triangulation growth method is to search for the edge formed by the two closest endpoints in the point concentration, and then search for the endpoints at the other end of the side according to the judgment method of the triangulation network, and then proceed to the next round of edges. The determination of this work is carried out continuously until the connection of all endpoints is completed. The specific steps are: first select the first endpoint, generally take the point at the geometric midline position to facilitate growth around, and then, search for the closest point near this endpoint, connect the starting endpoint with the searched point; and use the triangle division criterion to find the third point and connect it to get the initial division unit; the second step is based on the initial unit, with the three edges of the unit as the initial edge, search for the corresponding endpoint and connect it according to the circumscribed hollow sphere criterion. Finally, repeat the above steps until the points in the point set are processed.(2)
The segmentation method mainly divides the point set into small and disjoint subsets by dividing and categorizing the point set, minimizing it to the set containing only three points, and then connecting them by these three points. The Delaunay triangulation unit is formed, and then, the adjacent subsets are gradually merged, and finally the Delaunay triangulation of the entire point set is formed. The specific steps are: first divide the point set on the plane intox and y directions through the coordinate system, and sort on the y coordinate based on the x coordinate; then set the threshold so that the number of subsets will pass the judgment If the number is greater than the subthreshold, divide left and right until it reaches below the threshold, and then, the point set in the area can be connected to form a Delaunay triangulation.(3)
The idea of point-by-point interpolation is exactly different from that of the triangulation growth method. The triangulation growth method is from the inside to the outside, and the initial Delaunay triangulation unit is formed by points and expands outward, while the point-by-point interpolation method first forms a triangulation mesh in the polygon containing all the point sets and then establishes a triangulation mesh for the internal points. The specific steps are as follows: first, set a polygon to contain all the points in the point concentration, then construct the initial triangulation unit in the set polygon, and then iterate the next step; finally, insert any data point in the polygon, and place it in the triangulation. Search for the triangle containing the point and connect it with the triangle [40].


### 4.2. Analysis of the Tunnel Surface Calculation Example

Relying on the point cloud data of a certain tunnel, this paper constructs the surface based on MATLAB. The point cloud data is sorted at a certain interval after processing, and has a certain topological relationship. The central axis of the tunnel data is used as thex-axis, then x is arranged on the x-axis according to the order and stored in the point cloud text according to the corresponding order. When the points are arranged at a fixed distance, the triangulation method is the same as the quadrilateral surface. However, the formation of triangles increases the amount of calculation, and the boundary formed on the surface will affect the effect of the curved surface model, so based on the Delaunay triangulation method, the growth method is used to expand to adapt to the ordered point cloud. The specific algorithm flow chart is shown in Figure 1.Figure 1
Flow chart of meshing algorithm.Specific steps are as follows:(1)
import the point cloud data, extract the coordinate values of the points (x, y, z), sort the order of the point sets, and sort them according to the numerical value from large-to-small to form a three-column matrix P;(2)
We calculate the closest distance to another point according to the spatial value of each point and record and save it, connect the two points and find the third point closest to the midpoint through line segment calculation; and determine whether the third point meets the maximum internal angle, and the maximum is connected. The point gets the triangle bag as the unit point set;(3)
look for the fourth point on the opposite side of the third point of the unit point set, and connect the first, second, and fourth points if they are satisfied;(4)
We determine whether the fourth point satisfies the triangle circumcircle according to the maximum empty circumcircle criterion, if yes, reverse the diagonals, and the optimization ends.(5)
We determine whether the distance between the triangular nets is greater than 5 cm, if it is satisfied, it will be rejected, and finally, the above steps will be looped until the calculation is completed.According to the triangulation method, the sampling interval is 5. The point cloud data of cm is constructed. Because the tunnel surface is not a polyhedron, the triangulation that exceeds the point cloud interval is eliminated, and the final results are shown in Figures2 and 3.Figure 2
Tunnel surface constructed by meshing method.Figure 3
Local area of tunnel surface constructed by meshing method.From the results in Figures2 and 3, it can be seen that after the tunnel surface is constructed by meshing, 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. This result shows that the curved surface constructed by the meshing method has the shortcomings of poor smoothness and insufficient local accuracy.


### 4.3. Application Scope Analysis

In this method, each point in the point cloud is on the curved surface. The author scans the same tunnel twice to obtain a point cloud with a distance of 1 mm. The measuring points are arranged in a longitudinal direction of 1 cm and a circumferential direction of one degrees angle. The two fitting surfaces are compared. The mean is 0.0241 mm and the variance is 0.0002057 mm2. It can be seen that the height of the fitted surface is consistent with the actual surface, and the fitting accuracy is within 0.1 mm. However, the high-order surface of the surface is discontinuous, and the normal direction of the surface will change abruptly. Therefore, the meshing method can be used as a surface fitting method for tunnel contour detection. This method should not be used as a surface fitting method for tunnel monitoring measurement and flatness detection. However, if the point cloud density is large enough or the problem of sudden change in the normal direction of the surface is alleviated, it can also be used as a surface fitting method for tunnel monitoring and measurement.


## 4.1. Theory of Delaunay Meshing Method

The meshing method is to construct a mesh close to the original surface after the known data control points are connected through a certain relationship method. It is a common method for processing 3D models. Most of the existing meshing methods are as follows. It is improved on the theoretical basis of the Delaunay triangulation method. The Delaunay triangulation method usually converts a point cloud into a point set first, and then maps it from space to a plane, and then makes the points on the plane form a closed line segment, and then, according to the formed line segment, extends outwards, so it is suitable for point cloud data with a small amount of point cloud data and a tighter topological relationship, such as the generation of finite element meshes [5]. Therefore, for different types of point cloud data forms, the meshing method needs to be changed accordingly.The Delaunay triangulation method needs to grow the grid unit to form a complete surface. There are usually three widely used growth algorithms, the triangulation growth method, the segmented area method, and the point-by-point insertion method [39].(1)
The principle of the triangulation growth method is to search for the edge formed by the two closest endpoints in the point concentration, and then search for the endpoints at the other end of the side according to the judgment method of the triangulation network, and then proceed to the next round of edges. The determination of this work is carried out continuously until the connection of all endpoints is completed. The specific steps are: first select the first endpoint, generally take the point at the geometric midline position to facilitate growth around, and then, search for the closest point near this endpoint, connect the starting endpoint with the searched point; and use the triangle division criterion to find the third point and connect it to get the initial division unit; the second step is based on the initial unit, with the three edges of the unit as the initial edge, search for the corresponding endpoint and connect it according to the circumscribed hollow sphere criterion. Finally, repeat the above steps until the points in the point set are processed.(2)
The segmentation method mainly divides the point set into small and disjoint subsets by dividing and categorizing the point set, minimizing it to the set containing only three points, and then connecting them by these three points. The Delaunay triangulation unit is formed, and then, the adjacent subsets are gradually merged, and finally the Delaunay triangulation of the entire point set is formed. The specific steps are: first divide the point set on the plane intox and y directions through the coordinate system, and sort on the y coordinate based on the x coordinate; then set the threshold so that the number of subsets will pass the judgment If the number is greater than the subthreshold, divide left and right until it reaches below the threshold, and then, the point set in the area can be connected to form a Delaunay triangulation.(3)
The idea of point-by-point interpolation is exactly different from that of the triangulation growth method. The triangulation growth method is from the inside to the outside, and the initial Delaunay triangulation unit is formed by points and expands outward, while the point-by-point interpolation method first forms a triangulation mesh in the polygon containing all the point sets and then establishes a triangulation mesh for the internal points. The specific steps are as follows: first, set a polygon to contain all the points in the point concentration, then construct the initial triangulation unit in the set polygon, and then iterate the next step; finally, insert any data point in the polygon, and place it in the triangulation. Search for the triangle containing the point and connect it with the triangle [40].


## 4.2. Analysis of the Tunnel Surface Calculation Example

Relying on the point cloud data of a certain tunnel, this paper constructs the surface based on MATLAB. The point cloud data is sorted at a certain interval after processing, and has a certain topological relationship. The central axis of the tunnel data is used as thex-axis, then x is arranged on the x-axis according to the order and stored in the point cloud text according to the corresponding order. When the points are arranged at a fixed distance, the triangulation method is the same as the quadrilateral surface. However, the formation of triangles increases the amount of calculation, and the boundary formed on the surface will affect the effect of the curved surface model, so based on the Delaunay triangulation method, the growth method is used to expand to adapt to the ordered point cloud. The specific algorithm flow chart is shown in Figure 1.Figure 1
Flow chart of meshing algorithm.Specific steps are as follows:(1)
import the point cloud data, extract the coordinate values of the points (x, y, z), sort the order of the point sets, and sort them according to the numerical value from large-to-small to form a three-column matrix P;(2)
We calculate the closest distance to another point according to the spatial value of each point and record and save it, connect the two points and find the third point closest to the midpoint through line segment calculation; and determine whether the third point meets the maximum internal angle, and the maximum is connected. The point gets the triangle bag as the unit point set;(3)
look for the fourth point on the opposite side of the third point of the unit point set, and connect the first, second, and fourth points if they are satisfied;(4)
We determine whether the fourth point satisfies the triangle circumcircle according to the maximum empty circumcircle criterion, if yes, reverse the diagonals, and the optimization ends.(5)
We determine whether the distance between the triangular nets is greater than 5 cm, if it is satisfied, it will be rejected, and finally, the above steps will be looped until the calculation is completed.According to the triangulation method, the sampling interval is 5. The point cloud data of cm is constructed. Because the tunnel surface is not a polyhedron, the triangulation that exceeds the point cloud interval is eliminated, and the final results are shown in Figures2 and 3.Figure 2
Tunnel surface constructed by meshing method.Figure 3
Local area of tunnel surface constructed by meshing method.From the results in Figures2 and 3, it can be seen that after the tunnel surface is constructed by meshing, 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. This result shows that the curved surface constructed by the meshing method has the shortcomings of poor smoothness and insufficient local accuracy.


## 4.3. Application Scope Analysis

In this method, each point in the point cloud is on the curved surface. The author scans the same tunnel twice to obtain a point cloud with a distance of 1 mm. The measuring points are arranged in a longitudinal direction of 1 cm and a circumferential direction of one degrees angle. The two fitting surfaces are compared. The mean is 0.0241 mm and the variance is 0.0002057 mm2. It can be seen that the height of the fitted surface is consistent with the actual surface, and the fitting accuracy is within 0.1 mm. However, the high-order surface of the surface is discontinuous, and the normal direction of the surface will change abruptly. Therefore, the meshing method can be used as a surface fitting method for tunnel contour detection. This method should not be used as a surface fitting method for tunnel monitoring measurement and flatness detection. However, if the point cloud density is large enough or the problem of sudden change in the normal direction of the surface is alleviated, it can also be used as a surface fitting method for tunnel monitoring and measurement.


## 5. Poisson Surface Construction

### 5.1. Poisson Equation

The Poisson equation is a partial differential equation commonly used in mechanical engineering and theoretical physics in mathematics, and it is derived from the Laplace equation [7]. It is precisely because of its extremely extensive application foundation in the fields of electromagnetics and astronomy, coupled with its advantages in image processing, so modern architectural design, survey photography, film, and television animation, and other fields that require three-dimensional modeling are all adopted. The loose algorithm builds a three-dimensional model.The Poisson equation is first expressed as the Laplace equation without other terms. The△φ=0 Laplace equation is also called the harmonic equation, which is a kind of partial differential equation.The expression of the Laplace equation is as follows:(1)∇2u=∂2u∂x2+∂2u∂y2+∂2u∂z2=0.In the formula:∇2 known as the Laplace operator, the Laplace equation here is a second-order partial differential equation, in the case of three dimensions, the Laplace equation boils down to solving the second-order differentiable real independent variable x,y,z real function.Also, the Poisson equation is(2)△φ=f,f and φ can be equations of real or complex numbers on a manifold, so when a manifold is a member of Euclidean space, Poisson’s equation is ∇2φ=f, so based on the fact that the right-hand side of the Laplace equation is 0, the Poisson’s equation in three-dimensional coordinates can be written.(3)∂2∂x2+∂2∂y2+∂2∂z2φx,y,z=fx,y,z.


### 5.2. Poisson’s Construction Principle

The Poisson reconstruction is a very intuitive method. Its core idea is that the point cloud represents the position of the object surface, and its normal vector represents the direction inside and outside. By implicitly fitting an indicator function derived from an object, an estimate of a smooth object surface can be given. In this paper, due to the uniform distribution of data, the surface is constructed based on the uniform distribution of point cloud samples.The surface construction method based on the Poisson algorithm is a kind of mesh reconstruction algorithm based on implicit function, mainly, through the optimization interpolation to fit the surface, the surface obtained is an approximate surface, and there is a certain error. The main idea is to take the point cloud data containing outer normals as input data. According to our sample setS each S contains multiple Pi points with internal and external normals Ni. If all data points are on or near the surface of an unknown model M, Poisson equation is constructed through the internal relationship between the indicator function of M and the gradient, then, based on the indicator function and M isosurface extraction, the final output surface model of a process. Define an indicator function to describe the entity model, this function is defined as 1 inside the model and 0 outside the model, the indicator function χM is as follows:(4)χM=0,M∈themodeloutside,1,M∈themodelinternal.Letz=fx,y be a function in the region D of the plane, the function has the first continuous partial derivative, any point on the function x,y∈D, The vector ∂f/∂xi+∂f/∂yj is the gradient of the function z=fx,y at Px,y, Do remember gradfx,y, Is the gradfx,y=∂f/∂xi+∂f/∂yj.gradfx,y is a vector. It is in the same direction as the maximum derivative at that point, and the gradient is equal to the derivative in that direction. The gradient value of the indicator function χ can be calculated according to the above formula, since the indicator function does not change almost anywhere on the surface of the object and the gradient is a vector field of 0, the gradient of the indicator function χM is equal to the inner and outer normals of the model surface. The inner relation between a set of directed sampling points and the indicator function is constructed. The gradient of the indicator function can be used to represent the set of directed sampling points [41], That is,(5)∇χ=V⟶,∇ stands for differential operator and the above equation is Poisson’s equation.Figure4 represents the set of image points, the gradient of indicator function, indicator function, and model surface, respectively.Figure 4
Poisson building process.Thus, the problem of indicating function can be transformed into that of calculating gradient. The goal is to find an indicator function that minimizes the vector region whose gradient is equal to the point cloud data that isminχ=∇χ−V⟶. Then, using the divergence operator, find an indicator function χ, so that the Laplace operator (the divergence of the gradient) is equal to the divergence of V⟶, which stands for(6)Δχ=∇·∇χ=∇·V⟶.


### 5.3. Poisson Construction Process

#### 5.3.1. Define Space

First, the point cloud data must select the space of the function to discretize the data. However, such a unified structure becomes impractical to reconstruct for fine detail. Because the dimensions of space are cubic in resolution, the number of triangles increases quadratic, so it needs to make the implicit function represent only the reconstructed surface, point cloud is stored by octree [42], and improve the efficiency of data extraction. The vector V is then introduced, where V is the set of vector fields and Fo is the node function.(7)Foq=Fq−o⋅co⋅w1o⋅w3,o is the node, o⋅c indicates the node center, and o⋅w indicates the width of node o; at this point, the vector field V can be precisely a linear Fo set; the Poisson equation is expressed by matrix and the indicator function Fo can be calculated effectively.


#### 5.3.2. Selecting the Basis Function

When choosing the basis functionF, the objective is to select a function such that the vector field V defined in the equation can be accurately represented as a linear sum of node functions Fo.When selecting a function, you need the linear sum ofFo to represent the value of the function at a point in the set, let q be that point, so we just have to sum over the node o that is close to q. Therefore, we set F as the n-th convolution of the Box filter, resulting in the basis function F:(8)Fx,y,z=BxByBz∗n.Type:∗ stands for convolution symbol, ∗n represents n identical functions convolving the B (a) calculation formula of are as follows:(9)Ba=1,a༜0.5,0,otherwise.


#### 5.3.3. Defining the Vector Field

Because the target data is uniformly sampled, it is generally considered that the divided module is a constant. The vector fieldV needs to be made close to the gradient of the indicator function, so the gradient field expression defining the indicator function is approximately(10)V⟶q=∑s∈S∑o∈NgbrDsαo,Foq⋅s⋅N⟶.Type:S stands for convergence of points, s is the region adjacent to K at any point in the point set, NgbrD represents the eight nearest nodes of the current node o, and s⋅N⟶ is the vertex normal vector of node o.


#### 5.3.4. Solving Poisson’s Equation

After obtaining the vector fieldV, Poisson equation Δχ=∇⋅V⟶ is established. It can be calculated by solving the minimum function χ, vector V is the dimension ϑ. Therefore, for the coordinate o of node vo=∇⋅V⟶,Fo, for the solution χ, we need to approach v by the vector formed by the Laplace operator of the projection χ, the function, and Fo in the function space. This translates to ∑o∈ϑv−vo, that is, minx∈RϑLx−v is solved [43].


#### 5.3.5. Isosurface Extraction

BecauseL is relatively symmetric, it has a fixed structure in space. Therefore, the heavy grid method [44] can be used to solve the problem, and in order to obtain a reconstructed surface ∂M, one must first select a standard value and then compute from the indicator function. The standard values are chosen so that the extracted surface is very close to the position of the input sample. We do this by evaluating the sample location and using the average of the standard value surface extraction values. After the standard value is obtained, the surface can be obtained by the cube method [45], and finally, the modules in each region are spliced to form a complete surface.(11)∂M=q∈R3χq=γ,(12)V⟶q=∑s∈S∑o∈NgbrDsαo,Foq⋅s⋅N⟶.


### 5.4. Example Analysis of Tunnel Surface Construction

In this section, based on MATLAB, Poisson construction method is adopted to realize the surface construction of point cloud data of Thunder Grass Mountain tunnel. The programming process is shown in Figure5.Figure 5
Algorithm flow chart of constructing surface by Poisson.Point cloud data with an actual interval of 5 cm was used to construct the surface. In order to more intuitively reflect the surface effect, point cloud data of the arch with large curvature were selected to construct the surface. The tunnel surface constructed by the Poisson construction method is shown in Figure6.Figure 6
Tunnel surface constructed by Poisson.It can be seen from Figure6 that a complete surface can be obtained after the point cloud is constructed by Poisson, and the concave and convex contour of the surface is relatively clear. On this basis, the Y-axis area from −0.5 m∼3.0 m is enlarged to obtain Figure 7, and the results are shown in the figure.Figure 7
Local area of tunnel surface constructed by Poisson.Figure7 show the tunnel surface obtained by Poisson reconstruction method is continuous and complete, which can reflect the concave and convex condition of the tunnel surface.Figure8 shows that the Poisson construction surface (grey) is highly coincident with the point cloud (red). When the point cloud is further thinned, the surface can still be further smoothed. When the full-section tunnel surface is fitted, the range of the Poisson construction surface will be larger than the range of the point cloud surface, and manual intervention is required to set the boundary of the surface fitting. The Poisson-constructed surfaces of the point cloud obtained after scanning the same section of the tunnel twice are also almost completely coincident.Figure 8
The full-section tunnel surface obtained by Poisson construction.


### 5.5. Application Scope Analysis

The points in the point cloud of this method are all on the surface and the high order is continuous, and the normal direction of the surface will not change abruptly. Therefore, the Poisson-constructed surface can be used as a fitting surface for tunnel monitoring and measurement. However, this method requires manual intervention and the setting of surface fitting boundaries, which is not suitable for the surface fitting of tunnel contour detection with low requirements on accuracy and high efficiency. If the point cloud is thinned, the Poisson-constructed surface can also achieve a smoother effect, meeting the requirements of the reference surface for flatness detection.


## 5.1. Poisson Equation

The Poisson equation is a partial differential equation commonly used in mechanical engineering and theoretical physics in mathematics, and it is derived from the Laplace equation [7]. It is precisely because of its extremely extensive application foundation in the fields of electromagnetics and astronomy, coupled with its advantages in image processing, so modern architectural design, survey photography, film, and television animation, and other fields that require three-dimensional modeling are all adopted. The loose algorithm builds a three-dimensional model.The Poisson equation is first expressed as the Laplace equation without other terms. The△φ=0 Laplace equation is also called the harmonic equation, which is a kind of partial differential equation.The expression of the Laplace equation is as follows:(1)∇2u=∂2u∂x2+∂2u∂y2+∂2u∂z2=0.In the formula:∇2 known as the Laplace operator, the Laplace equation here is a second-order partial differential equation, in the case of three dimensions, the Laplace equation boils down to solving the second-order differentiable real independent variable x,y,z real function.Also, the Poisson equation is(2)△φ=f,f and φ can be equations of real or complex numbers on a manifold, so when a manifold is a member of Euclidean space, Poisson’s equation is ∇2φ=f, so based on the fact that the right-hand side of the Laplace equation is 0, the Poisson’s equation in three-dimensional coordinates can be written.(3)∂2∂x2+∂2∂y2+∂2∂z2φx,y,z=fx,y,z.


## 5.2. Poisson’s Construction Principle

The Poisson reconstruction is a very intuitive method. Its core idea is that the point cloud represents the position of the object surface, and its normal vector represents the direction inside and outside. By implicitly fitting an indicator function derived from an object, an estimate of a smooth object surface can be given. In this paper, due to the uniform distribution of data, the surface is constructed based on the uniform distribution of point cloud samples.The surface construction method based on the Poisson algorithm is a kind of mesh reconstruction algorithm based on implicit function, mainly, through the optimization interpolation to fit the surface, the surface obtained is an approximate surface, and there is a certain error. The main idea is to take the point cloud data containing outer normals as input data. According to our sample setS each S contains multiple Pi points with internal and external normals Ni. If all data points are on or near the surface of an unknown model M, Poisson equation is constructed through the internal relationship between the indicator function of M and the gradient, then, based on the indicator function and M isosurface extraction, the final output surface model of a process. Define an indicator function to describe the entity model, this function is defined as 1 inside the model and 0 outside the model, the indicator function χM is as follows:(4)χM=0,M∈themodeloutside,1,M∈themodelinternal.Letz=fx,y be a function in the region D of the plane, the function has the first continuous partial derivative, any point on the function x,y∈D, The vector ∂f/∂xi+∂f/∂yj is the gradient of the function z=fx,y at Px,y, Do remember gradfx,y, Is the gradfx,y=∂f/∂xi+∂f/∂yj.gradfx,y is a vector. It is in the same direction as the maximum derivative at that point, and the gradient is equal to the derivative in that direction. The gradient value of the indicator function χ can be calculated according to the above formula, since the indicator function does not change almost anywhere on the surface of the object and the gradient is a vector field of 0, the gradient of the indicator function χM is equal to the inner and outer normals of the model surface. The inner relation between a set of directed sampling points and the indicator function is constructed. The gradient of the indicator function can be used to represent the set of directed sampling points [41], That is,(5)∇χ=V⟶,∇ stands for differential operator and the above equation is Poisson’s equation.Figure4 represents the set of image points, the gradient of indicator function, indicator function, and model surface, respectively.Figure 4
Poisson building process.Thus, the problem of indicating function can be transformed into that of calculating gradient. The goal is to find an indicator function that minimizes the vector region whose gradient is equal to the point cloud data that isminχ=∇χ−V⟶. Then, using the divergence operator, find an indicator function χ, so that the Laplace operator (the divergence of the gradient) is equal to the divergence of V⟶, which stands for(6)Δχ=∇·∇χ=∇·V⟶.


## 5.3. Poisson Construction Process

### 5.3.1. Define Space

First, the point cloud data must select the space of the function to discretize the data. However, such a unified structure becomes impractical to reconstruct for fine detail. Because the dimensions of space are cubic in resolution, the number of triangles increases quadratic, so it needs to make the implicit function represent only the reconstructed surface, point cloud is stored by octree [42], and improve the efficiency of data extraction. The vector V is then introduced, where V is the set of vector fields and Fo is the node function.(7)Foq=Fq−o⋅co⋅w1o⋅w3,o is the node, o⋅c indicates the node center, and o⋅w indicates the width of node o; at this point, the vector field V can be precisely a linear Fo set; the Poisson equation is expressed by matrix and the indicator function Fo can be calculated effectively.


### 5.3.2. Selecting the Basis Function

When choosing the basis functionF, the objective is to select a function such that the vector field V defined in the equation can be accurately represented as a linear sum of node functions Fo.When selecting a function, you need the linear sum ofFo to represent the value of the function at a point in the set, let q be that point, so we just have to sum over the node o that is close to q. Therefore, we set F as the n-th convolution of the Box filter, resulting in the basis function F:(8)Fx,y,z=BxByBz∗n.Type:∗ stands for convolution symbol, ∗n represents n identical functions convolving the B (a) calculation formula of are as follows:(9)Ba=1,a༜0.5,0,otherwise.


### 5.3.3. Defining the Vector Field

Because the target data is uniformly sampled, it is generally considered that the divided module is a constant. The vector fieldV needs to be made close to the gradient of the indicator function, so the gradient field expression defining the indicator function is approximately(10)V⟶q=∑s∈S∑o∈NgbrDsαo,Foq⋅s⋅N⟶.Type:S stands for convergence of points, s is the region adjacent to K at any point in the point set, NgbrD represents the eight nearest nodes of the current node o, and s⋅N⟶ is the vertex normal vector of node o.


### 5.3.4. Solving Poisson’s Equation

After obtaining the vector fieldV, Poisson equation Δχ=∇⋅V⟶ is established. It can be calculated by solving the minimum function χ, vector V is the dimension ϑ. Therefore, for the coordinate o of node vo=∇⋅V⟶,Fo, for the solution χ, we need to approach v by the vector formed by the Laplace operator of the projection χ, the function, and Fo in the function space. This translates to ∑o∈ϑv−vo, that is, minx∈RϑLx−v is solved [43].


### 5.3.5. Isosurface Extraction

BecauseL is relatively symmetric, it has a fixed structure in space. Therefore, the heavy grid method [44] can be used to solve the problem, and in order to obtain a reconstructed surface ∂M, one must first select a standard value and then compute from the indicator function. The standard values are chosen so that the extracted surface is very close to the position of the input sample. We do this by evaluating the sample location and using the average of the standard value surface extraction values. After the standard value is obtained, the surface can be obtained by the cube method [45], and finally, the modules in each region are spliced to form a complete surface.(11)∂M=q∈R3χq=γ,(12)V⟶q=∑s∈S∑o∈NgbrDsαo,Foq⋅s⋅N⟶.


## 5.3.1. Define Space

First, the point cloud data must select the space of the function to discretize the data. However, such a unified structure becomes impractical to reconstruct for fine detail. Because the dimensions of space are cubic in resolution, the number of triangles increases quadratic, so it needs to make the implicit function represent only the reconstructed surface, point cloud is stored by octree [42], and improve the efficiency of data extraction. The vector V is then introduced, where V is the set of vector fields and Fo is the node function.(7)Foq=Fq−o⋅co⋅w1o⋅w3,o is the node, o⋅c indicates the node center, and o⋅w indicates the width of node o; at this point, the vector field V can be precisely a linear Fo set; the Poisson equation is expressed by matrix and the indicator function Fo can be calculated effectively.


## 5.3.2. Selecting the Basis Function

When choosing the basis functionF, the objective is to select a function such that the vector field V defined in the equation can be accurately represented as a linear sum of node functions Fo.When selecting a function, you need the linear sum ofFo to represent the value of the function at a point in the set, let q be that point, so we just have to sum over the node o that is close to q. Therefore, we set F as the n-th convolution of the Box filter, resulting in the basis function F:(8)Fx,y,z=BxByBz∗n.Type:∗ stands for convolution symbol, ∗n represents n identical functions convolving the B (a) calculation formula of are as follows:(9)Ba=1,a༜0.5,0,otherwise.


## 5.3.3. Defining the Vector Field

Because the target data is uniformly sampled, it is generally considered that the divided module is a constant. The vector fieldV needs to be made close to the gradient of the indicator function, so the gradient field expression defining the indicator function is approximately(10)V⟶q=∑s∈S∑o∈NgbrDsαo,Foq⋅s⋅N⟶.Type:S stands for convergence of points, s is the region adjacent to K at any point in the point set, NgbrD represents the eight nearest nodes of the current node o, and s⋅N⟶ is the vertex normal vector of node o.


## 5.3.4. Solving Poisson’s Equation

After obtaining the vector fieldV, Poisson equation Δχ=∇⋅V⟶ is established. It can be calculated by solving the minimum function χ, vector V is the dimension ϑ. Therefore, for the coordinate o of node vo=∇⋅V⟶,Fo, for the solution χ, we need to approach v by the vector formed by the Laplace operator of the projection χ, the function, and Fo in the function space. This translates to ∑o∈ϑv−vo, that is, minx∈RϑLx−v is solved [43].


## 5.3.5. Isosurface Extraction

BecauseL is relatively symmetric, it has a fixed structure in space. Therefore, the heavy grid method [44] can be used to solve the problem, and in order to obtain a reconstructed surface ∂M, one must first select a standard value and then compute from the indicator function. The standard values are chosen so that the extracted surface is very close to the position of the input sample. We do this by evaluating the sample location and using the average of the standard value surface extraction values. After the standard value is obtained, the surface can be obtained by the cube method [45], and finally, the modules in each region are spliced to form a complete surface.(11)∂M=q∈R3χq=γ,(12)V⟶q=∑s∈S∑o∈NgbrDsαo,Foq⋅s⋅N⟶.


## 5.4. Example Analysis of Tunnel Surface Construction

In this section, based on MATLAB, Poisson construction method is adopted to realize the surface construction of point cloud data of Thunder Grass Mountain tunnel. The programming process is shown in Figure5.Figure 5
Algorithm flow chart of constructing surface by Poisson.Point cloud data with an actual interval of 5 cm was used to construct the surface. In order to more intuitively reflect the surface effect, point cloud data of the arch with large curvature were selected to construct the surface. The tunnel surface constructed by the Poisson construction method is shown in Figure6.Figure 6
Tunnel surface constructed by Poisson.It can be seen from Figure6 that a complete surface can be obtained after the point cloud is constructed by Poisson, and the concave and convex contour of the surface is relatively clear. On this basis, the Y-axis area from −0.5 m∼3.0 m is enlarged to obtain Figure 7, and the results are shown in the figure.Figure 7
Local area of tunnel surface constructed by Poisson.Figure7 show the tunnel surface obtained by Poisson reconstruction method is continuous and complete, which can reflect the concave and convex condition of the tunnel surface.Figure8 shows that the Poisson construction surface (grey) is highly coincident with the point cloud (red). When the point cloud is further thinned, the surface can still be further smoothed. When the full-section tunnel surface is fitted, the range of the Poisson construction surface will be larger than the range of the point cloud surface, and manual intervention is required to set the boundary of the surface fitting. The Poisson-constructed surfaces of the point cloud obtained after scanning the same section of the tunnel twice are also almost completely coincident.Figure 8
The full-section tunnel surface obtained by Poisson construction.


## 5.5. Application Scope Analysis

The points in the point cloud of this method are all on the surface and the high order is continuous, and the normal direction of the surface will not change abruptly. Therefore, the Poisson-constructed surface can be used as a fitting surface for tunnel monitoring and measurement. However, this method requires manual intervention and the setting of surface fitting boundaries, which is not suitable for the surface fitting of tunnel contour detection with low requirements on accuracy and high efficiency. If the point cloud is thinned, the Poisson-constructed surface can also achieve a smoother effect, meeting the requirements of the reference surface for flatness detection.


## 6. Lagrange Interpolation

### 6.1. Lagrange Interpolation

Lagrange interpolation is usually used to fit curves. In the coordinate system, if we know the value ofy=fx at n plus 1 different points x0,x1…xn, the values of the function are y0,y1…yn, then the function x0,y0,x1,y1…xn,yn passes rough the n + 1 points. Then, it can be considered to construct a polynomial y=Pnx whose degree does not exceed n over the n + 1 points so that it satisfies(13)Pnxk=yk,k=0,1,…n∗.To estimate any pointξ,ξ≠xi,i=0,1,…n, we can use the value of Pnξ as an approximation of the exact value fξ. Equation ∗ is the interpolation condition and contains the xi (I = =0, 1, ..., n), the minimum interval formed. Among them, polynomials of less than n degrees that meet interpolation conditions exist and are unique.So, there are a total ofn points x0,y0,x1,y1…xn−1,yn−1 on the plane. We are going to draw a function fx that takes its graph through these n points; practice the following: let the set Dn be the set of indices of points, x,yDn=0,1,…n−1. Let us do n polynomials, Pjx,j∈Dn. For any of them, k∈Dn, we have Pkx,Bk=ii≠k,i∈Dn.(14)Pkx=∏i∈Bkx−xixk−xi.In the calculation formula,pkx is a degree N − 1 polynomial, and we meet the conditions ∀m∈Ik,pkxm=0 and pkxk=1.This is what we get(15)Lnx=Σj=0n−1yjpjx.


### 6.2. Lagrange Interpolation Surface Construction Process

Through its one-dimensional definition method, it can be deduced that when there is a known nodexi,yi, its two-dimensional value is fxi,yi, among i=1,2,…,m;j=1,2,…,n. In this case, the Lagrange interpolation polynomial is(16)Lx,y=∑1≤i≤m1≤j≤nfxi,yiLijx,y=∑j=1n∑i=1mMixfxi,yiNj.In the calculation formula,(17)Mix=∏r=1r≠imx−xrxi−xr,i=1,2,…,m;Njy=∏t=1t≠imy−ytyi−yt,j=1,2,…,n.In order to facilitate calculation in practical application, we can write it as a matrix.(18)Lx,y=N1yN2y…NnyTfx1,y1fx2,y1…fxm,y1fx1,y2fx2,y2…fxm,y2…………fx1,ynfx2,yn…fxm,ynM1xM2x…Mmx=NTy·F·Mx.Suppose there arek × s points xi′,yj′i=1,2,…,k;j=1,2,…,s to be interpolated, the formula to be reconstructed by interpolation can be expressed in the form of matrix.(19)L=Lx1′,y1′Lx2′,y1′…Lxk′,y1′Lx1′,y2′Lx2′,y2′…Lxk′,y2′……⋯…Lx1′,yS′Lx2′,yS′…Lxk′,yS′.In the calculation formula,Nj=Nyj′j=1,2,…,s;Mi=Mxi′i=1,2,…k.At the same time in Lagrange interpolation reconstruction using known points to interpolate a control point because polynomial times will increase more calculation, therefore, local interpolation can be performed by using the known points adjacent to known points to be constructed, so this gives us the formula [46]:(20)Lx,y=∑i∈Φpx∑j∈ΦqyMixfxi,yiNjy.In the calculation formula,Φpx represents the abscissa set of p points in the node whose abscissa is closest to X, and Φqy represents the abscissa set of q points in the node whose ordinate is closest to y. Therefore, the above formula represents the interpolation construction of the p×q points by using the known nodes closest to (x, y), where, because lkx=∏j=1j≠knx−xj/xk−xj,k=1,…,n is differentiable to order n minus 1 with respect to x, so equation (13) has a continuous derivative of order n-1, because Mix=∏r=1r≠imx−xr/xi−xr has partial derivative with respect to m minus 1 of x, so Njy=∏t=1t≠imx−xt/xi−xt has n minus 1 continuous partial derivatives with respect to y, so Lx,y has n − 1 and m − 1 partial derivatives with respect to x and y, respectively. If Φpx is the abscissa set of p points closest to the abscissa x of the interpolation point, and if Φqy is the set of x coordinates of q points closest to y of the interpolation point, then Mix and Njy in formulas (3)–(19) have continuous partial derivatives of order P−1 and q−1 with respect to x and y, respectively. So, the surfaces obtained during interpolation have continuous partial derivatives of order P−1 and q−1 with respect to x and y, respectively. It follows that when p,q≥2, it can guarantee the continuity of the surface.


### 6.3. Tunnel Surface Construction and Effect Analysis

The calculation is realized based on MATLAB. The specific method is as follows:(1)
First, we import point cloud data to obtain coordinate values and form two-dimensional functions of control points (X, Y) from X;(2)
The transpose matrix is obtained according to the calculation of control points;(3)
determine the number of interpolation points to form nodes and convert the nodes into matrices;(4)
Interpolating function is obtained by transposing matrix and two-dimensional function.According to the properties of Lagrange interpolation method, Lagrange interpolation is polynomial in nature and can fit a certain amount of point cloud data into a polynomial function surface. According to this property, multiple polynomial functions can be constructed fromX and Y. However, because of its characteristics, although the surface fitting effect can be achieved in a certain local range, the number of polynomial functions of multiple times fluctuates greatly, which is easy to lead to excessive calculation.This paper adopts Leica mountain tunnel point cloud data to carry out surface interpolation for local range point cloud. The data in this area have typical characteristics of tunnels and is located in tunnel vault with large curvature. The surface can be observed after interpolation. After the data of the experimental point cloud are processed, the surface is constructed by Lagrange interpolation method, as shown in Figures9 and 10.Figure 9
Tunnel surface constructed by Lagrange interpolation.Figure 10
Local area of tunnel surface constructed by the Lagrange interpolation.Analysis of Figures9 and 10 shows that the tunnel surface constructed by the Lagrange interpolation method can reflect the overall contour of the tunnel surface and show basic details of the tunnel surface in local areas. However, the surface constructed by this method has poor smoothness and a convex hull. The results show that the Lagrange interpolation method has some shortcomings in tunnel surface construction, such as low smoothness and partial convex hull.


### 6.4. Application Scope Analysis

The points in the Lagrange interpolation surface point cloud are all on the surface, and the high order is continuous, and the normal direction of the surface will not change abruptly. Therefore, the Lagrange interpolation surface can be used as a fitting surface for tunnel monitoring and measurement. However, considering that it takes a long time to calculate, it should not be used as a method for tunnel contour detection. If the point cloud is thinned, the Lagrange interpolation surface can also achieve a smoother effect because of its high-order continuous characteristics, which can meet the requirements of the reference surface for flatness detection.


## 6.1. Lagrange Interpolation

Lagrange interpolation is usually used to fit curves. In the coordinate system, if we know the value ofy=fx at n plus 1 different points x0,x1…xn, the values of the function are y0,y1…yn, then the function x0,y0,x1,y1…xn,yn passes rough the n + 1 points. Then, it can be considered to construct a polynomial y=Pnx whose degree does not exceed n over the n + 1 points so that it satisfies(13)Pnxk=yk,k=0,1,…n∗.To estimate any pointξ,ξ≠xi,i=0,1,…n, we can use the value of Pnξ as an approximation of the exact value fξ. Equation ∗ is the interpolation condition and contains the xi (I = =0, 1, ..., n), the minimum interval formed. Among them, polynomials of less than n degrees that meet interpolation conditions exist and are unique.So, there are a total ofn points x0,y0,x1,y1…xn−1,yn−1 on the plane. We are going to draw a function fx that takes its graph through these n points; practice the following: let the set Dn be the set of indices of points, x,yDn=0,1,…n−1. Let us do n polynomials, Pjx,j∈Dn. For any of them, k∈Dn, we have Pkx,Bk=ii≠k,i∈Dn.(14)Pkx=∏i∈Bkx−xixk−xi.In the calculation formula,pkx is a degree N − 1 polynomial, and we meet the conditions ∀m∈Ik,pkxm=0 and pkxk=1.This is what we get(15)Lnx=Σj=0n−1yjpjx.


## 6.2. Lagrange Interpolation Surface Construction Process

Through its one-dimensional definition method, it can be deduced that when there is a known nodexi,yi, its two-dimensional value is fxi,yi, among i=1,2,…,m;j=1,2,…,n. In this case, the Lagrange interpolation polynomial is(16)Lx,y=∑1≤i≤m1≤j≤nfxi,yiLijx,y=∑j=1n∑i=1mMixfxi,yiNj.In the calculation formula,(17)Mix=∏r=1r≠imx−xrxi−xr,i=1,2,…,m;Njy=∏t=1t≠imy−ytyi−yt,j=1,2,…,n.In order to facilitate calculation in practical application, we can write it as a matrix.(18)Lx,y=N1yN2y…NnyTfx1,y1fx2,y1…fxm,y1fx1,y2fx2,y2…fxm,y2…………fx1,ynfx2,yn…fxm,ynM1xM2x…Mmx=NTy·F·Mx.Suppose there arek × s points xi′,yj′i=1,2,…,k;j=1,2,…,s to be interpolated, the formula to be reconstructed by interpolation can be expressed in the form of matrix.(19)L=Lx1′,y1′Lx2′,y1′…Lxk′,y1′Lx1′,y2′Lx2′,y2′…Lxk′,y2′……⋯…Lx1′,yS′Lx2′,yS′…Lxk′,yS′.In the calculation formula,Nj=Nyj′j=1,2,…,s;Mi=Mxi′i=1,2,…k.At the same time in Lagrange interpolation reconstruction using known points to interpolate a control point because polynomial times will increase more calculation, therefore, local interpolation can be performed by using the known points adjacent to known points to be constructed, so this gives us the formula [46]:(20)Lx,y=∑i∈Φpx∑j∈ΦqyMixfxi,yiNjy.In the calculation formula,Φpx represents the abscissa set of p points in the node whose abscissa is closest to X, and Φqy represents the abscissa set of q points in the node whose ordinate is closest to y. Therefore, the above formula represents the interpolation construction of the p×q points by using the known nodes closest to (x, y), where, because lkx=∏j=1j≠knx−xj/xk−xj,k=1,…,n is differentiable to order n minus 1 with respect to x, so equation (13) has a continuous derivative of order n-1, because Mix=∏r=1r≠imx−xr/xi−xr has partial derivative with respect to m minus 1 of x, so Njy=∏t=1t≠imx−xt/xi−xt has n minus 1 continuous partial derivatives with respect to y, so Lx,y has n − 1 and m − 1 partial derivatives with respect to x and y, respectively. If Φpx is the abscissa set of p points closest to the abscissa x of the interpolation point, and if Φqy is the set of x coordinates of q points closest to y of the interpolation point, then Mix and Njy in formulas (3)–(19) have continuous partial derivatives of order P−1 and q−1 with respect to x and y, respectively. So, the surfaces obtained during interpolation have continuous partial derivatives of order P−1 and q−1 with respect to x and y, respectively. It follows that when p,q≥2, it can guarantee the continuity of the surface.


## 6.3. Tunnel Surface Construction and Effect Analysis

The calculation is realized based on MATLAB. The specific method is as follows:(1)
First, we import point cloud data to obtain coordinate values and form two-dimensional functions of control points (X, Y) from X;(2)
The transpose matrix is obtained according to the calculation of control points;(3)
determine the number of interpolation points to form nodes and convert the nodes into matrices;(4)
Interpolating function is obtained by transposing matrix and two-dimensional function.According to the properties of Lagrange interpolation method, Lagrange interpolation is polynomial in nature and can fit a certain amount of point cloud data into a polynomial function surface. According to this property, multiple polynomial functions can be constructed fromX and Y. However, because of its characteristics, although the surface fitting effect can be achieved in a certain local range, the number of polynomial functions of multiple times fluctuates greatly, which is easy to lead to excessive calculation.This paper adopts Leica mountain tunnel point cloud data to carry out surface interpolation for local range point cloud. The data in this area have typical characteristics of tunnels and is located in tunnel vault with large curvature. The surface can be observed after interpolation. After the data of the experimental point cloud are processed, the surface is constructed by Lagrange interpolation method, as shown in Figures9 and 10.Figure 9
Tunnel surface constructed by Lagrange interpolation.Figure 10
Local area of tunnel surface constructed by the Lagrange interpolation.Analysis of Figures9 and 10 shows that the tunnel surface constructed by the Lagrange interpolation method can reflect the overall contour of the tunnel surface and show basic details of the tunnel surface in local areas. However, the surface constructed by this method has poor smoothness and a convex hull. The results show that the Lagrange interpolation method has some shortcomings in tunnel surface construction, such as low smoothness and partial convex hull.


## 6.4. Application Scope Analysis

The points in the Lagrange interpolation surface point cloud are all on the surface, and the high order is continuous, and the normal direction of the surface will not change abruptly. Therefore, the Lagrange interpolation surface can be used as a fitting surface for tunnel monitoring and measurement. However, considering that it takes a long time to calculate, it should not be used as a method for tunnel contour detection. If the point cloud is thinned, the Lagrange interpolation surface can also achieve a smoother effect because of its high-order continuous characteristics, which can meet the requirements of the reference surface for flatness detection.


## 7. Conclusion and Prospect

As an information sensing device, the 3D laser scanning device can be used to build a 3D model after fitting the measured point cloud. The 3D model can be combined with 5 G, BIM, cloud computing, the Internet of things, and other technologies to achieve precise positioning and intelligent decision-making. In this way, precise tunnel excavation, intelligent construction, and intelligent operation and maintenance can be realized. Tunnel monitoring measurement, contour detection, and flatness detection require different fitting methods, and the abovementioned fitting methods can be applied to tunnel 3D laser scanning detection and monitoring to varying degrees.(1)
perform surface fitting first. In the data processing of tunnel detection and monitoring, the information in the three-dimensional laser scanning point cloud of the tunnel can be extracted as much as possible to overcome the shortcomings of scattered point clouds and improve the detection and monitoring accuracy.(2)
The meshing method fits the curved surface and the point cloud completely and can be used as a curved surface fitting method for tunnel contour detection. Because of its high-order discontinuity and the sudden change of the normal line, it is generally not suitable as a surface fitting method for tunnel monitoring and measurement, but it can still be considered as a surface fitting method for tunnel monitoring and measurement after the point cloud is encrypted. Considering the high-order discontinuity of this method, this method is not used as a fitting method for the data plane for tunnel flatness detection.(3)
The fitting surface and the point cloud of the Poisson reconstruction method are completely coincident and can be used as a surface fitting method for tunnel contour detection. However, this method will appear redundant surfaces in the full-section fitting of the tunnel, so the calculation efficiency is low, and it is not recommended. It is used as a fitting surface method for tunnel contour detection. Because the method is highly continuous and the fitting surface is stable, it is recommended to be used as the surface fitting method for tunnel monitoring and measurement. Considering the high-order continuity of this method, it can be used as a datum surface fitting method for tunnel flatness detection after the point cloud is thinned.(4)
The Lagrange’s construction surface is completely coincident with the point cloud and can be used as a surface fitting method for tunnel contour detection. However, this method has many iterations and takes a long time, so it is not recommended as a surface fitting method for tunnel contour detection. Because the method is highly continuous and the fitting surface is stable, it is recommended to be used as the surface fitting method for tunnel monitoring and measurement. Considering the high-order continuity of this method, it can be used as a datum surface fitting method for tunnel flatness detection after the point cloud is thinned.



---

*Source: 8432308-2022-09-28.xml*