POINT CLOUD MODELLING BASED ON THE TUNNEL AXIS AND BLOCK ESTIMATION FOR MONITORING THE BADALING TUNNEL , CHINA

Recent years have witnessed a growing investigation of terrestrial laser scanning (TLS) for monitoring the deformation of tunnels. TLS provides the ability to obtain a more accurate and complete description of the tunnel surfaces, allowing the determination of the mechanism and magnitude of tunnel deformation, because the entire surface of the tunnel is more concretely modelled rather than being represented by a number of points. This paper models and analyses the point clouds from TLS to detect the possible deformation of a newly built tunnel. In the application of monitoring the Badaling Tunnel for the Winter Olympics 2020 in Beijing, China, the proposed method includes the following components: the tunnel axis is automatically estimated based on a 3D quadratic form estimation; all of the point clouds are segmented into axis-based blocks; and representative points, solved by a singular value decomposition (SVD) method, are estimated to describe the tunnel surface and establish the correspondence of data between days. The deformations are detected in the form of the distance discrepancies of representative points and verified by the measurements using total station.


INTRODUCTION 1.1 Background
Terrestrial laser scanning (TLS) can collect millions of 3D point cloud data per second, as well as intensity information.Multiple scans taken at different times and positions from an advancing tunnel face can be processed to describe the profiles of a tunnel and how they change over time.This quantified information can be used to determine the deforming orientations and magnitudes of the tunnel.
The magnitude of deformation is related to the rock characteristics, the excavation method, and the type and location of installed support.Deformations serve as early warnings of more severe failures.Monitoring tunnel deformation is an indispensable part of tunnel construction and operation.
The advantages of using TLS in tunnels rather than other surveying methods include the following: 1) 3D coordinates can be collected on the geometry of the tunnel along its entire length rather than on specific sections, usually at intervals of several metres; 2) TLS collects data without the influence of illumination; and 3) any sections/areas of a tunnel can be extracted for study purposes.

Previous work
Early work in Lam (2006) described the basic methods of TLS used for tunnel applications both in the field and for computational algorithms.In geometric tolerance analysis, the author proposed three datum axes for assessment of horizontal alignment, vertical alignment and profiles in slicing.This work opened the door to actual TLS implementation.
An overview of TLS in tunnel excavation is first discussed and summarized.Gikas (2012) validated the smoothness/thickness of shotcrete layers at an excavation stage and compared them against the dimensional tolerances in the concrete lining formwork.The quantitative findings resulting from the case studies could not be directly generalized.Wang et al. (2014) also overviewed several applications of TLS in tunnels and emphasized how TLS can be used to study aspects of tunnels.Similar applications refer to Moisan et al. (2015) and Lague et al. (2013).
The optimization of scanning parameters, scan registration, georeference and survey network has been investigated for tunnel geometry inspection (Pejic 2013).The author noted that the noise should be considered because of the physical limits, such as the incidence angle of TLS in tunnel surveying.Subsequently, Roca-Pardinas et al. (2014) analysed and modelled an error model with the consideration of the influence of ranges and incidence angles by means of Monte-Carlo simulation.Taking account of point density, incidence angle and footprint size, Argueelles-Fraga et al. ( 2013) scanned circular cross-section tunnels and determined the maximum scan distance and angular sampling interval to control data quality and working time.In general, the main factors influencing the point quality were as follows: incidence angle variations when scanning in a long and narrow corridor; occlusion by lines or machines; reflection from water and mud; and atmospheric related factors, such as heavy dust.These researches provided field guidance for distance measurements during data collection.
Eling (2009) used the iterative closest point (ICP) method to transform point clouds from one scan station to a reference scan station in the Oker dam application.Later, the ICP method was applied for fine registration (Wang et al. 2012).Chmelina et al. (2012) investigated translations by the ICP method between scans under the assumption that the transformation parameters refer to the centre of gravity of the respective patch.They suggested that the density of the distribution can be adjusted according to the requirements so that the tunnel surface comes close to a developable surface by varying the width of the mesh.In general, the ICP method requires iteration and does not necessary converge.Because each patch is processed individually, the displacement vectors are not checked for consistency with the overall deformation behaviour or for whether the entire vector field can be explained physically by the tunnels under construction.However, the explicit point correspondence is challenging to confirm, which may cause lower precision of the registration.Subsequently, Ji et al. (2015) registered adjacent scans by using closing conditions and external geometric constraints to eliminate accumulate errors and ensure the shape correctness of a high-speed railway viaduct.For high-precision scanning, the total error budget for change detection is dominated by the point cloud registration error and the surface roughness (Lague 2013).The authors proposed a cloud-to-mesh method whereby a point cloud is meshed and the distance to another is computed along the local normal of the point cloud to compare distance changes between point clouds.
Lacking an explicit geometric interpretation, Han et al. (2013) proposed a minimum-distance projection algorithm to establish point correspondences.The advantage of this approach is that deformation signals along any profiles can be immediately identified.Using an elliptical fit method, the levels of real convergence have been measured and statistically tested in tunnels and shafts (Diederichs et al. 2014).Even with relatively high degrees of occlusion, up to 40%, an elliptic fit analysis was used to analyse general trends of deformation, and profile analysis was used for anomalous movements additionally supported to maintain stability (Walton et al. 2014).However, discrete profiles cannot describe the entire deformation of tunnel surfaces.
There are two basic approaches to describe the shape of an object surface: point segmentation based on criteria such as the proximity of points, and direct estimation of surface parameters (Vosselman et al. 2004).The second approach produces a unique and relatively simple mathematical surface when the object surface has a regular shape, as in the case of tunnels.It provides a convenient way to automatically estimate the tunnel centre axis, which is necessary for segmenting raw point clouds in tunnel applications.
In general, geodetic engineering operations cover the entire lifecycle of tunnel construction and operation.The data processing in tunnel applications can be grouped as follows: 1) geodetic control network configuration and registration work, including calibration; 2) mapping of the tunnel corridor; and 3) monitoring tunnel deformation.Although many tunnel applications have used TLS data, a reliable algorithm that full use of huge amounts of point clouds is still to be developed.The overwhelming amount of data that can be extracted should not simply be used in raw form but should also be used in describing deformation.Inspired by these significance contributions and suggestions, this paper investigates a method to describe the possible deformation with TLS point clouds.

