Iterative automatic global registration algorithm for multi-view point cloud of underground tunnel space

Abstract

Aiming at the narrow and long tunnel structure, few internal features, and a large amount of point cloud data, the existing registration algorithms and commercial software registration results are not ideal, an iterative global registration algorithm is proposed for massive underground tunnel point cloud registration, which is composed of local initial pose acquisition and global adjustment. Firstly, the feature point coordinates in the point cloud are extracted, and then the station-by-station registration is performed according to the Rodrigues matrix. Finally, the registration result is considered as the initial value of the parameter, and the global adjustment of all observations is carried out. The observation values are weighted by the selection weight iteration method and the weights are constantly modified in the iteration process until the threshold conditions are met and the iteration stops. In this paper, the experimental data, made up of 85 stations of point cloud data, are from the Xiamen subway tunnel, which is about 1300 m long. When the accumulated error of station-to-station registration is too large, several stations are regarded as partial wholes, and the optimal registration is achieved through multiple global adjustments, and the registration accuracy is within 5 mm. Experimental results confirm the feasibility and effectiveness of the algorithm, which provides a new method for point cloud registration of underground space tunnel.

Introduction

Compared with other technologies, LIDAR (Light Detection and Ranging) technology has the advantage of not being effected of environmental conditions and can carry out accurate measurement even in a dark tunnel, which makes it a preferred means for data acquisition in underground tunnels. ¹ Since LiDAR technology can quickly obtain high-precision tunnel point cloud data, characterized by high efficiency and safety, it is widely used in line adjusting and slope fine-tuning, crack extraction, ² and tunnel monitoring. ³ When using LiDAR technology to collect a tunnel with a large-span linear structure, it is usually scanned from different angles and positions, and then obtains the complete tunnel point cloud through data registration, which is an essential step for subsequent point cloud processing and analysis. ⁴ Therefore, the research on automatic registration method of the underground tunnel point cloud is of great significant, especially in the rapid development period of tunnel engineering. ²

The current multi-view point cloud registration methods basically fall into point-set-based registration and geometric feature-based registration ones. ⁵ One of the classic point-set-based registration methods is the ICP (Iterative Closest Point) algorithm, which calculates the spatial transformation parameters of the point cloud by the nearest neighbor pairs in adjacent point cloud, and the search and iterative calculation are continued until the convergence condition is met. ⁶ However, the calculation speed of the ICP algorithm is slow and the requirement on the initial poses of the point cloud is also strict. ⁷ To overcome these limitations, a lot of optimization algorithms have been proposed, such as KD-tree (K-Dimensional tree), ⁸ point-to-plane, ⁹ point-to-projection, ¹⁰ and other algorithms to speed up the search of the nearest neighbors, as well as methods to simplify the point set of matching objects and reduce the amount of calculation data. ¹¹ Geometric feature-based registration is to perform rigid body transformation by using more than three pairs of points, lines, faces, and other features with the same name. The feature point-based registration method, a currently most commonly used one, extracts geometric features such as normal, curvature, and image features, or directly extracts the artificial marks from the point cloud. Generally, the registration process of two point sets can be divided into coarse registration and fine registration. ¹² Yin et al. ¹³ proposed that combines coarse point cloud registration using SIFT (scale-invariant feature transform) feature points and fine registration based on the KD tree improved ICP algorithm. In the coarse registration process, the initial value of the point cloud position was obtained, and then the fine registration was used to improve the registration accuracy. Huang et al. ¹⁴ presented an automatic registration algorithm on the basis of point feature histograms, using the difference of point feature histogram to screen out the transformation matrix with the smallest error, and then further fine registration was carried out based on the ICP algorithm. Although the accuracy is satisfactory, most registration methods based on ICP algorithm have the problem of local convergence, which is likely to trapped into the local minimum due to the poor initial parameters.

