AUTOMATIC MONITORING OF TUNNEL DEFORMATION BASED ON HIGH DENSITY POINT CLOUDS DATA

An automated method for tunnel deformation monitoring using high density point clouds data is presented. Firstly, the 3D point clouds data are converted to two dimensional surface by projection on the XOY plane, the projection point set of central axis on XOY plane named Uxoy is calculated by combining the Alpha Shape algorithm with RANSAC ( Random Sampling Consistency ) algorithm, and then the projection point set of central axis on YOZ plane named Uyoz is obtained by highest and lowest points which are extracted by intersecting straight lines that through each point of Uxoy and perpendicular to the two -dimensional surface with the tunnel point clouds, Uxoy and Uyoz together form the 3D center axis finally. Secondly, the buffer of each cross section is calculated by K-Nearest neighbor algorithm, and the initial cross-sectional point set is quickly constructed by projection method. Finally, the cross sections are denoised and the section lines are fitted using the method of iterative ellipse fitting. In order to improve the accuracy of the cross section, a fine adjustment method is proposed to rotate the initial sectional plane around the intercept point in the horizontal and vertical direction within the buffer. The proposed method is used in Shanghai subway tunnel, and the deformation of each section in the direction of 0 to 360 degrees is calculated. The result shows that the cross sections becomes flat circles from regular circles due to the great pressure at the top of the tunnel.


INTRODUCTION
In recent years, with the rapid expansion of the city, the subway has become one of the major way of transport, the scale of subway project construction is expanding.Carrying out precise, fast and automated deformation monitoring of the subway, especially the subway in operation is so necessary.Although the accuracy of traditional deformation monitoring method (such as total station, convergence meter and so on) for subway tunnel is based on point cloud data (Cheng et al., 2015;Ju et al., 2014;Li et al., 2015;Liu et al., 2013;Lu et al., 2016;Tuo et al., 2013;Wang et al., 2013;Xu et al., 2016).Based on translated threedimensional point cloud into a binary image, Han (Han et al., 2013) extracted the central axis of the tunnel using skeleton method, calculated and adjusted the cross section.Finally, the coordinates of the stations obtained by their method are compared with the coordinates of the stations measured by the total station, which shows the advantages of the method in terms of speed and precision.Through projected the three-dimensional point cloud on two 2D planes, Kang (Kang et al., 2014) calculated the central axis of the tunnel using three curve fitting models and used a global extraction algorithm to minimize the deviations in the overlap segments, which improved the accuracy of the tunnel axis.Then the method of Han is used to obtain the cross section plane, and the local least squares fitting method is used to search the cross section points of each section in the 360 degree direction.Cheng (Cheng et al., 2016) improved the speed of boundary line extraction based on the method of Kang, calculated the initial cross-sectional plane based on one of the two boundary lines, and obtained the final cross-sectional plane by adjusting the initial cross-sectional plane twice.Finally, they used a morphological denoising algorithm to remove the non-lining points in the cross-section.Comparing with Han's method, their method shows the high accuracy, reliability and versatility of the tunnel in different shapes.
This paper studied the method of automatic extraction of tunnel cross section based on the point clouds obtained by mobile laser scanning technology (MMS).Firstly, taking into account the MMS can obtain absolute coordinates of point cloud directly, it is usually need to rotate the tunnel point clouds and make it parallel to a certain axis before projecting them onto 2D plane (Kang et al., 2014), the process is eliminated in our method and the amount of computation is reduced.Secondly, the construction speed of cross sectional point set is increased by constructing the buffer of each section in the method of k-nearest neighbor.Thirdly, a denoising method of iterative ellipse fitting is proposed, combining the fitting of the sectional curve line with the denoising of the sectional point set.Thirdly, in order to reduce the error, we make the initial cross sectional planes rotate in the horizontal and vertical directions at a certain angular interval around their respective intercept points by iteration and choose the optimal cross section at every intersect point as the final result.Finally, the method is validated by the data obtained by reverse measurement for the same tunnel.The specific process of the method is as follows.

