Research on Point Cloud Registering Method of Tunneling Roadway Based on 3D NDT-ICP Algorithm

In order to meet the needs of intelligent perception of the driving environment, a point cloud registering method based on 3D NDT-ICP algorithm is proposed to improve the modeling accuracy of tunneling roadway environments. Firstly, Voxel Grid filtering method is used to preprocess the point cloud of tunneling roadways to maintain the overall structure of the point cloud and reduce the number of point clouds. After that, the 3D NDT algorithm is used to solve the coordinate transformation of the point cloud in the tunneling roadway and the cell resolution of the algorithm is optimized according to the environmental features of the tunneling roadway. Finally, a kd-tree is introduced into the ICP algorithm for point pair search, and the Gauss–Newton method is used to optimize the solution of nonlinear objective function of the algorithm to complete accurate registering of tunneling roadway point clouds. The experimental results show that the 3D NDT algorithm can meet the resolution requirement when the cell resolution is set to 0.5 m under the condition of processing the point cloud with the environmental features of tunneling roadways. At this time, the registering time is the shortest. Compared with the NDT algorithm, ICP algorithm and traditional 3D NDT-ICP algorithm, the registering speed of the 3D NDT-ICP algorithm proposed in this paper is obviously improved and the registering error is smaller.


Introduction
The environment of the tunneling face in coal mines is complex and the intelligence degree is low, making it difficult for tunneling equipment to perceive the work environment intelligently [1].Environmental modeling for underground roadways is an effective way to solve the above problems.Sensors can be used to acquire the environmental point cloud data of tunneling roadways and establish the environmental point cloud map, which is convenient for guiding and carrying out the excavation and mining work [2,3].However, due to the limitations of the special underground environment, the point cloud map generated directly has low accuracy and large errors in general.At present, improving the accuracy of the underground environment map of coal mines is the primary difficulty in the intelligent construction of the tunneling face [4].
Technologies such as 3D scanning and DT (Digital Twins) have received extensive attention in the industrial field.In terms of building inspection, the German Spacetec company developed the TS3 vehicle-mounted tunnel scanning system.The system adopts the principle of three-dimensional laser scanning measurement, which can perform crack detection, tunnel boundary detection and temperature collection [5].The British Balfour Beatty Raily company (London, UK), developed the LaserFleX TM system, which accurately detects the boundary of the tunnel based on the laser three-dimensional scanning method [6].Columbia University in the United States used a three-dimensional laser scanner to collect data on the collapsed area and calculate changes in the amount of earthwork and changes in contour boundaries [7].Shi Xinxiao et al. verified the feasibility of 3D laser technology in roadway modeling, which provides a more accurate digital model for digital mine construction [8].Zhu Haibin et al. used 3D laser scanning technology to obtain 3D information of the roadway and establish a 3D digital model of the roadway, which improved the efficiency of surveying and mapping work in the roadway engineering [9].In the protection of cultural relics, W. Neubauer Et al. used the Rigel 3D laser scanner to collect data on the Sphinx of Giza, provide data support for the protection and maintenance of the Sphinx and store a complete 3D digital model of the Sphinx [10,11].In terms of digital cities, Vosselman et al. combined 3D laser scanning technology with 2D cadastral vector data, and extracted the height of houses through airborne laser scanning to build a 3D city model [12].The above research methods all provide solutions for environmental modeling of tunneling roadways.However, in the process of generating the point cloud map of tunneling roadways, such solutions would be affected by the internal constraints of the roadway and cause large errors.Therefore, in order to further improve the accuracy of environmental modeling of tunneling roadways, it is also necessary to register the initial point cloud of tunneling roadways.
At present, the commonly used 3D point cloud registering algorithm is the Iterative Closest Points (ICP) [13].The ICP algorithm basically meets the registering requirements of most 3D point clouds, but its running speed and accuracy mainly depend on the given initial transformation estimation, the size of the point cloud and the initial position accuracy [14].Chen H et al. proposed a Two Step ICP (TICP) algorithm [15].The method firstly uses the ICP algorithm for coarse registering on initial point cloud data, and then applies the obtained transformation matrix as initial transformation to the ICP algorithm again, which improves the registering efficiency.SI Choi et al. accelerated the ICP algorithm based on CUDA [16], but it is still difficult to meet the requirements of fast registering of actual point clouds.In order to overcome the above ICP algorithm problems, He et al. proposed a principal component analysis method for point cloud registration.This method mainly provides a good initial value for the ICP calculation by transforming the dimensions of the data.However, when there are many noise points in the point cloud, it is difficult for this algorithm to obtain a good initial value [17].Yang et al. proposed a Scale-ICP algorithm based on seven-dimensional space iteration.Compared with the traditional ICP algorithm, this method has a certain degree of improvement in iterative speed and registration accuracy, but it still relies on the iterative process, which leads to the problem of slow convergence of the algorithm [18].Sharp et al. proposed the ICPIF (Iterative Closest Points using Invariant Features) algorithm to improve the accuracy of finding point pairs [19].In recent years, another registering algorithm proposed by Biber P and Strasser W, Normal Distribution Transform (NDT), has gradually attracted the attention of scholars [20].The NDT algorithm uses standard optimization technology to determine the best registering of two point clouds [21].In the registering process, the features of corresponding points are not used to calculate the registering of point clouds.For a large number of point cloud data, the registering speed is faster than ICP, but its registering accuracy is not as high as ICP.After that, relevant scholars started the fusion research of NDT and ICP algorithms, aiming at realizing fast and accurate registering of point clouds.Sobreira H compared the NDT and ICP algorithms, and added a new algorithm to PCL, which can improve the matching efficiency of the ICP algorithm [22].Pang S fully tested the NDT and ICP algorithms through experiments, and the experimental results show that the NDT algorithm has better ability to deal with real adversity and higher computational efficiency [23].Attia M uses a combination of NDT and ICP algorithms.The experimental results show that using the ICP algorithm to coarsely register the point cloud and then using the NDT algorithm for refinement is better than directly using the ICP algorithm for refinement [24].Wang Qingshan et al. verified the reliability and accuracy of the combination of NDT and ICP algorithms in the campus environment [25].Zhang Guiyang et al. used NDT and ICP algorithms [26], and their experiments verified that an NDT-ICP algorithm can improve the speed and accuracy of point cloud registering.
Wang Yanming et al. applied an NDT-ICP point cloud registering algorithm on the basis of traditional robot 3D topography flexibility measurement data [27], which improves the speed and accuracy of robot topography flexibility measurement.Based on the above research, this paper optimizes and improves 3D NDT and ICP algorithms, and proposes a 3D NDT-ICP method suitable for rapid registering of 3D point clouds in tunneling roadways, which provides a technical scheme for environmental modeling of tunneling roadways.

