Global Registration of Subway Tunnel Point Clouds Using an Augmented Extended Kalman Filter and Central-Axis Constraint

Because tunnels generally have tubular shapes, the distribution of tie points between adjacent scans is usually limited to a narrow region, which makes the problem of registration error accumulation inevitable. In this paper, a global registration method is proposed based on an augmented extended Kalman filter and a central-axis constraint. The point cloud registration is regarded as a stochastic system, and the global registration is considered to be a process that recursively estimates the rigid transformation parameters between each pair of adjacent scans. Therefore, the augmented extended Kalman filter (AEKF) is used to accurately estimate the rigid transformation parameters by eliminating the error accumulation caused by the pair-wise registration. Moreover, because the scanning range of a terrestrial laser scanner can reach hundreds of meters, a single scan can cover a tunnel segment with a length of more than one hundred meters, which means that the central axis extracted from the scan can be employed to control the registration of multiple scans. Therefore, the central axis of the subway tunnel is first determined through the 2D projection of the tunnel point cloud and curve fitting using the RANSAC (RANdom SAmple Consensus) algorithm. Because the extraction of the central axis by quadratic curve fitting may suffer from noise in the tunnel points and from variations in the tunnel, we present a global extraction algorithm that is based on segment-wise quadratic curve fitting. We then derive the central-axis constraint as an additional observation model of AEKF to optimize the registration parameters between each pair of adjacent scans. The proposed approach is tested on terrestrial point clouds that were acquired in a subway tunnel. The results show that the proposed algorithm is capable of improving the accuracy of aligning multiple scans by 48%.


Introduction
Because underground structures, such as tunnels, require routine inspections and maintenance for their optimal use, efficient and accurate tunnel inspections are mandatory. The construction of 3D models of tunnels is important for such applications. Applications of laser technology are rapidly expanding with decreased costs and increased accuracy. Therefore, 3D laser scanners make it possible to obtain point clouds with high accuracy and high spatial resolution for 3D tunnel model construction. However, one of the biggest problems that are encountered when processing these scans is terrestrial point cloud registration, in which rigid transformation parameters (RTPs) are used to align one dataset with another. Because tunnels generally have tubular shapes, the distribution of the tie points between adjacent scans is usually limited to a narrow region. Therefore, the problem of registration error accumulation becomes inevitable when multiple scans must be registered. Bergevin et al. [1] presented an algorithm that considers the network of views as a whole and minimizes the registration errors of all views simultaneously. Inspired by that work, Benjemaa and Schmitt [2] extended pair-wise registration based on a multiz-buffer technique to a global registration. They applied rigid transformations to transform each moving surface immediately after its rigid transformation had been estimated. Similarly, Sharp et al. [3] proposed an analytical method to solve for global registration parameters that involves building a graph to describe the relationship between neighboring views. This approach then decomposes the graph into basis cycles so the nonlinear optimization problem can be solved over each basis cycle in a closed form. Hu et al. [4] built a topological graph to determine the best registration path for all range scans. Stoddart and Hilton [5] identified pair-wise correspondences between points in all views and then iteratively minimized the correspondence errors over all views using a descent algorithm. This basic technique was extended by Neugebauer [6] and Eggert et al. [7] using a multiresolution framework, surface normal processing, and boundary point processing. Williams and Bennamoun [8] suggested a further refinement by including individual covariance weights for each point. There is currently no consensus as to the best approach for solving the global registration problem. Kang et al. [9] proposed a global registration method that minimizes the self-closure errors across all scans through simultaneous least-squares adjustments. Additional sensors, such as Global Navigation Satellite System (GNSS), compasses, and tilt sensors, are often combined with TLS instruments to help solve or reduce the global registration problem.
In recent years, optimization algorithms that use a series of measurements observed over time have been introduced in point cloud registration. Ma and Ellis [10] proposed the unscented particle filter (UPF) algorithm [11] to register two point data sets in the presence of isotropic Gaussian noise. Moghari and Abolmaesumi [12] proposed a registration algorithm, the unscented Kalman filter (UKF) [13], which is based on the continuous assessment of point cloud registration parameters of two rigid bodies. In this paper, we regard the point cloud registration as a stochastic system and the global registration as a process that recursively estimates the rigid transformation parameters of each scan. The augmented extended Kalman filter (AEKF) is utilized to produce accurate estimates of rigid transformation parameters by eliminating the error accumulation that is caused by the pair-wise registration. Because the central axis extracted from a single scan predictably controls more than the tie points, it is employed to reduce the error accumulation for the registration of multiple scans. Therefore, the central-axis constraints are derived as the control condition of the AEKF, so the registration parameters between each pair of adjacent scans will be globally optimized.
We begin by proposing the global registration method based on the augmented extended Kalman filter in Section 2. Section 3 optimizes the AEKF by introducing the central-axis constraint as an additional observation model. Section 4 discusses the test results, and we offer conclusions and suggestions for further research in Section 5.