Purpose
The aim of this study is to develop a method that quantifies tunnel deformation in an underground railway station using highdensity and time series data acquired from TLS.
Specifically, this research focuses on: 1) automatically estimating the tunnel axis with quadratic forms of the tunnel surface; 2) segmenting all the point clouds into small rectangular blocks; and 3) finding the corresponding representative points for day comparisons by plane surface estimation in blocks.

Overview of the proposed method
Raw point clouds can be used for visualization, but to perform deformation analysis, data processing is required to analyse and establish the corresponding relations of random points.How one establishes a suitable geometric model to fit the real tunnel surface is greatly important for further deformation comparisons.It is also the focus of this research.A basic workflow of the proposed method includes the following components (Figure 1): pre-process data to eliminate outliers and register all scans to a uniform coordinate system; estimate the tunnel axis with an automatic method; segment all the point cloud into blocks; solve representative points with plane surface estimation; and compute and compare possible deformation between data from different days.
Figure 1.Workflow of the proposed method

Estimate tunnel axis with quadratic forms
A quadratic form can be used to evaluate the determinants and then test the form parameters to describe the object surface (Kutterer and Schoen, 1999;Hesse and Kutterer, 2006).The equation for a quadratic form estimation can be written as where   is the coordinate vector of a single point, M is the symmetric coefficient matrix,  is the coefficient vector,  is the scalar factor, and k = 1,2, , … , n, where n denotes the total number of points scanned by TLS.The parameters   ( = 1,2, … ,10) can also be expressed as The parameters  1 −  10 can be estimated bya means of singular value decomposition (SVD) (Felus and Burtch, 2009) without the steps of iteration.A determinant method is constructed to estimate four motion invariant parameters.Then, the shape of the object surface can be determined by looking for the test tree of automatic form recognition.The determinant method and the test tree are described elsewhere (Kutterer and Schoen, 1999).The advantage of this method is that only several parameters need to be estimated to describe the shape of the object surface.Thus, the coordinates   = [  ,   ,   ] of the tunnel axis (perpendicular to the profile direction) are derived by the estimated geometric parameters, where   and   are determined by the automatic form recognition and   is obtained from the maximum and minimum z-values of point clouds.The estimated tunnel axis is employed as the reference for further segmenting of point clouds.

Segment point clouds
Based on the determined tunnel axis, a segmentation method is used to identify the object surface and to establish a uniform framework for different days.Based on this framework, all of the point clouds are divided into small segmented blocks, and the definition of the blocks is explained as follows.
According to the directions of the plane surfaces, for example, we set the  direction as the tunnel axis direction and the  direction as the profile direction.The point clouds and corresponding coordinates of the tunnel axis are grouped into profiles with a certain thickness   of the tunnel axis.To conveniently perform segmentation, the tunnel axis is shifted to the left of the tunnel according to the point clouds.The moved tunnel axis is termed the moved-tunnel-axis, and the point coordinates of this axis are sorted by an ascending order.In the profile with a certain thickness, the horizontal values (   ,    ) of the four corner points of the profile are computed as follows: where  denotes the -th column number of the block, ∆ is the pre-defined block width along the profile direction, and (  1 ,   1 ) and (   ,  )  ) are the first and last coordinates of the movedtunnel-axis in the profile, respectively.The total number of profiles equals the number of rows of segmenting, and the columns of the segmentation are determined by the pre-defined block width.Thus, the four corners of the blocks are obtained to describe the framework of the tunnel surface.
The centre horizontal values of a block are averaged by their corresponding four corner points, and the vertical value of the block is computed by averaging all the vertical values of the point clouds in the block.Thus, the geometric centre point (  ) of the block is set as a reference point for estimating a representative point in the next step.

