Automatic Processing Method for Deformation Monitoring of Circle Tunnels Based on 3 D LiDAR Data

The application of 3D LiDAR technology has become increasingly extensive in tunnel monitoring due to the large density and high accuracy of the acquired spatial data. The proposed processing method aims at circle tunnels and provides a clear workflow to automatically process raw point data and easily interpretable results to analyze tunnel health state. The proposed automatic processing method employs a series of algorithms to extract point cloud of a single tunnel segment without obvious noise from entire raw tunnel point cloud mainly by three steps: axis acquisition, segments extraction and denoising. Tunnel axis is extracted by fitting boundaries of the tunnel point cloud rejection in plane with RANSAC algorithm. With guidance of axis, the entire preprocessed tunnel point cloud is segmented by equal division to get a section of tunnel point cloud which corresponds to a single tunnel segment. Then the noise in every single point cloud segment is removed by clustering algorithm twice, based on distance and intensity. Finally, clean point clouds of tunnel segments are processed by effective deformation extraction processor to get ovality and three-dimensional deformation nephogram.


Introduction
Full automatic deformation monitoring technology based on LiDAR can monitor the structure in a non-contact way and is expected to be one of the most important directions in the field of SHM (structure health monitoring) in the future [1][2][3][4] [5] [6].LiDAR techniques can measure several million points in a short time with a high accuracy and has been used by many researchers for generating three-dimensional models in real-estate industry for visualization, in building renovation projects for surveying, and monitoring structure health [7][8] [9][10] [11][12] [13][14] [15][16] [17].
Tunnel monitoring based on 3D LiDAR data has developed rapidly but insufficient attention has been paid to practical engineering applications and the degree of automation in point cloud processing stage too has apparently been ignored.Present tunnel monitoring system based on LiDAR is hard to put into use because of its time-consuming point clouds processing system especially in long tunnels where the time cost in point cloud processing can be five to eight times more than data acquisition.There are few data processors designed specifically for circle tunnel point clouds data.By the end of 2017, the subway tunnel with circle profile in operation had reached 617 kilometers in Shanghai, the longest in the world and there are 218 kilometers of subway tunnels under construction or in the planning stages.The massive tunnels with higher departure frequency and longer operation time put forward higher requirements to the mode and efficiency of monitoring system.
The techniques presented in recent literature on point cloud data processing are as shown below.Timothy Nuttens et al. [18] observed the differences between average radius values during the first three months after construction of a circular tunnel structure.But the authors manually processed the data; it wasn't automatic.Jen-Yu Han et al. [19] generated tunnel profiles at multiple epochs and proposed the minimum-distance projection (MDP) algorithm to establish point correspondences after which deformation signals can be identified.But the authors still extracted the centerline of the tunnel manually.Before long, Han and his colleagues [20] improved the approach as real 3D approach and estimated the MDP algorithm using directly the 3D dispersed point clouds.But the method is hardly put into use due to the huge computational costs it involves.Rodriguez-Cuenca et al. [21] focused on extracting street boundaries and proposed a method to detect edges of roads automatically, which have some reference for axis extraction.
In the field of 3D modeling, Andrey Dimitrov et al. [22] present a region growing algorithm based on geometrical continuity for robust context-free segmentation of unordered point clouds.The thought of region growing can be used for reference.There is no doubt that building modeling techniques plays an important role in BIM [23][24] [25].P. Tang et al. [26] considered the impact of related factors, for instance scanning distance, density of data and incidence angle and then proposed a model to estimate edge loss in scanned data.
Recent developments about feature identification which have some significant reference for denoising from 3D LiDAR data are mainly related to road surface features and ground structures.As road surface features become more and more complicated and critical to traffic safety, Jenny Guo at el. [27] proposed a fully automatic approach for reconstructing road surface features based on the LiDAR technology.Jaselskis et al. [28] applied laser scanning to improve transportation projects by measuring the volume of soil and rock and determining road elevations.Cabo et al. [29] proposed an algorithm to detect planar or quasi-planar surfaces from point clouds by using line clouds and identifying surfaces by grouping lines.Besides, the mobile laser scanner arose from 2000s and made a significant contribution to the feature identification field [30][31] [32][33][34] [35] [36].
Intensity collected by scanner contains extra information and can expand the applied range of the LiDAR data.Alireza G. Kashani et al. [37] proposed a cluster-based method to detect wind-induced roof covering damage automatically.In conclusion, the Lidar intensity "I" led to the highest accuracy among clustering features tested in damage detection.A.J. Kashani et al. [38] also studied the application of intensity in the wind damage detection.The proposed processing method applies the intensity to remove the noise points inside the tunnel, which are reflected by different materials.
Our aim was to achieve full-automation in LiDAR data processing to significantly improve the efficiency compared to the conventional manual way which is obviously time-consuming.The proposed automatic processing method contains three steps: axis acquisition, segments extraction and denoising to extract clean point cloud of a single tunnel segment from the entire raw tunnel point cloud.Tunnel axis is extracted by fitting boundaries of the tunnel point cloud rejection in plane with RANSAC algorithm.In the segments extraction step, the axis of tunnel point cloud with known number of tunnel segments is equally divided to get normal planes.Then the point cloud between two adjacent normal planes is extracted.Denoising is implemented in every point cloud segment by clustering algorithm based on distance and intensity, separately.

