Landslide Deformation Extraction from Terrestrial Laser Scanning Data with Weighted Least Squares Regularization Iteration Solution

Abstract

The extraction of landslide deformation using terrestrial laser scanning (TLS) has many important applications. The landslide deformation can be extracted based on a digital terrain model (DTM). However, such methods usually suffer from the ill-posed problem of a multiplicative error model as illustrated in previous studies. Moreover, the edge drift of commonly used spherical targets for point cloud registration (PCR) is ignored in the existing method, which will result in the unstable precision of the PCR. In response to these problems, we propose a method for extracting landslide deformations from TLS data. To archive the PCR of different period point clouds, a new triangular pyramid target is designed to eliminate the edge drift. If a fixed target is inconvenient, we also propose a PCR method based on total station orientation. Then, the use of the Tikhonov regularization method to derive the weighted least squares regularization solution is presented. Finally, the landslide deformation is extracted by DTM deference. The experiments are conducted on two datasets with more than 1.5 billion points. The first dataset takes Lashagou NO. 3 landslide in Gansu Province, China, as the research object; the point cloud data were collected on 26 February 2021 and 3 May 2021. The registration accuracy was 0.003 m based on the permanent triangular pyramid target and 0.005 m based on the total station orientation. The landslide deforms within 3 cm due to the ablation of the frozen soil. The second dataset is TLS data from the Lihua landslide in Chongqing, China, collected on 20 April 2021 and 1 May 2021. The overall deformation of the Lihua landslide is small, with a maximum value of 0.011 m. The result shows that the proposed method achieves a better performance than previous sphere-based registration and that the weighted least square regularization iterative solution can effectively reduce the ill-condition of the model.

1. Introduction

The extraction of landslide deformation is important for many applications [1,2,3,4], such as disaster management and deformation detection. Laser scanning can be used to acquire accurate and dense 3D points from a target surface and has unique advantages when it comes to landslide monitoring, deformation extraction, and disaster management.

As the most direct manifestation of landslide stability, landslide surface deformation has received special attention from scholars [5,6,7]. Conventional landslide monitoring methods, such as using the GNSS (Global Navigation Satellite System) [8,9,10], leveling [11,12], crack gauges [13,14], etc., can only obtain sparse measurements at a few locations or within the affected slopes. Therefore, it is difficult to interpret the overall deformation characteristics of a landslide [15,16]. Terrestrial laser scanning (TLS) is a ground-based active imaging method that rapidly acquires precise and dense 3D point clouds on the surface of objects through laser ranging. Guo et al. [17] used point cloud data and the digital elevation model visualization method capable of a sky view to carry out geological disaster identification research and verified the reliability of airborne LiDAR identification results through field investigation. Abellan et al. [18] used TLS technology to monitor a dangerous rock mass in Spain, and discussed the feasibility of millimeter-level high-precision monitoring. Kayen et al. [19] used TLS to monitor nearly 400 large landslides induced by the Chuetsu earthquake in Niigata Prefecture, Japan, and greatly improved the efficiency of post-earthquake disaster assessment. Liu et al. [20] proposed a landslide displacement monitoring method based on point cloud density characteristics; the method identified the slope variation area and directly reflected the landslide surface deformation.

PCR can merge these individually scanned period point clouds. The basic idea of PCR is to seek the best transformation parameters to transform a point cloud with a local coordinate system to the same reference system [21]. PCR is divided into real-time registration and accurate registration of point clouds in different periods. Real-time registration refers to the detection of the surrounding environment while scanning and registering the real-time scanning point cloud using mobile laser scanning with the acquired point cloud. It is widely used for fast modeling, indoor navigation, and simultaneous localization and mapping (SLAM). It is a low-precision PCR method. The point cloud data of different periods is typically obtained using fixed scanners (such as TLS). The main PCR methods for different period data are marker-based registration and data-based registration.

However, there are still many problems with the application of landslide monitoring using TLS technology [22,23]. When using TLS to monitor landslide deformation, it is necessary to collect a multi-period landslide point cloud and calculate the landslide deformation by comparing the spatial position of the point cloud. Unfortunately, each period point cloud is based on an independent coordinate system, and the reference data are not unified, resulting in the low accuracy and unreliable results of landslide deformation detection [24].

The widely used marker in TLS is the standard spherical target [25,26]. Spherical targets have high precision and are easy to carry, but they are expensive, require high precision in manufacturing and maintenance, and cause edge drift because of their excessive laser reflection angle. The spherical target is shown in Figure 1, and the spherical target cloud is shown in Figure 2. Since the edge drift is concentrated away from the scanning side, it does not receive uniformly distributed random noise, which will cause non-negligible calculation errors in the spherical center coordinates. An ICP algorithm is the most used data-based PCR method due to its simple iteration and fast convergence [27]. The principle of ICP is to find the nearest correspondence in the source and the target point cloud and use the distance matrix of the nearest correspondence to estimate the rotation and translation parameters [28].

The landslide deformation can be extracted based on DTM [29,30]. Unfortunately, the DTM adjustment model is a typical ill-posed multiplicative random error model [31]. The ill-posed error model widely exists in the processing of remote sensing causing the instability of the parameter solution and even serious deviations from the true value. Xu et al. [32] first proposed the least square solution for the multiplicative error model. Shi et al. [33] summarized the least square, weighted the least square, and bias-corrected the weighted least square method of the multiplicative error model, and deduced the accuracy formula for the parameter estimation. If these three methods are used to deal with the ill-posed multiplicative error model without considering the ill-posed feature of the coefficient matrix, the parameter estimation will be biased and unstable. To avoid the ill-posed coefficient matrix, Wang et al. [34] built a DTM using the points at the peaks in the simulated data. This model can only be used to process simulation data and cannot be applied to practical observational data. Due to the ill-posed coefficient matrix, the complex collinearity between the column vectors of the coefficient matrix causes the condition number of the normal equation to be too large and thus the solution is unstable. Therefore, research on ill-posed multiplicative random error model adjustment theory needs to be further studied. At present, the methods for dealing with ill-posed multiplicative error models mainly include truncated singular-value decomposition (TSVD) [35], the Tikhonov regularization method [36], the ridge estimation method [37], and the virtual observation value solution [38]. The Tikhonov regularization method is a general adjustment method that strictly adheres to the theory. It is commonly used to solve the ill-posed multiplicative random error model.

The works mentioned above have carried out detailed studies on the landslide deformation extraction from TLS. The PCR of different period point cloud and the DTM construction are the vital procedures for deformation extraction. However, the number of point clouds on a marker sphere is sparse in the scanning process of a landslide. Edge drift caused by an excessively large laser reflection angle inevitably occurs due to the characteristics of the sphere. Thus, the calculation error of the spherical center coordinates cannot be disregarded. In this research, a new triangular pyramid target was developed which can unify the data of the multi-period point cloud by three non-parallel surfaces. In addition, aiming to address the problems of poor site access conditions and inconvenient traffic, we also propose a new scheme of data acquisition and registration based on total station orientation. The DTM adjustment model is an ill-posed multiplicative error model. If an ill-posed coefficient matrix is not considered, the parameter estimation will be biased and unstable. Hence, we utilize the Tikhonov regularization method to derive the weighted least squares regularization solution.

2. Methodology

The method proposed in this paper mainly consists of three steps. First, the PCR of the multi-period point cloud data is performed. Then, considering the ill-multiplicative error model, a weighted least square regularization iterative solution is proposed for DTM construction. Finally, the landslide deformation is obtained by differential DTMs.

2.1. PCR of Multi-Period Point Cloud Data

2.1.1. PCR Based on Triangular Pyramid Target

Multi-phase point cloud is characterized by three fixed triangular pyramid vertices of concrete in the field. The source point set $X\triangleq {\left\{{\overrightarrow{x}}_i\right\}}_{i=1}^{N_x}(N_x\in \mathbb{N})$ and the target point set $Y\triangleq {\left\{{\overrightarrow{y}}_i\right\}}_{i=1}^{N_y}(N_y\in \mathbb{N})$ are assumed in the n-dimensional space ${ℝ}^n$. The objective of the rigid-body registration algorithm is to find a rotation matrix and translational matrix so that the source point set and target point set can optimally correspond in space. The optimization criterion based on Euclidean distance can be written as follows [39],