Solve representative points with plane estimation
The point clouds within a block are modelled as a representative point.According to the point clouds in a block, a total leastsquares method based on an SVD procedure (Horn and Johnson 1990, p.153;Felus and Burtch 2009) is introduced to estimate the representative point of the block, under the assumption that both the observations and coefficient matrix contain by random errors.
Because the selected point clouds in the block normally can be deemed as a planar surface, the normal vector  ̂ and the distance  ̂ of the block are estimated to evaluate the planar surface when there are more than    points (40 points for this dataset): where  ̂  is the coordinate vector of the point and  is the number of the points in a block.The covariance matrix   ̂ includes the stochastic values  n ̂n ̂ and   ̂ ̂ of the estimated parameters: where  ̂ 2 is the variance factor.According to the computational method of Drixler (1993), the correlations between the normal vector and the distance are not considered.The representative point  ̂ is estimated using the geometric centre point   of the block: The estimated covariance matrix  ̂ ̂ ̂ of the representative points is obtained by means of variance propagation: where  is the differential function with respect to  ̂ and  ̂.Thus, the coordinates and the covariance of the representative point are derived to estimate and evaluate the precision of the point clouds.

Study area
The Badaling Tunnel station lies approximately 80 km northwest of Beijing.This underground tunnel is part of the intercity railway built for the Winter Olympics 2020 from Beijing to Zhangjiakou.The designed speed varies from 250 to 350 km/h.A Leica P40 scanner, with a nominal 3D position accuracy of 3 mm at 50 m and angular accuracy of 8 ′′ , was applied to scan part of the tunnel between Dec. 2 nd -11 th , 2016.A Leica TS02plus2 total station with a nominal range accuracy of 1.5 mm + 2 ppm and angular accuracy of 2 ′′ was used to establish the monitoring network.Setting the origin of scan 3 as a starting point and scans 3 and 4 as a reference, five scans of point clouds each day were transformed into a uniform coordinate system in the Cyclone © software (Figure 2).Three days (Dec. 2 nd , 4 th and 9 th ) of point cloud data on the top surface of the tunnel were processed and compared to detect possible deformation.The proposed method was implemented in the Matlab © program.

Experimental results and analysis
In the reference coordinate system, the point clouds from the top tunnel surface are modelled as an elliptical cylinder according to estimates of the quadratic surface and looking for the test tree.
The four estimated parameters of the elliptical cylinder include the centre coordinates ( ̂ and ê  ) and the axes (a ̂ and b ̂) of the elliptical cylinder to express the shape of the top tunnel surface (Table 1), where  ̂ is the mean value of the point clouds in the direction of the y−axis.It is noted that the estimated parameters only approximately describe the coordinates of the tunnel axis.To segment the point clouds conveniently, according to the manually selected two points from the left bottom and right bottom of the point clouds, all point clouds were transformed along the z-axis with a rotation angle of 2. 23.Using this rotation angle, all the day data were transformed into a regular reference coordinate system in which the x -axis is along the profile direction and the y-axis is perpendicular to the profile direction.
Using the estimated centre axis and the pre-defined block size (approximately 23.80 cm × 40 cm ), the point clouds of the tunnel surface were segmented into 3195 blocks with 213 rows along the profile direction and 15 columns perpendicular to the profile direction.
Solved by the SVD algorithm in a block, a representative point and its standard deviation were obtained to represent the coordinates and precision of the block.All representative points are shown in Figure 3.  Setting the total station measurements on Dec. 2 nd as the reference point [0, 0, 0], the movements (  = [   ,    ,    ]) relative to the reference station at other days are computed: where  0  is the total station measurements from Dec. 2 nd and    ( = 1 and 2) denotes the total station measurements from Dec. 4 th and 9 th , respectively.The results are shown in Table 2, where the subscript  is ignored to simplify the expression.The referenced scan 3 shifted approximately 8.91 mm and 7.73 mm relative to the position on Dec. 4 th and 9 th , respectively.Considering the nominal accuracy of 1.5 mm + 2 ppm provided by the manufacturer of the total station, we determined that the network moved approximately 7. 41 mm and 6.23 mm relative to the position on Dec. 4 th and 9 th , respectively.
Dec. 2 nd 0 0 0 0 Dec. 4 th 7.30 -5.10 0 8.91 Dec. 9 th 4.90 -5.90 1.00 7.73 Table 2. Total station measured movements of scan 3 The distance discrepancies computed from the point clouds with the colour bar changing from 0 to 4 cm are shown in Figure 5 and Figure 6, which describe the distance discrepancies from Dec. 2 nd -4 th and Dec. 2 nd -9 th , respectively.The mean distance discrepancies estimated from TLS are approximately 15.89 mm and 13.46 mm from Dec. 2 nd -4 th and Dec. 2 nd -9 th , respectively.In many tunnel monitoring tasks, the deformation on the top surface is larger than that of the sides and the bottom of the tunnel.
The TLS results are estimated based on the point clouds on the top surface of the tunnel; while the total station observations are collected on the ground surface.This inequality caused the distance discrepancies of the TLS to be larger than that of the movements from the total station.The distance discrepancies on the right part are larger than those of the left part in both of the comparisons.Because the network was established on the left part, the point clouds on the right part contain larger errors according to error propagation.Data collection was influenced by construction work on Dec. 9 th , and therefore there is no data on the right part of the tunnel in Figure 6.

