Robust and automatic modeling of tunnel structures based on terrestrial laser scanning measurement

Abstract

The terrestrial laser scanning technology is increasingly applied in the deformation monitoring of tunnel structures. However, outliers and data gaps in the terrestrial laser scanning point cloud data have a deteriorating effect on the model reconstruction. A traditional remedy is to delete the outliers in advance of the approximation, which could be time- and labor-consuming for large-scale structures. This research focuses on an outlier-resistant and intelligent method for B-spline approximation with a rank (R)-based estimator, and applies to tunnel measurements. The control points of the B-spline model are estimated specifically by means of the R-estimator based on Wilcoxon scores. A comparative study is carried out on rank-based and ordinary least squares methods, where the Hausdorff distance is adopted to analyze quantitatively for the different settings of control point number of B-spline approximation. It is concluded that the proposed method for tunnel profile modeling is robust against outliers and data gaps, computationally convenient, and it does not need to determine extra tuning constants.

Introduction

A variety of research has been carried out in the engineering field to recognize and reconstruct three-dimensional (3D) objects. ¹ Terrestrial laser scanning (TLS) has been increasingly applied in the deformation monitoring of tunnel structures, since it can acquire full-field point cloud quickly instead of sparing single points, which is a typical traditional measurement manner. ² Deformation analysis usually requires the estimation of an analytical model, which could be, for example, a mesh model, cylinder or free-form curve and surface in the case of tunnels. ³ It is sometimes difficult to avoid the existence of outliers and data gaps in the point cloud data obtained by TLS measurement. ⁴ When gaps or blank exists in the measurement data, the filtering process will not fill the gap area since there is no data obtained with the scanning measurement. ⁵ The decreased resolution and density can result in a higher uncertainty of deformation analysis. ⁶ Therefore, a precise approximation is essential for the reliable deformation analysis in the field of engineering applications.

Surface approximation has been employed to establish a surface model and to extract deformation information for the analysis of structures. ⁷ The benefits of surface approximation lie mainly in the global parametrization and avoidance of noisy and scattered data points. ⁸ Since free-form B-spline is an indispensable method for global parametrization, the influence of outliers and data gaps on B-splines should be a concern. Outliers can deteriorate the result of B-spline approximation, with the consequence that the result of deformation analysis is inaccurate and unreliable. In adition, data gaps in a point cloud are harmful in B-spline fitting, because the missing data brings about oscillation near the gap area. ⁹

A traditional way to assure the precision is to delete outliers in advance of the approximation. ⁶ Some methods used frequently to eliminate point cloud outliers for tunnel deformation monitoring are listed here. Obvious outliers could be removed manually in certain software, and others could be reduced according to the distribution analysis of, for example, longitudinal, surrounding, and random errors of the measurement data. ¹⁰ Smoothing techniques were also used to filter the outliers, where segmentations were smoothed by averaging adjacent points. ⁵ Panholzer and Prokop ¹¹ filtered disturbing objects in the point cloud data using the wedge absolute method, where adjacent points are connected to a laser source and the angle formed was checked with a previously defined angular threshold value to determine disturbing objects. Moreover, a stability test was employed to distinguish outliers and the real deformation, which was performed on a local approximation of the overall analytical model. ¹¹ Segmented planes were computed through principle component analysis and Hough transformation to preserve the main shape of some rocky surfaces. ¹² However, these procedures were mostly local, time- and labor-consuming, and relied on the expertise of operators, which could be a problem when processing data of a large-scale, kilometers-long structure. ¹³

The performance of B-spline approximation with contaminated data was studied in the literature. Predefined parameters of B-splines, such as the number of control points, can affect the performance of B-splines against noise,^14–22 and the optimal parameters can be obtained as the values associated with the extremum of the likelihood ratio. ²³ Decreasing the model complexity using fewer parameters can filter some noise, but this sacrificed the accuracy of the description of the local shape. A suitable selection of B-spline parameters balances noise resistance and local shape information preservation better compared with polynomials. ²⁴ Combining data from other sensors helped to validate the precision of B-splines constructed from complex point cloud data. ²⁵