Automatic data processing
Data processing extracts applicable data from raw point clouds scanned by TLS (Terrestrial Laser Scanner) in field work and provides valuable information for engineers to analyze the structure's health state.The proposed processing method is significantly different from the conventional manual way in terms of the automaticity.The new processing method employs a series of algorithms which are efficient and stable enough to finish steps from tunnel axis acquisition to denoising.Fig. 1 illustrates the flow chart of the proposed automatic data processing method.

Registration of point clouds
In general, a single scan station can cover only a finite range and is hardly qualified for field work of an entire structure like a tunnel [39][40].Thus the entire tunnel point cloud is usually registered by different point clouds of multiple scan stations.There are mainly four methodologies for registration of point clouds: direct location [41], ICP methodology [42][43] [44], methodologies based on geometric features of point clouds [45] and methodologies based on targets.The mathematical essence of point clouds registration is to solve the coordinate conversion matrix [46] and the mathematical model of multiple registrations can be denoted as: (2,3) . . .
The proposed monitoring system chooses the methodology based on targets and the registration is completed by business software Cyclone which corresponds to the Leica Geosystem.The registration methodology based on target chooses the center of the target as homonymy points to calculate coordinate conversion matrix and then achieves the registration.At least 3 targets need to be scanned from two adjacent scan stations in order to finish registration in the field work.