(1)$\underset{R,\overrightarrow{t},c(i)\in \left\{1,2,\cdots ,N_y\right\}}{{\displaystyle \min }}{\displaystyle \sum _{i=1}^{N_x}{‖(R{\overrightarrow{x}}_i+\overrightarrow{t})-{\overrightarrow{y}}_{c(i)}‖}_2^2} \mathrm{s}.\mathrm{t}. R^{T}R=I_3,\det (R)=1$

ICP algorithm is widely used in rigid body registration due to its simple iteration and fast convergence. The ICP algorithm iteratively computes the rigid transformation parameters R and $\overrightarrow{t}$ until the objective function converges to a local minimum. In the kth iteration, the algorithm is implemented through the following two steps,

• (1). Setting up the matching point correspondence of two point sets,

  (2)$c_{k(i)}=\arg \min {‖(R_{k-1}{\overrightarrow{x}}_i+{\overrightarrow{t}}_{k-1})-{\overrightarrow{y}}_{c(i)}‖}_2^2$

  where $\mathrm{i}=1,\cdots ,N_x$, $R_{k-1}$, ${\overrightarrow{t}}_{k-1}$ is the argument of iteration k − 1.

• (2). The new rotation matrices and shift vectors are computed by minimizing the square distance [40],

  (3)$\left(R_k,{\overrightarrow{t}}_k\right)=\underset{R^{T}R=I_n,\det (R)=1,\overrightarrow{t}}{{\displaystyle \arg \min }}{\displaystyle \sum _{i=1}^{N_x}{‖R{\overrightarrow{x}}_i+\overrightarrow{t}-{\overrightarrow{y}}_{c_k(i)}‖}_2^2}$

The condition for the end of iteration for the ICP algorithm is that the registration error is small enough or the number of iterations reaches its maximum number, then the algorithm terminates, otherwise the iteration continues [41]. We use root mean square (RMS) $\epsilon$ to evaluate the registration accuracy,

(4)$\epsilon =\frac{1}{N}{{\displaystyle \sum _{i=1}^{N_s}‖R_k{\overrightarrow{x}}_i+{\overrightarrow{t}}_k-{\overrightarrow{y}}_{c_k(i)}‖}}_2^2$

where $R_k$, ${\overrightarrow{t}}_k$ is the argument of iteration k.

2.1.2. PCR Based on Total Station Orientation

The point cloud coordinates are converted to the independent coordinate system using the Bursa–Wolf transformation model. The Bursa–Wolf transformations are conformal 3D Cartesian coordinate transformations commonly used in surveying, photogrammetry, and geodesy [42].

2.2. Landslide DTM Using the Weighted Least Squares Regularization Solution

Unlike common Gauss–Markov models with additive random error [43,44,45], we suppose a set of observations with multiplicative random error $y_i\left(i=1,2,\cdots n\right)$, and the multiplicative error model can be defined as follows [46],

(5)$y_i=f_i(\beta )(1+{\epsilon }_i),i=1,2,\cdots ,n$

where ${\epsilon }_i$ is the zero mean random error, $f(\beta )$ and $\epsilon$ are the random variables, $\beta$ is the t-dimensional unknown parameter vector to be estimated, $\odot$ is the Hadamard product of matrices or vectors, and 1 is the n dimensional column vectors with all values set to 1.

Assuming that the variances of the elements in vector $\epsilon$ are equal and independent, if every $f_i(·)$ is a linear function of variable $\beta$, that is, $f_i(\beta )=x_i^T\beta$, $x_i^T$ is a t-dimensional row vector, the multiplicative error model can be written as follows,

(6)$y=(X\beta )\odot (1+\epsilon )$

The Formula (6) can be rewritten in the following form,

(7)$\overline{e}=y-X\beta =(X\beta )\odot {\epsilon }_m$

where $X={(x_1,x_2,\cdots ,x_n)}^T$.

It is evident that the accuracy of the observation value in the multiplicative error model is proportional to the signal intensity. The stronger the signal $X\beta$, the greater the error. The covariance matrix $Q_y$ is a function of the unknown parameter $\beta$. The commonly used methods for solving the multiplicative error model include the least squares method and the deviation correction weight least squares method. The weighted least square regularization algorithm is applied to the ill-posed multiplicative error model in this paper.

We utilize the least square (LS) method to calculate the parameter; the parameter estimation can be expressed as follows [47],

(8)${\widehat{\beta }}_{LS}=(X^{T}X{)}^{-1}{X}^{T}y$

where ${\widehat{\beta }}_{LS}$ is the least square estimation of the multiplicative error model.

According to the law of cofactor propagation, the cofactor matrix of $\overline{e}$ and $y$ is,

(9)$Q_{\overline{e}}=Q_y=\left[\begin{array}{cccc}{x}_1^T\beta & 0& 0& 0\\ 0& x_2^T\beta & 0& 0\\ 0& 0& \ddots & 0\\ 0& 0& 0& x_n^T\beta \end{array}\right]$

We also utilize an unbiased estimation of weighted least-square (bcWLS) method; the parameter estimation can be expressed as follows [48],

(10)${\widehat{\beta }}_{bc}{=(\mathrm{X}}^T{\mathrm{Q}}_y^{-1}{\mathrm{X})}^{-1}{\mathrm{X}}^T{\mathrm{Q}}_y^{-1}y$

where ${\widehat{\beta }}_{bc}$ is the parameter estimation using bcWLS method.

If the normal equation $N={(X^{T}X)}^{-1}$ is ill-posed, the LS solution is unreliable. Hence, to overcome the ill-posed state of the normal equation, we use the weighted least squares regularization iterative solution (RWLS). First, a regularization factor $\alpha$ is introduced to construct the regularization criterion of the ill-posed multiplicative error model,

(11)$\varphi ={\overline{e}}^T{Q}_{\overline{e}}^{-1}\overline{e}+\alpha {\beta }^T\beta =\min$

According to the definition of RWLS, the parameter estimation can be expressed as follows [49],

(12)${\widehat{\beta }}_{RWLS}={(A^T{Q}_{\overline{e}}^{-1}A+\alpha I_n)}^{-1}{A}^T{Q}_{\overline{e}}^{-1}y$

where ${\widehat{\beta }}_{RWLS}$ is the RWLS solution of ill-posed multiplicative error model, and $I_n$ is identity matrix. The term $\alpha I_n$ is added to the inversion of the normal equation, which effectively weakens the ill-posed state of the normal equation. Therefore, we can obtain a reliable and stable solution.

As evident from Equation (11), $‖y-A\widehat{\beta }‖$ and $‖\widehat{\beta }‖$ are the function of the regularization parameter $\alpha$, and the L curve method is used to calculate the regularization parameters $\alpha$. The main procedures are as follows [50,51],

• (1). Taking $‖y-A\widehat{\beta }‖$ as the x-axis coordinate and $‖\widehat{\beta }‖$ as the y-axis coordinate, we get multiple sets of coordinate points ($‖y-A\widehat{\beta }‖$, $‖\widehat{\beta }‖$),

• (2). These coordinate points are fitted to a curve similar to L shape, and the $\alpha$ value corresponding to the point with the largest curvature is used as the estimation of the regularization parameter.

As evident from Equation (12), $Q_{\overline{e}}$ is the nonlinear function of ${\widehat{\beta }}_{RWLS}$ on the left side of the Equation (12), and ${\widehat{\beta }}_{RWLS}$ is the nonlinear function of $Q_{\overline{e}}$ on the right side of Equation (12). Hence, there is no analytical solution, only a numerical solution [52]. We can only use approximation methods to iteratively calculate the numerical solution.

The detailed procedure of RWLS solution is as follows,

• (1). The LS solution is used as the RWLS initial value of iteration;

• (2). Compute the cofactor matrix $Q_{{\overline{e}}_i}$;

  (13)$Q_{{\overline{e}}_i}=\left[\begin{array}{cccc}{x}_1^T{\widehat{\beta }}_{RWLS}^i& 0& 0& 0\\ 0& x_2^T{\widehat{\beta }}_{RWLS}^i& 0& 0\\ 0& 0& \ddots & 0\\ 0& 0& 0& x_n^T{\widehat{\beta }}_{RWLS}^i\end{array}\right]$

