Trans. Nonferrous Met. Soc. China 25(2015) 3165-3172

Reconstructing DEM using TLS point cloud data and NURBS surface

Bing SONG^1,2, Nan-shan ZHENG^1,2, Dong-wei LI¹, Ran-li CHEN², Liang LI^1,2

1. School of Environment Science and Spatial Informatics, China University of Mining and Technology, Xuzhou 221116, China;

2. Key Laboratory for Land Environment and Disaster Monitoring, National Administration of Surveying, Mapping and Geoinformation (NASG),Xuzhou 221116, China

Received 14 December 2014; accepted 28 July 2015

Abstract: Underground coal mining inevitably results in land surface subsidence. Acquiring information on land surface subsidence is important in the detection of surface change. However, conventional data acquisition techniques cannot always retrieve information on whole subsidence area. This study focuses on the reconstruction of a digital elevation model (DEM) with terrestrial laser scanning (TLS) point cloud data. Firstly, the methodology of the DEM with terrestrial 3-dimensional laser scanning is introduced. Then, a DEM modeling approach that involves the application of curved non-uniform rational B-splines (NURBS) surface is put forward. Finally, the performance of the DEM modeling approach with different surface inverse methods is demonstrated. The results indicate that the DEM based on the point cloud data and curved NURBS surface can achieve satisfactory accuracy. In addition, the performance of the hyperbolic paraboloid appears to be better than that of the elliptic paraboloid. The reconstructed DEM is continuous and can easily be integrated into other programs. Such features are of great importance in monitoring dynamic ground surface subsidence.

Key words: terrestrial 3D laser scanning; NURBS surface; DEM modeling; performance analysis

1 Introduction

The terrestrial laser scanner (TLS), which is also referred to as LiDAR, uses a time-of-flight distance measurement of an infrared laser pulse that is reflected on a scanned object. A point cloud is then obtained and holds information such as the x, y, and z coordinates as well as the intensity of the reflected pulse for each reflected point [1]. The TLS offers the possibility of acquiring 3-dimensional terrain information with high accuracy and high spatial resolution. It can also contribute greatly to the knowledge on both the surface shape and the kinematics of land subsidence and provide data that geologists, geomorphologists, and geo-technics can use to interpret the phenomenon. The application of the TLS technique in land subsidence analysis is particularly interesting because it can produce a digital elevation model (DEM) of land surfaces with high accuracy and precision. Moreover, TLS devices have the advantage of providing a large amount of high-resolution data within a short period and allowing the development of precise and detailed descriptions of scanned areas [2,3]. Thus, such devices have been increasingly used to track the evolution of natural surfaces. The existing applications mainly include landslide and rockfall dynamics as well as coastal cliff erosion [4-7].

Land subsidence induced by coal mining is a complex spatio-temporal process. Monitoring land surface subsidence is known to provide a basis for studying such a complex process. Thus, traditional leveling, total stations, GPS-RTK, and GNSS precise point positioning have been used to monitor mining subsidence, but these methods can only acquire point coordinates or height information [9,10]. Although several approaches for deformation monitoring have been found to yield highly accurate results, the point-based information they acquire is still considered as a limitation [11,12]. Both DInSAR and the persistent scatterer interferometry (PSI) technique with programmed high-resolution TerraSAR-X data set have been utilized in a coal mining area. However, the typically high rate of subsidence and low density of detected persistent scatterers limit the application of the PSI technique, and the scarcity of sufficient SAR data sets presents another dilemma [13].

As a new measurement technique developed in recent years, the TLS can quickly achieve the whole basin monitoring of surface subsidence caused by coal mining. This device can also overcome the defects of the point-based observation mode. However, a weak topological relationship exists among the large amount of point clouds obtained by the TLS because of the lack of manual intervention, which is a challenge for the reconstruction of surface DEMs based on discrete points [14,15].

