An Accurate TLS and UAV Image Point Clouds Registration Method for Deformation Detection of Chaotic Hillside Areas

Abstract

Deformation detection determines the quantified change of a scene’s geometric state, which is of great importance for the mitigation of hazards and property loss from earth observation. Terrestrial laser scanning (TLS) provides an efficient and flexible solution to rapidly capture high precision three-dimensional (3D) point clouds of hillside areas. Most existing methods apply multi-temporal TLS surveys to detect deformations depending on a variety of ground control points (GCPs). However, on the one hand, the deployment of various GCPs is time-consuming and labor-intensive, particularly for difficult terrain areas. On the other hand, in most cases, TLS stations do not form a closed loop, such that cumulative errors cannot be corrected effectively by the existing methods. To overcome these drawbacks, this paper proposes a deformation detection method with limited GCPs based on a novel registration algorithm that accurately registers TLS stations to the UAV (Unmanned Aerial Vehicle) dense image points. First, the proposed method extracts patch primitives from smoothed hillside points, and adjacent TLS scans are pairwise registered by comparing the geometric and topological information of or between patches. Second, a new multi-station adjustment algorithm is proposed, which makes full use of locally closed loops to reach the global optimal registration. Finally, digital elevation models (DEMs, a DEM is a numerical representation of the terrain surface, formed by height points to represent the topography), slope and aspect maps, and vertical sections are generated from multi-temporal TLS surveys to detect and analyze the deformations. Comprehensive experiments demonstrate that the proposed deformation detection method obtains good performance for the hillside areas with limited (few) GCPs.

1. Introduction

Deformation detection has been defined as the surveying of a region in different epochs and the identification of geometric differences based on the captured data [1]. Its purpose is to determine if the geometric state of the area of interest has changed as well as the quantified change. The detection of the deformations of hillside areas has been an active and important task in mitigating hazards and property loss in the field of earth observation. The continuous occurrence of hazard events increased the demand for new techniques for further measurement and analysis. From the past to the present, traditional surveying techniques (e.g., static GPS, Real-time kinematic: RTK, total station survey) and remote sensing based techniques, such as close-range photogrammetry, unmanned aviation vehicles (UAVs), and laser scanning, are widely used to capture the spatial data of hillside areas by integrating them with one another or applying them individually. TLS is gaining attention due to the high point density and spatial resolution data that can be acquired in a short time. Its application in unstable hillside areas may provide important information for achieving a more complete knowledge of the active processes that could forecast potential geo-disasters [2].

In the literature of deformation detection based on laser scanning points or dense photogrammetric image points, both the registration and deformation detection techniques have been covered from different angles. Therefore, these two aspects are reviewed in the following two sub-sections.

1.1. Registration Techniques of Point Clouds

The registration technique of 3D points has been applied in various applications (e.g., mobile robotics, object recognition, and 3D reconstruction). The registration methods have been widely classified into two categories: Coarse and fine registration. The coarse registration provides initial values for fine registration to avoid falling into the local minimum. The widely used coarse registration methods include: Point signature [3], spin image [4], shape context [5], fast point feature histograms (FPFH) [6], normal vectors/principal curvature/curvature change [7,8,9], geometric primitives [10,11,12], etc. The most popular fine registration method is the ICP algorithm [13,14] and its variants [15,16,17,18]. ICP is the iterative closest point algorithm. It optimizes the transformation between two point clouds, by minimizing the square errors between correspondences. Furthermore, there are also many other recent methods, including 4PCS (4-points Congruent Sets) [19], SLAM (Simultaneous Localization and Mapping) [20], probabilistic registration [21], genetic algorithm [22], etc.

Both coarse and fine registration algorithms can be further categorized into pairwise registration and multi-station registration methods. Pairwise registration determines the transformation between two adjacent point clouds while the multi-station registration takes multiple point clouds into account to correct the cumulative errors introduced by pairwise registration. Most of the above methods belong to pairwise registration. More recently, Tombari et al. proposed a 3D descriptor that encodes histograms of the normal vectors of points within a support, describing the local reference frame of a feature point [23]. Similarly, Guo et al. constructed a unique local reference frame and applied rotational projection statistics to describe it [24]. Contrastingly, Geng et al. defined a new function to measure feature similarity by calculating the distance function based on the neighborhood constraints of geometrical features [25].

The above pairwise registration methods have demonstrated their superior performance for two adjacent point clouds. However, they are ineffective for multiple stations since cumulative errors are not considered adequately. To address this issue, various multi-station registration methods have been proposed. Williams et al. applied a weighted matrix encoding the correspondence information to calculate the optimal transformation iteratively [26]. Huber et al. constructed a connected sub-graph to order the procedures of operation to reduce the cumulative errors [27]. Similarly, Zhu et al. sequentially registered each station to a model constructed by other stations [28]. Their transformations are sequentially refined by traversing all stations. Guo et al. proposed a weighted motion method to provide more reliable pairwise registration results, improving the accuracy of multi-station registration [29]. Moreover, the multi-station adjustment algorithm is another efficient method to correct cumulative errors. However, the existing multi-station adjustment methods mainly solve this problem based on a global closed loop constructed by multiple stations [30,31,32].