• (3). The regularization parameter ${\alpha }_i$ is obtained by L curve method;

• (4). The iterative formula is as follows:

  (14)${\widehat{\beta }}_{RWLS}^{i+1}={(A^T{Q}_{{\overline{e}}_i}^{-1}A+{\alpha }_i{I}_n)}^{-1}{A}^T{Q}_{{\overline{e}}_i}^{-1}y,i=0,1,\cdots$

• (5). Repeat steps (2)–(4) until the value of the parameter subtraction is smaller than the threshold,

  (15)$‖{\widehat{\beta }}_{RWLS}^{i+1}-{\widehat{\beta }}_{RWLS}^i‖<{\epsilon }_{th}$

  where ${\epsilon }_{th}$ is defined as 10⁻¹⁰ in this paper.

The unit weight variance is always an important quality index for evaluating observation results and data processing results. It is defined as follows by using the LS method,

(16)${\widehat{\sigma }}_{LS}^2=\frac{V_{LS}^T{V}_{LS}}{tr\{{\widehat{Q}}_y-{(X^{T}X)}^{-1}{X}^T{\widehat{Q}}_{y}X\}}$

where ${\widehat{\sigma }}_{LS}^2$ = unit weight variance.

The unit weight variance using bcWLS can be expressed as follows [53],

(17)${\widehat{\sigma }}_{bc}^2=\frac{V_{bc}^T{\widehat{Q}}_y^{-1}{V}_{bc}}{r}$

where ${\widehat{\sigma }}_{bc}^2$ = unit weight variance, r = freedom degree, and its value correlates with the number of observation value and parameter estimation.

The unit weight variance of the ill-posed multiplicative random error model can be expressed as follows [54],

(18)${\widehat{\sigma }}_0^2=\frac{{\widehat{\overline{e}}}^T{Q}_{\overline{e}}^{-1}\widehat{\overline{e}}}{r}$

where $\widehat{\overline{e}}=A{\widehat{\beta }}_{RWLS}-y=[A{(A^T{Q}_{\overline{e}}^{-1}A+\alpha I_n)}^{-1}{A}^T{Q}_{\overline{e}}^{-1}-I]y$.

The unit weight variance using RWLS can be expressed as follows [52],

(19)${\widehat{\sigma }}_{RWLS}^2={\widehat{\sigma }}_0^{2}tr(Q_{{\widehat{\beta }}_{RWLS}\times {\widehat{\beta }}_{RWLS}})+tr(b_{{\widehat{\beta }}_{RWLS}}{b}_{{\widehat{\beta }}_{RWLS}}^T)$

where ${\widehat{\sigma }}_{RWLS}$ is the unit weight variance. In Equation (19), $b_{{\widehat{\beta }}_{RWLS}}$ is defined as follows,

(20)$b_{{\widehat{\beta }}_{RWLS}}=-\alpha {(A^T{Q}_{\overline{e}}^{-1}A+\alpha I_n)}^{-1}{\widehat{\beta }}_{RWLS}$

In this paper, we use LS, bcWLS, and the RWLS method to estimate the parameters and calculate the accuracy of parameters.

In addition, the noise of point cloud data has been shown with a nature of multiplicative random error [55,56,57]. Therefore, the multiplicative error model of DTM with 6 unknown parameters is as follows [33],

(21)$H({\mathrm{X}}_{\mathrm{i}},y_i)=F({\mathrm{X}}_{\mathrm{i}},y_i)\odot (1+{\epsilon }_m)$

where $({\mathrm{X}}_{\mathrm{i}},y_i)$ represents the x and y coordinates of the ground point, $H({\mathrm{X}}_{\mathrm{i}},y_i)$ represents the corresponding elevation of the point, and ${\epsilon }_m$ represents the multiplicative error vector,

(22)$F({\mathrm{X}}_{\mathrm{i}},y_i)={\beta }_1+{\beta }_2{x}_i+{\beta }_3{y}_i+{\beta }_4{x}_i{y}_i+{\beta }_5{x}_i^2+{\beta }_6{y}_i^2$

Then, we apply LS, bcWLS, and RWLS method to calculate the parameter estimation and unit weight variance. Finally, we extract the landslide deformation by the difference of DTM.

3. Experiments Data

3.1. Simulation Experiment

3.1.1. Simulation Experiment I: GNSS Elevation Point Disturbed by Multiplicative Error

GNSS is used to measure the elevation of the ground point of a road center line in a certain area. It is assumed that the ground point elevation and distance conform to the following functional model,

(23)$y=10+4x+2x^2+x^3+0.5x^4+2x^5$

In this example, the value range of x is 0–300 m, and 31 points are equally spaced within the value range of x. The truth values of the parameters are 10, 4, 2, 1, 0.5, and 2.

Assuming that y is disturbed by the multiplicative random error, where the multiplicative random errors are independent of each other and follow a normal distribution with a mean value of 0 and a standard deviation of 0.1, and the corresponding observation equation of the multiplicative error model is,

(24)$Y(x)=y(x)\odot (1+{\epsilon }_m)$

where Y is the vector disturbed by multiplicative error, 1 is the 31-dimensional column vector with all elements being 1, and ${\epsilon }_m$ is the 31-dimensional multiplicative random error vector. The simulation values of the ground elevation points obtained from Equation (23) are shown in Table 1.

To show that the elevation points are affected by multiplicative error, the data of the ground elevation points before and after interference from multiplicative error are plotted in Figure 3.

3.1.2. Simulation Experiment II

This simulation uses DTM data. This paper mainly considers that the model is ill-posed, so it simulates the ill-conditioned DTM model and generates the DTM using the interpolation function method [53]

(25)${\overline{f}}_{xy}={\displaystyle \sum _{i=1}^4{\beta }_i{f}_i(x,y)}$

where ${\overline{f}}_{xy}$ represents the elevation obtained by the interpolation function, x and y range from 0 m to 80 m, the interval is 2 m, and the truth values of the four parameters are −1.5, 20, 5, and −4, respectively.

The function $f_i(x,y)$ is as follows,

(26)$f_1(x,y)=\exp \{-({(x-22)}^2+{(y-22)}^2)/500\}$

(27)$f_2(x,y)=\exp \{-({(x-28)}^2+{(y-28)}^2)/500\}$

(28)$f_3(x,y)=\exp \{-({(x-25)}^2+{(y-25)}^2)/500\}$

(29)$f_4(x,y)=\exp \{-({(x-20)}^2+{(y-20)}^2)/500\}$

The multiplicative error model is as follows,

(30)$h(x,y)=\overline{f}(x,y)\odot (1+{\epsilon }_m)$

where $h(x,y)$ is the vector of the observation values disturbed by multiplicative errors, 1 is the 1681-dimensional column vector whose elements are all 1, and ${\epsilon }_m$ is the 1681-dimensional multiplicative error vector values which are independent of each other. To illustrate that the DTM model is affected by multiplicative error, the DTM without the errors are plotted in Figure 4 and the DTM influenced by multiplicative random error is plotted in Figure 5. As can be seen from Figure 4 and Figure 5, although the DTM with the standard deviation has a 0.1 multiplicative error, it has a great influence on the elevation, and the condition number of the normal equation $N=A^T{Q}_{\overline{e}}^{-1}A$ is 4.456 × 10⁴; therefore, the equation is seriously ill.

3.2. Actual Experiment

3.2.1. Actual Experiment I

• (I). Experiment 1