The above-mentioned methods are exclusive for registration between two-view point clouds. When it is used to register multi-view point clouds, it repeatedly registers two adjacent point sets and integrates them into one point set until the registration is completed. Considering that the accumulative error in the sequential registration of adjacent point sets, researchers have proposed a strategy of multi-view point clouds global registration. Chen et al. ¹⁵ proposed an accurate global registration method. Firstly, the initial transformation matrix was determined by the regional descriptor based on OBB (oriented bounding box), and then it was refined by the ICP algorithm, and the scanning was accurately registered in the object model. In Li et al. ¹⁶ a global optimization registration algorithm of multi-view point cloud was reported by Li et al. After coarse registration, a KD tree was used to search for the nearest pair of overlapping regions, and then the adjustment model was used to calculate the optimal spatial transformation parameters of each view point cloud. In registration, overlapping area and approximate point pair of the same name need to be searched, so the calculation efficiency was slightly lower. Evangelidis and Horaud ¹⁷ proposed to transform the multi-view registration problem into a clustering problem, and estimated the Gaussian mixture model by using Expectation Maximization algorithm to realize the calculation of clustering and registration parameters, but there were many parameters to be solved in this method. Wang et al. ¹⁸ developed a global registration algorithm by a low-rank sparse decomposition matrix, which was verified with point cloud data from experimental data sets. Sanchez et al. ¹⁹ suggested a method for global registration in a structured scene indoors, using Gaussian images of point cloud and histograms to calculate rotation and translation parameters, respectively. Nevertheless, this method requires rich geometric features in the indoor scene. By analyzing the error accumulation problem of the multi-view point cloud registration, it is concluded that the accuracy of global adjustment results will be significantly improved. ²⁰

The above global registration methods have made great progress in the field of data registration, but the registration of large-scale tunnel point cloud is unsatisfactory due to the special spatial structure of the long and narrow tunnel. Moreover, there are no obvious features inside the tunnel, and the two point sets are tremendously similar, as shown in Figure 1, extracting features from the inside of the tunnel is a key problem. Additionally, there are tens of millions of scanning points in a single scanning station, while the complete tunnel often has dozens or even hundreds of stations, so the registration of massive point cloud of multiple stations is also a key problem. Regarding the registration of tunnel point cloud, Gosllga et al. ²¹ fitted the best cylinder and extracted its axis to the circular cross-section tunnel data, and rotated the multi-period point cloud to the same reference frame by the coincidence axis. However, it was difficult to extract the tunnel axis, therefore the target for registration was recommended, for example, Yue et al. ²² arranged a black and white board target on the scene, extracted its center point as the control point of the registration coordinate system. As the planar target was affected by the scanning angle, the center point cannot be extracted when scanning from the side. In view of the narrow and long tunnel structure, few internal features, and a large amount of point cloud data, an iterative global registration algorithm is proposed by using target ball which is not limited by any scanning angle. The center coordinates of the target ball are used as features point for registration, and then the registration results are taken as the initial values for global adjustment. The posterior variance selection weight iterative method is used to assign different weights to the feature constraints, and the global optimal registration value is obtained iteratively.

Figure 1.

Point cloud of the tunnel: (a) outside the tunnel and (b) inside the tunnel (both the internal and external structure are very similar).

The algorithm takes the center of the target ball as the geometric feature for registration, rather than searching the nearest point. Its advantage is that it does not require a small initial position between the two point sets, and avoids massive tunnel point clouds participating in the operation, which is the biggest difference from the iterative registration method based on points in one set to the closest points set. For example, the Trimmed iterative closest point algorithm based on Least Trimmed Squares approach, ²³ the curve registration algorithm combining €-reciprocal correspondence with the robust Least Median of Squares motion estimation, ²⁴ and the nearest point iterative registration algorithm based on the distance distribution statistics. ²⁵ In contrast, the proposed strategy has higher registration efficiency, but the cost is that the accumulative error has a significant impact on the registration results, so the global adjustment is carried out, and the proper weights of these features are calculated in each iteration to decrease the error accumulative. In the iteration process, a posterior variance selection weight function is incorporated to detect the gross error. It is considered that the gross error is a random model with abnormally large variance, and the outliers are dynamically weighted to reduce the weight continuously. Luck et al. ²⁶ used a similar strategy in the simulated annealing algorithm, weighting based on the median squared distance from each data point to the surface model, and choosing a good starting point for ICP through SA. However, the limitation is that the error surface is relatively flat with many local minima. In contrast, the proposed method has strong ability of gross error detection, can effectively reduce the accumulative error, obtain reliable registration results, and has high efficiency and high accuracy in massive tunnel point cloud registration. When the error accumulated of point cloud registration in large-scale underground tunnel scene is too large, the local stations of the tunnel are treated as a whole and registered into the reference coordinate system segment by segment. High precision tunnel registration results provide basic data for tunnel construction and monitoring applications and contribute to better serve the urban underground space data management. Consequently, it is of great practical significance to study the registration method of the tunnel point cloud to provide a reliable and feasible solution for tunnel data registration in underground space.