For mountainous areas, the existing pairwise registration methods do not work well for the following reasons: (1) Registration primitives encoding distinctive information are difficult to select; (2) the similarity of different points’ local information affects the determination of correspondences. Besides, in practice, it is difficult for TLS stations to form an efficient closed loop. The application of existing multi-station adjustment methods leads to compromised results that restrict their applications (e.g., deformation detection). Therefore, it is necessary that an efficient registration algorithm is developed.

1.2. Deformation Detection Methods

Various remote sensing techniques have been applied in the field of deformation detection. Passive space-borne imaging based solutions can collect optical images with a centimeter ground resolution. Then, traditional photogrammetry based data processing pipelines are adopted to generate ortho-images, DSMs (digital surface models), and DTMs (digital terrain models) for further deformation measurement. Compared with passive space-borne imaging based solutions, airborne based and terrestrial laser scanning provide a direct way to rapidly capture high density 3D surface coordinates, which can be processed using filtering operations, such as a slope-based filter (e.g., Vosselman, [33]), to separate ground points and non-ground points for DTM generation. Although photogrammetric DTMs generated from passive space-borne images are not usually as accurate and precise as LiDAR-based (Light Detection and Ranging) DTMs [34], they can still achieve a centimeter accuracy that meets the requirements of various applications [35].

Generally, traditional airborne- and satellite-based remote sensing techniques, such as differential InSAR (Interferometric Synthetic Aperture Radar), are suitable for the detection of deformations over large areas (e.g., multiple square kilometers, [36,37,38]). On the one hand, since each scene image of InSAR has wide coverage, the average cost of a large area is low. However, if used for a small specific area, the cost is relatively high. On the other hand, the temporal resolution (collection frequency) of InSAR data may affect its application of deformation detection. Because many InSAR datasets have a long or fixed revisit time, if fast-changing or large deformations occurred during the acquisition period, the problem of decorrelation may be caused [39,40,41,42].

The use of radio controlled mini-UAVs to collect the data of small specific areas is convenient, flexible, and cheap. In the last five years, they have been widely used to collect high geometric and temporal resolution data (e.g., images and point clouds) with a consumer grade optical camera or a light weight laser scanner (e.g., Velodyne VLP-16, 0.83 kg; RIEGL VUX-1LR, 3.5 kg). Thus, the use of mini-UAVs for the detection of deformations in mountainous areas is possible and cost affordable. Related researches have been published [43,44,45,46]. However, mini-UAVs apply a top-down acquisition method to collect data, making it difficult to collect the data from some steep areas. Although an oblique camera can be mounted on mini-UAVs, it generally has a lower accuracy since more parameters will be introduced in the photogrammetric procedures [47].

Plus, TLS can provide huge amounts of measurements in a short time. It enables flexible collection of dense 3D points in areas of interest, generating a more detailed description of the slope topography, thus improving our understanding of deformations in hillside areas. Most existing methods set up various retro-reflective targets (GCPs) in the environment to ensure geo-referencing and the alignment of different views’ stations [48,49,50]. However, the GCPs should be deployed in each TLS survey. The deployment is generally time-consuming and labor-intensive, and is sometimes very difficult when dealing with steep terrain. Few studies have investigated the detection of deformations using limited GCPs based on TLS points.

Overall, this paper utilizes the UAV dense image points as a reference frame work to register multi-temporal TLS surveys automatically and accurately. By reducing the GCPs of each TLS survey, more efficient and economic deformation measurements of slope topography are obtained. The main contributions of the proposed method are as follows:

(1)

We propose an efficient pairwise registration algorithm based on patch primitives to register adjacent point clouds of hillside areas coarsely. The main feature of the method is that the information of the trend and topographic undulation of the mountains together with multi-scale information of each patch are applied to determine the correspondences robustly.

(2)

We propose a novel multi-station adjustment algorithm to accurately register TLS and UAV dense image points, which corrects the cumulative errors based on locally closed loops formed by adjacent stations. The introduction of virtual points makes the linearization of the condition equation group possible, reaching the global optimal alignment of all TLS stations.

(3)

Based on the registration techniques, we generate DEMs, slope and aspect maps, and vertical sections of multi-temporal TLS surveys for deformation detection and terrain morphological analysis, thus demonstrating the effectiveness of the high resolution deformation detection method with limited GCPs.

The rest of this paper is structured as follows. Following this introduction, Section 2.1 and Section 2.2 describe the study area and data acquisition. Section 2.4 and Section 2.6 give a detailed description of the proposed pairwise registration and multi-station adjustment method. Then, the registration algorithms and deformation detection results are validated and discussed in Section 3 and Section 4. Finally, the conclusions and future research directions are presented in Section 5.

2. Materials and Methods

2.1. Experimental Area

The experimental area is located in a mountainous area near the East Sea, China (Figure 1). The area covers about 0.15 km² and is distributed along the coastal zones with elevations ranging from 9.1 m to 109.0 m. The slope facing the East Sea is steep and harsh, and the top topography is relatively flat. A quad-rotor helicopter md4-1000 (Microdrones GmbH, Siegen, Germany) equipped with a digital camera consisting of an affordable remote sensing platform was adopted to collect the optical images over the study area. Besides, we set up TLS (Leica C10) near the coastal zones to collect 3D laser scanning points. The parameter descriptions of md4-1000, the digital camera, and the laser scanner are listed in