CALCULATION OF CENTRAL AXIS
The central axis is a spatial curve line which can expresses the posture and trend of the tunnel.The cross sectional plane orthogonal to tunnel can be calculated through the tangent vector of each point at central axis.Because of the complicated calculation of spatial curve line, the most commonly used method is to reduce the dimension, that is, project threedimensional point cloud data to two-dimensional plane or convert it into a two-dimensional image to extract the central axis.
As the MMS can obtain absolute coordinates of point cloud directly and the amount of data is very large.In order to avoid the process of rotating point clouds to make it parallel to a certain axis and then project them to two-dimensional plane in traditional method.The extraction of central axis is divided into two parts in our method, namely extraction of central axis points of XOY plane (Uxoy) and extraction of central axis points of YOZ plane (Uyoz).When the point sets Uxoy and Uyoz are all calculated, quadratic curve equation is used to fit the two sets respectively and central axis of tunnel can be obtained using equation 1.

Extraction of central axis points on XOY plane
Firstly, the 3D point cloud data is projected onto XOY plane and the boundary point set is detected through Alpha shape algorithm.
And then, the method of Kang (Kang et al., 2014) is used to calculate the projection point set of center axis in XOY plane named Uxoy.

Boundary point set detection:
The alpha shape algorithm (Edelsbrunner et al., 1983) As the tunnel lining contains bolt holes and other appendages, leading to the unused points which named as noise in original point cloud data.Therefore, we use a random sampling consistency (RANSAC) algorithm (Fischler et al, 1981) with good anti-noise effect to estimate the parameters of equation.
In order to improve the accuracy of the central axis, a method of Kang (Kang et al., 2014)     and  2 , and the K-nearest neighbor method (Wu et al., 2007) is used to build the buffer of the cross section by searching point clouds for each point in the point set  1 using a certain number n (the size of the n is determined experimentally to ensure that the radius of buffer is greater than or equal to ε).The buffer will increase with the increment of ε, and the corresponding calculation is also greater.Therefore, we should select a appropriate ε value based on the actual data in order to take into account the and accuracy simultaneously.
Figure 3 Construction of section buffer

Denoising of sectional points and fitting of sectional lines
The presence of the ancillary facilities in tunnel resulted in the presence of noise in tunnel sectional points, so it is necessary to remove unused points.There are a lot of mature filtering algorithm in the processing of airborne point cloud data (Dufour et al, 2013;Vosselman, 2000), and vehicle point cloud filter method is usually do some improvements on the basis of airborne point cloud filter algorithm.But many airborne filtering algorithms do not suitable for the tunnel for its special shape, and the tunnel cross section could be considered as an ellipse which is close to the standard circle after deformation (Walton et al, 2014).So an iterative elliptical fitting method is proposed to filter the noise in cross sectional points.
Because of the complicate of 3D curve line fitting, it is necessary to rotate the normal of each cross section and make it parallel to y-axis, and then fit the cross sectional curve line according to the X coordinates and Z coordinates.In the filtering method, the noise can be removed according to the distance from the points to the fitted ellipse, and the RANSAC algorithm is used to estimate the parameters of ellipse fitting.The specific filtering process is as follows: First, the initial ellipse of cross-sectional points with noise is fitted using RANSAC algorithm, and the initial parameters (Including the coefficients of elliptic equation, center coordinates, long axis, short axis, etc.) of ellipse are calculated.
Second, the shortest distance   from each cross-sectional points to initial ellipse are calculated and form the distance point Finally, the mean and standard deviation of d are calculated according to the following equation: If |  −   | > 2, the cross section points are considered to be the noise and removed, meanwhile, using the denoised cross sectional points of first iteration to fit the cross sectional line in the second iteration until all the noise is filtered.