The cloud point model of terrestrial laser scanning is a geometric model comprising elements of discrete points. This model cannot be used directly, and the land surface must be reconstructed through modeling such as numeric description and graphic representation. The modeling methods for land surface DEMs are generally based on the point, Delaunay triangulation, grid, or mixed approaches. The triangulated irregular network methods in particular have been widely used in DEM reconstruction because of their superior performance in terrain fitting [16,17]. The triangulation network methods are characterized with good flexibility and perfect boundary adaptability. However, the curved surface created by a triangulation network method can hardly be described by mathematical formulas. In addition, such curved surface is difficult to operate, thus making DEMs inconvenient to be utilized when extracting land subsidence or deformation information. Curved non-uniform rational B-splines (NURBS) surfaces have been widely applied in industrial design, scene modeling, and reverse engineering because of their favorable continuity, smoothness, affine invariance, perspective invariance, and local control characters [18]. NURBS is a parametric representation method that requires lesser storage space than traditional methods and thus allows the possibility to improve the management of extensive terrain data.

Coal is expected to remain as the primary energy source of China in the coming decades. However, underground coal mining inevitably causes land subsidence. Thus, studying the damage and deformation law of strata and surface caused by coal mining within the context of maintaining coal mining resources is significant. To analyze mining subsidence and the resulting ecological damage, DEM reconstruction within coal mining subsidence areas is essential.

This work aims to present an elaborate DEM reconstruction methodology based on TLS point cloud data with a curved NURBS surface. The methodology is presented on the basis of quadrilateral surface domain parameters for land subsidence induced by underground coal mining.

2 Preprocessing point cloud data acquired with TLS

The point cloud data were acquired with a Trimble GX200 3D laser scanner with a wavelength of 1540 nm. The laser scanner featured a 40°× 40° field of view, and the beam divergence was 170 μrad. One centrally located scanner position was used in each field plot (Fig. 1).

To obtain the topographic information of a large area, data acquisition was performed according to the following procedures: 1) topographic investigation and control point arrangement; 2) control measurement for georeferencing; 3) target setting; 4) laser scanning of terrain affected by underground coal mining; 5) target measurement and extraction [14].

After obtaining the surface point cloud data via terrestrial 3D laser scanning, pre-processing steps for point clouds, such as noise mitigation, point cloud registration, data compaction and surface reconstruction, must be conducted.

Fig. 1 Coal mining subsidence monitoring with terrestrial laser scanner

An important consideration in maximizing the use of TLS point clouds is data processing, which involves the removal of erroneous points, point cloud registration, data segmentation and thinning, mesh generation, and post processing. To obtain the points corresponding to the DEM, raw point clouds must be cleaned. In the current research, an effective mean shift algorithm based on statistical iteration was used to remove the noise within the point cloud data [19,20]. Moreover, the weights of the sampling points were assumed to be equal. The iteration was performed by Eq. (1).

                      (1)

where x_i is the value of the ith point within circle O and h denotes the specified radius of circle O. The kernel function g(·) can be set as the Epanechnikov kernel, a uniform kernel, a normal kernel, or other kernels.

A mean shift algorithm was adopted to identify the point with the maximum density in circle O. After eliminating the noise within the point cloud data, the resultant data were chosen as input for generating the land surface DEM.

3 DEM modeling with NURBS surface

Several approaches to the reconstruction of real models using 3D laser scanning data have been proposed in the fields of computer vision and computer graphics over the last two decades. A number of strategies for reducing data volume have also been suggested. Modeling itself, especially model reconstruction with the aid of analytical surfaces, can be considered as a means to reduce data volume, and fitting surfaces can be regarded as a means to reduce noise. Although the Bezier curve is very flexible and adaptive when applied to 3D modeling, it must be transformed to the standard NURBS, in which case all precision is lost.

NURBS is a mathematical model in computer graphics that is commonly used to generate and represent curves and surfaces. This model offers great flexibility and precision in the handling of both analytic (surfaces defined by common mathematical formulas) and modeled shapes. The modeled surface can be efficiently handled by computer programs while allowing easy human interaction. NURBS surfaces are functions of two parameters mapping to a surface in 3D space. The shape of the surface is determined by control points. NURBS surfaces can represent simple geometrical shapes in compact form [21,22].

