TOWARDS DESCRIBING FULL-SECTION DEFORMATIONS USING TERRESTRIAL LASER SCANNING IN THE BADALING TUNNEL (CHINA)

This paper focuses on the analysis of point clouds from terrestrial laser scanning to interpret possible deformations of the new Badaling Tunnel that was built for the Winter Olympics 2022 in the nearby of Beijing, China. A reference framework is established to compare data corresponding to various days with blocks of uniform columns and rows from an estimated tunnel axis. Filling holes and detecting outliers are performed for quasi-planar estimation, and refinement transformation is used to adjust the data errors between different days. Finally, the full-section deformations are detected in the form of distance discrepancies of representative points and are verified against total station measurements. * j.wang@bjut.edu.com; Tel.: +86-10-67396062 The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XLII-3/W10, 2020 International Conference on Geomatics in the Big Data Era (ICGBD), 15–17 November 2019, Guilin, Guangxi, China This contribution has been peer-reviewed. https://doi.org/10.5194/isprs-archives-XLII-3-W10-235-2020 | © Authors 2020. CC BY 4.0 License. 235


INTRODUCTION
The application of terrestrial laser scanning (TLS) techniques in tunnels is a topic of interest in surveying engineering. In such projects, the advantages of using TLS for geometric surveying, surface condition assessment (Chmelina et al., 2012;Bartolo and Salvini, 2019) or deformation measurement (Lindenbergh and Pietrzyk, 2015;Wang et al., 2017) are described in the following. Three-dimensional point coordinates can be collected to describe the complete surface of the tunnel vault and pavement along its overall length. This implementation is not limited to specific discrete sections typically spaced at intervals of several metres, as with total stations and profilometers or using image-based techniques (Scaioni et al., 2014;Attard et al., 2018). The independence from illumination and object textures (as in the case of photogrammetric techniques) is applied. A short time is necessary for data acquisition, which is allowed by last-generation phase-shift laser scanners (up to 10 6 Hz). In addition, TLS provides the same properties that other instruments may do, such as the simplicity of the adopted georeferencing procedures, high degree of automation, and a completely digital operational workflow.
A notable application of TLS concerns monitoring tunnel deformations within times by comparing point clouds or their by-products collected at various epochs. Several applications are reported in the literature (Wang et al., 2014), but one of the major remaining problems is the evaluation of the quality of the results. Indeed, typical deformations that may occur in tunnels may result in point displacements on the order of a few millimetres, and these displacements may be comparable to the measurement noise of TLS technology. Statistical analyses are then required to distinguish real deformations (Fey and Wichmann, 2017) from background noise and to establish quality measures of the observed changes.
The use of TLS point clouds for tunnel deformation monitoring can be considered a particular case of a more general problem, the Areal-based deformation measurement, see Capra et al. (2015) and Wunderlich et al. (2016). A thorough critical review of this topic was reported in Lindenbergh and Pietrzyk (2015), to whom we refer for the main background. In general, a distinction should be made between two categories of objects or surfaces. The first category addresses those objects that can be modelled using a mathematically known geometry at the local level (see, e.g., Lindenbergh and Pfeifer, 2005) or global level (see, e.g., Van Gosliga et al., 2006) to be used for interpolating the point clouds to be compared and to detect deformations. After interpolation, the analysis is conducted by comparing the parameters of the estimated regular surfaces. The second group addresses those surfaces that cannot be modelled with mathematic functions but have to be locally analysed (see, e.g., Scaioni et al., 2013). This second category entails predominantly natural surfaces, such as mountain slopes and coastal areas. In this second group, Lague et al. (2013) proposed a flexible method (M3C2) in which the distances from points of the second point cloud with respect to the first one are computed along the local normal vectors in the nearby of each point. This approach does not use a unique reference surface that limits the detection of deformations in only one direction and can be used to analyse objects that may present nonisotropic deformations. In addition, the uncertainty of scan coregistration and the local surface roughness is considered when discriminating between significant deformations and noise.
The basic methods for the application of TLS in tunnel deformation analysis are described by Lam (2006). TLS was proposed to effectively assess the various geometric tolerances of tunnel structures in as-built surveying. Gikas (2012) validated the smoothness/thickness of shotcrete layers at an excavation stage and compared these parameters against the dimensional tolerances in the concrete lining formwork. The quantitative findings resulting from the case studies were used to over-break volume computation. Considering the point density, incidence angle and footprint size, Argüelles-Fraga et al. (2013) scanned circular cross-section tunnels and determined the maximum scan distance and angular sampling interval to control data quality and data acquisition time. Roca-Pardiñas et al. (2014) analysed and formulated an error model that considered the influences of ranges and incidence angles by means of Monte-Carlo simulation.
Some studies have followed a strategy based on the extraction of multiple cross sections from the complete point clouds. This approach may sometimes be limited since the whole data set is not completely exploited due to subsampling (Wunderlich et al., 2016). In reality, in tunnel deformation analysis, certain dense cross sections may suffice to understand the presence of structural problems. Unlike other surveying techniques, TLS allows the selection of those cross sections for investigation at the post-processing level and not at the data acquisition time. Han et al. (2013) presented a solution in which a minimum-distance projection algorithm was applied to establish point correspondences when an explicit geometric interpretation of the cross section was missing. Puente et al. (2016) also used point clouds to extract vertical clearances in the tunnel environment. This method was not limited to extract cross sections. Walton et al. (2014) applied an elliptical-fit analysis to study the general deformation trend, while the analysis of local profiles was used to look for anomalous movements. The advantage of this approach is that deformations along any profile can be immediately identified. However, extracting profiles for deformation analysis abandons most of the point clouds collected in the field, and the discrete profiles cannot completely describe the deformation of tunnel surfaces. As proposed by Kang et al. (2013), continuous cross sections should be investigated to detect the deformation of tunnels. Continuous profiles with a pre-defined thickness save enriched coordinate information and provide a full-section description of the tunnel surface.
Occlusions and accessibility limitations may lead to the incomplete reconstruction of the point cloud, thus introducing holes. To model the complete tunnel surface, it is necessary to reconstruct the missing parts using the available point clouds. Algebraic, computational geometry and implicit methods are generally investigated for surface reconstruction and hole filling (Wang and Oliveira, 2007). Algebraic methods (e.g., Terzopoulos and Metaxas, 1991) recover a surface by fitting a smooth function and cannot be used to reconstruct complex geometry. Computational geometry methods may leave holes in under-sampled regions and require longer computation time. Implicit functions cannot be directly used for the reconstruction of surfaces with boundaries (e.g., Carr et al., 2001). Sharf et al. (2004) used an octree to copy and paste samples from a set of 'template' regions to fill those missing areas. This method can produce impressive results but does not guarantee that the true surface is reconstructed. These drawbacks become crucial in the presence of deformations as discussed in Vichitvejpaisal and Kanongchaiyos (2014).
Outlier detection should also be addressed in high-density point clouds obtained in tunnel surveying (Pejić, 2013;Cao et al, 2019). For example, outliers may be due to varying incidence angles in which some parts of the point cloud may be relatively small. The RANSAC algorithm (Fischler and Bolles, 1981) and its derivatives has been widely applied to remove outliers, even though the ratio with respect to inliers is large (see, e.g., Schnabel et al., 2007). Indeed, the RANSAC algorithm can robustly accommodate data containing up to 50% outliers (Roth and Levine, 1993). For this reason, RANSAC was selected to detect and remove outliers in the application presented here.
Inspired by these important contributions and suggestions (particularly when using completely raw point clouds), we developed a method that quantifies tunnel deformation by searching for corresponding points in a time series of highdensity point clouds collected during the construction stage of an underground railway line. Indeed, at this stage the deformations are typically larger than the ones occurring during serviceability stage and may vary from millimetres to centimetres.
The next Section 2 describes the proposed method. In Section 3 a case study is presented. Data collected in a tunnel were processed using the proposed method and validated (see discussion in Sect. 4). Finally, conclusions are drawn in Section 5.