Acquisition of circle tunnel axis
The acquisition of tunnel axis is the most primary task in automatic data processing [47] [48].Accurate axis acquisition can help locate the tunnel posture which is applied as the fundament to serve point clouds denoising, data filtering and structure point clouds segment extraction to ensure smooth execution of the proposed processing method.The proposed monitoring system employs the algorithm based on RANSAC to acquire tunnel axis, which can improve the accuracy of axis acquisition by dealing with the effect of errors and noise.The origin of the coordinate system is positioned on the scanner center, the x-axis and y-axis are perpendicular to each other in the horizontal plane and z-axis is straight up.The spatial point clouds are projected to plane XOY and plane YOZ and then the boundary points of the projection are acquired to fit quadratic curve based on RANSAC as the contour line.Finally, the tunnel axis is calculated through the upper and lower contour lines.The procedures of axis acquisition are as follows: 1. Compress point data if necessary to reduce computation.Take all of the z coordinates of point clouds as zero to obtain the projection of the tunnel point clouds in plane XOY; 2. Assume the axis of tunnel forms an acute angle with x-axis and extract the maximum x coordinate value as x max and minimum x coordinate value as x min among all the points coordinates; 3. The boundary points are extracted piecewise with interval ∆t due to the discreteness of point clouds.The factors which influence the interval ∆t are mainly the accuracy of data acquisition, accuracy of the contour line and the available time assigned to data post-processing.The proposed algorithm take ∆t as 10 mm.In every interval from x min to x max , extract the maximum and minimum of the y coordinates of all points coordinates as the upper and lower boundaries points; 4. Considering that the length of tunnel section is usually shorter than 1000 m with curvature which lies in the same side of point zero, the proposed algorithm employs the quadratic curve to fit the upper and lower boundaries of tunnel point clouds: Where the corner mark "ub" and "lb" represent the upper and lower boundaries respectively.The parameter estimation of fitting function employs the widely used robust RANSAC algorithm [49]. 5.The upper and lower boundaries and tunnel axis are supposed to be of the same shape considering the short distance and little curvature of tunnel section.In other words, the function graphs of upper and lower boundaries in plane XOY can be transformed by the tunnel axis graph through translation only and the two translation vectors are of the same size but in opposite directions.Assume the function of axis in plane XOY is: Transform the graph of axis function shown in Formula (3) to get the graphs of upper and lower boundaries function of point clouds projection and the translation vectors are (p, q) and (−p, −q) respectively: Simplify Formula (4) and contrast the coefficients: The parameters of tunnel axis projection function in plane XOY can be calculated by Formula (5).
The function of tunnel axis in plane YOZ can be calculated by a method almost the same as the one used above.In general, the bottom of the tunnel can hardly be scanned by the TLS as it is covered by the road inside the tunnel.Instead of fitting two boundaries, the algorithm used in plane YOZ can only extract and fit the upper boundary of the projection and obtain the axis by translating the upper boundary distance r along z-axis negative direction.In order to make the methodology clear, a sample of tunnel point cloud is cut out and the steps are sketched and denoted in Fig. 2.  The proposed algorithm employs the RANSAC algorithm to estimate the parameters for the fitting function.By contrast, the classic parameter estimation algorithm least square method is based on smoothness assumption and cannot detect and eliminate the abnormal data.However, the smoothness assumption is not available in most cases including 3D point clouds data which contain noise which cannot be compensated.Thus the RANSAC algorithm is the key in the process of axis acquisition.

Tunnel point clouds segments extraction
The tunnel structures in cities, shield tunnels, pipe-jacking tunnel or other circular tunnels, are mainly constituted of concrete tunnel linings or tunnel tubes.On the other hand, the radial deformation of tunnel primarily occurs in tunnel segment joints with the result that the stiffness of joints is smaller than that of the tunnel segment.So the monitoring based on TLS of radial deformation should mainly focus on the tunnel segments [50].Therefore, it's meaningful to segment tunnel point clouds and extract the segment to analyze the radial deformation of tunnel.
In order to simply the process and improve the efficiency on the premise of sufficient accuracy, the segment extraction module applies the method that divides the tunnel point clouds into small parts equally to extract point cloud segments.The procedures of tunnel point clouds segment extraction are as follows: 1. Cut out a section of tunnel point clouds with n linings and axis with equal length axis acquired in Section 2.2.m; 2. Divide the axis of the selected tunnel point clouds into n segments equally and get n-1 equal diversion points in the axis; 3. Divide the selected tunnel point clouds with known number of linings by the normal planes in equal diversion points into n segments then the tunnel point clouds segments extraction finish; 4. Extract the tunnel point cloud segments between two adjacent normal planes and then the tunnel point clouds segments extraction finish.