Mathematical Analysis of Tunneling Roadway Environment
A tunneling roadway is an independent confined space under the coal mine.Environmental modeling and twin model construction are convenient to realize environmentally intelligent perceptions of a tunneling face and create a safe "man-machine-environment" interaction foundation, as shown in Figure 1.

Selection of Point Cloud Acquisition Method for Tunneling Roadways
At present, the research on twin models of tunneling roadway environments is continuously developing [28].According to sensor selection, the commonly used methods for obtaining point clouds of tunneling roadways are mainly divided based on vision camera, based on millimeter wave radar, based on laser radar or laser scanner, etc. Figure 2 shows the use features of the existing point cloud acquisition methods in the tunneling environment.
As can be seen from Figure 2, Lidar has shown unique advantages in tunneling: high data accuracy and strong reliability.At the same time, Lidar obtains less point cloud data in an open and open environment, and obtains more point clouds in a confined space such as roadways, which provides good data support for subsequent tunneling roadway environment modeling.

Feature Analysis of Several Point Cloud Data in Tunneling Roadways
In tunneling roadways, most of the point cloud information can be scanned directly by Lidar, but some special point cloud information belongs to "discrete point cloud", which is easily confused with wrong point cloud, such as "anchor protection point cloud" (point cloud composed of anchor end points), as shown in the red circle in Figure 3.There are differences between the distribution of this kind of point cloud and the overall point cloud in the roadway, which is easy to be regarded as wrong point cloud and screened out by filtering algorithm in the preprocessing process.Firstly, the point cloud of a tunneling roadway is defined as R (Roadway Point Cloud).As the tunneling roadway needs to be supported by an anchor rod, the end points of the anchor rod form a group of relatively discrete point clouds R ar (i.e., anchor protection point clouds) inside the roadway.In the point cloud R, a neighborhood U is selected for any point in the anchor protection point cloud R ar .In the neighborhood U, the mean and variance of the mean distance between the point and all point clouds are calculated.
The points in the anchor protection point cloud R ar mainly have the following mathematical features: In Formula (1), F depends on the number of points in the neighborhood.