The rest of this paper is arranged as follows: firstly, the iterative global registration algorithm is introduced. Secondly, the registration experiment with actual tunnel data is described, and the results are discussed. Finally, a brief summary and prospect for future research are given.

The iterative global registration algorithm

Iterative global registration consists of three processes: data preprocessing, initial value parameters solution and global adjustment. The available features in the point cloud are first extracted to prepare for registration. When the registration base station is determined, the rotation matrix of each site and the coordinates of the same point are calculated as the initial parameters of the overall adjustment. Starting from the base station, search for adjacent feature points with the same name, register each site cloud to the base station through the Rodrigues matrix registration, and gradually expand the base station outwards. When more registration sites are involved, it is more likely to cause error accumulation, ²⁷ and the accuracy of station-to-station registration will correspondingly become lower. In this case, the global adjustment method is used to reduce the accumulative error, and all the spatial transformation parameters and the unknown points adjustment values are calculated. When the error margin is less than the specified threshold, the registration result is output. If the error is too large, the weight of each constraint is recalculated by the weight function, and the gross error is weakened or eliminated in the iterative calculation until the accuracy requirements are met, then the iteration is stopped and the registration point cloud is output the process as shown Figure 2.

Figure 2.

Iterative global registration flow chart, the whole process can be divided into three parts: data preprocessing, initial value solution, and iterative global adjustment (the red color in the second image represents the feature sphere and the green point cloud in the third image represents the base station).

Local coarse registration

Local coarse registration uses the feature constraints between the base station and the registration station and performs the station-to-station registration through the Rodrigues matrix. It uses three anti-symmetric elements instead of Euler angle to establish coordinate transformation model and then solves parameters, scale parameters, rotation parameters, and translation parameters respectively. By centering the coordinates, the translation parameters and rotation parameters are solved separately to reduce the ill conditioned severity of the normal equation. ²⁸ The anti-symmetric matrix, which is composed of three independent parameters, constructs the Rodrigues matrix as follows ²⁹ :

$$R={(I-S)}^{-1}(I+S)$$

(1)

Among them, I is an identity matrix, S is an anti-symmetric matrix composed of parameters a, b, and c.

$$S=\left[\begin{array}{ccc}0\hfill & -c\hfill & -b\hfill \\ \multicolumn{1}{c}{c}& 0\hfill & -a\hfill \\ \multicolumn{1}{c}{b}& a\hfill & 0\hfill \end{array}\right]$$

According to the principle of coordinate transformation, three pairs of features with the same name that are not on a straight line in space can be geometric features such as points, lines, and faces, and then the space transformation parameters can be solved.

In the experiment, we take the target sphere center as the feature, and the feature points of the same name $X_0=(x_0,y_0,z_0)$ and $X=(x,y,z)$ in the two point sets have the following relationship:

$$X0-(\lambda \mathrm{RX}+\mathrm{\Delta }X)=0$$

(2)

Among them, $\lambda$ is the scale parameter, $\mathrm{\Delta }X$ is the offset, $X\mathrm{m}$ and $X\mathrm{n}$ are the center of gravity coordinates, the scale is unchanged in the point cloud transformation, so $\lambda$  = 1.

Order,

$$\overline{X0}=X0-\mathrm{Xm},\overline{X}=X-\mathrm{Xn}$$

After centering, the coordinate transformation between the points with the same name is as follows (3):

$$\overline{X0}-R\overline{X}=0$$

(3)

The formula (3) is linearized to obtain the error formula of rotation parameters, as shown in formula (4):

$$V=\mathrm{At}-L$$

(4)

Among,

$$\begin{array}{c}\multicolumn{1}{c}{A=\left[\begin{array}{ccc}0\hfill & -\overline{z0}-z\hfill & -\overline{y0}-y\hfill \\ \multicolumn{1}{c}{-\overline{z0}-z}& 0\hfill & \overline{x0}+x\hfill \\ \multicolumn{1}{c}{\overline{y0}+y}& \overline{x0}+x\hfill & 0\hfill \end{array}\right]t}\\ =\left[\begin{array}{c}a\hfill \\ \multicolumn{1}{c}{b}\\ \multicolumn{1}{c}{c}\end{array}\right]L=\left[\begin{array}{c}\overline{x0}-x\hfill \\ \multicolumn{1}{c}{\overline{y0}-y}\\ \multicolumn{1}{c}{\overline{z0}-z}\end{array}\right]\hfill \end{array}$$