In this research, the NURBS surface was adopted for land surface modeling, i.e., terrain DEM. The performances of this model are provided in the following section.

3.1 NURBS surface

A k×l NURBS surface can be described by three kinds of equivalent expressions, namely, rational fraction, rational base function, and homogeneous coordinates [21-23].

The rational fraction is defined by Eq. (2). Its control vertices, weight factors, and knot vectors are computed by explicit formulas [24].

           (2)

where the control points d_i,j (i=0, 1, …, m; j=0, 1, …, n) present a topological mesh and ω_i,j represents the corresponding weight. The NURBS surface is the tensor product of two NURBS curves. Hence, two independent parameters u and v (with indices i and j, respectively) are used. N_i,k(u) and N_j,l(v) are the B-spline basis functions derived from the Cox–de Boor recursion formula in terms of knot vectors  and , respectively.

The rational basis function is defined as follows:

                  (3)

where  is a bivariate rational basis function defined by Eq. (4):

         (4)

The rational bivariate basis function  is the extension of the variable B-spline function. It has the same function graph and geometric properties as those of the irrational B-spline function. Such properties include local support, normativeness, and differentiability. Meanwhile, ω_i,j=1, and  is equal to .

The homogeneous coordinates are described in Eq. (5).

    (5)

where  and represents the control points with the corresponding weights or the homogeneous coordinates of control vertices d_i,j.

To generate a NURBS surface, a few parameters such as control vertices d_i,j, weight factors ω_i,j, and knot vectors U, V must be determined, along with their orders k and l. The weight factor ω_i,j is a shape parameter that only affects the curved surface within the rectangular domain , . When ω_i,j increases, the curved surface approaches the control points, and vice versa.

Similar to the curved non-rational B-spline surface, the NURBS surface can be classified into different types according to the knot vector along the direction of each parameter. The knots with multiplicity k+1 are generally selected to be the endpoints of each knot vector. Thus, the four corner points of the curved NURBS surface are exactly located at the corner points of the control mesh. Moreover, the unilateral partial derivative in the corner points of the curved surface is the same as that in the endpoints of the boundary.

3.2 Surface inverse algorithm

The inverse of the NURBS surface is the key for DEM reconstruction. It mainly involves the construction of a k×l cubic B-spline interpolation surface according to the provided data points P_i,j (i=0, 1, …, m; j=0, 1, …, n) with a topology matrix, which is important in the search for unknown interpolation surface control points d_i,j (i=0, 1, …, m+k-1; j=0, 1, …, n+l-1). Despite the fact that the surface inverse generally adopts the reverse solution for tensor product surfaces, the surface inverse is simplified into two phases of curve inverse. Thus, the B-spline interpolation surface equation is adopted as follows [23,24]:

            (6)

Suppose that Then, Eq. (6)

can be transformed into

                     (7)

The control points within the curve equation can be replaced with the m+k control curves c_i(v). If parameter v is fixed, then c_i(v) from m+k can be computed accordingly. Given that these points are regarded as control points, the interpolation curves with parameter u are determined. These curves are also called isoparametric curves. When the traversal of v within its domain is performed, the interpolation surface can be depicted with a large amount of isoparametric curves. The control points  within the n+l section curves can be inversed via B-spline interpolation, in which each curve has m+k control points. The control points are expressed as Eq. (8).

(i=0, 1, …, m; j=0, 1, …, n)                (8)

Other control curves must be determined to define the control points  within the section curves. These control curves can be calculated with Eq. (9) on the basis of the m+k interpolation curves.

                       (9)

Within the B-spline interpolation curve, (m+k)× (n+l) control points d_i,j (i=0, 1, …, m+k-1; j=0, 1, …, n+l-1) exist. These control points can be calculated with Eqs. (8) and (9). Finally, the surface inverse achieved via B-spline interpolation can be performed by substituting the control points d_i,j and the parameters u and v into the B-spline interpolation formula derived for curved surfaces. Figure 2 shows the surface inverse process in detail. The result is shown in Fig. 3.