Discussion
Due to the influence of tunnel surface roughness, the 3D quadratic form only approximately describes the shape of the tunnel surface.However, the coordinates of the estimated tunnel axis used in this research need not be in the real tunnel centre axis; they are only references used to segment all of the point clouds over days.
Because the two sides of the tunnel are straight, only top surface data were extracted and processed in this application.To better describe the tunnel deformation, full-section modelling of the tunnel is necessary in the future.If the tunnel surface is sheltered by obstacles or heavy dust, outlier detection methods need to be considered.
Considering the point density and deformation magnitude, the definition of a block size should be considered for different applications.A block with a large size can include sufficient point clouds but corresponds to a loss of the deformation description, whereas a block with a small size may lose computing efficiency.In applications, the block size should consider both the objects and the need for monitoring accuracy.
In this application, artificial targets were mounted on the ground surface and the side wall.It would be better if more targets were mounted on the top and side walls of the tunnel, so that the estimated representative points from the TLS and corresponding targets measured by the total station can be compared to obtain more reasonable monitoring results.
Excluding the possible movements between days from the results of TLS, the distance discrepancies still contain some errors from the TLS field work plane estimation of blocks and surveyors (such as the incidence angle).To fully optimize for detection of minimal deformation, proper geodetic error propagation should be investigated further to fully determine the data quality of TLS in structural monitoring applications.Both GPS and TLS results should be combined to describe and compare the deformation in global coordinate systems.For larger projects, optimal point density should be considered to save field work time and data processing.
The observations from the total station demonstrate the ground surface moved approximately 7.41 mm and 6.23 mm, while the mean distance discrepancies estimated from the TLS were approximately 15.89 mm and 13.46 mm at Dec. 2 nd -4 th and Dec. 2 nd -9 th , respectively.The ground surface movements could be deemed as a reference to compare the top surface movements, however, real deformations on the top surface need to be further separated from the derived distance discrepancies, because these distance discrepancies include both the real deformations and the errors from the TLS.
Except for the manually selected corner points from the point clouds and the pre-defined block size, the entire procedure is automatic and takes approximately 18.82 minutes using an Intel(R) Core(TM) i7-4790 CPU with a 64-bit operating system and 16 GB RAM, which consists primarily of computing time rather than data input time.

CONCLUSIONS AND FUTURE WORK
In this paper we present a deformation detection method in the construction stage of a tunnel based on point clouds from TLS.
To compare the data on different days, a uniform reference tunnel model is established using estimated parameters of the 3D quadratic form.A segmentation is then conducted to divide the object surface into axis-based blocks.The point clouds are estimated as a representative point in a block, which is solved by the SVD method in a quasi-planar surface.These estimated points establish explicit correspondences of the data from different days.The proposed method was applied to monitor the Badaling Tunnel, which was built for the Winter Olympics 2020 in Beijing, China.The network measured from the total station observations moved approximately 8.91 mm and 7.73 mm, while the mean distance discrepancies estimated from the TLS data were approximately 15.89 mm and 13.46 mm relative to the data for Dec. 2 nd -4 th and Dec. 2 nd -9 th , respectively.
This work is the first step towards the application of TLS in the Badaling Tunnel.In future work, we will explore and develop an outlier detection method for better reconstruction.Furthermore, more comprehensive profiles will be considered to give a more concrete 3D model.

Figure 2 .
Figure 2. Experiment (top) and top surface point clouds of tunnel (bottom)

Figure 3 .
Figure 3.Estimated representative points on the top tunnel surface on Dec. 2 nd

Figure 4 .
Figure 4. Standard deviations of representative points on Dec.2 nd with colour bar changing from 0 to 0.3 mm

Table 1 .
Estimated parameters of the elliptical cylinder[m]