The error equation of the translation parameter is expressed by the following (5):

$$V=\mathrm{\Delta }X-(X0-\mathrm{RX})$$

(5)

Global adjustment

The global adjustment is solved by the beam method adjustment model, and the result of station-to-station registration is used as the initial value to set up the global error equation. In the equation, the characteristic points $X\mathrm{t}0=(\mathrm{xt}0,\mathrm{yt}0,\mathrm{zt}0)$ and their observation values $X\mathrm{t}=(\mathrm{xt},\mathrm{yt},\mathrm{zt})$ are related in the following manner:

$$X\mathrm{t}0-(\lambda \mathrm{RX}\mathrm{t}+\mathrm{\Delta }X)=0$$

(6)

Among them, $\lambda$  = 1 is the scale parameter and $\mathrm{\Delta }X$ is the offset. It is linearized and expanded to get the error equation (7)

$$V1=A1t+\mathrm{BX}-L1$$

(7)

Among, $t={\left[\begin{array}{cccc}\begin{array}{cc}\begin{array}{cc}\mathrm{d}\mathrm{\Delta }x\hfill & \mathrm{d}\mathrm{\Delta }\mathrm{y}\hfill \end{array}\hfill & \mathrm{d}\mathrm{\Delta }\mathrm{z}\hfill \end{array}\hfill & \mathrm{d}\mathrm{\Delta }a\hfill & \mathrm{d}\mathrm{\Delta }\mathrm{b}\hfill & \mathrm{d}\mathrm{\Delta }\mathrm{c}\hfill \end{array}\right]}^T$ is the space transformation parameter correction number, $A1$ is the transformation parameter coefficient matrix, $X={\left[\begin{array}{ccc}\mathrm{d}x0\hfill & \mathrm{d}y0\hfill & \mathrm{dz}0\hfill \end{array}\right]}^T$ is the number of corrections to be fixed, $B$ is the undetermined coefficient matrix, and $L1$ is the observation residual. The matrix method equation is shown in the following formula (8):

$$\left[\begin{array}{cc}A{1}^T\mathrm{PA}1\hfill & A{1}^T\mathrm{PB}\hfill \\ \multicolumn{1}{c}{B^T\mathrm{PA}1}& B^T\mathrm{PB}\hfill \end{array}\right]\left[\begin{array}{c}t\hfill \\ \multicolumn{1}{c}{X}\end{array}\right]-\left[\begin{array}{c}A{1}^T\mathrm{PL}1\hfill \\ \multicolumn{1}{c}{B^T\mathrm{PL}1}\end{array}\right]=0$$

(8)

It is abbreviated as follows:

$$\left[\begin{array}{cc}N11\hfill & N12\hfill \\ \multicolumn{1}{c}{N21}& N22\hfill \end{array}\right]\left[\begin{array}{c}t\hfill \\ \multicolumn{1}{c}{X}\end{array}\right]-\left[\begin{array}{c}L1\hfill \\ \multicolumn{1}{c}{L2}\end{array}\right]=0$$

By solving the algorithm equation, the space transformation parameters, and the correction value of the fixed point can be obtained.

$$t={(N11-N12N{22}^{-1}N21)}^{-1}(L1-N12N{22}^{-1}L2)$$

$$X=N{22}^{-1}(L2-N21t)$$

Through the least-square adjustment, the posterior variance is obtained as follows, where r is the redundant observation.

$${\stackrel{\mathrm{\Lambda }}{\sigma }}_0^2=\frac{V^T\mathrm{PV}}{r}$$

(9)

The iterative global registration algorithm

Once a large observation error occurs in the feature constraints in the registration station, the accuracy of the adjustment will be insufficient. The selection weight iteration method is used for gross error detection, and the weight of large error constraint is gradually reduced through iteration to reduce or eliminate the influence of gross error. ³⁰ The iterative global registration algorithm constrains the weighting and solving of the observations continuously, and controls the error of the observation corrections within a certain threshold range until the registration is completed. The iterative process is as follows:

• (1) Initial value of the local coarse registration solution parameter: the initial value is the approximate value of each constraint observation in the base station coordinate system, which is provided by Rodrigues registration from station to station;

• (2) Perform global adjustment;

• (3) Modify the weight matrix P: when there is a significant difference between the variance of the observation value and the variance of the observation value of this group, there may be gross errors. According to the weight function (10) of the posterior variance selection weight iteration method, the observation value exceeding the threshold value is re-weighted, and the weight matrix P is replaced by P1.

$$P_{i,j}^{v+1}=\{\begin{array}{c}{p}_i^{v+1}=\frac{{\stackrel{\mathrm{\Lambda }}{\sigma }}_0^2}{{\stackrel{\mathrm{\Lambda }}{\sigma }}_i^2},T_{i,j}<F_{a,1,\mathrm{ri}}\hfill \\ \frac{{\stackrel{\mathrm{\Lambda }}{\sigma }}_0^2}{{\stackrel{\mathrm{\Lambda }}{\sigma }}_{i,j}^2}=\frac{{\stackrel{\mathrm{\Lambda }}{\sigma }}_0^2\mathrm{ri},j}{V_{i,j}^2},T_{i,j}\ge F_{a,1,\mathrm{ri}}\hfill \end{array}$$

(10)

Among them, $T_{i,j}=\frac{{\stackrel{\mathrm{\Lambda }}{\sigma }}_{i,j}^2}{{\stackrel{\mathrm{\Lambda }}{\sigma }}_i^2}$ and $F_{a,1,\mathrm{ri}}$ are inspection volume, $F_{a,1,\mathrm{ri}}$ is generally taken as 4.13, which is equivalent to the situation of significance level $\alpha =0.1\mathrm{\%}$ and efficacy $\beta =80\mathrm{\%}$ .

• (4) Repeat the previous two steps: return the solution results from t, “X” and weight matrix “P” to formula (8) and re-calculate the global adjustment, and use the weight function to modify the weights until the maximum correction number of the observed values is less than the threshold, that is, then stop the iteration immediately as shown in Figure 3.

Figure 3.

Iterative global registration flow chart.

Comparison of methods

Take the point cloud of the three stations tunnel as an experimental data for registration, as shown in Figure 4, showing the spatial location relationship before registration. Among them, the three site clouds are represented by different colors of light blue (LB), red (R), and dark blue (DB) respectively, and the light blue point cloud is taken as the reference station.

Figure 4.

Tunnel point cloud before registration.

Through feature registration, global adjustment, iterative global adjustment, and ICP algorithm, the results are shown in Figure 5, showing the point cloud after registration and the cross-section point cloud at a certain position, which is convenient for viewing the stratification. In addition, Table 1 shows the detailed data of the control point coordinates and registration errors in each station. The results of feature registration are shown in Figure 5(a), the RMS of the adjacent scanning stations are 10 mm and 3 mm respectively, with slight stratification. The registration results of global adjustment are shown in Figure 5(b), the error of global adjustment is obviously lower than that station-to-station feature registration, the RMS is 4 mm and 2 mm respectively, and there is almost no stratification in the section. The iterative global adjustment algorithm is used to register the above three stations, and the large error constraints are reweighted. As shown in Figure 5(c), the RMS is 2 mm and the registration effect is excellent. The results of the coarse registration of feature and fine registration of the ICP algorithm are shown in Figure 5(d). The RMS is 2 mm and there is no stratification in the point cloud section. The registration accuracy of the ICP algorithm is equivalent to that of the global iterative registration algorithm, but the real-time calculation time of the ICP algorithm is relatively long. A similar conclusion can be drawn from the comparison, the proposed method has advantages of registration accuracy and registration time in processing massive tunnel point cloud data.

Figure 5.

Registration results of tunnel point cloud: (a) feature registration results, with slight stratification, (b) global adjustment registration results, basically no stratification, (c) iterative global adjustment registration results, without stratification, and (d) ICP registration results, without stratification.

Table 1.

Original coordinates of target control points and station by station registration error.

	ID	X, Y, Z	X, Y, Z	Feature registration/mm			Global registration/mm			Iterative global registration/mm			ICP registration/mm
				W	Error	RMS	W	Error	RMS	W	Error	RMS	W	Error	RMS
LB-R	1	−9.336, 1.441, 0.259	1.359, −9.987, −0.126	1	8	10	1	7	4	0	7	2	1	4	2
	2	−15.82, 4.512, 0.148	1.355, −2.799, −0.244	1	13		1	4		0.36	4		1	2
	3	5.045, −1.173, 0.399	5.149, −24.096, 0.023	1	10		1	2		1	2		1	2
	4	−1.306, 2.310, 0.934	5.580, −16.863, 0.556	1	7		1	2		1	2		1	1
R-DB	5	−1.716, 4.956, −0.300	16.728, −9.505, 1.409	1	3	3	1	2	2	1	2	2	1	2	2
	6	−1.220, 12.523, 0.194	11.554, −3.957, 1.898	1	4		1	2		1	2		1	3
	7	−6.649, 21.794, −0.096	1.086, −1.552, 1.602	1	2		1	1		1	1		1	2

LB, R, and DB represent the scanning station of the corresponding color, ID represents target ball serial number, W represents the weight, Error represents registration error of each point, RMS represents the total root mean square error of station, and the error values are all accurate to the millimeter level.

Experiment and discussion

In this paper, the experimental data, made up of 85 stations of point cloud data by terrestrial laser scanner, are from the Xiamen subway tunnel, which is about 1300 m long and passes a turning with a radius of 650 m. According to the principle of iterative global registration, the algorithm is implemented on the visual studio platform by the C # programming language. The program is used to read the spherical coordinates of all stations for registration. The experiment prevents a large number of points from being involved in the calculation, improves the calculation speed and yields better registration results, which testifies to the feasibility of this method.

Data preprocessing

During acquisition, target balls are set as feature markers. On two different stations, there are at least three common targets that are not on a straight line, and four common targets that are not on the same plane. Through the improved RANSAC algorithm combined with the linear least square method, point clouds are extracted and spherical surfaces are fitted. The selection method of the initial iteration point in the RANSAC algorithm is improved. The star structure radiating from one point to eight directions is taken as the arrangement shape of seed points, which are used for fitting. After the spherical parameters are solved for the first time, according to the scanning density of space points, the region segmentation of depth image data is automatically carried out through iteration. Based on the distance threshold and qualified points, the linear least square method is used to fit the segmented quadric surface accurately.^31,32 The tunnel interior is shown in Figure 6, there are almost no other features except track and box. From this point of view, the result of algorithm extraction is reliable, and there is almost no interference from other spherical features. The following Figure 7 shows the process of automatically extracting the target ball from three different scenes: rails, wooden box, and tunnel ground. The spherical point (b) is detected from the tunnel point cloud (a) and then fitted to the target ball (c). With five to seven target balls arranged in each station, after the target balls are extracted from each station, the center coordinates and labels of the balls are recorded to prepare for registration.

Figure 6.

Scanning scene and scanning results.

Figure 7.

Automatic extraction of target balls: The figure shows the extraction process of target balls in three different scenes in the tunnel. The target balls are located on the rails, wooden boxes, and the ground, where (a) represents the original point cloud in the tunnel. The red dots in (b) indicate the detected spherical point cloud, and (c) represents the fitted sphere.

Registration results and accuracy analysis

After reading the spherical center coordinates of all stations at one time through the program, the RMS error of the local coarse registration by station-to-station is shown in Figure 8(a). The coarse registration RMS error of the former part station is less than 1 cm, and the maximum float is no more than 5 cm, but the accumulated RMS error from the 45th station is grave, whose maximum RMS error reaches 6.4 m. Given that the coarse registration RMS error of the latter part of the station is too large, and the result of the global adjustment directly does not converge, all the stations are divided into two parts. For the first 44 stations, the iterative global adjustment is performed directly and the result is shown in Figure 8(b). After adjustment, the RMS error of each station is effectively reduced, being no more than 5 mm. For other stations, the method of segmented adjustment is adopted. The tunnel is segmented according to the station sequence number, and the two segments are connected by one overlapping station point cloud. With the overlapping station selected as a temporary base station, the small segments are first coarsely registered, and the error is controlled at the centimeter level, as shown in Figure 8(c). Then the local small segment station is adjusted as a whole, and adopted is the optimal solution for each segment iterative global registration. According to the station serial number, the later segment point cloud is registered to the previous segment point cloud one by one, and then converted to the base station coordinate system. The iterative global registration RMS error of all stations is shown in Figure 8(d).

Figure 8.

Line chart of registration RMS error: (a) coarse registration RMS error of all stations, (b) global registration RMS error of former part stations, (c) coarse registration RMS error of latter part stations, and (d) iterative global registration RMS error of all stations. In the figure, the abscissa is the station ID, the ordinate is the registration RMS error, the polyline indicates the change trend of the RMS error, and the dots represent the corresponding registration RMS error of each station cloud.

After the iterative global registration, the RMS error of each station is controlled below 5 mm, and the final registration result is shown in Figure 9. In order to check the registration quality, the section data of different positions are cut from the point cloud and checked for stratification. The registration results have not been delaminated, indicating that the registration quality is excellent, which further attests to the accuracy of the method and its practicability of the point cloud registration in the tunnel.

Figure 9.

Iterative global registration results: (a) global point cloud and (b) section of the circular tunnel and rectangular tunnel. The short line in (a) represents the random position of the cut section, and the registration result has high accuracy without any stratification.

Comparison of registration results

We use different commercial software to process the tunnel point cloud data and show the results in Figure 10. Where (a) is the registration result of the feature point registration method in Cloud Compare, obviously, the RMS of the two stations data registration by this method is not high, only 2.8 cm. Because some points are blocked by intersecting tunnels, it is difficult to accurately select feature points, which is not suitable for large-scale tunnel point cloud registration. Among them, (b) is the registration result of 50 stations data by feature registration method and ICP algorithm in the SCENE, and the RMS is 1.3 cm. This method can automatically find out the target balls in the tunnel for coarse registration, and then use ICP for fine registration, which takes about 30 min. However, the automatic registration sometimes will be mismatch into multiple segments, which requires manual operation and will take a longer operation time. In addition, with the increase of data volume, the registration error is difficult to control, even cannot complete the registration of all point cloud data. Among them (c) is the result of cyclone registration using feature registration method and ICP algorithm, the RMS is 1.2 cm. This method has good accuracy but complex operation. It needs manual fitting of the target ball for registration, which takes a long time, about 60 min. Moreover, when it is used to process massive tunnel data, the accumulated error is difficult to eliminate. In contrast, the proposed algorithm eliminates the accumulated error through global adjustment and iterative weighting calculation to further improving the registration accuracy for massive tunnel point cloud data within 30 min. In conclusion, the method proposed in this paper can effectively make up for the shortcomings of commercial software in tunnel point cloud registration and has advantages in both efficiency and accuracy.

Figure 10.

The registration results of different commercial software: (a) is cloud compare, (b) is scene, and (c) is cyclone. Among them, (a) shown the result of two stations, the registration precision is 2.8 cm, which is not very ideal, (b) is the 50 stations registration result, after automatic registration is divided into five segments, the final registration accuracy within 1.3 cm, and (c) show the registration with 85 stations, and section at a certain location, the registration result is about 1.2 cm.

Conclusion

Aiming at the problem of massive point cloud registration in underground tunnel space, this paper proposes an iterative global registration method, which uses the target ball as the feature points for registration. First, the approximate values of each constraint in the base station coordinate system are solved by local Rodrigues coarse registration, and then the iterative global adjustment algorithm is used to adjust all constraint values to reduce error accumulation. For values with large errors, the weights are re-weighted by the posterior variance selection weight iteration method and perform multiple iteration adjustments until the parameter corrections value meets the threshold requirements. With the tunnel point cloud of 1300 m as the research object and the registration, results are compared with other methods, the experiment demonstrates the feasibility, efficiency, and accuracy of the algorithm that the iterative global registration algorithm proposed in this paper and provides a new idea for the registration of multi-view point cloud registration in underground space, which is of tremendous value in practice.

However, due to the spatial limitations of the tunnel, where a large number of control points are distributed in the narrow and long space, and their spatial distribution is lacking in variety, it is impossible to guarantee there are collinear or coplanar multiple control points. In order to prevent the divergence of matrix calculation during the global adjustment caused by collinear or coplanar conditions, the spatial distribution of control points should be detected in future research, and the collinear or coplanar control points should be screened out to further optimize the global registration model.