Global Registration Using an Augmented Extended Kalman Filter
The Kalman filter, which is also known as linear quadratic estimation (LQE), is an algorithm that uses a series of measurements observed over time, which contain noise (random variations) and other inaccuracies, and produces estimates of unknown variables that tend to be more precise than those that are based on a single measurement.
In this paper, the point cloud registration is regarded as a stochastic system, and the global registration is the process that recursively estimates the rigid transformation parameters of each scans. Therefore, a Kalman filter is used to produce accurate estimates of the rigid transformation parameters by eliminating the error accumulation that is caused by the pairwise registration.
Because the rigid transformation model is nonlinear, we utilize the extended Kalman filter (EKF) [14], which is the nonlinear version of the Kalman filter and is the de facto standard in the theory of nonlinear state estimation, navigation systems and GPS, to estimate the six rigid transformation parameters (three for the translation and three for the rotation). As a global registration process, the RTPs that are acquired by pair-wise registrations should be globally optimized. Therefore, the system state is augmented to contain the RTPs of all pair-wise registrations that have been completed, so the optimized RTPs in the global reference frame are estimated in terms of the RTPs of the new registration and its preceding registration. This paper presents a design for an augmented extended Kalman filter (AEKF) for the global registration of tunnel point clouds.
Suppose that N scans of 3D point cloud data are expressed as: {V 1 , V 2 ,. . ., V N }. The i-th scan V i is represented as {v ij = (x ij , y ij , z ij ) | j = 1,2,. . .,M i }, and M i represents the point number of the i-th scan.
For a scale factor of 1, the rigid transformation between adjacent scans is parameterized as follows: The system covariance matrix is a symmetric matrix that is expressed as where p ij (k) is a block matrix that represents the covariance matrix between X i (k) and X j (k). System status model. Because the RTPs of each pair-wise registration are static, the system state transition equation becomes where f(.) is the state-transition model. The state-transition model is a unit matrix I, which can be ignored.

System augmented model
During the global registration process, when a new pair-wise registration is considered at time tk, its RTPs are added to the system state vector. The RTPs in the global reference frame are estimated in terms of the RTPs of the new registration and its preceding registration using Eq (5).
where X new (k) represents the augmented RTPs of the currently considered registration, X r (k −1) denotes the RTPs of the preceding registration, g(.) is a system-augmented function, the RTPs of the new pair-wise registration are T M = [ΔT X ΔT Y ΔT Z Δφ Δω Δκ], and ω(k) describes a variety of uncertainties in the pair-wise registration and modeling process, which is assumed to comply with the Gaussian distribution and is thus expressed as a white noise vector N (0, Q).

Observation model
An observation model is established to optimize the RTPs of all of the pair-wise registrations by minimizing the differences between the corresponding 3D point pairs transformed into the common reference frame. Therefore, the observation value Z is the difference between a corresponding point pair, so the observation model is derived as

Status augmentation
When a new pair-wise registration is completed, X new (k) is added to the system state vector.
The system state vector and the system covariance matrices are then augmented as follows: where X new denotes the augmented RTPs of the new pair-wise registration, X + and P + represent the a posteriori system state estimation, Xand Prepresent the a priori system state estimation, rg is the Jacobian matrix of the system augmented model, and Q is the covariance matrix of the system noise.

Observation updates
In terms of the observation model and the augmented system state vector, the new information _(k), the new information variance S(k) and the Kalman gain W are computed as follows: The state vector and covariance matrix of the system are updated as X þ ðkÞ ¼ X À ðkÞ þ W Ã _ðkÞ P þ ðkÞ ¼ P À ðkÞ À W Ã SðkÞ Ã W T g ð12Þ where X + and P + represent the a posteriori system state estimation, and Xand Prepresent the a priori system state estimation.