The Lashagou NO. 3 landslide was a loess landslide located at Lashagou village, Jishishan county, Gansu province, China. Its geographical coordinates are 35° 35′ 23″E, 102° 57′ 12″N. It is in the transition zone between the Tibetan Plateau and the Loess Plateau, which is typical of an excavation-related shallow loess engineering landslide [58,59]. The main stratum of the landslide is Holocene residual slope silty clay and the main stratum of the slide bed is middle and lower Pleistocene fluvial and lacustrine maroon silty clay. The central width of the slide body is 60–100 m, the axial length is 80–130 m, the area is about 25,000 m², the thickness is 3–8 m, and the volume is about 150,000 cubic meters, belonging to the middle-shallow medium and small-sized landslide groups. The rear wall of the landslide was 0.2–0.4 m high, the natural slope was 15°, and the main slide direction was 85°. The front slope of this landslide group is mainly damaged by slime in the surface 1–3 m, but the lower soil is relatively stable, and the back edge of the landslide group is in the shape of an armchair and a gentle arc, which belongs to the traction cohesive soil landslide group. The optical remote sensing imagery of the landslide is shown in Figure 6a. As shown in Figure 6b, the toe of the Lashagou NO. 3 landslide is close to the G310 highway, which is activated or generated by the slope excavation during the construction of the G310 mountain highway. The main effects of the Lashagou NO. 3 landslide are shown in Figure 7. Figure 7b shows a local diagram of a crack 10 cm wide. Figure 7c shows the local displacement diagram of a 0.6 m staggered platform in the vertical direction of road surface formed by landslide. Figure 7d shows the local map of the road cracks formed by the landslide.

In this study, a Leica P50 laser scanner was used for data acquisition. The scanner is a ground-based 3D laser scanner with medium and long range, which has the advantages of obtaining high-quality point cloud data in harsh environments and a wide scanning range with minimal noise. The parameters of the Leica P50 laser scanner are shown in Table 2. In addition, A Leica Nova TS30 total station and some prism spheres are also required to collect data. The parameters of the Leica TS30 total station are shown in Table 3.

• (I). Triangular pyramid target point cloud acquisition

In this paper, a triangular pyramid target was developed, and three concrete triangular pyramids were poured in the study area, among which any triangular pyramid target is shown in Figure 8. The size of the triangular pyramid template is an equilateral triangle with a side length of 50 cm and a height of 15 cm. The point cloud of the triangular pyramid is shown in Figure 9.

• (II). Data collection based on total station orientation

In the monitoring of the landslide and other deformation forms, points in the stable zone should be taken as the monitoring benchmark. The coordinates of the four control points in the stable zone in this paper are shown in Table 4.

The point cloud of the prism sphere was collected using a Leica P50 scanner. We used the Leica Nova TS30 total station to collect the 3-D coordinates of the prism sphere simultaneously. The experiment was divided into three steps,

• (1). Leica TS30 total station was used to collect the coordinates of the prism balls near the landslide area.

The total station was set up at any position, and the coordinates of any two control points in the stability area were intersected to obtain the coordinates of the current measuring station. Then the central coordinates of the prism balls near the landslide area were measured. The photos of the prism balls are shown in Figure 10.

• (2). The coordinates of the prism balls near the landslide area were obtained by fitting point cloud data collected by TLS simultaneously.

The Leica P50 scanner was set up at any position, and the point cloud data of the prism balls were obtained. Then the center coordinates of the prism balls were calculated by fitting the point cloud data of the prism balls. The transformation parameters from the point cloud coordinate system to the total station coordinate system can be obtained by using the Bursa–Wolf model.

• (3). The unification of the data.

The point cloud coordinates of each period were converted to the control point coordinate system of the stable area by using the transformation parameters to realize the unification of the data.

The creep deformation of a landslide is a slow dynamic process. In this paper, a one-station measurement method was adopted to avoid the error caused by multi-station data PCR, which improved the accuracy of the landslide deformation monitoring. The Leica P50 scanner was set on the highest point of the slope opposite the landslide, and the scanning accuracy was set to 3 mm@ 10 m. Considering the influence of permafrost ablation and other factors on the landslide, two period point cloud data were collected on 26 February 2021 (temperature 3 degrees Celsius, before permafrost ablation) and 3 May 2021 (temperature 28 degrees Celsius, after permafrost ablation) to analyze the landslide deformation. The data collection photo from period I (May 3, 2021) is shown in Figure 11, and the point cloud of Lashagou NO.3 landslide is shown in Figure 12.

3.2.2. Actual Experiment II

The Lihua landslide was on the left side of the Lihua road at Nan’an District, Chongqing, China. The relative height difference was about 10 m. It was a typical excavation-induced shallow loess engineering landslide. Once a large-scale instability occurs, it may block the Lihua road and threaten the lives and property of the nearby residents and pedestrians. We collected a two-phase point cloud of the landslide with a Leica P50 3D laser scanner on 20 April 2021 (phase 1) and 1 May 2021 (phase 2). The Lihua landslide point cloud of phase 1 is shown in Figure 13. The UAV image of the Lihua landslide acquired on 4 June 2022 is shown in Figure 14.

4. Analysis of Result and Discussion

4.1. The Results of Simulation Experiments

4.1.1. The Result of Simulation Experiment I

In simulation experiment 1, the square root of unit weight variance is 0.3. The parameter values obtained using the three methods are shown in Table 5.

The comparison between the RWLS method and the classical LS method is shown in Figure 15.

The condition number of the normal equation $N=A^T{Q}_{\overline{e}}^{-1}A$ is 8.3098 × 10⁵, the equation is ill-posed. The regularization parameter values in this algorithm change dynamically, and the change of the regularization parameters with the number of iterations is shown in Figure 16.

It is evident from Figure 15 that ground elevation points are influenced by multiplicative error, resulting in serious deviation. It is evident from Table 5 that the ${‖\Delta \beta ‖}_2$ obtained by LS method is the largest. The ${‖\Delta \beta ‖}_2$ calculated by the bcWLS method is 13.61, which still deviates greatly from the truth value because the ill-posed state of the model is not considered. The ${‖\Delta \beta ‖}_2$ calculated by the RWLS algorithm in this paper is 3.28, and the parameter estimation is closer to the true value, indicating that the RWLS algorithm has a certain effect on reducing ill-condition. It is evident from Table 5 that the bcWLS method does not consider the ill-posed nature of the model, which led to ${\widehat{\sigma }}_0$ deviating from the true value. However, the RWLS algorithm takes these factors into account and the ${\widehat{\sigma }}_0$ is closer to the true value, further verifying the advantages of the algorithm in this paper.

4.1.2. The Result of Simulation Experiment II

Following Xu et al. (2013), it is here assumed that ${\widehat{\sigma }}_0$ is 0.3, and that the multiplicative random error vectors are independent from each other. The parameter estimation, ${‖\Delta \beta ‖}_2$ and ${\widehat{\sigma }}_0$, calculated by LS, bcWLS, and RWLS, are listed in Table 6. The DTM made using the RWLS method is shown in Figure 17. The regularization parameters with the number of iterations are shown in Figure 18.

As can be seen from Table 6, the parameter estimation obtained by bcWLS deviates seriously from the true value because the ill-posed state of the model is not considered, while the result obtained by the RWLS algorithm is closer to the true value. The parameter ${‖\Delta \beta ‖}_2$ calculated by RWLS is 1.88. It is less than that of the LS method and bcWLS method, which further indicates that the algorithm in this paper has a certain effect on reducing the ill-posed state. It is evident from Table 6 that the ${\widehat{\sigma }}_0$ obtained by the LS method is 4.52, and that ${\widehat{\sigma }}_0$ seriously deviates from the true value. The ${\widehat{\sigma }}_0$ calculated by the bcWLS method is 4.36, which is closer to the truth value, because it considers the weight of the observation value. The RWLS method considered the influence of ill-condition on the ${\widehat{\sigma }}_0$, hence the ${\widehat{\sigma }}_0$ is the closest to the true value.

4.2. The Result of the Actual Experiment

• (I). The result of experiment I

4.2.1. The result of PCR Based on Triangular Pyramid Target

Three planes can be fitted using the point cloud of a triangular pyramid, and then triangular pyramid vertices can be intersected through the three planes. The coordinates of three triangular pyramid vertices in phase 1 and phase 2 are shown in Table 7 and Table 8, respectively.

The initial landslide point cloud of phase 1 and phase 2 is shown in Figure 19. The translation vector obtained by the fine algorithm is (0.134, 0.080, −0.038). In this paper, the rotation matrix is transformed into the rotation angle around the coordinate axis, in which the rotation angle around the X axis is 0.164°, the rotation angle around the Y axis is −0.180°, and the rotation angle around the Z axis is −46.886°.