Fig. 2 Surface inverse algorithm

Fig. 3 Surface inverse through B-spline curved surface

3.3 DEM reconstruction

To effectively reconstruct the DEM of topographic surfaces, data segmentation must be conducted. In this work, the NURBS surface algorithm is adopted to extract characteristic lines from the topographic data obtained from terrestrial laser scanning. For a continuous topographic surface with minor curvature change, the terrain surface curvature detection technique is employed for data segmentation. Considering the occurrence of fractures or large curvatures in the terrain, the contour detection method is used. The terrain is also divided into a few surfaces. In other words, the terrain can be subdivided. For any subdivision, the terrain can be meshed. Figure 4 shows the surface patches with meshes.

After data segmentation and mesh processing, the surface DEM can be ultimately generated with the NURBS interpolation (Fig. 5).

Fig. 4 Details of meshes of surface paches

Fig. 5 DEM reconstruction with NURBS surface

4 Performance analysis of DEM modeling

Given the few errors within the data acquired by the TLS, an experimental research was conducted to evaluate the performance of the DEM. The simulated elevation errors were used as inputs for DEM modeling for reconstruction considering instrument specification.

Firstly, the DEM model was generated with the application of the TLS point cloud data after being preprocessed. In this model, the projection range was defined in [-3, 3]×[-3, 3] within the horizontal plane. Both the elliptic paraboloid z=(x²+y²)/18 and the hyperbolic paraboloid z=(x²-y²)/18 were used to shape the valley and ridge of the topographic terrain, with the absolute elevation differences being less than 1 (Figs. 6 and 7). The aim was to assess modeling precision with clean, error-free data.

Fig. 6 Elliptic paraboloid for shaping topographic terrain

Fig. 7 Hyperbolic paraboloid for shaping topographic terrain

Secondly, after uniformly sampling the digital terrain model with an interval distance of 1, 7×7 sampling points were obtained. These points were randomly added with certain errors that showed a normal distribution elevation. Given that the nominal accuracy of a single point obtained by the Trimble GX200 laser scanner is ±6 mm/50 m, the error range could be within ±12 mm while the sampling range is lower than 50 m. A curved quasi-uniform double cubic B-spline surface was chosen to shape the terrain. Firstly, 9×9 control points were inverted, and then both elliptic paraboloid and hyperbolic paraboloid were interpolated according to the curved B-spline surface. The results show no obvious differences after the incorporation of the normal error (Figs. 8 and 9). Although observation and simulated errors are included, the interpolated curved surfaces are acceptable in terms of their shapes.

Fig. 8 Interpolation of elliptic paraboloid

Fig. 9 Interpolation of hyperbolic paraboloid

For performance comparison purposes, 19×19 check points were uniformly selected from the fitted terrain surface and their surface elevations were computed base on the fitted surface model. The differences between the fitted surface and the original elliptic paraboloid as well as those between the fitted surface and the hyperbolic paraboloid were viewed as true errors. The distributions of errors from the elliptic paraboloid and hyperbolic paraboloid are shown in Figs. 10 and 11, respectively.

Fig. 10 DEM modeling error with elliptic paraboloid

Fig. 11 DEM modeling error with hyperbolic paraboloid

The comparisons in terms of two scenarios, namely, error within the whole domain and error of the domain unaffected by vertices with multiplicity, are shown in Tables 1 and 2, respectively. The results indicate that both paraboloids are effective when applied in terrain surface reconstruction for processing terrestrial laser scanning data. Notably, the performance of the hyperbolic paraboloid appears to be better than that of the elliptic paraboloid.

Table 1 Error of whole domain (mm)

Table 2 Error of domain unaffected by vertices with multiplicity (mm)

The consideration of data errors increases the total errors of the DEM. However, the modeling errors within the whole domain are found to be equivalent to the data source errors. Such errors do not greatly affect the DEM. Although the accuracy of modeling at the center of the domain is high, given that this region is unaffected by points with multiplicity, the precision of input data appears to be particularly important.

5 Conclusions