AEKF with Central-Axis Constraint
Because tunnels generally have tubular shapes, the distribution of tie points between adjacent scans is usually limited to a narrow region. As a result, the registration of point clouds is prone to error accumulation. The central axis that is extracted from a single scan can be more than 100 meters long, while the shift between adjacent scans in a tunnel has to be short (e.g., 20 meters) because long shifts between adjacent scans will lead to large differences between the corresponding scanning angles of incidence to the same object point, which will affect the accuracy of the identification of tie points for the registration. Therefore, the overlap between the fitted axes is expected to be much larger than that between the tie points (Fig 1), from which a constraint can be derived to control the error accumulation.
First, we propose a global extraction algorithm that is based on segment-wise quadratic curve fitting to extract the central axis of the subway tunnel. The central-axis constraint is then derived as an additional observation model of the AEKF to optimize the registration parameters between each pair of adjacent scans.

Determination of the central axis of a tunnel based on 2D projections
The central axis of a tunnel is continuously extracted using a 2D projection of the point cloud and curve fitting using the RANSAC algorithm, and the axis is optimized using a global extraction strategy that is based on segment-wise fitting.
Estimation of the boundary points. The tunnel point clouds are projected onto the XOY plane, from which we extract the boundary points of both sides of the tunnel. Therefore, the shape of the tunnel does not influence the extraction of the two parallel bounding lines on the XOY plane. An algorithm for boundary point extraction is proposed using a moving window. Fig 2 shows a circular window with a predefined radius that is centered at the point of interest P. All of the points within the window are considered as the neighboring points of point P. The polar angles of the neighboring points are computed relative to point P (e.g., α 1 ). We then calculate the differences between consecutive polar angles. If point P is a boundary point, the difference Δα i + 1,i between boundary points P i and P i + 1 is much larger than the difference Δα i,i − 1 between boundary point P i and interior point P i − 1 . Therefore, once the difference is greater than a predefined threshold, point P is labeled as a boundary point.
Fitting of the bounding lines. The bounding lines of a tunnel usually contain segments that are straight lines, curves and transition curves, which are parameterized as follows: Straight line model: Transition curve model: Curve model: where a and b are the parameters of a straight line, c, d, e, and f are the parameters of a transition curve, and g, h, and k are the parameters of a curve. The bounding line fitting process includes the estimation of multiple models. To ensure the robustness of the fit, the RANSAC algorithm [15] is used to estimate the parameters of the three models. Instead of using as much data as possible to obtain an initial solution and attempting to eliminate the invalid data points, RANSAC uses as small an initial data set as is feasible and enlarges this set using consistent data when possible. The RANSAC paradigm contains three unspecified parameters: (1) the error tolerance, which is used to determine  Extraction of boundary points using a moving window. All of the points within the moving window are considered to be neighboring points of point P. The polar angles of the neighboring points are computed relative to point P (e.g., α1). If point P is a boundary point, the difference between the consecutive polar angles of boundary points P i and P i + 1 is much greater than the difference between the polar angles of point P i and interior point P i − 1 . Therefore, once the difference is greater than a predefined threshold, point P is labeled as a boundary point. whether a point is compatible with the model; (2) the number of subsets to attempt; and (3) the threshold t, which is the number of compatible points and is used to determine that the correct model has been found. The determination of these three parameters is discussed in the introduction to RANSAC [15].
A statistical testing algorithm is proposed to implement a hypothesis testing process and automatically detect the initial models from the extracted boundary points to ensure that the proper model will be selected to fit each segment of the bounding line. The statistical test is implemented using the straight line, transition curve and curve models.
Because the mathematical models of the bounding line segments to be fitted are determinate, their parameters that are calculated by consecutive inlier sets are expected to converge in the direction of the optimal solution. Therefore, our strategy is based on a histogram to dynamically evaluate the convergence of the hypothesis models during the hypothesis testing process. Different convergent clusters of the hypothesis models will be presented in the histogram. We select the oldest model parameter set as the reference point for each convergent cluster of model parameter sets. The more the hypothesis models converge to a cluster, the more possible it is that the reference model of the cluster is correct. Therefore, we use the degree of convergence of a cluster to detect the initial models. The degree of convergence is a percentage that describes the number of models that converge to a model cluster. It is calculated by dividing the number of hypothesis models in the cluster by the total number of hypothesis models.
To determine the convergence of a newly computed model to a model cluster, we need to evaluate the deviation between the new model and the reference model of the cluster. We construct vectors with two (straight line), four (transition curve) and three (curve) dimensions for each set of model parameters. The Euclidian distances between different vectors are computed to describe the deviation between the new model and the reference model of the cluster. If the deviation is smaller than the predefined threshold, the number of FMs in the cluster is increased by one, and the degree of convergence is updated accordingly. If the new hypothesis model does not converge to any existing cluster, it will be regarded as a new cluster in the histogram.
In this method, the histogram of the candidate model parameter sets is updated during each iteration using the newly calculated hypothesis model parameter set. When the degree of convergence of a candidate model parameter set reaches a predefined threshold, the candidate model parameter set is detected as an initial model to fit the bounding line segment. If the degree of convergence fails to reach the threshold after a predefined number of iterations, we interpret that there is no such model.
To visualize the statistical results, we illustrate them as a histogram (Fig 3). The horizontal axis denotes the mean value of the model parameters, and the vertical axis represents the degree of convergence of each cell. A high degree of convergence for a parameter reflects a high probability of finding the initial model. After the initial model is detected, RANSAC is used to robustly estimate the optimized model parameters. Two, four, and three points are used to estimate the model parameters to fit a straight line, a transition curve and a curve, respectively. The criterion that is used to identify outliers is based on the deviations of the tested points from the fitted model. The inlier bounding points of a certain model are classified as a segment that is used in the following global optimization. The final optimal parameters are computed by a least-squares adjustment using the obtained inlier points.
Fitting of the central axis. After fitting the bounding lines, the boundary points are evenly resampled. To extract the central axis of the tunnel, the normal vector V l of the left bounding curve at boundary point P l is determined (Fig 4). A straight line orthogonal to the normal vector reaches the right bounding curve from P l and generates point P l '. Theoretically, the radial line from point P l ' that is orthogonal to V l ' reaches the left bounding curve at point P l , so the extracted central-axis point is the midpoint of the line P l P l '. However, because the bounding curves are subject to errors that are generated from the fitting processes, the radial orthogonal to V l ' produces point P l " instead of point P l . Fig 4 shows that M l ' and M l " are the midpoints of P l P l ' and P l ' P l ", respectively. The extracted central-axis point is determined as M l , which is the average of points M l ' and M l ". The same process is implemented from boundary point P r on the right bounding curve to extract the point on the central axis as point M r . The presented strategy to fit a bounding line is used to generate the central axis based on the extracted central-axis points. Because the extraction of the central axis is implemented on the XOY plane, the height of the central axis is determined as the mid-height of the tunnel points.

