Automatic Extraction of Tunnel Lining Cross-Sections from Terrestrial Laser Scanning Point Clouds

Tunnel lining (bare-lining) cross-sections play an important role in analyzing deformations of tunnel linings. The goal of this paper is to develop an automatic method for extracting bare-lining cross-sections from terrestrial laser scanning (TLS) point clouds. First, the combination of a 2D projection strategy and angle criterion is used for tunnel boundary point detection, from which we estimate the two boundary lines in the X-Y plane. The initial direction of the cross-sectional plane is defined to be orthogonal to one of the two boundary lines. In order to compute the final cross-sectional plane, the direction is adjusted twice with the total least squares method and Rodrigues’ rotation formula, respectively. The projection of nearby points is made onto the adjusted plane to generate tunnel cross-sections. Finally, we present a filtering algorithm (similar to the idea of the morphological erosion) to remove the non-lining points in the cross-section. The proposed method was implemented on railway tunnel data collected in Sichuan, China. Compared with an existing method of cross-sectional extraction, the proposed method can offer high accuracy and more reliable cross-sectional modeling. We also evaluated Type I and Type II errors of the proposed filter, at the same time, which gave suggestions on the parameter selection of the filter.


Introduction
Over the past years, terrestrial laser scanning (TLS) has been used in inverse engineering [1][2][3], for tracking the changes of natural surfaces via the comparison of different point clouds [4][5][6][7][8] and in the estimation of forest attributes [9][10][11]. The use of the TLS technique has also become popular in tunnel engineering due to the various advantages over conventional geodetic devices, such as laser beam profilers and total stations, which take more time to acquire data [12] and cannot offer high-density 3D datasets. The existing applications of TLS for tunnels contain geological feature detection [13], deformation analysis [14][15][16], and cross-sectional extraction [12,17], and the automatic processing of tunnel point clouds has received increasing research attention [12,18].
The checking of tunnel cross-sections is the primary method for deformation monitoring and clearance inspection. In addition to conventional geodetic surveys, several methods have been developed to extract tunnel cross-sections based on digital photogrammetry or the TLS technique. Combining photogrammetry and laser-lit spots, Wang et al. [19] improved a profile-image method for measuring cross-sections. Their method overcomes the limit of the number of points that conventional geodetic surveying has, but it is difficult to provide sufficient lighting conditions in an actual tunnel. TLS can offer active measurement even in a tunnel where there is no light, which makes it a preferred technique for the extraction of tunnel cross-sections. Several scholars used standard geometric models, such as an ellipse [15] or a circular cylinder [14], to approximate tunnel point clouds, and then estimated cross-sections from the fitted models. However, the non-lining points (i.e., pipes and equipment

Estimation of the Tunnel Boundary Lines in the X-Y Plane
An entire tunnel is generally a long tube, so it can be scanned along its centerline station by station. The tunnel point clouds are registered in a user-specified coordinate system using sphere reference targets. In this system, the origin is near the tunnel entrance, and the Y axis is oriented along the direction of the initial segment of the tunnel. Several algorithms [22][23][24] have been developed to quickly extract features from LIDAR point clouds based on projection and gridding. Projecting the scanned data onto the X-Y plane can simplify the 3D tube to a long and narrow 2D object, from which the two boundary points' groups are extracted from the both sides of the 2D object. In order to improve the speed of extraction, an algorithm for extracting the boundary point groups is proposed using a fixed grid.
The projections of the entire tunnel point clouds in the X-Y plane are discretized using a square grid. A grid size that is too large or too small will decrease the computational efficiency or the extraction precision of boundary cells, respectively. The appropriate size of the grid is about one-twentieth of the width of the tunnel. The value of is used to determine whether points exist or not in the cell , has a value of 1 if points exist and 0 if there are no points. The empty cell ( = 0) is obviously a non-boundary cell; in Figure 2, a 9 × 9 sub-gird consisting of the cell and eight neighboring cells is used to determine if the cell is a boundary cell or not when = 1, which is formulated as:

Estimation of the Tunnel Boundary Lines in the X-Y Plane
An entire tunnel is generally a long tube, so it can be scanned along its centerline station by station. The tunnel point clouds are registered in a user-specified coordinate system using sphere reference targets. In this system, the origin is near the tunnel entrance, and the Y axis is oriented along the direction of the initial segment of the tunnel. Several algorithms [22][23][24] have been developed to quickly extract features from LIDAR point clouds based on projection and gridding. Projecting the scanned data onto the X-Y plane can simplify the 3D tube to a long and narrow 2D object, from which the two boundary points' groups are extracted from the both sides of the 2D object. In order to improve the speed of extraction, an algorithm for extracting the boundary point groups is proposed using a fixed grid.
The projections of the entire tunnel point clouds in the X-Y plane are discretized using a square grid. A grid size that is too large or too small will decrease the computational efficiency or the extraction precision of boundary cells, respectively. The appropriate size of the grid is about one-twentieth of the width of the tunnel. The value of N ij is used to determine whether points exist or not in the cell ij, N ij has a value of 1 if points exist and 0 if there are no points. The empty cell (N ij = 0) is obviously a non-boundary cell; in Figure 2, a 9 × 9 sub-gird consisting of the cell ij and eight neighboring cells is used to determine if the cell ij is a boundary cell or not when N ij = 1, which is formulated as:  The center points of non-boundary cells will be used instead of all points in them for the further extraction of boundary points ( Figure 3). As shown in Figure 4, the point of interest is an arbitrary point in a boundary cell. We use its eight neighbor cells for avoiding the incorrect extraction of the pseudo-boundary points near the bounding rectangle . An angle criterion is proposed based on the distribution of neighboring points (the center points of non-boundary cells and all points in the boundary cells) of point in the 9 × 9 sub-gird. Cartesian coordinates of the neighbor points are converted to polar coordinates, with the pole set at point and polar axis L oriented along the positive direction of the X axis. The angular coordinates (e.g., 1 ) of all neighboring points are sorted by value, and then the differences (e.g., ∆α −1, ) between two consecutive neighboring coordinates are computed. Point is a boundary points, if the maximum difference exceeds a pre-specified threshold ( ), and a non-boundary point, otherwise. This angle threshold ( ) is set to 175°, and can work well even in a curved tunnel, because the radius of curvature of the curve segment is large enough (generally greater than 200 m) to consider that the curved boundary in the area of 9 × 9 sub-gird is nearly straight. The center points of non-boundary cells will be used instead of all points in them for the further extraction of boundary points ( Figure 3).  The center points of non-boundary cells will be used instead of all points in them for the further extraction of boundary points ( Figure 3). As shown in Figure 4, the point of interest is an arbitrary point in a boundary cell. We use its eight neighbor cells for avoiding the incorrect extraction of the pseudo-boundary points near the bounding rectangle . An angle criterion is proposed based on the distribution of neighboring points (the center points of non-boundary cells and all points in the boundary cells) of point in the 9 × 9 sub-gird. Cartesian coordinates of the neighbor points are converted to polar coordinates, with the pole set at point and polar axis L oriented along the positive direction of the X axis. The angular coordinates (e.g., 1 ) of all neighboring points are sorted by value, and then the differences (e.g., ∆α −1, ) between two consecutive neighboring coordinates are computed. Point is a boundary points, if the maximum difference exceeds a pre-specified threshold ( ), and a non-boundary point, otherwise. This angle threshold ( ) is set to 175°, and can work well even in a curved tunnel, because the radius of curvature of the curve segment is large enough (generally greater than 200 m) to consider that the curved boundary in the area of 9 × 9 sub-gird is nearly straight. As shown in Figure 4, the point of interest P is an arbitrary point in a boundary cell. We use its eight neighbor cells for avoiding the incorrect extraction of the pseudo-boundary points near the bounding rectangle ABCDE. An angle criterion is proposed based on the distribution of neighboring points (the center points of non-boundary cells and all points in the boundary cells) of point P in the 9 × 9 sub-gird. Cartesian coordinates of the neighbor points are converted to polar coordinates, with the pole set at point P and polar axis L oriented along the positive direction of the X axis. The angular coordinates (e.g., α 1 ) of all neighboring points are sorted by value, and then the differences (e.g., ∆α i−1,i ) between two consecutive neighboring coordinates are computed. Point P is a boundary points, if the maximum difference exceeds a pre-specified threshold (T), and a non-boundary point, otherwise. This angle threshold (T) is set to 175 • , and can work well even in a curved tunnel, because the radius of curvature of the curve segment is large enough (generally greater than 200 m) to consider that the curved boundary in the area of 9 × 9 sub-gird is nearly straight.  Since the tunnel boundary line can be used to determine the initial directions of cross-sectional planes, a cubic polynomial function is chosen to smooth and represent tunnel boundary points, which is parameterized as follows: where b , b , b , and b are the parameters of a boundary line. The points belonging to the measurement errors or different structures (i.e., boundary points of refuge recesses) will also be extracted if they meet the angle criterion. The RANSAC algorithm was first proposed by Fischler and Bolles [25] and it is an iterative method to estimate parameters of a mathematical model from a set of observed data which contains outliers. Hence, we adopt this algorithm to find the real boundary points and estimate the parameters of Equation (2).

Extraction of Cross-Sections
After estimation of the two tunnel boundary lines, we define the initial direction of cross-sectional planes to be orthogonal to one of the two boundary lines. Due to the construction and measurement errors and rough concrete linings, even a fine estimation of tunnel boundary lines cannot ensure that the plane orthogonal to it is the real cross-sectional plane. In order to find the real cross-sectional plane, we make adjustments twice to the direction of the initial cross-sectional plane based on the estimations of the local centerline (line in Figure 5) and upper boundary line (line in Figure 6).  Since the tunnel boundary line can be used to determine the initial directions of cross-sectional planes, a cubic polynomial function is chosen to smooth and represent tunnel boundary points, which is parameterized as follows: where a b , b b , c b , and d b are the parameters of a boundary line. The points belonging to the measurement errors or different structures (i.e., boundary points of refuge recesses) will also be extracted if they meet the angle criterion. The RANSAC algorithm was first proposed by Fischler and Bolles [25] and it is an iterative method to estimate parameters of a mathematical model from a set of observed data which contains outliers. Hence, we adopt this algorithm to find the real boundary points and estimate the parameters of Equation (2).

Extraction of Cross-Sections
After estimation of the two tunnel boundary lines, we define the initial direction of cross-sectional planes to be orthogonal to one of the two boundary lines. Due to the construction and measurement errors and rough concrete linings, even a fine estimation of tunnel boundary lines cannot ensure that the plane orthogonal to it is the real cross-sectional plane. In order to find the real cross-sectional plane, we make adjustments twice to the direction of the initial cross-sectional plane based on the estimations of the local centerline (line l c in Figure 5) and upper boundary line (line l u in Figure 6).  Since the tunnel boundary line can be used to determine the initial directions of cross-sectional planes, a cubic polynomial function is chosen to smooth and represent tunnel boundary points, which is parameterized as follows: , and b are the parameters of a boundary line. The points belonging to the measurement errors or different structures (i.e., boundary points of refuge recesses) will also be extracted if they meet the angle criterion. The RANSAC algorithm was first proposed by Fischler and Bolles [25] and it is an iterative method to estimate parameters of a mathematical model from a set of observed data which contains outliers. Hence, we adopt this algorithm to find the real boundary points and estimate the parameters of Equation (2).

Extraction of Cross-Sections
After estimation of the two tunnel boundary lines, we define the initial direction of cross-sectional planes to be orthogonal to one of the two boundary lines. Due to the construction and measurement errors and rough concrete linings, even a fine estimation of tunnel boundary lines cannot ensure that the plane orthogonal to it is the real cross-sectional plane. In order to find the real cross-sectional plane, we make adjustments twice to the direction of the initial cross-sectional plane based on the estimations of the local centerline (line in Figure 5) and upper boundary line (line in Figure 6).  As shown in Figure 5, point S n (x S n , y S n ) is selected along one boundary line that was estimated in Section 2.1, depending on the location of the cross-section of interest. The initial cross-sectional plane l is a vertical plane determined from point S n (x S n , y S n ) orthogonal to the boundary line. In the first adjustment, the subsets G 1 and G 2 (black points) of the two boundary point groups in the X-Y plane are extracted from between the two planes that are parallel to plane l at a distance d. As illustrated in Section 2.1, the two boundary lines of the tunnel can be considered as two straight lines on a small scale, so the value of d should be small, but at the same time it must be large enough to provide enough data for the estimation of the centerline l c (the intersection line of the real cross-sectional plane and X-Y plane should be orthogonal to this centerline). We propose an algorithm to directly estimate the centerline l c from the subsets G 1 and G 2 with the total least squares method [26]. The least squares and total least squares methods assess the fitting accuracy in different ways [27]: the least squares method minimizes the sum of the squared vertical distances from the data points to the fitting line, while the total least squares method minimizes the sum of the squared orthogonal distances from the data points to the fitting line.
Since two local boundary lines l 1 and l 2 are parallel and have the same distance from the centerline l c , the boundary lines are formulated as: where a and b are the parameters of centerline l c , and c is the parameter of boundary lines l 1 and l 2 . By using subsets G 1 and G 2 , the constraint equation for fitting the local centerline l c is derived from Equation (3): where: where n 1 and n 2 , respectively, denote the numbers of points in subsets G 1 and G 2 .
We define U ∑ V T to be the singular value decomposition of the augmented matrix [B L], where ∑ = diag (σ 1 , σ 2 , σ 3 , σ 4 ) and σ 1 > σ 2 > σ 3 > σ 4 . Since B is a full rank matrix (all points in subsets G 1 and G 2 are different), the total least squares approximate solutionX for X is given by: After the estimation of the centerline l c , plane l is adjusted to be plane l whose direction from point S n (x S n , y S n ) is orthogonal to centerline l c , so the cross-sectional plane l is formulated as: In the second adjustment, a point group G p is extracted from between the two planes that are parallel to plane l at a distance d (see Figure 5), and then projected onto the vertical plane traversing centerline l c (see Figure 6). The upper boundary points are extracted using the method proposed in Section 2.1. To find the angle θ of the rotation, the upper boundary line l u is fitted by using those boundary points with the total least squares method. Vector is the normal vector of plane ′ , and vector is the unit vector whose direction from point is along the projection of plane ′ onto the X-Y plane. According to Rodrigues' rotation formula [28], vector is rotated by an angle of about the axis in the direction of , from which we can obtain the normal vector ′ of plane ′′ : with = (1, , 0) and = ( √1+ 2 , − 1 √1+ 2 , 0). Point is given a height value by the average height of subsets 1 and 2 . Combining point and vector ′ , the final cross-sectional plane ′′ is represented by: The final cross-section is extracted using the projection of the nearby points onto plane ′′ , where the nearby points whose orthogonal distances to plane ′′ are less than ′ /2 are extracted from the raw tunnel point clouds.

A Filtering Algorithm for Non-Lining Points Removal
Many morphological filtering methods that are frequently employed in signal processing, image analysis, and bare-earth extraction fields [29][30][31], are applicable for noise removal. It is difficult to find a simple function to approximate a non-circular lining cross-section, especially when a tunnel lining has deformed. Similar to the idea of the morphological erosion, we propose an angle-based filter for removing non-lining points without the limit of the shape of tunnel lining cross-section.
As shown in Figure 7, the theoretical cross-section of the tunnel consists of three tangent circles, and it can be absolutely positioned onto the final cross-sectional plane ′′ by using the central axis (the intersection of plane ′′ and the vertical plane traversing centerline ) of the extracted cross-section and vertex (the intersection of plane ′′ and line are illustrated in Figure 6). After the location of the theoretical cross-section, the operator of our filter is defined as an angle of degrees ( < 180°). The vertex of the angle is positioned at each cross-sectional point , and the two sides of the angle make an angle ∕ 2 with the positive axis, where the axis is taken to be along the normal line to the theoretical cross-section at point , and the axis is drawn through Vector v is the normal vector of plane l , and vector u is the unit vector whose direction from point S n is along the projection of plane l onto the X-Y plane. According to Rodrigues' rotation formula [28], vector v is rotated by an angle of θ about the axis in the direction of u, from which we can obtain the normal vector v of plane l : with v = (1, a, 0) and u = a √ 1+a 2 , − 1 √ 1+a 2 , 0 . Point S n is given a height value z S n by the average height of subsets G 1 and G 2 . Combining point S n and vector v , the final cross-sectional plane l is represented by: The final cross-section is extracted using the projection of the nearby points onto plane l , where the nearby points whose orthogonal distances to plane l are less than d /2 are extracted from the raw tunnel point clouds.

A Filtering Algorithm for Non-Lining Points Removal
Many morphological filtering methods that are frequently employed in signal processing, image analysis, and bare-earth extraction fields [29][30][31], are applicable for noise removal. It is difficult to find a simple function to approximate a non-circular lining cross-section, especially when a tunnel lining has deformed. Similar to the idea of the morphological erosion, we propose an angle-based filter for removing non-lining points without the limit of the shape of tunnel lining cross-section.
As shown in Figure 7, the theoretical cross-section of the tunnel consists of three tangent circles, and it can be absolutely positioned onto the final cross-sectional plane l by using the central axis (the intersection of plane l and the vertical plane traversing centerline l c ) of the extracted cross-section and vertex p u (the intersection of plane l and line l u are illustrated in Figure 6). After the location of the theoretical cross-section, the operator of our filter is defined as an angle of α degrees (α < 180 • ). The vertex of the angle is positioned at each cross-sectional point p i , and the two sides of the angle make an angle α/2 with the positive y i axis, where the y i axis is taken to be along the normal line to the theoretical cross-section at point p i , and the x i axis is drawn through p i perpendicular to the y i axis. The designed filter searches for cross-sectional points inside it. Since the surface of the tunnel lining is normally rough, we add a confidence interval d v (point p i is shifted a distance of d v along the positive y i axis to be point p v i ) for our filter. A point p i (p v i ) at the vertex of the angle is accepted as a lining point if there are no other cross-sectional points p j inside the angle. The filter function for defining the set Lp of lining points is mathematically represented as follows: where CSp is the set of all cross-sectional points. perpendicular to the axis. The designed filter searches for cross-sectional points inside it. Since the surface of the tunnel lining is normally rough, we add a confidence interval (point is shifted a distance of along the positive axis to be point ) for our filter. A point ( ) at the vertex of the angle is accepted as a lining point if there are no other cross-sectional points inside the angle. The filter function for defining the set of lining points is mathematically represented as follows: where is the set of all cross-sectional points.

Data Acquisition
The proposed method was tested in a double-track railway tunnel with the length of 619 m in Sichuan, China. As shown in Figure 8, the point cloud dataset was captured by the Faro X130 terrestrial laser scanner (Lake Mary, FL, United States) with 31 scans, and the distance of adjacent scans was about 20 m. Three sphere reference targets were laid between the two adjacent scans, which ensures the stability of the scanning position registration. All scans were registered together in a userspecified coordinate system using Faro Scene software (Lake Mary, FL, United States). The details of the point cloud dataset are listed in Table 1. MATLAB (Natick, MA, USA) was used to implement the data processing and analysis, as well as the visual representation in the following subsections.

Data Acquisition
The proposed method was tested in a double-track railway tunnel with the length of 619 m in Sichuan, China. As shown in Figure 8, the point cloud dataset was captured by the Faro X130 terrestrial laser scanner (Lake Mary, FL, USA) with 31 scans, and the distance of adjacent scans was about 20 m. Three sphere reference targets were laid between the two adjacent scans, which ensures the stability of the scanning position registration. All scans were registered together in a user-specified coordinate system using Faro Scene software (Lake Mary, FL, USA). The details of the point cloud dataset are listed in Table 1. MATLAB (Natick, MA, USA) was used to implement the data processing and analysis, as well as the visual representation in the following subsections.

Detection of Tunnel Boundary Points in the X-Y Plane
As described in Section 2.1, the tunnel point cloud dataset was projected onto the X-Y plane and then discretized to improve the speed of extraction using a 0.5 m resolution grid (the width of this tunnel is about 10 m). As shown in Figure 9a, the dark dots are the extracted boundary points, but some of them are the boundary points of refuge recesses and will affect the determination of the final cross-sectional planes. The polynomial fitting was used for smoothing of tunnel boundary lines and the elimination of outliers, and the RMSE of the discrepancies is 46.5 mm. It is possible to eliminate these points of refuge recesses for a better extraction of tunnel boundary points (Figure 9b), because the heights of them are lower than the nearby lining-boundary points and the projections of them in the X-Y plane are outside the fitted boundary lines.

Detection of Tunnel Boundary Points in the X-Y Plane
As described in Section 2.1, the tunnel point cloud dataset was projected onto the X-Y plane and then discretized to improve the speed of extraction using a 0.5 m resolution grid (the width of this tunnel is about 10 m). As shown in Figure 9a, the dark dots are the extracted boundary points, but some of them are the boundary points of refuge recesses and will affect the determination of the final cross-sectional planes. The polynomial fitting was used for smoothing of tunnel boundary lines and the elimination of outliers, and the RMSE of the discrepancies is 46.5 mm. It is possible to eliminate these points of refuge recesses for a better extraction of tunnel boundary points (Figure 9b), because the heights of them are lower than the nearby lining-boundary points and the projections of them in the X-Y plane are outside the fitted boundary lines.

Extracting Results
One of the two fitted boundary lines was used to determine the initial cross-sectional planes. Since the discrepancies of the fitted boundary line are slightly large, the initial directions are inaccurate. To optimize the directions of cross-sectional planes, the proposed method in Section 2.2.1 was implemented to obtain the final directions by using = 40 cm. Based on the final cross-sectional planes, and by using ′ = 5 mm, ten cross-sections were extracted to test the accuracy of our method and shown in Figure 10.

Extracting Results
One of the two fitted boundary lines was used to determine the initial cross-sectional planes. Since the discrepancies of the fitted boundary line are slightly large, the initial directions are inaccurate. To optimize the directions of cross-sectional planes, the proposed method in Section 2.2.1 was implemented to obtain the final directions by using d = 40 cm. Based on the final cross-sectional planes, and by using d = 5 mm, ten cross-sections were extracted to test the accuracy of our method and shown in Figure 10.

Extracting Results
One of the two fitted boundary lines was used to determine the initial cross-sectional planes. Since the discrepancies of the fitted boundary line are slightly large, the initial directions are inaccurate. To optimize the directions of cross-sectional planes, the proposed method in Section 2.2.1 was implemented to obtain the final directions by using = 40 cm. Based on the final cross-sectional planes, and by using ′ = 5 mm, ten cross-sections were extracted to test the accuracy of our method and shown in Figure 10. Figure 10. Extraction of ten cross-sections. Figure 10. Extraction of ten cross-sections.

Assessment of Extracting Accuracy
The proposed method was compared to the method of cross-sectional estimation as specified by Han et al. [12] since Han's method has achieved a high accuracy in comparison to total station surveying. As shown in Figure 11, two comparisons were made for the cross-sections extracted from the same location: (1) the clear width W c of the cross-section; and (2) the height H vr from the vertex of cross-section to the top of the inner rail. As shown in Figure 6, the theoretical discrepancy (∆h) between the height (h) of our cross-section and the height (h ) of Han's (cross-sections are extracted using a vertical plane) can be represented as follows:

Assessment of Extracting Accuracy
The proposed method was compared to the method of cross-sectional estimation as specified by Han et al. [12] since Han's method has achieved a high accuracy in comparison to total station surveying. As shown in Figure 11, two comparisons were made for the cross-sections extracted from the same location: (1) the clear width of the cross-section; and (2) the height from the vertex of cross-section to the top of the inner rail. As shown in Figure 6, the theoretical discrepancy (∆ℎ) between the height (ℎ) of our cross-section and the height (ℎ ′ ) of Han's (cross-sections are extracted using a vertical plane) can be represented as follows:  Table 2). The clear widths of cross-sections extracted by the proposed method are very close to Han's, and the height discrepancies are very close to the theoretical discrepancy (∆ℎ). Specifically, the height discrepancy will increase with an increase in tunnel grade. The discussions above indicates that our method is able to offer high accuracy and more reliable tunnel cross-sections. In another word, the horizontal coordinates of cross-sectional points achieve high accuracy, and the vertical coordinates are more reliable.  Let h = 8.7 m (the theoretical value of H vr ) and tanθ = 0.021 (ten cross-sections located where the theoretical tunnel grade is 21‰) then ∆h = −1.9 mm. The average discrepancies of W c and H vr are −0.4 mm and −2 mm, respectively, and the RMSEs are 0.8 mm and 2.1 mm ( Table 2). The clear widths of cross-sections extracted by the proposed method are very close to Han's, and the height discrepancies are very close to the theoretical discrepancy (∆h). Specifically, the height discrepancy will increase with an increase in tunnel grade. The discussions above indicates that our method is able to offer high accuracy and more reliable tunnel cross-sections. In another word, the horizontal coordinates of cross-sectional points achieve high accuracy, and the vertical coordinates are more reliable.  Figure 10 shows that there are a lot of non-lining points in the cross-sections. These non-lining points belong to pipes and catenary equipment ( Figure 12). To eliminate non-lining points' interference in the safety assessment of tunnel linings, the filtering algorithm mentioned in Section 2.2.2 was used for the removal of non-lining points. The angle α of the operator must be large enough to ensure that almost all non-lining points are removed while, at the same time, the confidence interval d v cannot be set to a value that is too small; otherwise the operator will lose a lot of real lining points. The example in Figure 13a shows that 29.3% of lining points were removed from the cross-section (ID = 1) by using α = 165 • and d v = 0. Hence, the confidence interval d v was set to 1 cm, which is slightly larger than the undulation of this tunnel surface, and only 0.275% of lining points were removed (Figure 13b).  Figure 10 shows that there are a lot of non-lining points in the cross-sections. These non-lining points belong to pipes and catenary equipment ( Figure 12). To eliminate non-lining points' interference in the safety assessment of tunnel linings, the filtering algorithm mentioned in Section 2.2.2 was used for the removal of non-lining points. The angle of the operator must be large enough to ensure that almost all non-lining points are removed while, at the same time, the confidence interval cannot be set to a value that is too small; otherwise the operator will lose a lot of real lining points. The example in Figure 13a shows that 29.3% of lining points were removed from the cross-section (ID = 1) by using = 165° and = 0. Hence, the confidence interval was set to 1 cm, which is slightly larger than the undulation of this tunnel surface, and only 0.275% of lining points were removed (Figure 13b).     Figure 10 shows that there are a lot of non-lining points in the cross-sections. These non-lining points belong to pipes and catenary equipment ( Figure 12). To eliminate non-lining points' interference in the safety assessment of tunnel linings, the filtering algorithm mentioned in Section 2.2.2 was used for the removal of non-lining points. The angle of the operator must be large enough to ensure that almost all non-lining points are removed while, at the same time, the confidence interval cannot be set to a value that is too small; otherwise the operator will lose a lot of real lining points. The example in Figure 13a shows that 29.3% of lining points were removed from the cross-section (ID = 1) by using = 165° and = 0. Hence, the confidence interval was set to 1 cm, which is slightly larger than the undulation of this tunnel surface, and only 0.275% of lining points were removed (Figure 13b).  To assess the performance of the angle-based filter with different angles, the quantitative evaluations of the ten cross-sections resulted in Table 3

Conclusions
We presented an automated and effective method for extracting tunnel lining cross-sections from terrestrial laser scanning (TLS) point clouds. This can be applicable in a tunnel with an arbitrarily-shaped lining cross-section. In our method, the tunnel point cloud dataset was projected onto the X-Y plane to extract the boundary points of both sides. By using these tunnel boundary points, the initial direction of the cross-sectional plane was determined, and then adjusted with the total least squares method. The cross-sectional plane, using Rodrigues' rotation formula, was adjusted again for capturing the final cross-sectional points. To generate the bare-lining cross-section, an angle-based filter algorithm was developed for removing non-lining points based on the morphological erosion.
The proposed method was validated on the point cloud dataset of a real railway tunnel. The results of cross-sectional extraction were compared with an existing method, which showed that the clear widths of cross-sections achieved high accuracy (RMSE of 0.8 mm) and the cross-sectional heights were more reliable. The results of no-lining point removal indicated that the proposed filter was able to offer a good classification for cross-sectional points. The performance of the filter will deteriorate with the outermost outliers of the tunnel point clouds increasing; how to reduce these outliers is among our planned future work.