1) Terrestrial 3D laser scanning is a new approach to DEM reconstruction. It can be employed to extract information on land surface deformation in underground coal mining areas.

2) The curved NURBS surface algorithm based on quadrilateral surface domain parameters was studied for DEM modeling. The performance of the DEM reconstruction was analyzed using errors from both the data source and the modeling algorithm. The results indicate that the precision can reach a couple of millimeters and that high precision occurs at the domain unaffected by knots with multiplicity.

3) The application of cubic quasi-uniform B-spline interpolation reduces the modeling accuracy in the knots with certain multiplicity.

4) The approach does not pay significant attention to specific topographic terrains such as discontinuous terrains and to the efficiency of control point searching. Further study in this area thus must be conducted.

References

[1] DIMITRI L, NICOLAS B,  L. Accurate 3D comparison of complex topography with terrestrial laser scanner: Application to the Rangitikei canyon (N-Z) [J]. ISPRS Journal of Photogrammetry and Remote Sensing, 2013, 82: 10-26.

[2] SLOB S, HACK R. 3D terrestrial laser scanning as a new field measurement and monitoring technique [C]//Engineering Geology for Infrastructure Planning in Europe: A European Perspective. Heidelberg: Springer, 2004: 179-189.

[3] PIROTTI F, GUARNIERI A, VETTORE A. State of the art of ground and aerial laser scanning technologies for high-resolution topography of the earth surface [J]. European Journal of Remote Sensing, 2013, 46: 66-78.

[4] WAWRZYNIEC T F, MCFADDEN L D, ELLWEIN A, MEYER G, SCUDERI L, MCAULIFFE J, FAWCETT P. Chronotopographic analysis directly from point-cloud data: A method for detecting small, seasonal hill slope change, Black Mesa Escarpment, NE Arizona [J]. Geosphere, 2007, 3(6): 550-567.

[5]  A, CALVET J, VILAPLANA J M, BLANCHARD J. Detection and spatial prediction of rockfalls by means of terrestrial laser scanner monitoring [J]. Geomorphology, 2010, 119(3-4): 162-171.

[6] GUZZETTI F, MONDINI A C, CARDINALI M, FIORUCCI F, SANTANGELO M, CHANG K T. Landslide inventory maps: New tools for an old problem [J]. Earth-Science Reviews, 2012, 112: 42-66.

[7] BARBARELLA M, FIANI M. Monitoring of large landslides by terrestrial laser scanning techniques: Field data collection and processing [J]. European Journal of Remote Sensing, 2013, 46: 126-151.

[8] OLSEN M J, JOHNSTONE E, KUESTER F, DRISCOLL N, ASHFORD S A. New automated point-cloud alignment for ground-based light detection and ranging data of long coastal sections [J]. Journal of Surveying Engineering, 2011, 137(1): 14-25.

[9] HE Guo-qing, YANG Lun, LING Geng-di, JIA Cai-feng, HONG Du. Mining subsidence sciences [M]. Xuzhou: China University of Mining and Technology Press, 1991. (in Chinese)

[10] GAO Jing-xiang, LIU Chao, WANG Jian, LI Zeng-ke. A new method for mining deformation monitoring with GPS-RTK [J]. Transactions of Nonferrous Metals Society of China, 2011, 21(S3): s659-s664.

[11] CAI Chang-sheng, LUO Xiao-min, ZHU Jian-jun. Modified algorithm of combined GPS/GLONASS precise point positioning for applications in open-pit mines [J]. Transactions of Nonferrous Metals Society of China, 2014, 24(5): 1547-1553.

[12] ZHENG Nan-shan, CAI Liang-shi, BIAN He-fang, LIN Cong. Hybrid particle filtering algorithm for GPS multipath mitigation [J]. Transactions of Nonferrous Metals Society of China. 2014, 24(5): 1554-1561.

[13] LIU Zhen-guo, BIAN Zheng-fu, LEI Shao-gang, LIU Dong-lie, SOWTER Andrew. Evaluation of PS-DInSAR technology for subsidence monitoring caused by repeated mining in mountainous area [J]. Transactions of Nonferrous Metals Society of China, 2014, 24(10): 3309-3315.