The point cloud after unifying the data through the obtained rotation and translation parameters is shown in Figure 20. The image of the Lashagou NO.3 landslide is shown in Figure 21. It is evident that the two-phase point cloud has been accurately registered based on the triangular pyramid target. The point cloud data of the two periods are unified, which lays a foundation for landslide deformation detection.

The standard sphere target has often been used in previous studies [50,51,52,53]. Figure 1 shows that edge drift caused by an excessively large laser reflection angle will inevitably occur. Thus, the calculation error of the spherical center coordinates cannot be disregarded. Compared with previous studies [52,53], our methods handle this problem by designing a new triangular pyramid marker. As shown in Figure 9, three planes can be fitted through the triangular pyramid point cloud, and then the triangular pyramid vertices can be intersected through the three planes. The registration accuracy of the point cloud based on the triangular pyramid target is 0.003 m.

4.2.2. The result of PCR Based on Total Station Orientation

Figure 22 shows the point cloud after fine registration using the total station orientation theory. The two-phase point cloud has been accurately transformed into the control point coordinate system in the stable area, realizing the unification of the multi-phase point cloud data. The method is rigorous in theory, easy to operate in practice, free from the restriction of landslide terrain conditions, and has a wide applicability. The registration accuracy based on the total station orientation was 0.005 m.

By comparing the two PCR method with other PCR methods, it is evident that:

• (1). In terms of data acquisition, to ensure the accuracy of registration, the data acquisition mode which is based on a fixed triangular pyramid target is simple and is capable of fast operation. It is suitable for good visibility and convenient transportation and is easy to cast in a concrete environment. In addition, based on the total station orientation data acquisition mode, it is not limited by field conditions and has a wide applicability. Especially, in a situation where there are some control points in the stable area, the point cloud can be directly and effectively transformed into the existing coordinate system, which is of great significance to the continuous dynamic monitoring of landslide disasters based on point cloud data.

• (2). In terms of PCR, the PCR theory based on a fixed triangular pyramid target and total station orientation is strict. The PCR based on a fixed triangular pyramid target can overcome the problem of low precision in the sphere-based PCR method. The advantage of PCR based on total station orientation is that the point cloud coordinates can be converted to the existing control point coordinate system, which is favorable for the utilization of survey area engineering.

The proposed approach strives at addressing the unstable precision of different period TLS data. Figure 23 shows the TLS data of the 1 m wide detailed image shown in Figure 20.

As shown in Figure 23, the result after applying the target-based (triangular pyramid) PCR method suggests that the precision is stable. This may also indicate that in TLS data such methods tend to retain point cloud detail features. Compared with the spherical target [25,26], the triangular target can overcome the edge drift. After applying the PCR method based on total station backward orientation, the precision of PCR is high. Compared with a target-based PCR method [27], the station setup is flexible. In addition, it is suitable for the situations that are not easily accessible.

Our PCR method is easy to use. The center of sphere has often been used for PCR in previous studies [25]. Our conclusion about PCR based on spherical target is similar to [27], which is that it is difficult to calculate the spherical center precisely because of edge drift. Thus, a new PCR target is necessary. Compared with previous studies, our method handles this problem by designing a triangular target. The vertices can be precisely intersected by three triangular planes.

4.2.3. Landslide Based on DTM Difference

Parameter estimation values and their accuracy values are shown in Table 9. The coefficient matrix of the LS method and bcWLS method is singular and the parameter estimation does not converge due to the influence of the ill-condition. The algorithm in this paper takes ill-condition into account and has advantages for processing real data.

The DTM of the Lashagou NO.3 landslide is shown in Figure 24. The deformation of the Lashagou NO.3 landslide is shown in Figure 25. It is evident from Figure 25 that from 26 February 2021 to 3 May 2021, the overall deformation of the landslide was small, and its maximum value was 0.031 m. The landslide will continue to maintain a stable state without special factors, such as earthquakes, heavy rainfall, and artificial excavation.

• (II). The result of experiment II

The point cloud of the Lihua landslide (phase I) is shown in Figure 26. The DTM disturbed by multiplicative error obtained in this paper is shown in Figure 27.

The estimation and accuracy of the parameters are shown in Table 10.

The ${\widehat{\sigma }}_0$ obtained by the LS method is too large, which does not conform to reality, further indicating that the LS method cannot process actual data. It can also be seen from Table 10 that the error of the parameter estimation obtained by the LS method is the largest, and that the error obtained by the RWLS method is much smaller than LS and bcWLS method.

The deformation of the Lihua landslide is shown in Figure 28. It is evident from Figure 28 that from 20 April 2021 to 1 May 2021, the overall deformation of the Lihua landslide was small, with a maximum value of 0.011 m.

As shown in Figure 26, the DTM has a huge deviation because of the multiplicative random errors. As indicated in Table 9 and Table 10, the ${\widehat{\sigma }}_0$ obtained by LS and bcWLS deviates from the true value, but RWLS can obtain a relatively reasonable result. The results show that the RWLS method is applicable to solve the ill-posed state of DTM in practical applications.

The ill-posed problem of the coefficient matrix has been ignored in previous studies, which results in an unstable or non-convergent solution [62]. We used the Tikhonov regularization method to derive the RWLS solution for an ill-posed multiplicative error model. The results of simulation I and simulation II suggest that RWLS’s solution has a better performance than other methods, such as LS [47] and bcWLS [48]. In the test on two actual landslide data sets, the LS and bcWLS methods dealt poorly with DTM construction [29,30,54], because the DTM has a huge deviation due to the random errors. Our method takes ill-condition into account. Hence, the results of actual experiment I and actual experiment II are consistent with the results of the simulation data. Thus, the RWLS method provides a new solution for the ill-posed multiplicative error model.

5. Conclusions

In this paper, TLS technology was used to extract the deformation of a landslide. We discussed the major factors affecting the landslide deformation, and attempted to reveal the relevant deformation mechanism. A new measurement and data processing strategy is proposed in this paper, aiming at extracting the deformation of a landslide. Finally, the landslide deformation was extracted by using two datasets. The first dataset is the point cloud data of the Lashagou NO.3 landslide collected on 26 February 2021 and 3 May 2021. The second dataset is the point cloud data of the Lihua landslide collected on 20 April 2021 and 1 May 2021. The main conclusions drawn from this study are as follows.

In this study, we designed a new triangular pyramid PCR marker, and a PCR experiment for different period TLS datasets based on the new marker was carried out. Compared with the commonly used spherical target, the edge drift of the sphere is overcome by the new marker-based PCR proposed in this paper. The registration accuracy based on the permanent triangular pyramid is 0.003 m. If it is not convenient to make a fixed target, we also propose a new PCR method based on total station orientation. The registration accuracy based on the unified theory of backward orientation using the total station method is 0.005 m. The results show that the new method precisely aligns the point cloud system with the engineering independent coordinate system.

The weighted least square regularization iterative method can effectively solve the ill-posed multiplicative error problem. In the tests of two simulation data set, the accuracy and reliability of the RWLS algorithm outperformed previous methods. Then, two actual landslide experiment were conducted. The results show that the RWLS solution can effectively reduce the ill-posed state of the model. The deformation of the Lashagou NO.3 landslide was small, with a maximum value of 0.031 m. The maximum deformation of the Lihua landslide was 0.011 m.

In general, the monitoring results of the DTM subtraction based on TLS data technology show that the Lashagou NO.3 landslide and Lihua landslide will continue to maintain a stable state without special factors, such as earthquakes, heavy rainfall, or artificial excavations. In cases where the point cloud system and engineering independent coordinate system are different and the DTM is disturbed by the multiplicative random error, our method is a more practical choice with stabler precision compared with previous methods.

Author Contributions

Methodology, L.Z., X.M., Z.X., S.Z., C.H. and G.C.; software, L.Z., Z.X. and G.C.; validation, Y.Z.; writing—original draft preparation, L.Z.; writing—review and editing, C.H.; visualization, C.H.; funding acquisition, S.Z. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

Not applicable.

Conflicts of Interest

The authors declare no conflict of interest.

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures and Tables

Figure 1. Spherical target.

Figure 2. Point cloud of spherical target.