Overview of the proposed method
The processing pipeline starts after certain preliminary stages. First, multiple scans are registered into a common coordinate system using appropriate techniques that result in a registration error lower than measurement noise. Second, those points in the data set that are located on the pavement of the tunnel are manually eliminated. The workflow of the proposed method is described in Figure 1.

Estimate of the tunnel axis and setup of the reference frame
To conduct a fast segmentation of the registered point clouds, these clouds are transformed from the "raw" reference system into a right-handed Cartesian coordinate system. The transversal profile direction is defined as the -axis, while the -axis is towards the heading of the tunnel, and the vertical -axis is towards the top of the tunnel (Fig. 2).  ) are the coordinates of points in the 'new' coordinate system; the yellow block is the example of a segmented block).
The transformation is done based on three manually selected points from the raw point cloud data ( where ∆ = [0.5 * , 0, 0] is the shift vectors, is the minimum -value of the point cloud, is the rotation matrix with the angle computed from the corner points under the assumption that a scale factor has no real influence in a short distance, and is the vector with coordinates of the point clouds in the regular coordinate system. Therefore, the new coordinate system can be used to establish a uniform reference framework. In the case where the heading direction of the tunnel varies, the transformation parameters must be changed. In the 'new' coordinate system, the tunnel surface in the point cloud can be segmented into columns and rows. The row data set is processed using profiles with a certain thickness ( ) along the heading direction of the tunnel, and the columns are defined along the profile direction. In each profile, the -value of the tunnel axis is searched from the maximum -value in the point clouds, and the y-values of the tunnel axis are searched based on the minimum y-value of the point clouds. Although the estimated tunnel axis may not coincide with the physical tunnel axis, this axis is suitable to be used to compare multiepoch data.
In the segmentation, the profile numbers are set as the row numbers of the blocks. Then, parallel with the tunnel axis, the number of columns is computed with a pre-defined block width ∆ . The coordinates ( , , ) of the -th point from the point clouds are grouped into the corresponding block with the number − ∆ , . Thus, all the point clouds are divided into small segmented blocks (see the yellow example block in Figure 2) of equal size that are similar to the concept of octrees (Samet, 2006).

Quasi-Planar Estimation and Outlier Detection
Since in a tunnel surface the neighbouring geometry often remains invariant within each block, points in neighbouring blocks can be used to reconstruct the missing surface due to an occlusion. Thus, to fully describe a tunnel surface, holes are filled by processing of their neighbouring points. In the -th block of -th column and the -th row without enough points, the points' coordinates ( = , , ) in the -th block are recovered on the basis of the points ([ , , ]) belonging to the neighbouring -th block: Due to the small size of each block after segmentation, points within each block generally feature a quasi-planar surface. Given this hypothesis of local planarity, in the case that the number of points in the block is over a minimum threshold points (in general 100 points), the RANSAC algorithm is applied to exclude those points that have not been scanned from the tunnel surface. To detect local deformations between different days, an approach similar to the one proposed in Lindenbergh and Pfeifer (2005) is applied to data sets that fall inside a given block. The central horizontal coordinates of a block are averaged by their corresponding four corner points, while the vertical central value is averaged on the basis of all points of the block. The geometric centre point ( ) of the block is set as an initial reference point for estimating a representative point. The normal vector = , , and the distance of the block are estimated to evaluate the planar surface when there are more than points (100 points for this data set). According to the computational method of Drixler (1993), the correlations between the normal vector and the distance are not considered. This model is solved using the singular value decomposition (SVD) procedure (Horn and Johnson, 1990;Felus and Burtch, 2008).
The representative point is estimated using the geometric centre point of the block: where = 1 . The estimated covariance matrix Σ of the representative points is obtained by means of variancecovariance propagation. The standard deviations (std. devs.) of the representative points are computed by the variances of the representative points in the directions of the -axis, -axis and -axis, respectively. Thus, the representative points are used to depict the tunnel's surface instead of the raw point clouds.

Compensation of Systematic Misalignments and Multi-Epoch Comparison
In reality, the quality of the surveyed point clouds may be influenced by many error sources, such as the monitoring network's geometry, the registration process and environmental factors. Since, in practice, it is challenging to detect individual error components, a local transformation is accomplished to compensate for the possible presence of systematic errors. The representative points at each epoch are used as identical points for the estimation of the refinement similarity transformation (three translations ∆ , one scale and three rotation parameters ): ) T denotes the coordinates of the identical representative points from the epoch that are moved, and =( , , ) T denotes the coordinates of the identical representative points from the reference epoch. The SVD algorithm is used to compute the transformation parameters within an errors-in-variables (EIV) model (Felus and Burtch, 2008) where coordinates observed at both epochs are treated as stochastic variables (see more details on the adopted algorithm in Wang, 2013). All corresponding representative points from corresponding blocks are processed as identical points. The application of the estimated transformation allows us to mitigate the effects of systematic misalignments between multi-epoch data.
Ideally, if there are no deformations, the average coordinates of the corresponding representative points should be identically the same, but the random and residual systematic errors could make them different. Since the coordinate differences of these representative points may reveal possible deformations, the 3D distances between the corresponding representative points between two epochs are used to describe the local and global tunnel deformations: where (Σ ) and (Σ ) are the estimated covariance matrix Σ of the representative point coordinates in the reference and compared epochs, respectively.

Description
The Badaling Railway Tunnel lies approximately 80 km northwest of Beijing. This tunnel is part of the intercity railway built from Beijing to Zhangjiakou for the Winter Olympics 2022 (Fig. 3a). The total length of this tunnel is 12.01 km, and its width and height are approximately 7 m and 4.3 m, respectively. The rock structure of the tunnel contains weak weathering monzogranites stratum. Rock cracks have developed, and the rock mass are fractured. Therefore, deformation monitoring is essential during construction works of the tunnel.

Data Acquisition
To establish the geodetic reference frame for deformation monitoring activities, a geodetic network was setup and measured inside the tunnel. This network was also linked to some GNSS benchmarks located outside the tunnel for the purpose of geo-referencing in the mapping grid. A highprecision Leica Geosystems TS02plus2 total station with a nominal range accuracy of (1.5+2 p.p.m.) mm and an angular accuracy of ±2 was used to establish the geodetic network. The scan positions were measured daily by total station.
A Leica Geosystems P40 TLS (Figure 3d) with a nominal 3D position accuracy of ±2 mm at 50 m and an angular accuracy of ±8 was applied to scan the tunnel between the 2 nd and 13 th December 2016. All scan standpoints were measured with the total station and their 3D coordinates determined with respect to the geodetic network. Five scans (Figures 3b and 3c) were surveyed in each day (epoch) and registered into the geodetic coordinate system using Cyclone© software. The registration accuracy was better than ±4 mm using the identical artificial targets. Since the other scan positions were damaged or blocked due to construction, only Scan Stations 3 and 4 were set as references for registration. The movements between Scan Stations 3 and 4 in a day were not considered. The point clouds close to the ground were filtered out because of the rough surface and because they were partially sheltered by instruments (Fig. 3d). The point density on the tunnel's surface after registration varied from 2 mm to 1.5 cm (Fig. 3e). Four observation epochs were recorded on the 3 rd , 9 th , 11 th and 13 th December 2016 (we did not scan every day). All the algorithms adopted for this analysis were implemented in MathWorks MATLAB © code. Since the point cloud scanned on the 9 th December is more complete than the others, this data set was selected as the reference for comparing the other three epochs. The results from the comparison of the 3 rd vs 9 th December epochs, the 9 th vs 11 th Decemer epochs, and the 9 th vs 13 th December epochs are named Subset 1, Subset 2 and Subset 3, respectively.
Considering each subset, first, the displacements of the positions of each scan stations' ( ) were computed: where 0 = [ 0 , 0 , 0 ] are the coordinates in the reference epoch, and ] are the coordinates from the other epoch in the subset. As shown in Table 1, Scan Station 3 moved approximately 7.72 mm, 0.28 mm and 2.34 mm in Subsets 1, 2 and 3, respectively. Scan Station 4 moved approximately 13.00 mm, 4.30 mm and 1.87 mm in Subsets 1, 2 and 3, respectively. The local 3D displacements of Subset 1 (7.72 mm and 10.36 mm) in the case of both Scan Stations 3 and 4 were significantly larger than the ones in Subsets 2 and 3.

Deformation Analysis Results
The first operation before comparing both point clouds in each subset was the pairwise alignment aimed at the definition of the same block structures. Since the vertical alignment was already conducted during data acquisition by levelling, a 2D transformation should be estimated on the basis of two manually selected corresponding points. These points were chosen from two opposite corners of the raw point clouds. A horizontal rotation angle of 0.036°and a translation of 2.29 m along the -axis were found. Thus, a series of profiles with respective widths of 20 cm have been directly segmented. The tunnel axis has been searched by the maximum z-value in each profile minus the height of the tunnel. Using the estimated tunnel axis and a pre-defined block size (20 cm x 40 cm), the point clouds were then segmented into 3526 blocks, with 86 rows along the profile direction and 41 columns perpendicular to it. The tunnel axis and the segmented blocks were unchanged within all the subsets, and all the subsets were geo-referenced in the geodetic reference framework.
In the cases of blocks without a sufficient number of points (N=100), the hole-filling algorithm was applied by searching the neighbouring points of each incomplete block. The numbers of blocks that have required the hole-filling procedure are shown in Table 2. The data set acquired on 3 rd December has the largest number of holes because of the lowest point density relative to all the other epochs.
To remove outliers, RANSAC algorithm was individually applied to each block. Then, by processing the remaining inlier points where there were more than 100 points within a block, coordinates of the representative point of each block and their std. devs. were estimated using the SVD algorithm. As an example, in Figure 4 the representative point coordinates and their std. devs. are reported for the data set collected on 9 th December. Here, the colour pattern changes from 0 mm to 1 mm depending on the std. devs. of the representative points. Figure 4 shows that point density influences the precision of the data, particularly in the areas with an insufficient number of points in the point cloud. Those points on the tunnel's top show larger std. devs. since the surface here was obstructed by ventilation pipes (Figure 3d). In Table 3, means and std. devs. of representative points computed after applying RANSAC filtering algorithm are compared to the ones obtained without using RANSAC.
Local 3D displacements measured by total stations 1 (3 rd vs 9 th ) 2 (9 th vs 11 th ) 3 (9 th vs 13 th ) 3D displacements of Scan Std. devs. of the representative points have generally improved after the removal of outliers based on RANSAC. Since the point density on 3 rd December was much lower than the one of other epochs, we obtain a larger number of holes to be filled. Std. devs. of the representative points relevant to that epoch were substantially larger. However, these std. devs. (in the order of millimetres) were smaller than the displacements measured by total station (in the order of centimetres) reported in Table 1. This fact indicates that the representative point can be used to interpret deformations if there are no systematic errors or network stability problems. RMSs and std. devs. of these points are shown in Table 4 for the three subsets. Considering the magnitudes of the distance discrepancies, we analysed these results, which include either deformations and systematic errors. Thus, we applied the refinement transformation to the effect of systematic errors in a subset.
After estimating and applying the refinement transformation, the RMSs and std. devs. of all the distances between the representative points were computed again (see Table 5). By comparing Tables 4 and 5, the RMSs decreased by 3.33 cm, 4.04 cm and 4.56 cm in Subsets 1, 2 and 3, respectively. However, by comparing the displacements of the scan stations measured by total station (see Table 1), the results in Table 5 still include both real deformations and residual errors.
The distances between representative points are displayed in Figure 5 for all three subsets using a colour pattern to represent the std. devs. The changes on the top of the tunnel are larger than the ones in the middle part. The most substantial variation occurs at the bottom of the tunnel walls since the point clouds collected at the bottom sector were influenced by the rougher The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XLII-3/W10, 2020 International Conference on Geomatics in the Big Data Era (ICGBD), 15-17 November 2019, Guilin, Guangxi, China surface resulting from some unremoved outliers due to the presence of some pieces of equipment and electrical wires during scanning. Indeed, results from total station observations denote that 3D displacements between scan stations are 10.36 mm, 2.01 mm and 2.11 mm in Subsets 1, 2 and 3, respectively. However, results from TLS show that the distances were 2.47 cm, 1.25 cm and 2.17 cm for the same three groups (see Table 6). Note that the TLS point clouds are measured from the top and two sides of the tunnel surface, while the total station observations are acquired on the ground surface. In tunnel monitoring tasks, the deformation on the top and side surfaces is often larger than that at the bottom of the tunnel due to pressure effects, implying that the distances between representative points in corresponding blocks as derived from TLS data analysis are often larger than the displacements between scan stations on the ground as obtained from total station measurements. Notwithstanding the application of RANSAC filtering and the refinement transformation estimation, approximately 1-2 cm systematic errors remain in the processed data.   Table 6. RMSs of the distances between representative points in corresponding blocks obtained from tls data analysis and displacements of the scan station obtained from total station observations.
Congruence tests are used to check out the significance of the derived displacements. The distances between epochs are tested to assess if the differences of the coordinates between two epochs are significant with respect to their variances. Under the assumption that the tunnel had deformations, an F-test (Wunderlich et al., 2016;Niemeier, 2008) is conducted to check the significance of the distance discrepancies at the 90% quantile with infinite degrees of freedom: If the empirical quantity exceeds the 90% quantile of the Fdistribution, the distances of the representative points significantly differ between two epochs. Otherwise, the differences are likely explained by random deviations in the measurements. The quantities of the three subsets indicate that the tunnel had significant movements. However, an adequate estimation of the covariance matrix has to be further investigated due to various influencing factors.

DISCUSSION
Limited by the operational environment, several total station measurements are processed to compare the full-section TLS results. It would be more reasonable if more targets were mounted on the top and sidewalls of the tunnel so that the estimated representative points from TLS and the artificial targets measured by total stations could be compared in a more comprehensive manner.
To fully optimise the detection of the smallest deformations, proper geodetic error propagation should be further investigated to fully describe the data quality of the TLS (see Wujanz et al., 2019). The point density and incidence angle that may influence data precision are not considered in this study. According to classical deformation analysis, within the time span of a day, the tunnel's geometry is considered stable. Strictly speaking, this hypothesis cannot hold due to the long observation time with several scans for a complete surface (Wunderlich et al., 2016).
The distances between representative corresponding blocks were used to analyse tunnel deformations, but these computed distances may include both real deformations and errors. These may be due to multiple causes, such as network geometry, the adopted registration process and changes of environmental factors. Real deformations need to be further separated from the estimated distances using statistical techniques.
The entire procedure is implemented in an automatic way except for the preliminary alignment of raw point clouds, an operation that is necessary to establish a common reference framework. To this end, two points are manually selected from the raw data. Moreover, the block size and the RANSAC input parameters must be input by the user.
The processing of the data set reported in the case study section takes approximately 18.82 minutes using an Intel(R) Core(TM) i7-4790 CPU with a 64-bit operating system and 16 GB RAM; data processing consists primarily of computing time rather than data input time.

CONCLUSIONS
We have presented a full-section deformation detection method to be used in the construction stage of a tunnel based on point clouds obtained from terrestrial laser scanning. To compare data from different days, a uniform tunnel block framework is estimated, and then segmentation is conducted to split the tunnel surface into blocks. To fill the holes and detect outliers, point clouds from neighbouring blocks and the RANSAC-based filtering method were used to estimate a representative point per each block. This point is obtained using the SVD method in a quasi-planar surface. Refinement transformation is estimated and applied to adjust systematic errors that may occur between different days. These estimated points establish explicit correspondences for comparisons between different epochs. The proposed method was applied to monitor the Badaling Tunnel, which was built for the Winter Olympics 2022 in Beijing, China. Deformations are analysed in the form of the distance discrepancies of representative points and verified by total station measurements.