Environmental Constraints for Accurate Modeling of Tunneling Roadways
The non-contact measurement accuracy of Lidar in a confined space is mainly affected by the size of the light spot.In the process of scanning and measuring tunneling roadways, the laser beam irradiates the surface of the measuring object to form a light spot, and the radar laser point is theoretically the central position of the light spot.However, the laser point may in fact be located at any position of the spot, so there is an influence of the spot size on the measurement accuracy of the laser point.There are many factors that affect the spot size.Among the many factors, scanning features (incident angle and distance) have the greatest effect [29].When the incident angle and measurement distance change, the spot area would be affected, which affects the laser point measurement accuracy of Lidar, as shown in Figure 4. M. Bitenc showed Equation ( 2) for the change in the long axis X of the elliptical spot [30].In Equation ( 2), S is the measurement distance, H is the laser beam linewidth, γ is the incident angle and R is the spot diameter.
As shown in Figure 5, combined with the tunneling roadway environment, it can be seen that Lidar has good applicability underground, but it is affected by the inevitable constraints of the roadway environment.For example, 1 meteorological conditions such as temperature, air pressure and humidity inside the roadway would affect the atmospheric refractive index, and thus affect the distance measurement accuracy of Lidar. 2 The temperature gradient and atmospheric vibration in the direction of the light path in the roadway would affect the direction of the light, and then increase the angle measurement error of Lidar. 3 The roadway equipment body shields part of the laser radar measurement field of vision.As a result, the tunnel point cloud information generated by a single scanning of Lidar is not perfect, and random errors and some outliers are easily superimposed-i.e., the wrong point cloud R e is generated.
Therefore, it is necessary to use a point cloud registering method to splice laser point clouds from different viewing angles to improve the modeling results of roadway environment.However, the wrong point cloud R e generated in the generation process of laser point cloud would have certain influence on the registering accuracy and efficiency of tunneling roadway point clouds.When the number of points in the wrong point cloud R e is large, it would interfere with the search of corresponding points of point clouds and generate wrong registering point pairs, which can easily cause large deviation in coarse registering of point clouds and local convergence in fine registering of point clouds.Therefore, it is necessary to preprocess the initial point cloud of a tunneling roadway and design an appropriate registering method.

Pretreatment of Point Cloud Data in Tunneling Roadways
At present, the filtering methods for laser scattered point clouds can be generally divided into two types.The first type is to convert the point cloud model into a grid model for filtering processing; the second type is to filter directly by a filtering method.Considering that anchor protection point cloud R ar is easily confused with wrong point cloud R e in tunneling roadways, it would screen out some point clouds of a tunneling roadway and destroy the overall structure of the point cloud of the tunneling roadway by preprocessing the initial point cloud of the tunneling roadway by the second type of direct filtering method, so the first type of filtering method is selected to preprocess the initial point cloud of the tunneling roadway.After the point cloud model of the tunneling roadway is converted into a grid model, the point cloud data in each grid is filtered and screened to ensure the integrity of the point cloud structure and remove a large number of unnecessary point clouds.In this paper, a Voxel Grid filtering method is used to preprocess point cloud data.The specific mathematical derivation is as follows: (1) According to the coordinate set of point cloud data, the maximum values x max , y max and z max the minimum values x min , y min and z min on the three coordinate axes of X, Y and Z are obtained.
(2) According to the maximum and minimum values on the three coordinate axes of X, Y and Z, the side length l x , l y and l z of the minimum bounding box of the point cloud are obtained.
(3) Set the side length of the voxel small grid cell and divide the X, Y and Z coordinate axes into M, N and L parts equally, then the minimum bounding box is divided into In Equation (4), • means rounding down and sum is the total number of voxel small grids.
(4) Number the small grid of each voxel-the numbers are (i, j, k)-and determine the voxel cell to which each data point belongs.
(5) Carry out point cloud reduction filtering.Calculate the center of gravity of each voxel small grid, and replace all points in the voxel small grid with the center of gravity.
If the center of gravity does not exist, all points in the small grid are replaced by the data point closest to the center of gravity.
where c ijk , p i and k are the center of gravity, data points and points of voxel small grids, respectively.