Figure 3. GNSS elevation point before and after being disturbed by multiplicative error.

Figure 4. DTM.

Figure 5. DTM with multiplicative error.

Figure 6. (a) Optical remote sensing imagery of Lashagou NO. 3 landslide acquired from Google Earth; (b) UAV image of Lashagou NO. 3 landslide.

Figure 7. (a) Study area and photos showing Lashagou NO. 3 landslide disease. (b) ground fracture (c) staggered platform (d) tensile crack caused by the landslide.

Figure 8. Triangular pyramid target.

Figure 9. The point cloud of triangular pyramid.

Figure 10. Prism ball target.

Figure 11. Prism ball target.

Figure 12. Point cloud of phase 1.

Figure 13. Lihua landslide point cloud of phase 1.

Figure 14. UAV image of Lihua landslide.

Figure 15. GNSS elevation obtained by the RWLS method and LS.

Figure 16. The regularization parameters with the number of iterations.

Figure 17. DTM by the RWLS method.

Figure 18. The regularization parameters with the number of iterations.

Figure 19. Initial point cloud of landslide in two phases.

Figure 20. Landslide point cloud after the PCR of two period data.

Figure 21. Photos of Lashagou NO.3 landslide.

Figure 22. Point cloud after accurate registration.

Figure 23. PCR based on triangular target (1 m wide detail in Figure 20).

Figure 24. Lashagou NO.3 landslide DTM.

Figure 25. Deformation of Lashagou NO.3 landslide.

Figure 26. Point cloud of Lihua landslide.

Figure 27. DTM disturbed by multiplicative error.

Figure 28. Deformation of Lihua landslide.

Table 1

Simulation Data.

Point Number	y	Point Number	y
1	8.8529	17	52.9943
2	9.3072	18	66.6330
3	10.0079	19	80.6002
4	8.0547	20	80.1413
5	13.7458	21	106.8200
6	13.1323	22	112.3601
7	12.5329	23	136.7752
8	16.5769	24	184.5394
9	13.1390	25	254.6235
10	17.2795	26	242.7289
11	19.0292	27	323.5583
12	22.8096	28	360.2978
13	26.2168	29	482.2399
14	27.0681	30	453.0284
15	34.8363	31	595.4323
16	40.9079

Table 2

Parameters of Leica P50 3D laser scanner [60].

Parameter	Value
Scan range mode	from 0.4 to120 mfrom 0.4 to 270 mfrom 0.4 to 570 m,>1 km
Scan Rate	up to 1,000,000 points per second
Vertical/horizontal field-of-view	360°/290°
Range noise *	0.4 mm rms at 10 m0.5 mm rms at 50 m
Operating temperature	−4° F to + 122° F
Dual-axis compensator	accuracy 1.5″

* at 78% albedo.

Table 3

Parameters of Leica TS30 total station [61].

Parameter		Value
Accuracy of angle	Horizontal and Vertical	0.5″
Distance Measurement Range	Round Prism (GPR1)	3500 m
Accuracy of distance	Standard (prism)	1 mm + 1 ppm

Table 4

The coordinates of control points.

Control Point	x/m	y/m	z/m
SCP1	531.375	533.190	761.551
SCP2	535.557	522.346	762.020
SCP3	479.548	473.771	766.402
SCP4	412.943	401.501	772.336

Table 5

Parameter estimation and accuracy.

Method	${\widehat{\mathit{\beta }}}_1$	${\widehat{\mathit{\beta }}}_2$	${\widehat{\mathit{\beta }}}_3$	${\widehat{\mathit{\beta }}}_4$	${\widehat{\mathit{\beta }}}_5$	${\widehat{\mathit{\beta }}}_6$	${\widehat{\mathit{\sigma }}}_0$	${‖\mathbf{\Delta }\widehat{\mathit{\beta }}‖}_2$
True value	10	4	2	1	0.5	2	0.3	—
LS	19.64	−147.21	442.38	−453.75	189.92	−25.37	12.0984	678.47
bcLS	11.45	−2.19	6.27	7.96	−8.11	3.99	20.98	13.61
RWLS	10.77	2.12	1.46	2.52	−1.46	2.46	1.26	3.28

Table 6

Parameter estimation and accuracy.

Method	${\widehat{\mathit{\beta }}}_1$	${\widehat{\mathit{\beta }}}_2$	${\widehat{\mathit{\beta }}}_3$	${\widehat{\mathit{\beta }}}_4$	${\widehat{\mathit{\sigma }}}_0$	${‖\mathbf{\Delta }\widehat{\mathit{\beta }}‖}_2$
true value	−3.5	15	8	−3	0.3	—
LS	21.68	21.00	−11.25	−14.87	4.52	34.38
bcWLS	−6.44	18.99	8.59	−1.69	4.36	3.18
RWLS	−1.94	14.91	7.46	−3.91	0.42	1.88

Table 7

Coordinates of three triangular pyramid vertices in phase 1.

Point Number	x/m	y/m	z/m
CP1	5.989	22.562	−4.090
CP2	12.648	15.101	−4.803
CP3	32.494	−6.676	−6.480

Table 8

Coordinates of three triangular pyramid vertices in phase 2.

Point Number	x/m	y/m	z/m
CP1	−12.412	19.621	−4.134
CP2	−2.413	19.380	−4.847
CP3	27.049	18.977	−6.523

Table 9

Parameters of DTM estimation and accuracy.

Method	${\widehat{\mathit{\beta }}}_1$	${\widehat{\mathit{\beta }}}_2$	${\widehat{\mathit{\beta }}}_3$	${\widehat{\mathit{\beta }}}_4$	${\widehat{\mathit{\beta }}}_5$	${\widehat{\mathit{\beta }}}_6$	${\widehat{\mathit{\sigma }}}_0$
LS	−0.368	0.340	0.128	−0.007	−0.015	0.0003	13.12
bcWLS	−0.265	0.135	0.313	−0.023	−0.062	0.0005	6.93
RWLS	−0.165	0.102	0.033	0.016	−0.002	0.0267	2.34

Table 10

The estimation and accuracy of DTM parameters.

Method	${\widehat{\mathit{\beta }}}_1$	${\widehat{\mathit{\beta }}}_2$	${\widehat{\mathit{\beta }}}_3$	${\widehat{\mathit{\beta }}}_4$	${\widehat{\mathit{\beta }}}_5$	${\widehat{\mathit{\beta }}}_6$	${\widehat{\mathit{\sigma }}}_0$
LS	−0.532	1.234	0.568	−0.146	−0.051	0.021	26.16
bcWLS	−0.389	0.657	0.425	−0.072	−0.081	0.0063	8.24
RWLS	−0.132	0.231	0.046	0.024	−0.006	0.001	2.34

References