Control points determination of B-spline approximation could be counted as one of the most vital steps for point cloud processing. Although plenty of methods of fitting B-splines to measurement data have been introduced,^26–28 only a few dealt with a robust method of control point placement which is immune to outliers and gaps in the measurement data, and which, moreover, satisfies a certain tolerance regarding the real shape. To date, the control point positions are obtained mostly through least squares (LS) estimation. ²⁹ Some researchers determined the optimal number of control points, which was usually a model selection problem, by means of the Akaike Information Criterion, the Bayesian Information Criterion, or statistical learning theory. ³⁰ The quality of the control point positions estimated was evaluated by the average estimated standard deviation, which focused on the quality of estimation procedure regarding datasets instead of the real shape of the object. ³¹ Minimizing energy function has recently been adopted to adjust control points which could be applied in smoothing and gap filling. ³² Xu et al. ³³ employed an adaptive robust estimator based on Student’s t-distribution to fit a B-spline curve to measurement data. This estimator involves the degree of freedom of the t-distribution as a parameter or tuning constant, which is adapted to the data given according to the outlier characteristics.

Statistical method based on ranking was first proposed for simple location problems by Wilcoxon ³⁴ and Hettmansperger and McKean. ³⁵ R-estimators have gained increasing attention in geodetic science.^36,37 Kargoll ³⁸ described an algorithm for computing R-estimates of the parameters of a general linear model with satellite gravity gradiometer data. Duchnowski ³⁹ adapted the theory of R-estimation to the adjustment of leveling networks with homoscedastic observations in the context of deformation analysis. To the best of our knowledge, the rank-based method has not been studied in the underground space reconstruction and TLS data processing.

Actually, there are not plenty of research works on robust 3D model reconstruction for tunnels. The geometries of tunnels are macroscopic range, from several meters to kilometers. The resolution of tunnel measurement by TLS is ca. several millimeters. This requires a specific B-spline model which could adapt the characteristics of such a structure scale and measurement resolution. Furthermore, measurement data of such structures contain occlusions and blanks, disturbing objects, and coexisting sparse and dense points. This may result in the instability of B-spline approximation, which is unacceptable for the automation process. In this article, we propose a robust solution of such problems based on the rank methods.

The contents of the article are arranged as follows: section “Introduction” introduces the background and the state of the art. Section “Method” is a detailed description of the proposed method, which includes B-spline algorithm, rank-based Wilcoxon method, and the Hausdorff distance method as a metric of the B-spline modeling. In section “Data analysis,” we analyze the tunnel profile based on the proposed method, and compare with LS B-spline approximation and the manually denoised data. The comparative results are presented in section “Results.” Finally, we draw the conclusion of the article in the final section.

Method

B-spline algorithm

The B-spline curve is adopted to reconstruct a 3D model describing the real geometry of the profile. A B-spline curve is defined after equation (1) ⁴⁰

$$C\left(u\right)=\sum _{i=0}^n{N}_{i,p}\left(u\right){\mathrm{P}}_i,0\le u\le 1$$

(1)

where C(u) is a B-spline curve, $\{{\mathrm{P}}_i\}$ are the n + 1 control points and $\{N_{i,p}(u)\}$ are the pth-degree B-spline basis functions ⁴⁰ defined on knot vector U

$$\mathbf{U}=\{0,\dots ,0,u_{p+j},\dots ,u_{m-p-1},1,\dots ,1\}$$

(2)

The number of knots is m + 1 and related to the number of control points n + 1 and the degree p by

$$m=n+p+1$$

(3)

The control points $\{{\mathrm{P}}_i\}$ in the B-spline curve approximation problem are unknown. Since the B-spline curve is a linear function, the control points can be obtained in a LS sense. Three main steps are required to compute the control points when the degree p and the control point number n + 1 are fixed. The first step is parametrization, which precomputes parameter value $\overline{u}$ for the N data points given $\{Q_k\},k=1,\dots ,N$ . The parametrization method adopted in this article is based on the chord length, ⁴¹ where the length of each chord is defined as $|Q_k-Q_{k-1}|,k〉1$ . The parameter value $\overline{u}$ is determined with the chord length after