Denoising of point clouds data
Various errors and noises exist in the acquired data and denoising is thus very essential.In general, the point clouds scanned by TLS in tunnel consist of parts of points reflected by the surface of the facility and equipment in the tunnel and this part of data is useless for subsequent data processing.Besides, noise is generated in point clouds data because of other human factors and environmental factors such as dust in tunnel, high reflective surface and miss operation.In a word, the noise should be removed to minimize its influence on the subsequent processing operations.
The proposed algorithm contains a special module based on clustering for noise removal with main focus on the redundant point clouds caused by the equipment and facilities in the tunnel.The raw point clouds of the structure contain mainly two kinds of noise: the first kind of noise is generated as a result of existence of facilities and equipment in the tunnel; the second kind of noise is generated due to the unusual reflectivity of scanned surface, generally light-reflecting and light-absorption surface.The first kind of noise is obvious and it is necessary to eliminate it or the data processing can't continue.The second kind of noise has a negative effect on the result accuracy and is difficult to locate and distinguish because the noise distributes in point clouds randomly.The cost of removal of the second kind of noise is fairly high.So the second kind of noise is eliminated in the follow processor which is based on RANSAC algorithm.
The proposed algorithm for denoising employs K-means algorithm, which is one of the simplest and commonly used algorithms, to do cluster analysis based on distance and intensity respectively.Firstly, the structure point clouds are clustered based on the difference between distance between the tunnel axis and the facilities in the tunnel and the distance between the tunnel axis and the surface of tunnel linings.However, some facilities are very close to the tunnel surface and the connecting pieces are in direct contact with the tunnel lining so that the difference of distance is too small to cluster.In conclusion, it's hard to eliminate the point clouds generated by the facilities efficiently and accurately.
On the other hand, the LiDAR intensity is based on the optical strength which highly depends on the reflectance properties of the scanned materials [51].So different materials may reflect laser beams with different intensity strengths even with an equal distance and incidence angle.The facilities and equipment in the tunnel are generally made of metal and rubber or the surface is covered with paint, which is totally different from the concrete tunnel lining surface with respect of reflectance properties.Set ω as dilution factor applied in uniform sampling of the massive data to reduce the computational time.
Where N i is the number of points in segment i and ω is the dilution factor defined by the amount of data and requirement of algorithm efficiency, N A is the actual point number in calculation.
The second clustering analysis is based on intensity and is implemented on the result of the first clustering analysis.Akca suggested that intensity values can supply decent additional information for identification [52].Where C 0 is the initial clustering center, I J is the intensity value, C represents class str or class equ and N c is the total point number of class str or class equ; 2. Calculate the distance between every residual point and two initial clustering centers and cluster the points to the nearest class.Then recalculate the clustering centers of the two classes: Where k is the iterations; 3. The iterative calculation ends when the difference between a new clustering center and last clustering center is equal to or less than the threshold.
Through the verification by experiments, the proposed denoising algorithm can remove the point clouds reflected by facilities and equipment efficiently and accurately.Fig. 7 illustrates the two-dimensional probability distribution of distance and intensity of the proposed algorithm.

Extraction of circle tunnel deformation
The proposed monitoring system based on TLS can make full use of massive spatial data acquired by the scanner and also deepen our cognition of structure state greatly.This chapter digs deeper into the tunnel point clouds segments to get more comprehensive information of tunnel's structural health.
According to the existing researches and engineering experiences, tunnel deformation is one of the most direct indicators used to analyze the tunnel condition [53] [54].So the proposed processing method chooses the radial deformation as the monitor indicator to reflect the safety condition of tunnel structure.Timothy et al. [18] calculated the algebraic deviations and mean absolute error and root mean square error between the measured radius value and the design radius to assess the condition of tunnels.
The radial deformation of circle tunnels which are usually constructed by assembling different segments is usually caused by compression and reversal of tunnel joints.In general, the circular tunnel structure in operation gets deformed before failure under the nonuniform stress from surrounding stratums and can be approximately treated as an ellipse with tiny difference between long axis and minor axis.Hereby, the ovality and three-dimensional deformation nephogram are employed to indicate radial deformation of tunnel structure [55].
The processor employs the RANSAC algorithm for curve-fitting to acquire the ellipse of the tunnel cross-section formed by scanned points and analyze the radial deformation.The ovality is illustrated by the following formula where T represents the tunnel ovality: Where F is the length of long axis and f is the length of minor axis of ellipse fitted, R is the design radial of tunnel.The ovality and three-dimensional deformation nephogram can make full use of the massive information contained in point clouds collected by TLS and demonstrate the deformation of almost every point in a tunnel to reflect its health condition.