1. Yao, L.; Sun, H.; Zhu, L.; Zhou, Y. Development and application of deformation monitoring system for lanslide at Funchunjiang Dam. Surv. Rev.; 2014; 46, pp. 444-452. [DOI: https://dx.doi.org/10.1179/1752270614Y.0000000109]

2. Gikas, V. 3D Terrestrial laser scanning for geometry documentation and construction management of highway tunnels during excavation. Sensors; 2012; 12, pp. 11249-11270. [DOI: https://dx.doi.org/10.3390/s120811249]

3. Vezocnik, R.; Ambrozic, T.; Sterle, O.; Bilban, G.; Pfeifer, N.; Stopar, B. Use of terrestrial laser scanning technology for long term high precision deformation monitoring. Sensors; 2009; 9, pp. 9873-9895. [DOI: https://dx.doi.org/10.3390/s91209873] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/22303152]

4. He, M.; Zhu, Q.; Du, Z.; Hu, H.; Ding, Y.; Chen, M. A 3D shape descriptor based on contour clusters for damaged roof detection using airborne LiDAR point clouds. Remote Sens.; 2016; 8, 189. [DOI: https://dx.doi.org/10.3390/rs8030189]

5. Kasperski, J.; Delacourt, C.; Allemand, P.; Potherat, P.; Jaud, M.; Varrel, E. Application of a terrestrial laser scanner (TLS) to the study of the Séchilienne Landslide (Isère, France). Remote Sens.; 2010; 2, pp. 2785-2802. [DOI: https://dx.doi.org/10.3390/rs122785]

6. Crepaldi, S.; Zhao, Y.; Lavy, M.; Amanzio, G.; Suozzi, E.; De Maio, M. Landslide analysis by multi-temporal terrestrial laser scanning (TLS) data: The Mont de la Saxe landslide. Rend. Online Della Soc. Geol. Ital.; 2015; 35, pp. 92-95. [DOI: https://dx.doi.org/10.3301/ROL.2015.72]

7. Wang, D.; Hollaus, M.; Schmaltz, E.; Wieser, M.; Reifeltshammer, D.; Pfeifer, N. Tree stem shapes derived from TLS data as an indicator for shallow landslides. Procedia Earth Planet. Sci.; 2016; 16, pp. 185-194. [DOI: https://dx.doi.org/10.1016/j.proeps.2016.10.020]

8. Wang, G.; Philips, D.; Joyce, J.; Rivera, F. The integration of TLS and continuous GPS to study landslide deformation: A case study in Puerto Rico. J. Geod. Sci.; 2011; 1, pp. 25-34. [DOI: https://dx.doi.org/10.2478/v10156-010-0004-5]

9. Cina, A.; Piras, M. Performance of low-cost GNSS receiver for landslides monitoring: Test and results. Geomat. Nat. Hazards Risk; 2015; 6, pp. 497-514. [DOI: https://dx.doi.org/10.1080/19475705.2014.889046]

10. Senkaya, M.; Babacan, A.E.; Karslı, H.; San, B.T. Origins of diverse present displacements in a paleo-landslide area (Isiklar, Trabzon, northeast Turkey). Environ. Earth Sci.; 2022; 81, pp. 1-24. [DOI: https://dx.doi.org/10.1007/s12665-022-10372-2]

11. Stiros, S.C.; Vichas, C.; Skourtis, C. Landslide monitoring based on geodetically derived distance changes. J. Surv. Eng.; 2004; 130, pp. 156-162. [DOI: https://dx.doi.org/10.1061/(ASCE)0733-9453(2004)130:4(156)]

12. Tzenkov, T.; Gospodinov, S. Geometric analysis of geodetic data for investigation of 3D landslide deformations. Nat. Hazards Rev.; 2003; 4, pp. 78-81. [DOI: https://dx.doi.org/10.1061/(ASCE)1527-6988(2003)4:2(78)]

13. Klimeš, J.; Rowberry, M.D.; Blahůt, J.; Briestenský, M.; Hartvich, F.; Košťák, B.; Rybář, J.; Stemberk, J.; Štěpančíkováet, P. The monitoring of slow-moving landslides and assessment of stabilisation measures using an optical–mechanical crack gauge. Landslides; 2012; 9, pp. 407-415. [DOI: https://dx.doi.org/10.1007/s10346-011-0306-4]

14. Baroň, I.; Supper, R. Application and reliability of techniques for landslide site investigation, monitoring and early warning–outcomes from a questionnaire study. Nat. Hazards Earth Syst. Sci.; 2013; 13, pp. 3157-3168. [DOI: https://dx.doi.org/10.5194/nhess-13-3157-2013]

15. Shi, X.; Hu, X.; Sitar, N.; Kayen, R.; Qi, S.; Jiang, H.; Wang, X.; Zhang, L. Hydrological control shift from river level to rainfall in the reactivated Guobu slope besides the Laxiwa hydropower station in China. Remote Sens. Environ.; 2021; 265, 112664. [DOI: https://dx.doi.org/10.1016/j.rse.2021.112664]

16. Guo, F.; Meng, X.; Qi, T.; Dijkstra, T.; Thorkildsen, J.K.; Yue, D.; Chen, G.; Zhang, Y.; Dou, X.; Shi, P. Rapid onset hazards, fault-controlled landslides and multi-method emergency decision-making. J. Mt. Sci.; 2022; 19, pp. 1357-1369. [DOI: https://dx.doi.org/10.1007/s11629-021-6941-x]

17. Guo, C.; Xu, Q.; Dong, X.; Liu, X.; Yu, J. Geohazard Recognition by Airborne LiDAR Technology in Complex Mountain Areas. Geomat. Inf. Sci. Wuhan Univ.; 2021; 46, pp. 1538-1547.

18. Abellán, A.; Vilaplana, J.; Calvet, J.; García-Sellés, D.; Asensio, E. Rockfall monitoring by terrestrial laser scanning-case study of the basaltic rock face at Castellfollit de la Roca (Catalonia, Spain). Nat. Hazards Earth Syst. Sci.; 2011; 11, pp. 829-841. [DOI: https://dx.doi.org/10.5194/nhess-11-829-2011]

19. Kayen, R.; Pack, R.; Bay, J.; Collins, B. Ground-Lidar visualization of surface and structural deformatios of the Niigata Ken Chuetsu. Earthq. Spectra; 2006; S1, pp. 147-162. [DOI: https://dx.doi.org/10.1193/1.2173020]

20. Liu, W.; Xie, M. Landslide monitoring based on point cloud density characteristics. Rock Soil Mech.; 2020; 41, pp. 3748-3756.

21. Zhao, J.; Du, S.; Liu, Y.; Saif, B.S.; Hou, Y.; Guo, Y.C. Evaluation of the stability of the palatal rugae using the three-dimensional superimposition technique following orthodontic treatment. J. Dent.; 2022; 119, 104055. [DOI: https://dx.doi.org/10.1016/j.jdent.2022.104055] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/35121138]

22. Yang, M.; Gan, S.; Yuan, X.; Gao, S.; Zhu, Z.; Yu, H. Study on data acquisition and registration experiment of terretrial laser scanning point clouds under complicated banded terrain condition. Bull. Surv. Mapp.; 2018; 5, pp. 35-40.

23. Zeybek, M.; Şanlıoğlu, İ. Accurate determination of the Taşkent (Konya, Turkey) landslide using a long-range terrestrial laser scanner. Bull. Eng. Geol. Environ.; 2015; 74, pp. 61-76. [DOI: https://dx.doi.org/10.1007/s10064-014-0592-x]

24. Hsieh, Y.C.; Chan, Y.C.; Hu, J.C.; Chen, Y.Z.; Chen, R.F.; Chen, M.M. Direct measurements of bedrock incision rates on the surface of a large dip-slope landslide by multi-period airborne laser scanning DEMs. Remote Sens.; 2016; 8, 900. [DOI: https://dx.doi.org/10.3390/rs8110900]

25. Xian, Y.; Xiao, J.; Wang, Y. A fast registration algorithm of rock point cloud based on spherical projection and feature extraction. Front. Comput. Sci.; 2019; 13, pp. 170-182. [DOI: https://dx.doi.org/10.1007/s11704-016-6191-1]

26. Kümmerle, J.; Kühner, T.; Lauer, M. Automatic calibration of multiple cameras and depth sensors with a spherical target. Proceedings of the 2018 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS); Madrid, Spain, 1–5 October 2018; pp. 1-8.

27. Zhou, Y.; Han, D.; Hu, K.; Qin, G.; Xiang, Z.; Ying, C.; Zhao, L.; Hu, X. Accurate virtual trial assembly method of prefabricated steel components using terrestrial laser scanning. Adv. Civ. Eng.; 2021; 2021, 9916859. [DOI: https://dx.doi.org/10.1155/2021/9916859]

28. Li, J.; Hu, Q.; Ai, M. Point cloud registration based on one-point ransac and scale-annealing biweight estimation. IEEE Trans. Geosci. Remote Sens.; 2021; 59, pp. 9716-9729. [DOI: https://dx.doi.org/10.1109/TGRS.2020.3045456]

29. Chen, Z.; Zhang, B.; Han, Y.; Zuo, Z.; Zhang, X. Modeling accumulated volume of landslides using remote sensing and DTM data. Remote Sens.; 2014; 6, pp. 1514-1537. [DOI: https://dx.doi.org/10.3390/rs6021514]

30. Artese, S.; Perrelli, M. Monitoring a landslide with high accuracy by total station: A DTM-based model to correct for the atmospheric effects. Geosciences; 2018; 8, 46. [DOI: https://dx.doi.org/10.3390/geosciences8020046]

31. Pawluszek, K.; Borkowski, A.; Tarolli, P. Sensitivity analysis of automatic landslide mapping: Numerical experiments towards the best solution. Landslides; 2018; 15, pp. 1851-1865. [DOI: https://dx.doi.org/10.1007/s10346-018-0986-0]

32. Xu, P.; Shimada, S. Least squares parameter estimation in multiplicative noise models. Commun. Stat. Simul. Comput.; 2000; 29, pp. 83-96. [DOI: https://dx.doi.org/10.1080/03610910008813603]

33. Shi, Y.; Xu, P.; Peng, J.; Shi, C.; Liu, J. Adjustment of measurements with multiplicative errors: Error analysis, estimates of the variance of unit weight, and effect on volume estimation from LiDAR-type digital elevation models. Sensors; 2014; 14, pp. 1249-1266. [DOI: https://dx.doi.org/10.3390/s140101249] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/24434880]

34. Wang, L.; Yu, D. Virtual observation method to ill-posed total least squares problem. Acta Geod. Cartogr. Sin.; 2014; 43, pp. 575-581.

35. Xu, P. Truncated SVD methods for discrete linear ill-posed problems. Geophys. J. Int.; 1998; 135, pp. 505-514. [DOI: https://dx.doi.org/10.1046/j.1365-246X.1998.00652.x]

36. Koch, K.R.; Kusche, J. Regularization of geopotential determination from satellite data by variance components. J. Geod.; 2002; 76, pp. 259-268. [DOI: https://dx.doi.org/10.1007/s00190-002-0245-x]

37. Hoerl, A.E.; Kennard, R.W. Ridge regression: Biased estimation for nonorthogonal problems. Technometrics; 1970; 12, pp. 55-67. [DOI: https://dx.doi.org/10.1080/00401706.1970.10488634]

38. Wang, L.; Chen, T. Virtual Observation Iteration Solution and A-Optimal Design Method for Ill-Posed Mixed Additive and Multiplicative Random Error Model in Geodetic Measurement. J. Surv. Eng.; 2021; 147, 04021016. [DOI: https://dx.doi.org/10.1061/(ASCE)SU.1943-5428.0000363]

39. Du, S.; Zheng, N.; Ying, S.; Liu, J. Affine iterative closest point algorithm for point set registration. Pattern Recognit. Lett.; 2010; 31, pp. 791-799. [DOI: https://dx.doi.org/10.1016/j.patrec.2010.01.020]

40. Li, P.; Wang, R.; Wang, Y.; Tao, W. Evaluation of the ICP algorithm in 3D point cloud registration. IEEE Access; 2020; 8, pp. 68030-68048. [DOI: https://dx.doi.org/10.1109/ACCESS.2020.2986470]

41. Pomerleau, F.; Colas, F.; Siegwart, R.; Magnenat, S. Comparing ICP variants on real-world data sets. Auton. Robot.; 2013; 34, pp. 133-148. [DOI: https://dx.doi.org/10.1007/s10514-013-9327-2]

42. Deakin, R.E. A note on the Bursa-Wolf and Molodensky-Badekas Transformations; School of Mathematical and Geospatial Sciences, RMIT University: Melbourne, VIC, Australia, 2006; pp. 1-21.

43. Golub, G.H.; Van Loan, C.F. An analysis of the total least squares problem. SIAM J. Numer. Anal.; 1980; 17, pp. 883-893. [DOI: https://dx.doi.org/10.1137/0717073]

44. Oktaba, W. Tests of Hypotheses for the General Gauss-Markov Model. Biom. J.; 1984; 26, pp. 415-424. [DOI: https://dx.doi.org/10.1002/bimj.4710260411]

45. Jeng, F.C.; Woods, J.W. Compound Gauss-Markov random fields for image estimation. IEEE Trans. Signal Processing; 1991; 39, pp. 683-697. [DOI: https://dx.doi.org/10.1109/78.80887]

46. Shi, Y.; Xu, P. Adjustment of measurements with multiplicative random errors and trends. IEEE Geosci. Remote Sens. Lett.; 2020; 18, pp. 1916-1920. [DOI: https://dx.doi.org/10.1109/LGRS.2020.3010827]

47. Davidon, W.C. New least-square algorithms. J. Optim. Theory Appl.; 1976; 18, pp. 187-197. [DOI: https://dx.doi.org/10.1007/BF00935703]

48. Heyde, C.C. Quasi-Likelihood and Its Application: A General Approach to Optimal Parameter Estimation; Springer: New York, NY, USA, 1997; pp. 32-69.

49. Gregor, J.; Fessler, J.A. Comparison of SIRT and SQS for regularized weighted least squares image reconstruction. IEEE Trans. Comput. Imaging; 2015; 1, pp. 44-55. [DOI: https://dx.doi.org/10.1109/TCI.2015.2442511]

50. Hansen, P.C. Analysis of discrete ill-posed problems by means of the L-curve. SIAM Rev.; 1992; 34, pp. 561-580. [DOI: https://dx.doi.org/10.1137/1034115]

51. Vogel, C.R. Non-convergence of the L-curve regularization parameter selection method. Inverse Probl.; 1996; 12, 535. [DOI: https://dx.doi.org/10.1088/0266-5611/12/4/013]

52. Wang, L.; Chen, T.; Zou, C. Weighted least squares regularization iteration solution and precision estimation for ill-posed multiplicative error mode. Acta Geod. Cartogr. Sin.; 2021; 50, pp. 589-599.

53. Xu, P.; Shi, Y.; Peng, J.; Liu, J.; Shi, C. Adjustment of geodetic measurements with mixed multiplicative and additive random errors. J. Geod.; 2013; 87, pp. 629-643. [DOI: https://dx.doi.org/10.1007/s00190-013-0635-2]

54. Shi, Y.; Xu, P. Comparing the estimates of the variance of unit weight in multiplicative error models. Acta Geod. Geophys.; 2015; 50, pp. 353-363. [DOI: https://dx.doi.org/10.1007/s40328-014-0096-y]

55. Hill, C.A.; Harris, M.; Ridley, K.D.; Jakeman, E.; Lutzmann, P. Lidar frequency modulation vibrometry in the presence of speckle. Appl. Opt.; 2003; 42, pp. 1091-1100. [DOI: https://dx.doi.org/10.1364/AO.42.001091] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/12617227]

56. Leigh, C.L.; Kidner, D.B.; Thomas, M.C. The use of LiDAR in digital surface modelling: Issues and errors. Trans. GIS.; 2009; 13, pp. 345-361. [DOI: https://dx.doi.org/10.1111/j.1467-9671.2009.01168.x]

57. Hladik, C.; Alber, M. Accuracy assessment and correction of a LIDAR-derived salt marsh digital elevation model. Remote Sens. Env.; 2012; 121, pp. 224-235. [DOI: https://dx.doi.org/10.1016/j.rse.2012.01.018]

58. Jia, C.; Jin, Z.; Yang, P.; Tang, Z. Stability analysis on landslide in section K181 + 840 ~ K182 + 040 of Lin − Da Highway. J. Lanzhou Petrochem. Polytech.; 2018; 18, pp. 26-28.

59. Jia, X.; Yin, X.; Ren, Y.; Fu, Z. Climate change of Linxia of Gansu Province in recent 43 years. J. Arid Meteorol.; 2012; 30, pp. 249-254.

60. Leica ScanStation P50—Long Range 3D Terrestrial Laser Scanner. Available online: https://leica-geosystems.com/products/laser-scanners/scanners/leica-scanstation-p50 (accessed on 3 June 2022).

61. Leica TS30 Champion’s League. Available online: https://www.bandwork.my/product/pdf/731201461738PMTS30_Brochure_en.pdf (accessed on 3 June 2022).

62. Kobler, A.; Pfeifer, N.; Ogrinc, P.; Todorovski, L.; Oštir, K.; Džeroski, S. Repetitive interpolation: A robust algorithm for DTM generation from Aerial Laser Scanner Data in forested terrain. Remote Sens. Environ.; 2007; 108, pp. 9-23. [DOI: https://dx.doi.org/10.1016/j.rse.2006.10.013]