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

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.


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 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.

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.
(N y ∈ N) are assumed in the n-dimensional space R 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], 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 → 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) The new rotation matrices and shift vectors are computed by minimizing the square distance [40], 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) ε to evaluate the registration accuracy, where R k , → t k is the argument of iteration k.

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].

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 (i = 1, 2, · · · n), and the multiplicative error model can be defined as follows [46], where ε i is the zero mean random error, f (β) and ε are the random variables, β is the t-dimensional unknown parameter vector to be estimated, 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 ε are equal and independent, if every f i (·) is a linear function of variable β, that is, f i (β) = x T i β, x T i is a t-dimensional row vector, the multiplicative error model can be written as follows, The Formula (6) can be rewritten in the following form, where X = (x 1 , x 2 , · · · , 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β, the greater the error. The covariance matrix Q y is a function of the unknown parameter β. 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], whereβ LS is the least square estimation of the multiplicative error model. According to the law of cofactor propagation, the cofactor matrix of e and y is, We also utilize an unbiased estimation of weighted least-square (bcWLS) method; the parameter estimation can be expressed as follows [48], whereβ 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 α is introduced to construct the regularization criterion of the ill-posed multiplicative error model, According to the definition of RWLS, the parameter estimation can be expressed as follows [49],β whereβ RW LS is the RWLS solution of ill-posed multiplicative error model, and I n is identity matrix. The term α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β and β are the function of the regularization parameter α, and the L curve method is used to calculate the regularization parameters α. The main procedures are as follows [50,51], (1) Taking y − Aβ as the x-axis coordinate and β as the y-axis coordinate, we get multiple sets of coordinate points ( y − Aβ , β ), (2) These coordinate points are fitted to a curve similar to L shape, and the α value corresponding to the point with the largest curvature is used as the estimation of the regularization parameter.
As evident from Equation (12), Q e is the nonlinear function ofβ RW LS on the left side of the Equation (12), andβ RW LS is the nonlinear function of Q 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 e i ; (3) The regularization parameter α i is obtained by L curve method; (4) The iterative formula is as follows: (5) Repeat steps (2)-(4) until the value of the parameter subtraction is smaller than the threshold, where ε th is defined as 10 −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,σ whereσ 2 LS = unit weight variance. The unit weight variance using bcWLS can be expressed as follows [53], whereσ 2 bc = 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],σ 2 0 =ê The unit weight variance using RWLS can be expressed as follows [52], whereσ RW LS is the unit weight variance. In Equation (19), bβ RW LS is defined as follows, 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], where (X i , y i ) represents the x and y coordinates of the ground point, H(X i , y i ) represents the corresponding elevation of the point, and ε m represents the multiplicative error vector, 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.

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, 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, where Y is the vector disturbed by multiplicative error, 1 is the 31-dimensional column vector with all elements being 1, and ε 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.

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

Simulation Experiment II
This simulation uses DTM data. This paper mainly considers that the model is illposed, so it simulates the ill-conditioned DTM model and generates the DTM using the interpolation function method [53]  where 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, The multiplicative error model is as follows, 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 ε 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

Actual Experiment
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 Plat eau, 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 stra tum 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 2 , 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 wal 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 soi landslide group. The optical remote sensing imagery of the landslide is shown in Figure   Figure 5. DTM with multiplicative error.

Actual Experiment
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 2 , 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.
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. 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.   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. 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.  (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. (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.
The creep deformation of a landslide is a slow dynamic process. In this paper, a onestation 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.

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 excavationinduced shallow loess engineering landslide. Once a large-scale instability occurs, it may The creep deformation of a landslide is a slow dynamic process. In this paper, a onestation 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.

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 excavationinduced shallow loess engineering landslide. Once a large-scale instability occurs, it may

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.  Figure 13. The UAV image of the Lihua landslide acquired on 4 June 2022 is shown in Figure 14.

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 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.  Figure 16. The condition number of the normal equation N = A T Q −1 e A is 8.3098 × 10 5 , 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.  Figure 16. It is evident from Figure 15 that ground elevation points are influenced by multipli cative error, resulting in serious deviation. It is evident from   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 ∆β 2 obtained by LS method is the largest. The ∆β 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 ∆β 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σ 0 deviating from the true value. However, the RWLS algorithm takes these factors into account and theσ 0 is closer to the true value, further verifying the advantages of the algorithm in this paper.

The Result of Simulation Experiment II
Following Xu et al. (2013), it is here assumed thatσ 0 is 0.3, and that the multiplicative random error vectors are independent from each other. The parameter estimation, ∆β 2 andσ 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. Following Xu et al. (2013), it is here assumed that 0 σ is 0.3, and that the multiplicative random error vectors are independent from each other. The parameter estimation, 2 β Δ and 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 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 0 σ obtained by the LS Following Xu et al. (2013), it is here assumed that 0 σ is 0.3, and that the multiplicative random error vectors are independent from each other. The parameter estimation, 2 β Δ and 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 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 0 σ obtained by the LS Figure 18. The regularization parameters with the number of iterations.
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 ∆β 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σ 0 obtained by the LS method is 4.52, and thatσ 0 seriously deviates from the true value. Theσ 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σ 0 , hence theσ 0 is the closest to the true value.  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 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 Tables 7 and 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 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. 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. 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. 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:

The result of PCR Based on Total Station Orientation
(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. 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.
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 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.

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.
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.   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.    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 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. The estimation and accuracy of the parameters are shown in Table 10. Theσ 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 Tables 9 and 10, the 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 As shown in Figure 26, the DTM has a huge deviation because of the multiplicative random errors. As indicated in Tables 9 and 10, theσ 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.

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.