Error analysis
Errors of the proposed processing method mainly consist of single point location error, measuring error and algorithm error [55].The single point location error is defined by the accuracy of scanner hardware.The TLS point clouds data were acquired from the Leica Scanstation C10 and measurements taken with this type of laser scanner resulted in range precision of 4 mm and range distance 134 m @ 90% reflectivity.The measuring error can be restricted within the permissible range efficiently by adjusting the distance between two adjacent stations and incidence angle.The algorithm error occurs as a result of some unavoidable limitations of algorithm can transmit and amplify the deviation and is difficult to evaluate [56][57][58] [59].
An operation is simulated to analyze the error produced in the proposed monitoring system, a conventional monitoring system based on total station is applied for subsidiary and the surveying result is compared to quantify the error and obtain the realistic accuracy of the proposed processing method.
Experiment was conducted in a pipe-jacking tunnel lining with a 4m diameter as shown in Fig. 10.Measurements of 55 points in the same cross-section of the tunnel lining were taken by conventional monitoring system based on total station twice and then the tunnel lining was scanned by TLS.The instrument parameters applied in two monitoring systems are shown in Tab. 1. Fig. 11 illustrates the tunnel lining profiles fitted by two sets of data collected by total station and point clouds scanned by TLS respectively.distance between the point with same central angle measured by total station and TLS in different measurements as measurement error of the proposed monitoring system and obtain the histograms (Fig. 12).Input the measurements illustrated in Fig. 12 for statistical analysis through Jarque-Bera test, the measurement error shown in Fig. 12 left follows the normal distribution with mean value of 0.1mm and standard deviation of 1.2mm and the right also follows normal distribution with mean value of -0.1mm and standard deviation of 1.2mm.

Application
The proposed monitoring system is aimed at achieving automation and high efficiency of tunnel structure health monitoring within the permissible engineering range.This section describes application of the proposed processing method to a practical project and compares the efficiency of point cloud processing method the automatic way and the conventional manual way.

Practical monitoring project
Shanghai metro line 7 was put into operation in 2009 with a total length of approximately 44.35km, which is constructed with two side-by-side tunnels with an inner diameter of 5.5m by slurry balance shield.The tunnel is formed by lining segments loop by loop with a length of 1.2m and every single loop is assembled by six lining segments with one top block, one bottom block and four standard blocks.
The LiDAR data is scanned by Leica Scanstation C10 with middle resolution in Shanghai metro line 7 and Fig. 13 shows the realistic picture and point cloud of the tunnel.The quality of the point clouds depends on the following factors: instrument mechanism, atmospheric conditions, object surface properties and scan geometry [60][61][62] [63].Javier Roca-Pardinas et al. [64] thought that distance to the object and angle of incidence mostly affect the accuracy of point clouds and mainly the angle of incidence.Pejić [65] thought the incidence angle shouldn't exceed 78 • theoretically while Lichti [8] defines the limit as 65 • or measurement error would rise sharply.Besides, Soudarissanane [66] conducted a series of experiments which indicated the incidence angle shouldn't exceed 70 • .Choose the distance between two adjacent stations 12 m, namely, 8 tunnel lining segments which are easy to locate in field work and maximum incidence angel 62less than the incidence angel limit 65 • .The whole monitoring zone is 150 m long and contains 100 tunnel lining segments.There are totally 12 TLS scan stations and every two adjacent point clouds are assembled by three HDS targets between them.

Monitoring project efficiency analysis
The proposed processing method mainly focuses on automation of point cloud processing.The proposed processing method processes the massive data in a fully automatic way to significantly improve the efficiency of point cloud processing and allow the engineers to do other more meaningful work.
The steps of manual point clouds processing are as follows: first label the tunnel linings in a fixed direction for convenience of subsequent processing; then segment the tunnel point clouds; remove the obvious noise, namely, the points reflected by the facilities and equipment inside the tunnel; implement ellipse-fitting for tunnel segment point clouds through business software to extract the axis; finally export the processed information.As to subway tunnel, the time consumed in processing a single tunnel lining point cloud segment is about 6 minutes according to experience.
The monitored zone contains 100 tunnel lining segments and the 3D point clouds data scanned in the zone is processed by conventional method and the proposed method separately, by the same computer.Tab. 2 denotes the difference between results of the two methods.According to the comparison from the above table, the largest component of computational time consumed in the conventional method is the manual noise removal, lining segment division and tunnel axis extraction.These manual operations are not only time-consuming but also affect accuracy adversely.
The conventional method of point cloud processing mainly depends on manual operations and business software, which is a time-consuming and labor intensive work with low accuracy.The business software is designed for all customers in various industries, so most of the functions are universal and not specially designed for tunnel point clouds processing.The manual extraction of tunnel lining segments cannot provide precise measurements due to manual partition of tunnel point clouds which totally depends on the difference in joint intensity with naked eye.The manual removal of obvious noise reflected by facilities inside the tunnel is very cursory because the selection tool in Cyclone (business software applied for raw point clouds data processing of Leica) is fixed and it does not suit the geometry of tunnel point clouds.On the other hand, the defects in manual removal of noise also result in over denoising or inadequate denoising.
According to the analysis above, the proposed monitoring system based on TLS which was successfully applied in practical engineering achieves monitoring of tunnel structure with automation and high efficiency within the engineering permissible range.