Adjustment of the cross section
The cross sections are extracted on the basis of central axis in our method.However, the central axis still exists fitting error even with a very fine method in extraction process, resulting in the error of initial cross section point set.In addition, the error in data acquisition process also impact the accuracy of initial cross section.Therefore, a fine adjustment method is proposed to make the final cross sections closer to the theoretical one, it make each initial cross section plane rotate at a certain angle in the horizontal and vertical directions within the buffer zone, as shown in figure 4, the black circle is assumed to be the initial cross section plane, the red circle is the result of initial cross sectional plane rotation around the intercept point in horizontal plane and the blue circle is the result of the initial cross sectional plane rotation around the intercept point in vertical plane.Then the method of 3.1 is used to obtain the sectional points set and the method of 3.2 is used to filter the noise and fit the section curve line after each rotation process.Finally, it is considered that the cross section which has the minimum area at every

Result of cross sections:
Based on the central axis, the initial 18 cross sectional planes are calculated in the intervals of 2 meters and the initial intercept point is set at 1 m from the tunnel port.In order to increase the density of cross sectional points, the method in Section 3.1 is used to calculate the initial cross sectional point sets.The buffer radius threshold ε in this paper is 10000 (figure 6) and the points within the distance range of ±5 mm (range accuracy of the data acquisition ) are projected onto each cross sectional plane to construct 18 initial cross sectional point sets finally (figure 7).direction within the scope of -2 0 and +2 0 at the interval of 0.1 0 in our experiment, that is to say, every cross sectional plane is rotated for 82 times in its respective buffer and 82 cross sectional point sets at each intercept point of the tunnel will be constructed.
Figure 8 is the comparison of initial cross sectional point sets (black dots) and the rotated cross sectional point sets (red dots) when rotation of the cross sectional plane are +2 0 .
Figure 8 Comparison of cross sectional point sets The area of initial cross sectional point sets and optimal cross sectional point sets are calculated using the parameters of each sectional curve line, the comparison result is in table 2. From the result we can see that the area of the optimal sections are all less than or equal to the area of initial sections, and the corresponding rotation angle of the optimal section is in the range of -1 0 to +1 0 , most of them are in the range of -0.5 0 to 0.5 0 .reverse measurement data and adjust the point set iteratively using the method in section 3.3.The optimal sectional curve line is compared with the optimal curve line calculated by direct measurement data at the same intercept point in section 4.2.2 and the proximity of the two curve lines is expressed by sectional points which are interpolated in the degree of 0 to 360 in their planes.As shown in figure 9, the differences between the two cross sections range from -1mm to 2mm and the closer to the horizontal and vertical direction, the greater the differences.The mean and the RMSE (Mean square error) of the differences are as small as 0.398mm and 1.072mm respectively, note that although the distributions of point cloud are different, but the cross-sectional results calculated at the same location using our method are very close.As in table 3, there are still some noises on the lining of tunnel have not been filtered when the number of iteration n is equal to 8, but the noises are better filtered and the sectional curve line is better fitted when the number of iteration n exceeds 10.
Experiments show that the number of cross sectional point set is basically stable when the number of iterations exceeds 10 and almost no noise can be delete.with minimum areas will be selected as the final ones. Figure 10 is the denoising result of the final cross sectional point sets.
Figure 10 Denoisig result of the final sectional point sets