Global adjustment of the central axis using segment-wise fitting
Because the extraction of the segments of the bounding lines and the central axis on the XOY plane using the three models may be affected by noise in the tunnel points, there may be deviations in the overlapping parts of adjacent fitted models (Fig 5). Therefore, we propose a global extraction algorithm to minimize the deviations.
To maintain consistency between adjacent fitted models, the divided segments overlap each other by a small amount, and a global least-squares adjustment is developed to implement the multiple model fitting of all of the segments together by minimizing the deviations in the overlapping parts of adjacent fitted models. Using Eqs (13)- (15), the constraints are derived between a straight line, a transition curve and a curve, respectively, and are added to the adjustment model. For example, Eq (16) parameterizes the constraint between a straight line and a transition curve: where a i and b i are the line parameters of segment i, c j , d j , e j , and f j are the transition curve parameters of segment j, and Y is the Y coordinate of a point in the overlap region between segments i and j.  To extract the central axis of the tunnel, the normal vector V l of the left bounding curve at boundary point P l is determined. A straight line orthogonal to the normal vector reaches the right bounding curve from P l and generates point P l '. Theoretically, the radial line from point P l ' that is orthogonal to V l ' reaches the left bounding curve at point P l , so the extracted central-axis point is the midpoint of line P l P l '. However, because the bounding curves are subject to errors that are generated from the fitting process, the radial that is orthogonal to V l ' produces point P l " instead of point P l . M l ' and M l " are the midpoints of P l P l ' and P l ' P l ", respectively. The extracted central-axis point is determined as M l , which is the average of points M l ' and M l ". The same process is implemented from boundary point P r on the right bounding curve to extract the point on the central axis as point M r . Based on the extracted central-axis points, the presented strategy to fit a bounding line is used to generate the central axis.
which are derived from Eq (17) for the overlap region between segment j and segment j + 1.
The form of C ij depends on the models of the two overlapping segments. In the proposed global least-squares adjustment system, Eq (16) is used as a constraint equation and is weighted with a large value (e.g., 10) instead of with 1, as occurs in an observation equation. Based on the coefficient matrix B, we calculate the optimized parameters of the bounding line segments by following the least-squares strategy. After the bounding lines are fitted, the method presented in Section 3.1 is implemented to extract the central-axis points, which we use to generate the globally optimized central axis using the proposed global least-squares adjustment system.