Accurate Registering Method of Point Cloud in Tunneling Roadways
After preprocessing the point cloud data, four groups of experimental point clouds are obtained, one of which is recorded as the reference point cloud M(x, y, z) and the other as the point cloud to be registered N(x, y, z).The two have overlapping areas, O(x, y, z) and O = M ∩ N. In order to reduce the influence of roadway internal environment constraints, reduce the registering error of tunneling roadway point clouds and improve the registering speed, this paper integrates the advantages of NDT and ICP algorithms, and proposes a registering method of tunneling roadway point clouds based on a 3D NDT-ICP algorithm.
Most of the traditional 3D NDT-ICP algorithms use two-step combination registering.First, the 3D NDT algorithm is used to correct the initial input point cloud pose, and then the ICP algorithm is used to register the point cloud after pose correction.The specific process is shown in Figure 6.In the traditional 3D NDT-ICP algorithm, the 3D NDT algorithm first performs coarse registration on the point cloud that needs to be registered, corrects the initial position of the point cloud to be registered and makes it approach the reference point cloud.Then, the ICP algorithm is used for fine registration, and the point cloud to be registered after the 3D NDT algorithm is registered.In the traditional two-step 3D NDT-ICP algorithm, the project files need to be compiled separately, a total of two times, the efficiency is low and the registration process is not coherent.
As shown in Figure 7, the flow of the 3D NDT-ICP algorithm proposed in this paper is firstly optimized according to the environmental features of tunneling roadways, and then the point cloud coordinate transformation matrix solved by the 3D NDT algorithm is taken as the initial matrix of the ICP algorithm to ensure the consistency of registering.Then, a kd-tree is used for point pair search in the registering process of ICP algorithm to improve the registering speed of point cloud.Finally, the Gauss-Newton method is used to iteratively optimize the objective function of ICP algorithm to solve the optimal coordinate transformation parameters between point clouds and complete the accurate registering of point clouds.

Solution of Point Cloud Coordinate Transformation Parameters in Tunneling Roadways
Firstly, the reference point cloud M(x, y, z) data sample is evenly divided into several regular 3D voxel units with the same size by a 3D NDT algorithm, and then the probability distribution of each 3D point cloud position in the 3D voxel unit is expressed by normal distribution; the expression is shown in Equation (7).The points of the corresponding point pair in the reference point cloud M are denoted as M i , M i = (x M , y M , z M ), M i ∈ M, i = 1, 2, 3, . .., and the points in the point cloud N to be registered are denoted as The parameter C is the covariance matrix of the 3D point clouds in each voxel unit and q is the mean value of the 3D point clouds in each voxel cell.The parameter c is a constant.Specific definitions of q and C are shown in Equations ( 8) and (9).
Equation ( 11) initializes the transformation parameters for the point cloud M and N: According to coordinate transformation parameters, each point cloud sample of the point cloud to be registered N(x, y, z) is mapped into a reference point cloud M(x, y, z) data sample coordinate system.The mapped point cloud is recorded as N , the probability distribution of each 3D point mapping in N is summed and the coordinate transformation parameters are evaluated: Through a Hessian matrix method, s(p) is optimized to maximize the value of s(p).The problem of solving the optimal transformation of the matrix is regarded as the process of minimizing s(p).The realization of this process is solved by the Newton method and the Hessian matrix.Let f = s(p); in order to minimize the function f , the following equations must be processed for each operation: In Equation ( 13), g is the transposed gradient of f , and the specific elements can be expressed as: H is the Hessian matrix of f , and its elements can be expressed as: After obtaining ∆p, we can update ∆p with the following formula:

Transfer of Point Cloud Coordinate Transformation Parameters in Tunneling Roadways
Any point P(x p , y p , z p ) in the reference point cloud is M, and any point Q(x q , y q , z q ) in the point cloud is to be registered N.Meanwhile, P and Q satisfy P ∈ (M ∩ O) and Q ∈ (N ∩ O), and the external conversion relationship of the two independent 3D coordinate systems is further described by using the seven-parameter spatial similarity transformation model through the ICP algorithm [31], as shown in Formula (17).
where t x , t y and t z are the three components along the coordinate axis direction, α, β and γ, and are the three angle parameters rotating around the coordinate axis; w is the scale transformation factor between coordinate systems, which generally defaults to 1.Then, the final result of Equation ( 11) is passed to the ICP algorithm to initialize the rotation matrix and translation matrix in Equation (17).
In order to speed up ICP registering, this paper uses a kd-tree to retrieve roadway point cloud in the ICP algorithm part.This data search method is widely used in the field of point cloud research, so it is not listed as the focus of introduction.Assuming that the two groups of point clouds in the tunneling roadway jointly generate m groups of corresponding point pairs, when Equation ( 19) obtains the optimal solution, the two groups of point clouds complete the solution of the spatial position conversion relationship.
In order to facilitate the derivation of the subsequent iterative optimization equation, in Equation ( 19) M i − (RN i + T) 2 = d 2 = D, d represents the Euclidean distance between M i and N i , as shown in Equation (20):

Iterative Optimization of Point Cloud Coordinate Transformation Parameters in Tunneling Roadways
The objective function E(R, T) in Equation ( 19) needs to be optimized.That is, Equation ( 23) is minimized.Due to the nonlinearity of the objective function, the estimated values of parameters cannot be obtained by solving the extreme values of linear least squares and other multivariable functions, so complex optimization algorithms are needed to solve this problem.In this paper, the Gauss-Newton method is used to optimize the objective function in ICP algorithm, so as to avoid local convergence of tunneling roadway point clouds in the process of precise registering.
The general form of the nonlinear system model is: Among them, The sum of squares of Euclidean distances at a certain point between two groups of point clouds in the roadway before and after registering is: We can substitute Equation ( 21) into Equation ( 22) for iterative operation processing: Among them, f (x) = ( f 1 (x), f 2 (x), . . ., f m (x)) T , x = (t x , t y , t z , w, α, β, γ) T .If the point of the k iteration is x k , it can be seen from the Taylor expansion equation that: A(x k ) is a multidimensional first-order partial derivative matrix, namely a Jacobian matrix: In the above equation, x j = (t x , t y , t z , w, α, β, γ) T .With A k = A(x k ), the Equation ( 25) is combined and brought into the Equation ( 23), the following results can be obtained: Among them, d k = x − x k .For the minimum problem of minS(x), the least square method can obtain: when A k T A k is reversible, there are: Because d k = x − x k , when x k+1 = x = x k + d k , in combination with Equation ( 29), there is an iterative equation: the Hessian matrix of minS(x) at point x k .
In sum, it can be concluded that ) is the Gauss-Newton equation and ) can be called the Gauss-Newton direction.

Experimental Verification
In the experiment, the operating system of the upper computer is Windows 10, the running memory is 8 G, and the software is Visual Studio 2017.

Point Cloud Data Preprocessing Experiment
First, make four sets of PCD files that simulate the point cloud of mine tunnels: 1 (test1-a.pcd,test1-b.pcd), 2 (test2-a.pcd,test2-b.pcd), 3 (test3-a.pcd,test3-b.pcd)and 4 (test4-a.pcd,test4-b.pcd). 1 test1-a.pcd is the reference point cloud, test1-b.pcd is the point cloud to be registered; 2 test2-a.pcd is the reference point cloud, and test2-b.pcd is the point cloud to be registered; 3 test3-a.pcdis the reference point cloud, test3-b.pcd is the point cloud to be registered; 4 test4-a.pcd is the reference point cloud, and test4-b.pcd is the point cloud to be registered.Then, the Voxel Grid filtering algorithm and passthrough filtering algorithm are used to preprocess the four groups of laser point clouds 1 , 2 , 3 and 4 respectively.Since the reference point cloud and the point cloud to be registered are essentially the same, this paper preprocesses the reference point cloud in the four groups of point clouds.The experimental results are shown in Figure 8, and the experimental data are shown in Table 1.
As can be seen from Table 1, the pass-through filtering reduces the number of laser point clouds more than the Voxel Grid filtering method.However, combined with Figure 8a,b, it can be seen that the pass-through filtering can easily destroy the point cloud structure of tunneling roadways, which results in a loss of integrity, while the Voxel Grid filtering method reduces the number of point clouds while ensuring the structural integrity of point clouds.It is more suitable for preprocessing point cloud data of tunneling roadways.The four groups of laser point clouds filtered by the Voxel Grid algorithm are recorded as experimental point clouds for subsequent experimental research.As can be seen from Table 1, the pass-through filtering reduces the number of laser point clouds more than the Voxel Grid filtering method.However, combined with Figure