Deformation Analysis
In order to analyze the deformation of the tunnel, the amount of deformation of each section in each direction is obtained and comparing with 2.75m which is the standard radius of tunnel lining.The deformation curve line of 18 cross sections in the 360 degree direction is drawn in figure 11, it shows that deformation of tunnel is distribute in the range of -21mm to +13mm, and the maximum deformation of each section is distributed near the degree of 0 (the right side of the section), 90 (the top of the section) and 180 (the left side of the section).From the whole of the curve we can see that the negative value of the deformation is more than the positive value significantly and the reduction is greater than the increase in the amount of deformation.These information indicating that vertical shrinkage deformation of the tunnel is greater than horizontal expansion deformation, the tunnel sections are deformed into flat circles from regular circles.axis through the intercept points can be calculated and the cross sectional point sets can be obtained by projected the points near the cross sectional plane within a threshold range onto the sectional plane.In order to increase the obtaining speed, the tunnel point cloud is blocked by calculating the buffer of each section before projection.In addition, an iterative denoising method based on RANSAC algorithm is proposed to remove the noise from the cross point cloud and fitted the section curve line.
In order to increase the accuracy of the cross section, we adopt a fine adjustment method to rotate the initial cross sectional plane in horizontal and vertical direction and choose the section with the minimum area as the optimal one.The reverse measurement data of the same tunnel is used to verify our method.Finally, the method in this paper is applied to the real tunnel and analyzed its deformation.Our method can be used for point cloud data obtained through exchanging station and later stitching, or obtained directly by moving laser scanning technology.However, there is still need for improvement in the method, such as the denoising process of the point cloud may remove some minor deformation points, and in the fine adjustment method of the cross-section, we use the area of the sectional curve line which obtained by fitting method to measure the size of the cross section may lead to some errors and we will study how to improve his accuracy in the future research.
is used to interpolate the projection point set of central axis on XOY surface.As shown in figure 1, L1 and L2 are the boundary curve lines of tunnel on XOY plane which are fitted by the above method.Firstly, L1 is evenly resampled.Select a random point   from point set L1 and the tangent vector   of curve line L1 is calculated at   , the perpendicular line of   is intersect to another boundary curve line L2 at point   ′ ; Secondly, the tangent vector   ′ of curve line L2 is calculated at   ′ , the perpendicular line of   ′ is intersect to the boundary curve line L1 at point   ′′ .Because of the fitting error in the pre-calculation process,   ′ and   ′′ must not coincide.So the final central axis point M1 is the average of points   ′ and   ′′ , which are the midpoint of   and   ′ ,   ′ and   ′′ respectively.The other central axis points are also calculated using the same method and construct point set Uxoy finally.

Figure 1
Figure 1 Interpolate of central axis points

Figure 2
Figure 2 Detection of Uyoz In addition, it's can't guarantee to search for highest and lowest points in all directions because of the discrete of point cloud.To solve this problem, the point clouds of tunnel are projected onto the XOY plane first, and search for the points within the range of ∆ (∆ is the maximum density of the point cloud data) for each point in point set Uxoy, and then the highest and lowest points (Zmax and Zmin) can be chosen from these points and considered to be the boundary points of tunnel in YOZ plane.

ε
The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XLII-2/W7, 2017 ISPRS Geospatial Week 2017, 18-22 September 2017, Wuhan, China intercept point is the optimal ones.

Figure
Figure 4 Rotation of the initial cross sectional plane

Figure 5
Figure 5 Results of Uxoy and Uyoz

Figure 9
Figure 9 Differences of cross sectional lines Result of different number of iterationsFinally, the noise are filtered and the sectional curve line are fitted of every cross sectional sets at every rotation angle by iteration and the areas of the sections are calculated using the parameters of sectional curve line.The cross sectional point sets

Figure 11
Figure 11 Deformation curve lines

Figure 12
Figure 12 Oblateness of the cross sections

2.1.2 Extraction of Uxoy:
After obtaining the boundary points, it is necessary to fit boundary curve line, and then the projection point set of central axis on XOY surface can be interpolated based on the boundary curve line.A quadratic curve equation is used to fit the boundary points, which is expressed by formula 2: The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XLII-2/W7, 2017 ISPRS Geospatial Week 2017, 18-22 September 2017, Wuhan, China distances calculation.Since only the points that near each sectional plane are projected onto it, so the tunnel point clouds should first be blocked in order to improve the efficiency.

Table 1 ,
(X1, Y1, Z1) is the point on the direct central axis, (X2, Y2, Z2) is the point on the reverse central axis, and D is the distances between them.We can see that the distances are about 0.15cm, the mean of the distance is 0.16cm and RMS is 0.007cm.From 4.1 we can know that the average distance of registration is 0.1cm.Therefore it is considered that the central axes calculated by the two groups of data are basically consistent.