Conclusions
This article proposes an automatic circle tunnel point cloud processing method based on 3D LiDAR data that extracts the tunnel deformation through ovality and a three-dimensional nephogram.
The proposed processing method puts forward a series of efficient and accurate algorithms to automatically process point clouds of circle tunnel.Firstly it extracts the tunnel axis using algorithms based on RANSAC; then extracts tunnel segments from the entire tunnel point clouds precisely by equally dividing the axis; next removes the obvious noise by twice clustering; finally obtains the ovality and tunnel deformation nephogram for evaluation of circle tunnel radial deformation.
The proposed monitoring system mainly focuses on circle tunnels.According to the analysis of the field contrast test and practical application in Shanghai metro line 7, the proposed processing method significantly improves the efficiency of point cloud processing within the engineering permissible range.The indicators employed in result analysis especially nephogram obtain more abundant reference information for the evaluation of tunnel structure and can fully develop the advantage of massive information from 3D LiDAR clouds.

Figure 2 .
Figure 2. The sketch of every single step in tunnel axis acquisition.(a) is the projection of compressed tunnel point cloud in the plane XOY.(b) illustrates the upper points extracted from every point cloud interval from xmin to xmax.(c) illustrates the lower boundary points extracted from every point cloud interval from xmin to xmax.(d) illustrates the point cloud projection and all boundaries including inline and outline upper boundaries and inline and outline lower boundaries fitted through the points extracted in (b) and (c) after RANSAC.(e) illustrates all boundaries fitted through boundary points after RANSAC.(f) extract the inline upper, lower boundaries and axis extracted from inline boundaries.

Fig. 3 Figure 3 .
Fig.3illustrates procedure of tunnel point clouds segments extraction.The tunnel axis and tunnel point clouds segments are processed by a series of efficient and robust automatic algorithms to remove the interference from obvious noise and is then applied to analyze the tunnel structure condition.

Preprints
(www.preprints.org)| NOT PEER-REVIEWED | Posted: 18 June 2018 doi:10.20944/preprints201806.0283.v1 The proposed algorithm employs the second clustering based on intensity at the outcome of clustering based on distance to assist the analysis and improve the accuracy.The flowchart of the proposed algorithm is illustrated by Fig.4.

Figure 4 .Figure 5 .
Figure 4. Flow chart of proposed algorithm based on twice clustering based on distance and intensity orderly to denoise

Fig. 6
denotes the intensities distribution of point cloud in experiments.It's obvious that most of points concentrate in the area around value -1500 to -1000.

Figure 6 .
Figure 6.Intensity distribution of raw point clouds in experiments

Figure 7 .
Figure 7. Two-dimensional probability distribution of distances and intensities

AmplifyFigure 8 .Figure 9 .
Figure 8.The ovality of a tunnel point cloud segment produced by the proposed monitoring system.(a) is the comparison between the ellipse fitted by point cloud and the standard circle which represents the designed schema.(b) is the unfolding drawing of (a) and illustrates a more intuitional comparison about deformation at different central angles.

Figure 12 .
Figure 12.Histograms of differences between measurements by TLS and two measurements of the total station.(a) denotes the quantitative distribution of measurement error of the proposed monitoring system which takes the first set of data from total station as standard; (b) denotes the quantitative distribution of measurement error of the proposed monitoring system which takes the second set of data from total station as standard.

PreprintsFigure 13 .Figure 14 .
Figure 13.The realistic picture and point cloud of Shanghai metro line 7 scanned by TLS in this monitoring project.There is a Leica C10 scanstation in tunnel

Figure 15 .
Figure 15.The overall deformation nephogram of the Shanghai metro line 7

Table 2 .
Comparison in post-processing time between two methods