Observation model derived from the central-axis constraint
The extracted central axis comprises multiple segments that can be parameterized using Eqs (13)- (15). We derive the central-axis constraint by adding an additional observation model into the AEKF system. The observation model is established to minimize the deviation that describes how the point on the axis segment in the analyzed scan does not fit the corresponding axis segment model that was determined in the fixed scan (Fig 6), when the point is transformed into the global reference frame. Therefore, the observation value Z is computed as (for the straight line model): where (a, b) represent the linear parameters, Z mid denotes the middle height of the tunnel points in the overlap, h(.) is the observation model, and υ(k) denotes a variety of uncertainties in the scanning measurement and the transformation of coordinates, which is supposed to comply with the Gaussian distribution and is thus expressed as a white noise vector N (0, R) The observation model is linearized as: where rh A is the Jacobian matrix derived from Eqs (21) and (22)). Because the observation model is established to minimize the deviation that describes how the point on the axis segment in the analyzed scan does not fit the corresponding axis segment model that was determined in the fixed scan, the selection of the number of observation equations is flexible in the proposed algorithm. The more observation equations that are added to the AEKF system, the greater the contribution that the constraint has on the registration. Thus, a constraint that is derived from a low-quality central-axis segment is assigned a small number of observation equations in our algorithm.
Therefore, the observation model of the AEKF system is modified as: The RTPs of all of the pair-wise registrations are then optimized by minimizing both the differences between the corresponding 3D point pairs and the point-to-model deviation, which is expected to improve the robustness and accuracy of the proposed global registration approach. The improvement can be attributed to the multiple overlaps of the central axis among the scans, while we are most likely able to find the tie points only from the overlaps between two consecutive scans.

Experimental Results
The proposed approach was tested on real datasets (Fig 7) that were acquired by a RIEGL LMS VZ-400 laser scanner in a subway tunnel in Shanghai, China. Twelve scans were captured with an average shift of 10 m between the scanning centers. The scans cover tunnel segments with different shapes for a distance of 450 m. The scan's angular resolution was 0.046°, and the range accuracy was ±5 mm.  Table 1 lists the accuracies of the pair-wise registrations; the average accuracy is 0.021 m.

Pair-wise registration
As shown in Fig 7, the tunnel has a tubular shape, so the average length of the areas covered by the tie points is only 30 m, while the mean length of the tunnel segment that is captured by one scan is 150 m. Therefore, the accumulation of registration error is inevitable (Fig 9a). To estimate the error accumulation, Fig 9b-9g illustrate cross sections that were extracted at the same position from scans 1 and 4, scans 2 and 5, scans 3 and 6, scans 7 and 10, scans 8 and 11, scans 8 and 12, respectively, which were transformed into the same coordinate system using the pair-wise registration results. The distinct deviations between the two cross sections in Fig  9b-9g are due to the accumulation of the pair-wise registration errors.
Twenty check point pairs were selected from the two cross sections in Fig 9b-9g. Table 2 list the differences. The average deviation is 0.0463 m, which is larger than the average value of 0.022 m in Table 1.
As proposed in Section 2, we implemented the augmented extended Kalman filter to eliminate the error accumulation that was caused by the pair-wise registration to accurately estimate the rigid transformation parameters.