Table 1. Parameters for the md4-1000, digital camera, and Leica C10.

Equipment Types	Helicopter md4-1000	SONY ILCE-7R	Leica C10
Critical Parameters	Cruising speed: 15.0 m/s. Rate of climb: 7.5 m/s. Weight: 2650 g. Maximum load: 2000 g. Recommended load: 800 g. Flying time: over 50 min Flight range: min. 500 m with remote control.	Focal length: 35 mm. Resolution: 36.4 Megapixel. Picture size: 7360 × 4912 Pixel. Shutter speed: 1/8000 s. Weight: 407 g. Camera size: 127 × 94 × 48 mm.	Range: 300 m. Scan rate: up to 50,000 points/s. Spot size: 4.5 mm. Point spacing: 1 mm. Filed-of-view: Horizontal 360°; Vertical 270°. Laser class: 2.

. Based on the data of the flat areas (both UAV and TLS collect well), this paper takes the geo-referenced UAV image points as the framework, and registers multi-temporal TLS surveys to the UAV image points.

2.2. Data Acquisition: UAV Flying and Terrestrial Laser Scanning

We collected the optical images with a md4-1000 UAV equipped with SONY ILCE-7R. The geo-referenced UAV image points were used as the framework to align multi-temporal TLS surveys. The UAV flight was performed in September 2014, and the flight task was executed with a predefined flying routine (four strips) by a remote controlled computer (shown in Figure 2). The flight elevation was about 150 m with a ground resolution of 1.6 cm and the flight range was up to 20 minutes. The wind scale was less than force-3. Longitudinal overlap was over 70% and side overlap was over 55%. To attain a smaller image blur and better illumination, the shutter speed of the camera was set to 1/1000 s. Seven GCPs (two for check points) were set up to geo-reference the UAV dense image points accurately. These conditions together ensured the quality of the UAV image points (plane accuracy 3.2 cm, elevation accuracy 4.6 cm). Five periods of TLS surveys were performed from March 2013 to January 2015. The TLS (Leica C10) was set up near the coastal zone to capture the 3D point clouds of slope areas with an average point span of 0.08 m. For each survey, the TLS captures data with several scans, ranging from 3 to 6 stations.

Three main procedures were scheduled to ensure the quality of DEM (plane and elevation accuracy were less than 5.0 cm) generated using the drone-acquired data: (1) Designing the flight plan, referring to the desired flight height, ground sampling distance, and image overlap; (2) selecting the location of GCPs (well-distributed in bare-land, no-change area) based on geospatial data; and (3) placing and measuring GCPs on site. Seven GCPs were surveyed by GPS-RTK, the horizontal and vertical accuracy was expected to be within 0.10 m and 0.15 m, respectively.

2.3. Method Workflow

Two independent pipelines were deployed to process TLS point clouds and optical images before accurately registering TLS data to the UAV dense image points. The proposed method encompasses three key components: Patch based pairwise registration of TLS stations, multi-station adjustment of TLS and UAV points, and terrain deformation analysis. The workflow for the detection of deformation with the UAV and TLS data is shown in Figure 3.

2.4. Pairwise Registration between TLS Stations

2.4.1. Patch Primitive Extraction and Multi-Scale Descriptor Construction

For a chaotic mountainous environment, mismatched correspondences are easily introduced since many points’ local geometric information at a single (predefined) scale is similar. Apart from the local geometric information, the trend and topographic undulation of mountains is also important information. To make full use of the information, we proposed the extraction of patch primitives to register mountainous data robustly. To ensure the completeness and consistency of patches, the operation of 3D smoothing is necessary. A 3D Gaussian kernel [51] was applied to smooth the 3D point clouds, written as:

$$G(i;{\delta }^2)=\frac{1}{2\pi {\delta }^2}\exp \left(-\frac{1}{2{\delta }^2}{d}_i^2\right)$$

(1)

where d_i denotes the distance from a given point to its neighborhood point, p_i, and $\delta$ is the mean square root error. The larger $\delta$ is, the larger the smoothing region should be, and larger size features will be smoothed. Then, 3D point clouds were convolved with Gaussian kernel to realize the 3D smoothing. Based on the smoothed point clouds, the method below was used to extract patch primitives.

Select a seed point randomly, and take a fixed length as the clustering radius. The clustering radius is determined with respect to the median scene size of the stations, which is determined by the bounding box of each point cloud. Cluster points according to the normal vector of each point. The normal vector of each point is estimated by the principal components analysis (PCA) algorithm and the directions of normal vectors are unified by the consistent propagation algorithm. For each cluster obtained in the previous iteration, select the point with the largest curvature change, M_c, as the new seed point. The value of curvature change is calculated by the normal vectors of neighborhood points, as shown in Formula (2). Based on the new seed point, cluster again to get new clustered points. Repeat the above procedures until the clustered points are unchanged from the points of last iteration. For each patch, select the newest seed point as the feature point of this patch, and calculate the average normal vector as its normal vector. Figure 4 illustrates the curvature change rendering of two adjacent stations after smoothing. It demonstrates that the overlap areas have similar geometric information, suggesting that 3D smoothing ensures the consistency of patches from adjacent stations.