[14] HUANG Cheng-liang. Application of 3D laser scanning technology in monitoring building deformation [D]. Xuzhou: China University of Mining and Technology, 2009. (in Chinese)

[15] SHI You-feng, GAO Xi-feng. Quasi-uniform B-spline function using in digital terrain model [D]. Chengdou: Southwest Jiaotong University, 2008. (in Chinese)

[16] TAN Ren-chun, DU Qing-yun, YANG Pin-fu, ZHANG Shan-shan. Optimized triangulation arithmetic in modeling terrain [J]. Geomatics and Information Science of Wuhan University, 2006, 31(5): 436-439. (in Chinese)

[17] HUANG Zheng-ge, CHEN Jian-jun, ZHENG Yao. Triangulated irregular network based chunk gridding algorithm for terrain rendering [J]. Journal of Zhejiang University: Engineering Science, 2009, 43(10): 1939-1945. (in Chinese)

[18] GERSHENFELD N. The nature of mathematical modeling [M].UK: Cambridge University Press, 1999.

[19] COMANICIU D, RAMESH V, MEER P. The variable bandwidth mean shift and data-driven scale selection [C]//Proceedings of the Eighth IEEE International Conference on Computer Vision. Vancouver, BC, Canada: IEEE, 2001: 438-445.

[20] NING Ji-feng, ZHANG Lei, ZHANG David, WU Cheng-ke. Robust mean-shift tracking with corrected background weighted histogram [J]. IET Computer Vision, 2012, 6: 62-69.

[21] PIEGL L. On NURBS: A survey [J]. IEEE Computer Graphic and Application, 1991, 11(1): 55-71.

[22] SCHAEFER S, GOLDMAN R. Non-uniform subdivision for B-splines of arbitrary degree [J]. Computer Aided Geometric Design, 2009, 26(1): 75-81.

[23] W. Cubic B-spline curves and surfaces in computer aided geometric design [J]. Computing, 1977, 19(1): 29-34.

[24] WANG Xiao, CHEN Xiu-wan. Three dimensional terrain visualization on the internet using NURBS method [C]//IEEE International Geoscience and Remote Sensing Symposium, IGARSS 2008. Boston: IEEE, 2008: 1341-1344.

基于地面激光扫描点云数据和NURBS曲面的DEM重建

宋 冰^1,2，郑南山^1,2，厉东伟¹，陈冉丽²，李 亮^1,2

1. 中国矿业大学 环境与测绘学院，徐州 221116；

2. 国家测绘地理信息局 国土环境与灾害监测重点实验室，徐州 221116

摘  要：煤炭资源的井工开采不可避免地引起地表沉陷，快速高精度获取地面沉降信息尤为重要。然而，传统的数据获取方法难以满足要求，本研究基于地面激光扫描点云数据重建数字高程模型（DEM）。在介绍地面三维激光扫描重建DEM原理的基础上，提出一种融合NURBS曲面的DEM重建方法，并分析不同曲面拟合模型重建DEM的性能特征。结果表明：基于点云数据和NURBS曲面重建的DEM能够较好地满足精度要求，而马鞍面拟合方法的精度高于椭圆抛物面拟合方法的精度。采用NURBS曲面算法重建的DEM表面连续，易于互操作，对监测分析煤矿开采引起的地表动态沉降过程具有重要意义。

关键词：地面三维激光扫描；NURBS曲面；DEM建模；精度分析

(Edited by Mu-lan QIN)

Foundation item: Project (51174206) supported by the National Natural Science Foundation of China; Project (2014ZDPY29) supported by the Fundamental Research Funds for the Central Universities; Project (SZBF 2011-6-B35) supported by the Priority Academic Program Development of Higher Education Institutions (PAPD) of Jiangsu Province, China

Corresponding author: Nan-shan ZHENG; Tel/Fax: +86-516-83591330; E-mail: znshcumt@163.com

DOI: 10.1016/S1003-6326(15)63947-4