Global registration using the augmented extended Kalman filter
The pair-wise registration results were used to construct the augmented extended Kalman filter system with which the global registration was implemented. To compare the registration accuracies, cross sections were also extracted at the same positions from scans 1 and 4, scans 2 and    5, scans 3 and 6, scans 7 and 10, scans 8 and 11, and scans 8 and 12 (Fig 10). Fig 10 shows that the deviations between the two cross sections are less than the deviations shown in Fig 9; the average deviation computed from the results in Table 3 decreases to 0.0301 m.
Although the accumulation of errors of the pair-wise registration was reduced somewhat, the AEKF system was still implemented based on tie points that were shown to have a limited distribution. Therefore, as presented in Section 3, the constraint derived from the overlapping axes was introduced into the AEKF system to globally optimize the transformation parameters. Augmented extended Kalman filter with the central-axis constraint To derive the central-axis constraint, the central-axis segments were first extracted from the twelve scans.
Fitting of the central axis based on the 2D projection. As proposed in Section 3.1, the tunnel points of the dataset (e.g., scans 1 and 12) were projected onto the XOY plane as shown The bounding lines of a tunnel may comprise segments of straight lines, curves and transition curves. The proposed statistical testing algorithm was implemented using the three models to automatically detect the corresponding initial model from the extracted boundary points. Fig 12 shows the peaks in the histograms for the straight lines, transition curves and curves of the extracted boundary points that were fit using the three models.
Based on the detected initial models (e.g., the straight line and transition curve), we fit the bounding lines using RANSAC (Fig 13a and 13c). In Fig 13c,  11a were eliminated (highlighted in green) by RANSAC. The boundary points were evenly resampled in terms of the fitted model parameters. The method presented in Section 3.1 was implemented using these fitted parameters to extract the central-axis points (Fig 13a and  13c). Fig 13b and 13d show 3D views of the central axes that were generated from the centralaxis points. As shown in Fig 14, the central axis that was fit from the tunnel points consists of three segments. The yellow box on the left highlights the overlap between the curve and the transition curve, while the box on the right shows the overlap between the transition curve and the straight line. To test the fitting accuracy, we set 24 Y coordinates along the central axis within the overlap zones highlighted in the two yellow boxes and computed their corresponding points on the two adjacent segments. The deviations between the corresponding points were then calculated and are shown in Table 4. Large deviations with an RMSE of 26 mm (detailed views) are present within the overlap zones between the segments due to the noise in the tunnel point dataset.

the noise points shown in Fig
Global extraction of the central axis using segment-wise fitting. To optimize the extraction results, the global least squares adjustment proposed in Section 3.2 was implemented to  To test the fitting accuracy, we set 24 Y coordinates along the central axis within the overlap zones highlighted in the two yellow boxes. Their corresponding points on the two adjacent segments were computed. The deviations between the corresponding points were then calculated and are shown in Table 4. The RMSE of the deviations was reduced from 26 mm to 2 mm by the global extraction process.  Augmented extended Kalman filter with segment-wise axes. As presented in Section 3.2, the central axes were extracted by a global extraction strategy using segment-wise fitting, which produced many axis segments. Those segments were employed to construct an additional observation model according to Section 3.3 to control the error accumulation. Fig 16 shows that the deviations between the two cross sections are less than those shown in Fig 10. The average deviation computed from the results in Table 5 decreases to 0.024 m, which is only half of the mean value of 0.0463 m that is calculated from Table 2.

Conclusions
In this paper, we proposed a global registration approach that is based on an augmented extended Kalman filter and central-axis constraints. The point cloud registration was regarded as a stochastic system, so we utilized AEKF to produce accurate estimates of the rigid transformation parameters by eliminating the error accumulation caused by the pair-wise registration. Moreover, the central axis that was extracted from the scan was used to control the registration of multiple scans. The central axis of the subway tunnel was determined through a global extraction algorithm based on segment-wise quadratic curve fitting. We then derived the centralaxis constraint as an additional observation model of AEKF to optimize the registration parameters between each pair of adjacent scans.
The proposed algorithm was implemented using a terrestrial laser scanning dataset that was acquired in a subway tunnel. The experimental results show that when multiple scans are aligned into a common coordinate frame, the consecutive implementation of pair-wise registration causes error accumulation (from 0.022 m to 0.046 m). The results also illustrate that the application of the global registration based on AEKF can reduce the error accumulation (from 0.046 m to 0.030 m). The proposed algorithm of global extraction based on segment-wise fitting was implemented on the experimental point clouds to impose the central-axis constraint and achieved an accuracy of 2 mm. The presented AEKF with the segment-wise axis constraint was shown to improve the accuracy of aligning multiple scans by 48% (0.024 m versus 0.046 m).
Because the extraction of the central axis is important to the global registration method, future work will focus on improving the robustness and applicability of our algorithm for the extraction of the tunnel's central axis. Moreover, tunnels with more complex shapes and sharper curves will be considered in future work.