$$\begin{array}{l}c=\frac{1}{n}{\displaystyle \sum _{j=1}^n\Vert v_i-v_{neighbor,j}\Vert }\\ M_c=\frac{1}{n}{\displaystyle \sum _{j=1}^n\Vert c_i-c_{neighbor,j}\Vert }\end{array}$$

(2)

Note that mountainous features have multi-level details, and this multi-scale information of different features has significant differences that can be used to accurately distinguish them. To describe this multi-scale information of a feature, multiple feature responses (e.g., the curvature, curvature change, or other geometric information) are calculated across a set of discrete scales. In this paper, for each patch, we calculated the curvature change values of its feature point within various neighborhood ranges to comprehensively describe its multi-scale information.

The specific steps are as follows: According to the corresponding relationship between the original points and the smoothed points, find the original points corresponding to each feature point of patches. Then, for each corresponding point, search the original points within various neighborhood ranges (e.g., the neighborhood radii of four scales are set as: 1.0 m, 1.5 m, 2.0 m, and 2.5 m, respectively) to calculate its curvature changes, forming the multi-scale descriptor of each patch.

2.4.2. Matching Strategy and Transformation

To improve the correct matching rate of patches, we applied the triple matching strategy. The detailed procedures are as follows:

Step 1:

Select three non-colinear patches from the left and right station, respectively. For each station, calculate the mixed product of the normal vectors of patches, and calculate three distances between them.

Step 2:

Judge whether the angles between the patches of the left station are identical with that of the right station by comparing their mixed products. If not, repeat Step 1.

Step 3:

Compare the multi-scale information of each patch feature point from two stations for equality. For each scale, the difference threshold is determined by the mean square error of all points within the patch. If three pairs of similar points cannot be obtained, repeat Step 1.

Step 4:

Compare the corresponding distances of two stations for equality. If three correspondences are obtained, repeat the above steps until enough matched points are obtained.

The pseudo-code of tripe matching strategy is listed on the next page.

We calculated the rotation and translation parameters separately to obtain the transformation parameters robustly. After the above matching, any four non-collinear pairs of feature points from patches are selected to form two lines. Inspired by the method of [52], the transformation parameters were calculated as follows. Let $u_1,v_1$ be the vectors of two lines from the left station. Then, we have:

$$\{\begin{array}{l}{a}_1=u_1\\ b_1^{\prime }=v_1-(v_1\cdot a_1)a_1,b_1=b_1^{\prime }/\left|\right|b_1^{\prime }\left|\right|\\ c_1=a_1\times b_1\end{array}$$

(3)

where $(v_1\cdot a_1)$ indicates a dot product, Cartesian coordinate frame, $M_1=\left[a_1,b_1,c_1\right]$. Similarly, $M_2$ can be derived from the right station. The rotation matrix is:

$$R=M_1{M}_2^T$$

(4)

Further, the voting principle was used to determine the final rotation parameters. Then, we calculated the translation parameters by the distance differences between the perpendicular midpoints of two pairs of lines (midpoints are unique and easily obtained). Then, the translation parameters were optimized with the voting principle again. Finally, the transformation parameters were obtained as the initial values for fine registration (multi-station adjustment).

Algorithm Matching of triple patches

1. lPts and rPts are feature point array of triple patches selected randomly

2. lMixedP = (lPts[0].normal, lPts[1].normal, lPts[2].normal);

3. rMixedP = (rPts[0].normal, rPts[1].normal, rPts[2].normal);

4. //lPts[i].normal, rPts[j].normal are the normal vectors of corresponding patch, i∈[0,3], j∈[0,3]

5. // lMixedP and rMixedP are the mixed products of normal vectors of triple patches.

6. if | lMixedP- rMixedP | > Threshold then

7.  return

8. end if

9. //Match patches by the multi-scale informaiton of feature points

10. for i from 0 to lPts.size()-1

11.    for j from 0 to rPts.size()-1

12.     if isimilar(lPts[i].multiscales, rPts[j].multiscales ) = = true then

13.      onePair.left = i; onePair.right = j; PairArray.push_back(onePair);

14.      break

15.     end if

16.    end for

17. end for

18. if PairArray.size() != 3 then

19.  return

20. end if

21. //Match patches by the distances between feature points

22. flag = 0

23. for i from 0 to 2

24.   for j from i+1 to 2

25.     if distance(lPts[PairArray[i].left], lPts[PairArray[j].left]) = =

26.      distance(rPts[PairArray[i].right], rPts[PairArray[j].right])

27.      then flag++

28.     end if

29.    end for

30. end for

31. if flag != 3 then

32.  PairArray.erase(); //PairArray is the array of corresponding points

33.  end if

2.5. Dense Image Points Generation from UAV Optical Images

Compared with traditional aerial photogrammetric surveying, UAV provides a flexible and cost-effective solution for capturing high spatial resolution images, however, the images may have orientation changes and distortions because of possible trembling of the light weight UAV. The 3D data generated from the images, if characterized by a proper accuracy and resolution, can be considered a faithful numerical description of the slope topography [53], which can be helpful to evaluate the deformations of the slope.

First, the blurred images were discarded, and the lens distortions of the captured images were corrected. Starting from at least two images with known EO (EO refers to the exterior orientation, it describes the position and orientation of the camera when the image is taken with six parameters: Projection center coordinates and the rotations around axes), image matching techniques were applied to generate linking points, whose spatial intersection gives general 3D coordinates in the control points’ coordinate system. The resulting set of 3D points constitutes the so-called point cloud. The remaining images were imported into the Pix4D software [54] for aerial triangulation adjustment (construction of free network) and multi-view stereo technology to obtain dense image points. Then, the GCPs were used to geo-reference the dense image points to the control points’ coordinate system. After processing, the check points show that the achieved dense image points were convincing (plane accuracy 3.2 cm, elevation accuracy 4.6 cm).

The TLS scans were pairwise registered considering Section 2.4. Then, the coarse registered TLS scans were accurately aligned with the dense image points by the multi-station adjustment algorithm in the next section.

2.6. Multi-Station Adjustment Based on a Locally Closed Loop

The registration of multi-temporal TLS surveys to the UAV dense image points can transform them into a uniform spatial reference framework. Particularly, using the UAV dense image points as the reference framework will also reduce the GCPs of the filed survey. Patch based pairwise registration in Section 2.4 was applied to provide initial values between the UAV dense image points and the pairwise registered TLS point clouds. Few patch primitives were extracted from motion-changed regions (e.g., collapse areas), ensuring the robustness of the patch based coarse registration.

On the other hand, cumulative errors were introduced since TLS point clouds were pairwise registered. TLS stations and UAV dense image points do not form an effective closed loop, such that traditional multi-station adjustment methods cannot be used to correct the errors. To register the TLS and UAV data of a hillside area accurately, we constructed a condition equation group that can be conveniently linearized for multi-station adjustment based on locally closed loops formed by adjacent stations.

Three-fold, or upwards of three-fold, overlap usually exists between TLS stations and UAV data. Taking three-fold overlap, for example, to minimize the cumulative errors, multi-station adjustment should be conducted under the constraints of the correspondences from three overlap regions. Suppose three stations are: T₁, T₂, and T₃ (T₁ is fixed), respectively. Provided by the initial values, the nearest neighbor search algorithm is used to obtain the corresponding points from each overlap region. Suppose $P_l^{12}$ and $P_r^{12}$ are the corresponding points of T₁ and T₂, and $P_l^{23}$ and $P_r^{23}$ are the corresponding points of T₂ and T₃, respectively. According to the correspondences, condition equations can be formulated:

$$\begin{array}{l}scan1~2:P_l^{12}=R_{12}\cdot P_r^{12}+T_{12}\\ scan2~3:P_l^{23}=R_{23}\cdot P_r^{23}+T_{23}\end{array}$$

(5)

Here, R₁₂ and T₁₂ are the transformation parameters between T₁ and T₂, and R₂₃ and T₂₃ are the transformation parameters between T₂ and T₃. According to the theory of vision, Formula (5) can be linearized. For instance, the linearized error equation of T₁ and T₂ can be written as:

$$\begin{array}{l}scan1~2:V_{12}=A_{12}\cdot X_{12}-L_{12}\\ \{\begin{array}{l}{A}_{12}=\left[\begin{array}{l}100-Z_r^{12}0-Y_r^{12}\\ 0100-Z_r^{12}{X}_r^{12}\\ 001X_r^{12}{Y}_r^{12}0\end{array}\right]\\ X_{12}={\left[d\mathsf{\Delta }{X}_{12}d\mathsf{\Delta }{Y}_{12}d\mathsf{\Delta }{Z}_{12}d{\phi }_{12}d{\omega }_{12}d{\kappa }_{12}\right]}^T\\ L_{12}=\left[\begin{array}{l}{X}_l^{12}\\ Y_l^{12}\\ Z_l^{12}\end{array}\right]-R_{12}^0\cdot \left[\begin{array}{l}{X}_r^{12}\\ Y_r^{12}\\ Z_r^{12}\end{array}\right]-\left[\begin{array}{l}\mathsf{\Delta }{X}_{12}^0\\ \mathsf{\Delta }{Y}_{12}^0\\ \mathsf{\Delta }{Z}_{12}^0\end{array}\right]\end{array}\end{array}$$

(6)

Here,$R_{12}^0,\mathsf{\Delta }{X}_{12}^0,\mathsf{\Delta }{Y}_{12}^0,\mathsf{\Delta }{Z}_{12}^0$ are the initial transformation parameters, $d\mathsf{\Delta }{X}_{12},d\mathsf{\Delta }{Y}_{12},d\mathsf{\Delta }{Z}_{12}$, and $d{\phi }_{12},d{\omega }_{12},d{\kappa }_{12}$ are the residuals of the initial translation parameters and rotation angles. $P_r^{12}={\left[X_r^{12},Y_r^{12},Z_r^{12}\right]}^T$ are the observations and $V_{12}={\left[V_{12}^X{V}_{12}^Y{V}_{12}^Z\right]}^T$ are their corrections to balance Formula 6.

Suppose $P_l^{13}$ and $P_r^{13}$ are the corresponding points of T₁ and T₃. However, two sets of transformation parameters between T₁ and T₃ make it difficult for its condition equation to be linearized. To form new equations that can be linearized, we introduced a virtual point to transform the original condition equation into two condition equations, written as:

$$scan1~3:\{\begin{array}{l}{P}_v^{13}=R_{23}\cdot P_r^{13}+T_{23}\\ P_l^{13}=R_{12}\cdot P_v^{13}+T_{12}\end{array}$$

(7)

Here, $P_v^{13}$ is a virtual point in T₂ that corresponds to $P_r^{13}$ and $P_l^{13}$.

According to Formula (6), we constructed the error equation group as follows:

$$\begin{array}{l}scan1~2:V_{12}=A_{12}\cdot X_{12}-L_{12}\\ scan2~3:V_{23}=A_{23}\cdot X_{23}-L_{23}\\ scan1~3:\{\begin{array}{l}{V}_{13}^v=A_{23}^v\cdot X_{23}-L_{23}^v\\ V_{13}=A_{12}^v\cdot X_{12}-L_{12}^v\end{array}\end{array}$$

(8)

The matrix expression can be written as:

$$\left[\begin{array}{l}{V}_{12}\\ V_{23}\\ V_{13}^v\\ V_{13}\end{array}\right]=\left[\begin{array}{l}{A}_{12}0\\ 0A_{23}\\ 0A_{23}^v\\ A_{12}^{v}0\end{array}\right]\cdot \left[\begin{array}{l}{X}_{12}\\ X_{23}\end{array}\right]-\left[\begin{array}{l}{L}_{12}\\ L_{23}\\ L_{23}^v\\ L_{12}^v\end{array}\right]$$

(9)

Here, similar to Formula (6), $V_{23}$ are the corrections of $P_r^{23}={\left[X_r^{23},Y_r^{23},Z_r^{23}\right]}^T$, and $V_{13}^v$ and $V_{13}$ are the corrections of $P_r^{13}$ and $P_v^{13}$, respectively. These corrections are used together to balance Formulas (8) and (9).

Based on the initial values of registration, we calculated the residuals of R₁₂, T₁₂ and R₂₃, T₂₃ iteratively. For upwards of a three-fold overlap, multi-station adjustment can also be operated according to the above method. The locally closed loops are considered sequentially (the adjusted TLS stations are fixed) to optimize each TLS station separately, correcting the cumulative errors, and reaching the global optimal alignment of all TLS stations. Figure 5 shows the registration result of the UAV data and one TLS survey. Edge areas (red circles) demonstrate that good fusion of the TLS stations and UAV data was achieved, reflecting the effectiveness of the registration method.

2.7. DEM Generation and Deformation Comparison between Observations

DEMs of multi-temporal TLS datasets are usually used to analyze the deformations of areas of interest. To focus on the ground data of a hillside, we applied the progressive triangulated irregular network filtering algorithm [55] to filter the isolated vegetation before the construction of DEMs. The filter threshold varies based on the terrain smoothness, using a coarse-to-fine scheme. Regular grids of DEMs were created since it is the simplest and the most efficient approach in terms of storage and manipulation. The resolution of the generated DEM was 0.2 m (horizontal resolution), providing a detailed description of the scanned area. Since all the surveys were transformed into the control points’ coordinate system, it is feasible to compare the multi-temporal DEMs to measure and analyze the terrain deformations, and analyze the terrain surface of the areas of interest by marker pegs measured previously. The comparison was implemented in ArcGIS 10.5, where the elevation difference of each grid of the two-point clouds was measured to calculate the volumetric deformation of the changed areas.

In deformation detection, the slope and aspect are two vital geomorphic parameters, as they efficiently describe the relief and structure of the terrain surface. The slope of one surface point denotes the degree of the terrain surface inclination at this point, and the aspect represents the direction of the largest elevation change. For a given point on a terrain surface, z = f(x, y), the slope, S, and aspect, A, are defined as the function of the gradients in the X and Y directions (West-East and North-South):

$$\begin{array}{l}S=\arctan \sqrt{f_x^2+f_y^2}\\ A={270}^{\circ }+\arctan \left(\frac{f_y}{f_x}\right)-{90}^{\circ }\frac{f_x}{\left|f_x\right|}\end{array}$$

(10)

where f_x and f_y are the gradients at the West-East and North-South directions, respectively. The key step in determining the slope and aspect of one surface point is the calculation of f_x and f_y. For one grid DEM, f_x and f_y are obtained by the finite differential of a local moving window (e.g., 7 × 7 window).

3. Results

3.1. Registration Experiments

3.1.1. Extraction of Patch Primitives

According to Section 2.4.1, to ensure the completeness and consistency of the patches extracted, the TLS point clouds were first smoothed by the 3D Gaussian kernel. Then, for the point clouds of each station, patch primitives were extracted, where the clustering radius of patch was set to 2.0 m according to the method mentioned above. The extraction results of six stations are shown in Figure 6. In Figure 6, clusters of points of each color represent one patch primitive (colors were rendered randomly).

Figure 6 demonstrates that the trend and topographic undulation of the hillside area is described effectively by numerous patches in different directions. It also illustrates that the iterative patch extraction method ensures the similarity of patches in overlapping regions. The extraction results of six TLS stations suggest that lots of patches were extracted, ensuring a sufficient number of corresponding patches in the overlapping regions.

3.1.2. Accuracy Evaluation of Pairwise Registration and Multi-Station Adjustment

To evaluate the accuracy of the coarse registration (patch based pairwise registration) and fine registration (multi-station adjustment), we used the TLS data collected in March 2013 as an example. The errors of the coarse registration between the TLS stations are shown in Figure 7 (left columns). Based on the initial values provided by the pairwise registration, the errors of the multi-station adjustment are shown in Figure 7 (right columns). Figure 8 shows the error distribution of the multi-station adjustment between UAV dense image points and TLS stations.

Table 2. Mean registration errors between UAV dense image points and TLS stations.

Stations	Pairwise Registration (m)	Multi-Station Adjustment (m)	Number of Overlap Points
T1 and T2	0.037	0.013	2,443,652
T2 and T3	0.075	0.024	598,762
T3 and T4	0.067	0.028	517,641
T4 and T5	0.035	0.017	874,186
T5 and T6	0.052	0.016	1,193,318
UAV and TLS	0.054	0.023	8,092,421

shows the mean errors of the pairwise registration and multi-station adjustment. The mean errors were calculated by the mean distance of the nearest points in the overlap areas.

The left columns in Figure 7 suggest that the patch based pairwise registration can realize the coarse alignments between TLS stations. Table 2 shows that the mean errors of pairwise registrations can be controlled within 0.10 m. It illustrates the good robustness of this patch based registration method. The right columns in Figure 7 illustrate that the proposed multi-station adjustment can improve the registration accuracy effectively. The mean errors after multi-station adjustment can be controlled within 0.03 m (as shown in Table 2). This suggests that the initial values provided by pairwise registration are valid, thus preventing the multi-station adjustment from interfering with the local optimum. It also illustrates the effectiveness of the proposed multi-station adjustment method, which obtained the corrections of the registration parameters of each station by making full use of all locally closed loops, thus realizing a high accurate alignment between the TLS points and UAV dense image points. Moreover, Figure 7 and Table 2 show that the registration errors between station 2 and station 3 are relatively high (about 0.075 m). This is because extensive vegetation exists in the overlapping regions, leading to differences between the TLS stations from different views. However, the proposed method can still realize the registration between the stations, demonstrating a reliable and stable solution for the accurate registration of complex scenes.

3.1.3. Registration Results of Multi-Temporal TLS Surveys

According to Section 2.4 and Section 2.6, to reduce the GCPs and ensure the robustness of the registration, patch primitives were used to realize the coarse registration of TLS stations of the same period. Then, the UAV dense image points from September 2014 were taken as the reference framework to transform multi-temporal TLS point clouds into the control points’ coordinate system (WGS-84) by the proposed multi-station adjustment method. Finally, the mean registration errors between the UAV dense image points and multi-temporal TLS surveys were calculated to check the transformation accuracies. The registration results of the multi-station adjustment for each TLS survey are shown in Figure 9. The transformation results of the multi-temporal TLS surveys are shown in Figure 10. The mean errors between the UAV dense image points and five TLS surveys after multi-station adjustment are shown in

Table 3. Mean errors between UAV points and five TLS surveys after multi-station adjustment.

Five TLS Surveys	TLS of March 2013	TLS of August 2013	TLS of November 2013	TLS of September 2014	TLS of January 2015
Mean errors (m)	0.023	0.043	0.016	0.037	0.013
Number of overlap points	8,092,421	5,507,327	11,334,672	4,710,162	10,263,734

, respectively.

The statistical results in Table 3 demonstrate that the mean errors in the overlapping regions between five TLS surveys and the UAV dense image points are 0.023 m, 0.043 m, 0.016 m, 0.037 m, and 0.013 m, respectively, indicating that the registration between the TLS point clouds and UAV image dense points was achieved. In particular, visual checking found that the registrations between the cliff areas were also achieved. Such good registration results lay a solid foundation for the generation of high quality DEMs.

3.2. Experiments of Deformation Evaluation

DEMs of Each Survey and Deformation Maps between Multi-Temporal DEMs

Figure 11 shows the multi-temporal DEMs generated with a 0.2 m horizontal resolution of the study coastal area. To analyze the evolution of the terrain morphology, the terrain surface deformations over five time periods from March of 2013 to January of 2015 were calculated and are illustrated in Figure 12. The two adjacent temporal DEMs were also compared to investigate the extent of the deformation and are illustrated in Figure 13.

As seen in Figure 11, the multi-temporal DEMs accurately reflect the elevations changes of each period time and describe the terrain relief and topography in detail. The terrain surface deformations for over half of the past two years were calculated in the whole study area from the five surveys, resulting in 14,024 m³ of total earthwork. Figure 12 shows that the terrain surfaces change obviously each year. The deformations are mainly located at the two sides of the mountain slopes, because these areas undergo large excavations every year, resulting in large terrain volume changes as illustrated in Figure 13a,c. Besides the influence of excavation activities, heavy rainfalls during the summer of 2013 also caused large topographic changes (Figure 13a).

4. Discussion

4.1. Slope and Aspect Maps of Multi-Temporal TLS Surveys

Slope and aspect maps for each survey were also made for the terrain morphological analysis. Both were derived from the above DEMs. The slopes were classified into nine classes ranging from 0° to 82°, with slopes greater than 82° merged into one class. The aspects were classified into 10 classes corresponding to the cardinal directions. Figure 14 illustrates the slope and aspect maps, respectively.

As shown in Figure 14, areas with slopes over 36° account for about 80% of the total study area. The aspect maps show that the largest areas are located in the northeast direction, contributing about 30% of the total area. Combining the slope and aspect maps, we found that the terrain relief is mostly along the northeast direction, which is also the main direction of the topographic changes. Besides, Figure 14 demonstrates that the slope and aspect maps also changed with the excavations or landslides. To demonstrate this clearly, we selected three areas (as the black circles show) from the areas of interest. Based on the slope and aspect information of the areas, the deformations of the slope and aspect are consistent with the topographic deformations shown in Section 3.2.

To demonstrate the topographic changes of the areas of interest further by combining the slope and aspect information, the aspect-slope scatter points’ maps of the first survey (March 2013) and last survey (January 2015) based on the slope and aspect maps of the three areas were provided. Figure 15 shows the aspect–slope scatter plot of these three areas.

Figure 15a shows that area A has a wide range of slope (from 10 to 80 degrees), and two main directions (northeast and northwest, 40 and 320 degrees). After the excavation or landslides, the aspect changed obviously. The direction of northwest (320 degrees) was changed to east (40 degrees). From Figure 15b, we can see that area B has a small change in the slope (from a range of 35–60 degrees to the range of 35–70 degrees), but little in the aspect. Figure 15c shows that both ranges of the slope and aspect are wide (0–65 degrees for slope, 0–360 degrees for aspect). It also suggests that area C is an isolated hill. From Figure 15c, it is shown that a small area has changed in both the slope and aspect (the top blue points disappeared).

4.2. Vertical Sections of Marker Pegs

The profiles of the filtered data from five surveys were further used for the terrain surface analysis of areas of interest. Five marker pegs were set at the areas of interest. These pegs were used to obtain vertical sections of each period of point clouds. The distribution of the marker pegs in March 2013 DEM is illustrated in Figure 16. For each peg position, the vertical section was interpolated along its aspect direction (as shown in Figure 17, each color represents each survey data).

Figure 17 illustrates that all the marker peg positions have large slopes, and the deformation areas are mainly located between 30 m and 70 m. The top portions of the vertical sections from the No.2 and No.3 marker pegs are missing. This is because dense vegetation exists there, and the filtered results do not affect the terrain analysis of the areas of interest. The green boxes in Figure 17a,b demonstrate that the No.1 and No.2 peg positions experienced obvious “erosion” that was mainly caused by excavation activities.

The red boxes in Figure 17a,e illustrate that the partial later surfaces are higher than the earlier surfaces, which suggests that serious landslides occurred there. Additionally, the red boxes in Figure 17b,d show that minor landslides occurred at the No.2 and No.4 marker peg positions. It is noteworthy that the portions above the four red boxes are very steep (slopes over 55°). Furthermore, comparison of the data shows that most of the landslides occurred between March and August 2013 (rainy season). Figure 17 demonstrates that the topographical changes of the marker pegs are consistent with the results in Figure 12 and Figure 13, demonstrating the accuracy and reliability of the proposed method.

5. Conclusions and Future Work

DEM analysis based on TLS datasets is a common method for deformation detection. To transform multi-temporal TLS surveys into the control points’ coordinate system, most existing methods depend on a variety of GCPs. However, the deployment of various GCPs is time-consuming and labor-intensive, and cumulative errors are not fully considered. Compared with the cumbersome deployment of GCPs, a mini-UAV with a consumer grade optical camera provides a more beneficial method. In this paper, we investigated a deformation method based on a new registration algorithm that accurately registers TLS scans to UAV dense image points with limited GCPs. To fulfill the task, we pairwise registered the TLS scans by extracting patch primitives, and introduced a locally closed loops based multi-station adjustment algorithm for the accurate alignment of TLS and UAV data. The captured data from TLS and UAV were then registered in a uniform spatial reference framework. Second, the deformations of the experimental area were detected and analyzed according to the results generated by the five period surveys. The statistical results demonstrated that the deformation areas and the terrain volumes changes were effectively located, showing the success and effectiveness of the proposed method. Compared with the control survey (error within 0.03 m), this method had less precision (considering the errors of dense image points and the registration errors, the error was about 0.10 m), but was adequate for the application of this paper.

Some issues are still worthy of attention. Such applications have only been focused on specific site locations, rather than allowinfgfor extensive use in a variety of study areas. This is due to the cumbersome nature of the TLS instruments, making it difficult to deploy them in harsh mountain environments.