$${\overline{u}}_k=\{\begin{array}{c}0,\mathrm{if}k=1\hfill \\ \multicolumn{1}{c}{{\overline{u}}_{k-1}+\frac{|Q_k-Q_{k-1}|}{\sum _{k=1}^N|Q_k-Q_{k-1}|}}\\ \multicolumn{1}{c}{1,\mathrm{if}k=N}\end{array},\mathrm{if}k=2,\dots ,N-1$$

(4)

The second step for B-spline curve approximation is the knot vector determination. We employ a knot placement technique ⁴¹ in this article because it can reflect the distribution of $\{{\overline{u}}_k\}$ and result in a positive, definite, and well-conditioned matrix for the control point solution. The knot vector $\mathbf{U}=\{u_0,\dots ,u_m\}$ is computed after equation (5), ⁴¹ where $t=\mathrm{int}\left(d\right)$ means the largest integer satisfying $t\le d$

$$u_{p+j}=(1-\alpha ){\overline{u}}_{t-1}+\alpha {\overline{u}}_t(j=1,\dots ,n-p)$$

(5)

With $d=(m+1)/(n-p+1),t=\mathrm{int}\left(\mathrm{jd}\right),\alpha =\mathrm{jd}-t$ .

Based on the parameterization and knot vector determination, the basic functions are computed through Cox-de Boor recursion formulas ⁴¹

$$\begin{array}{c}{N}_{i,0}\left(\overline{u}\right)=\{\begin{array}{c}1,\mathrm{if}{u}_i\le \overline{u}\le u_{i+1}\hfill \\ \multicolumn{1}{c}{0,\mathrm{otherwise}}\end{array}\hfill \\ N_{i,p}\left(\overline{u}\right)=\frac{\overline{u}-u_i}{u_{i+p}-u_i}{N}_{i,p-1}\left(\overline{u}\right)+\frac{u_{i+p+1}-\overline{u}}{u_{i+p+1}-u_{i+1}}{N}_{i+1,p-1}\left(\overline{u}\right)\hfill \end{array}$$

(6)

and assembled for design matrix A, which contains $N\times (n+1)$ elements

$${\mathbf{A}}_{N\times (n+1)}=\left[\begin{array}{ccc}{N}_{0,p}\left({\overline{u}}_1\right)\hfill & \dots \hfill & N_{n,p}\left({\overline{u}}_1\right)\hfill \\ \multicolumn{1}{c}{⋮}& \hfill & ⋮\hfill \\ \multicolumn{1}{c}{N_{0,p}\left({\overline{u}}_N\right)}& \dots \hfill & N_{n,p}\left({\overline{u}}_N\right)\hfill \end{array}\right]$$

(7)

The unknown control point matrix is obtained through the Gauss-Markov model

$$\mathbf{X}={\left({\mathbf{A}}^T\mathbf{A}\right)}^{-1}{\mathbf{A}}^T\mathbf{l}$$

(8)

where the observation matrix consists of x, y, and z coordinates of N data points

$${\mathbf{l}}_{N\times 3}=\left[\begin{array}{ccc}{x}_1\hfill & y_1\hfill & z_1\hfill \\ \multicolumn{1}{c}{\dots }& \dots \hfill & \dots \hfill \\ \multicolumn{1}{c}{x_N}& y_N\hfill & z_N\hfill \end{array}\right]$$

(9)

and the unknown control point matrix is written as

$${\mathbf{X}}_{(n+1)\times 3}=\left[\begin{array}{ccc}{P}_{x0}\hfill & P_{y0}\hfill & P_{z0}\hfill \\ \multicolumn{1}{c}{\dots }& \dots \hfill & \dots \hfill \\ \multicolumn{1}{c}{P_{\mathrm{xn}}}& P_{\mathrm{yn}}\hfill & P_{\mathrm{zn}}\hfill \end{array}\right]$$

(10)

The B-spline method is advantageous for deformation data processing since it is capable of modeling the shape of a structure flexibly. Control point determination is an important step in B-spline approximation, which is usually realized by means of the LS method. This article improves the control point adjustment by employing an R-estimator based on Wilcoxon scores in order to obtain robust B-splines. The main computational framework is shown in Figure 1, where iterative steps are carried out to optimize the control point positions. First of all, we compute the parametrization, knot vector, and control points according to standard methods introduced in Piegl and Tiller. ⁴¹ Thereafter, the results of the parametrization and knot vector determination are input into the computing package of the aforementioned R-estimator, which also requires the specification of the design matrix A. Initial values for the control point positions can be determined as the LS solution, which corresponds to the standard B-spline approximation method as shown in the literature, such as in the previous work. ⁴²

Figure 1.

Framework of rank-based B-spline modeling.

Rank-based Wilcoxon method

B-spline models can be expressed in the form of a linear model

$$\mathbf{l}+\mathbf{r}=\mathbf{A}\mathbf{\beta }$$

(11)

where l represents the observation vector (consisting of the measure point coordinates), r the vector of residuals, $\mathbf{A}$ the design matrix, and $\mathbf{\beta }$ the parameter vector (consisting of the unknown control point coordinates). The well-known LS solution

$${\hat{\mathbf{\beta }}}_{\mathrm{ols}}={\left({\mathbf{A}}^T\mathbf{A}\right)}^{-1}{\mathbf{A}}^T\mathbf{l}$$

(12)

is obtained by minimizing the L2 norm

$$Q\left(\mathbf{\beta }\right)=\sum _{i=1}^n{\mathbf{r}}_i^2$$

(13)

of the residuals $r_i={\mathbf{A}}_i\mathbf{\beta }-{\mathbf{l}}_i$ regarding the unknowns $\mathbf{\beta }$ (where ${\mathbf{A}}_i$ denotes the ith row of the design matrix). In this article, we employ the R-estimator which obtain the solution by minimizing the pseudo norm

$$D\left(\mathbf{\beta }\right)=\sum _{i=1}^{n}a\left(R\left({\mathbf{r}}_i\right)\right){\mathbf{r}}_i$$

(14)

The rank function R assigns to each residual $r_i$ its rank or position within the sorted residuals ${\mathbf{r}}_{\left(1\right)}\le \dots \le {\mathbf{r}}_{\left(n\right)}$ as the index occurring in the brackets, so that $R({\mathbf{r}}_1)=1,\dots ,R({\mathbf{r}}_n)=n$ , when multiple residuals take identical values, they also share the same, average rank. We use the Wilcoxon scores generated by

$$a\left(R\left({\mathbf{r}}_i\right)\right)=\sqrt{12}·(\frac{R\left(\mathbf{r}\right)}{n+1}-\frac{1}{2}),R\left({\mathbf{r}}_i\right)=1,\dots ,n$$

(15)

Thus, the ranks become centered and standardized to emulate the behavior of residuals for any given number n of residuals given in a practical situation. The point of minimizing $D(\mathbf{\beta })$ is to replace the quadratic residuals in the L2 norm $Q(\mathbf{\beta })$ by a rather linear form, since the LS estimator tends to level out large residuals (i.e. outliers) with a deleterious effect on the adjusted model function. Instead of squaring outliers, an outlier is multiplied by one of the values $a(i)$ , yielding a less extreme contribution to the target function to be minimized. Consequently, an outlier is not leveled out but more or less ignored in the minimization process, so that it can show up in the estimated residuals

$${\hat{\mathbf{r}}}_{\mathrm{rbm}}=\mathbf{A}{\hat{\mathbf{\beta }}}_{\mathrm{rbm}}-\mathbf{l}$$

(16)

regarding the rank-based method much more clearly than in the (leveled) estimated residuals

$${\hat{\mathbf{r}}}_{\mathrm{ls}}=A{\hat{\mathbf{\beta }}}_{\mathrm{ls}}-\mathbf{l}$$

(17)

regarding the LS estimates. Using initial values ${\mathbf{\beta }}_0$ , the solution ${\mathbf{\beta }}_{\mathrm{rbm}}$ can be computed in each iteration step $s=1,\dots ,.\mathrm{itermax}$ by

$${\hat{\mathbf{\beta }}}_{\mathrm{rbm}}^{\left(s\right)}={\mathbf{\beta }}_0+\mathrm{\Delta }{\mathbf{\beta }}_{\mathrm{rbm}}^{\left(s\right)}$$

(18)

with

$${\widehat{\mathrm{\Delta }\mathbf{\beta }}}_{\mathrm{rbm}}^{\left(s\right)}={\left({\mathbf{A}}^T{\mathbf{P}}^{\left(s\right)}\mathbf{A}\right)}^{-1}{\mathbf{A}}^{\mathrm{T}}{\mathbf{P}}^{\left(s\right)}{\hat{\mathbf{r}}}_{\mathrm{rbm}}^{\left(s\right)}$$

(19)

An iteratively reweighted LS algorithm is formed based on the weight matrix

$${\mathbf{P}}^{\left(s\right)}=\left[\begin{array}{ccc}{p}_1^{\left(s\right)}\hfill & \cdots \hfill & 0\hfill \\ \multicolumn{1}{c}{⋮}& \ddots \hfill & ⋮\hfill \\ \multicolumn{1}{c}{0}& \cdots \hfill & p_n^{\left(s\right)}\hfill \end{array}\right]$$

(20)

whose individual weights are determined by

$$p_i^{\left(s\right)}=\frac{\sqrt{12}·(\frac{R({\hat{\mathbf{r}}}_{i,\mathrm{rbm}}^{\left(s\right)}-{\alpha }_0)}{n+1}-\frac{1}{2})}{({\hat{\mathbf{r}}}_{i,\mathrm{rbm}}^{\left(s\right)}-{\alpha }_0)},(i=1,\dots ,n)$$

(21)

where ${\alpha }_0$ is the median of the current residuals ${\hat{r}}_{i,\mathrm{rbm}}^{\left(s\right)}$ . Here, the ranks obtained $R(r_1),\dots ,R(r_n)$ are transformed into order-preserving scores $a(R(r_1)),\dots ,a(R(r_n))$ , so that $a(R(r_1))\le \dots \le a(R(r_n))$ . These scores are defined by means of some score function $\phi$ via the relationship

$$a\left(i\right)=\phi \left(\frac{i}{n+1}\right)$$

(22)

The Wilcoxon scores used in this article are generated by $\phi (u)=\sqrt{12}·(u-1/2)$ . The individual weights of the weight matrix ${\mathbf{P}}^{\left(s\right)}$ are determined by

$$p_i^{\left(s\right)}=\frac{a\left(R\left({\hat{\mathbf{r}}}_{i,\mathrm{rbm}}^{\left(s\right)}\right)\right)}{{\hat{\mathbf{r}}}_{i,\mathrm{rbm}}^{\left(s\right)}-{\alpha }_0},i=1,\dots ,n$$

(23)

where the median ${\alpha }_0$ of the current residuals ${\hat{\mathbf{r}}}_{i,\mathrm{rbm}}^{\left(s\right)}$ is subtracted to ensure positive weights. ⁴³

B-spline approximation algorithm is investigated and developed by adopting the R-estimator. Figure 1 presents the framework of rank-based B-spline modeling which involves methods to determine the weight matrix (f_w) and the weighted LS fit (f_v). ³⁸ An additional weight matrix is estimated based on the initial control points, by executing the weight computing routine f_w. Then, the improvement matrix for control points is obtained through the improvement computing routine f_v. Subsequently, a decision is made based on the criteria values of improvement and iteration times. If the values pass the criteria, optimal control points will be gained, otherwise, the control points will be updated with the aforementioned improvement matrix.

The key step in the framework is the rank-based Wilcoxon method, which includes weight and improvement computation, as well as the decision-making criteria.

Hausdorff distance

The Hausdorff distance is suitable for many areas, especially for the problem of shape comparison and matching. ⁴⁴ Hossain et al. ⁴⁵ investigated a linear time algorithm for computing the Hausdorff distance in the context of content-based image analysis, which is a significant metric in computer vision including image matching, shape retrieval and moving-object detection.

A robust approximation is supposed to be able to eliminate the influence of outliers automatically. Therefore, the performance of robustness is assessed by agreement of the adjusted observations with sets of data which are denoised beforehand (i.e. from which the outliers are removed). The Hausdorff distance is employed to measure the distance between the denoised point cloud and the approximation results. It is defined by equation (24), where a and b are points in two different sets of points A and B, and $\mathrm{a},\mathrm{b}$ means a norm between points a and b, which is Euclidean norm in this article. A generalized form of Hausdorff distances is computed after equation (25)

$$\mathrm{HD}(\mathrm{A},\mathrm{B})=\underset{\mathrm{a}\in \mathrm{A}}{max}\underset{\mathrm{b}\in \mathrm{B}}{min}\Vert \mathrm{a},\mathrm{b}\Vert$$

(24)

$$\mathrm{H}(\mathrm{A},\mathrm{B})=max(\mathrm{HD}(\mathrm{A},\mathrm{B}),\mathrm{HD}(\mathrm{B},\mathrm{A}))$$

(25)

The different sequences of point sets A and B express different physical meanings and result in a different distance value. When assessing the B-spline approximation, we consider the denoised data as a reference and match the approximation to the reference.

Data analysis

The original point cloud data is approximated by using the rank-based and LS method, where a comparison is carried out between the two. The algorithm of the rank-based B-spline modeling is presented in Figure 2, where the input data is the measurement data named D_m and the denoised data D_d . Both sets of data contain X, Y and Z coordinates of profile points. A difference is that the D_d has been processed where the outliers have been deleted, and the number of points is less than that of the D_m . The methods involved are the rank-based method (rbm), the least squares (ls) method, weight matrix f_w and weighted LS fit f_v. The results obtained are estimated control point coordinates.

Figure 2.

Algorithm of the rank-based B-spline modeling.

The improvement of the rank-based method from the point of infrastructure measurement is shown in Figure 3, where the point cloud depicts the tunnel profiles. The dimensions of the tunnel is 2.5 m in height and 5.5 m in width. The tolerance value tol adopted for rank-based B-spline modeling for the tunnel profiles is 10e-6. In order to bring out the difference between the rank-based and the LS method visually, the B-spline approximations involving different numbers of control points (cp) are shown for the same point cloud in Figures 3(a)–(c), where the cp is varied from 7 to 9. In each of these figures, the blue points represent the point cloud data with noise (including outliers) and the red and green lines represent the B-spline fitted by means of the rank-based and the LS method, respectively.

Figure 3.

Comparison between the rank-based and LS method: (a) cp = 9, (b) cp = 8, and (c) cp = 7.

It can be observed in Figure 3(a) that the rank-based method produces a solution that is practically unaffected by the outliers, whereas the LS solution is clearly attracted by the outliers near the bottom of the region indicated by the dashed blue line. As the tunnel structure itself does not sag in reality, the former solution is much more realistic than the latter, especially in the virtual box of Figure 3(a). It is concluded that under the same conditions, the rank-based method is more robust than the LS method when cp = 9.

According to the virtual box of Figure 3(b), the LS method is more sensitive to noise than the rank-based method when cp = 8. However, the rank-based solution for cp = 8 is notably less smooth than the rank-based solution for cp = 9. It is, thus, revealed that the rank-based method is rather insensitive to the outliers but quite sensitive to the choice of the number of control points.

It can be discovered that the rank-based method is still more robust than the LS method when cp = 7, as indicated by the virtual box of Figure 3(c). In comparison to the Figures, the difference between the rank-based and the LS method is diminished. It means that if the number of control points are fewer, the under-fitting condition of the B-spline will affect the optimization effect of the rank-based method. Thus, it is qualitatively confirmed that the rank-based method has superior robustness than the mostly adopted ordinary least square (OLS) in current literature works, and under the same conditions, the optimization performance of the rank-based method is related to the number of control points.

Results

The Hausdorff distance, which describes how far away two subsets of a metric space are from each other,^46,47 is employed to quantitatively analyze the improvement of the rank-based over the LS solution. The numerical values of this distance measure are shown in Table 1 for three different point clouds (“Data1,”“Data2,” and “Data3”).

Table 1.

Comparison of the rank- and LS-based B-spline modeling.

	cp	Hausdorff distance (10e-2)
		rbm	ls	Improvement (%)
Data1	4	3.44	4.79	39.24
	5	3.78	5.79	53.17
	6	2.63	6.56	149.43
	7	3.51	7.91	125.36
	8	2.52	8.44	234.92
	9	3.65	11.52	215.62
	10	3.03	10.88	259.08
Data2	4	10.26	12.48	21.64
	5	7.21	7.62	5.69
	6	5.04	7.92	57.14
	7	4.01	9.6	139.40
	8	4.73	10.65	125.16
	9	4.84	14.68	203.31
	10	3.36	11.37	238.39
Data3	4	15.53	15.57	0.26
	5	16.28	17.28	6.14
	6	10.9	19.4	77.98
	7	17.13	19.96	16.52
	8	16.54	18.8	13.66
	9	13.74	18.4	33.92
	10	8.22	18.34	123.11

cp: control point; rbm: rank-based method; LS: least squares.

The improvements of the rank-based method over the LS method are calculated by means of the formula ( D_L − D_R )/ D_R . Here, D_L represents the Hausdorff distance between the LS B-spline approximation and denoised (i.e. outlier-free) point cloud data. D_R denotes the Hausdorff distance between the rank-based B-spline approximation and denoised point cloud data. According to Table 1, all of the Hausdorff distance values regarding the rank-based method are smaller than those concerning the LS method. This demonstrates quantitatively that the rank-based method has superior noise immunity compared to the LS method. Furthermore, it can be observed that the improvement takes the largest value which is marked in bold, when the largest of the numbers of control points (cp = 10) is used for B-spline approximation. Moreover, if the B-spline is neither over- nor under-fitted, the superiority of the rank-based method is proportional to the number of control points. Therefore, in the process of fitting mass point cloud data which frequently occur in large projects, such as surveys of long tunnels, the rank-based method, by virtue of its superior fitting efficiency, can significantly reduce the time and labor of data denoising required by standard LS methods.

Considering the accuracy of the B-spline modeling highly depends on the number of control points, the optimal parameters selection are significant for the deformation analysis based on TLS measurement. Table 1 reveals that these three data have reached the maximum value of improvement when cp = 10. Taking the fitting speed of large structures into consideration, the optimal number of control points for the tunnel cross section is 10. In future, we will carry out model selection method in order to check the optimal number of control points of the B-spline modeling. The benefit of the proposed method is the automation and robustness in the parametric surface modeling without point cloud preprocessing which could significantly reduce time- and labor-consuming for large-scale structures.

Conclusion

In this article, a rank-based algorithm, which is a robust and intelligent method immune to outliers and data gaps, was adopted for B-spline approximation and applied in the context of tunnel profile analysis. Outliers and data gaps of point cloud data generally lead to imprecision of the fitted B-spline model. The traditional method of deleting the outliers before the approximation is time- and labor-consuming for large scale structures. To alleviate these burdens, we use the rank-based Wilcoxon method to adjust B-spline control points, which can automatically be immune to the influence of outliers and data gaps, thus, reconstruct robustly the tunnel models. In summary:

1. To date, less consideration has been given to the robustness of control point determination. In this article, the rank-based Wilcoxon method is employed in the particular context of approximating tunnel profile measurements, and the results are compared with those obtained via the LS method.

2. Qualitative comparisons between the rank-based and LS methods are presented with various numbers of control points. The result imply that the rank-based method has a comparative advantage in the large projects, especially for high requirements of time-consuming and at the same time need to consider the robustness.

3. Hausdorff distance is involved to analyze the improvement of the rank-based to LS method quantitatively and the relation between the improvement and the number of control points was investigated.

4. The rank-based method can reduce the time and labor cost of denoising for mass point cloud data fitting, thereby this method can improve the fitting efficiency significantly, especially for many large projects, such as extra-long tunnels and bridges.

Consequently, the robustness of the rank-based Wilcoxon method contributes significantly to the intelligence of the free-form modeling to the computational efficiency of handling large-scale projects. In a future study, we will investigate the parameters of the rank-based method and compare them with other robust methods. Artificial intelligence technology will be combined for an automatic identification of parameters of the rank-based method for intelligent modeling.