Parameter Optimization Experiment of NDT Algorithm
Influenced by Equation ( 7), the size of cells in 3D NDT algorithm affects the registering accuracy.Previously, the scholar Magnusson made a systematic comment on the cell size setting in 3D NDT method.When the cell setting is too large, it cannot well represent the features of point cloud.However, when the cell setting is too small, it is vulnerable to the noise of radar scanning equipment and may not be able to calculate the Gaussian distribution due to insufficient data in the cell.At the same time, Magnusson concluded through system experiments that the cell resolution of a 3D NDT algorithm is usually between 0.5 m and 2 m [32], which is ideal for laser scanning equipment.Combined with the above conclusions, in order to ensure that the cell resolution interval is sufficient, this paper chooses to set the cell resolution of the 3D NDT algorithm to 0.25 m, 0.5 m, 1 m, 1.5 m and 2 m.Based on four groups of experimental point clouds, experiments are carried out to provide the basis for selecting algorithm parameters for tunneling roadway environment modeling.The experimental results are shown in Figure 9, and the algorithm registration time is shown in Figure 10.In Figure 9, the green display part is the reference point cloud, and the red display part is the point cloud to be registered.When designing the 3D NDT-ICP algorithm, this article aims to quickly obtain high-precision registration point clouds.Since the accuracy of the NDT algorithm has little effect on the overall 3D NDT-ICP algorithm in this paper, this paper focuses on the time spent by the NDT algorithm.It can be seen from Figure 10 that when using the 3D NDT algorithm to register the point cloud of the tunnel and roadway, setting the pixel resolution to 0.5 m can reduce the time spent by the algorithm to the shortest.and Xi is the true value of the Euclidean distance between the corresponding points.Under absolutely ideal conditions, the corresponding points are fully registered and the distance is zero.Therefore, the value of Xi here is 0. Figure 13 shows X, Y and Z three-axis root mean square error after the four sets of experimental point clouds are registered.As can be seen from Figure 13, when registering the experimental point cloud, the 3D NDT-ICP algorithm proposed in this paper has the smallest registering error compared with the NDT algorithm, ICP algorithm and traditional 3D NDT-ICP algorithm.
Combined with Figures 12 and 13, it can be seen that the 3D NDT-ICP algorithm proposed in this paper reduces the registering error of tunneling roadway point clouds and saves the registering time of point clouds when the experimental point clouds are matched.
As shown in Figure 14, the requirements for point cloud registering accuracy in different stages of the tunneling face are shown.Combined with Figure 13, it can be seen that the algorithm proposed in this paper meets the requirements of the quality control of tunneling roadways forming in the first stage and the constraint accuracy of roadheader motion space in the second stage.It is beneficial for the high-efficiency registration of largescale point cloud data of tunneling roadways, and enhances the real-time and accuracy of the construction of twin models of tunneling face environment.Subsequent research is still needed to further improve the accuracy of point cloud models of tunneling roadways.

Figure 1 .
Figure 1.Digital twin conceptual model of tunneling face.

Figure 2 .
Figure 2. Comparison of point cloud acquisition methods.

Figure 4 .
Figure 4. Influence of incident angle on long axis of light spot.

Figure 5 .
Figure 5. External features of a tunneling roadway.
(a) The first and second groups of point cloud pretreatment (b) The third and fourth groups of point cloud pretreatment

Figure 14 .
Figure 14.Accuracy requirements of an underground tunneling face.

Table 1 .
Experimental data of laser point cloud pretreatment.

Table 1 .
Experimental data of laser point cloud pretreatment.