AMETHOD FOR MEASURING LARGE-SCALE DEFORMATION OF LANDSLIDE BODIES BASED ON NAP-OF-THE-OBJECT PHOTOGRAMMETRY

Geological disasters such as landslides and debris flows pose a serious threat to human life and property. To mitigate this risk, monitoring and early warning systems are essential. However, monitoring high-angle landslide areas can be challenging due to the steep and complex terrain, making it difficult to carry out large-scale and refined deformation measurements using existing methods. This paper proposes a method to measure the large-scale deformation of landslide bodies based on nap-of-the-object photogrammetry. The method uses UAVs to acquire high-resolution 3D models of steep and high landslide areas at multiple time periods. Then, 3D model matching is employed to obtain accurate variation information for terrain deformation measurement. To obtain fine 3D models, a terrain-adaptive nap-of-the-object photogrammetric flight planning is applied to design the optimal photographic positions and directions for capturing ultra-high-resolution images. The images are processed using photogrammetric principles and technologies to produce fine 3D models. For terrain deformation measurement, an algorithm is proposed to obtain 3D correspondences by fusing DEM differential and 3D model texture matching. The 3D points variation vectors are then calculated, and the large-scale deformation measurement results of the landslide body can be derived after the vectors are aggregated. Experiments were conducted on the Lijiebeishan landslide in Gansu Province, western China. The results showed that the proposed deformation measurement method was highly effective in accurately detecting areas with displacement greater than 5 cm, and the large-scale deformation trend is consistent with GNSS predictions. In conclusion, the proposed method is an effective way to measure the large-scale deformation of landslide bodies in high-angle landslide areas, providing a valuable tool for monitoring and early warning systems. * Corresponding author


INTRODUCTION
China's high mountain terrains are prone to geological disasters such as landslides and debris flows, which pose a significant threat to human life and property.On average, over 20,000 landslides and debris flow disasters occur each year, resulting in thousands of casualties.These disasters can occur suddenly, as was the case in June 2022 when a train traveling from Guiyang to Guangzhou was hit by a debris flow, resulting in 9 injuries and deaths.To eliminate this threat, monitoring, early warning, and control measures have been adopted.Among these measures, monitoring is the most fundamental work.In the early stages of monitoring, manual observation methods were typically used, but they have limitations such as low efficiency, high construction difficulty, and occupational hazards (Yang, 2005).New technologies such as GNSS, close-range photogrammetry, Terrestrial Laser Scanning (TLS), and Synthetic Aperture Radar (SAR) have emerged since the 20th century, providing strong support for geological disaster monitoring (Yang et al., 2019).Placing GNSS observation points is the most direct and effective way to obtain landslide deformation (Li et al., 2021;Shen et al., 2021).However, due to the steep and complex terrain of high-angle landslides, finding suitable placement points is difficult.Even if placement points are found, the construction and maintenance is extremely difficult and costly, making it impossible to carry out largescale monitoring (Ma et al., 2022;Kuang et al., 2022).Synthetic Aperture Radar Interferometry (InSAR) can carry out largescale monitoring, but its observation accuracy is limited (Cenni et al., 2021;Liu et al., 2022).Close-range photogrammetry and TLS also suffer from accuracy limitations at long-range observations (Yang et al., 2018;Zhao et al., 2022).Therefore, further research is needed to develop more accurate and efficient monitoring methods for high mountain terrains to ensure the safety of human life and property.
The rapid development of unmanned aerial vehicle (UAV) technology has led to significant advancements in UAV photogrammetry, enabling increasingly convenient and fast data acquisition.This has resulted in a crucial role in conventional terrain mapping, urban 3D modeling, and geological disaster monitoring.Using UAV photogrammetry for geological disaster monitoring involves two primary tasks.The first is to capture multi-temporal photographs of the designated area using UAVs, which are then processed to generate time-series fine 3D model data.Various methods, including aerial vertical, terrain, oblique, and optimized views photography, are fully capable of terrain mapping and urban 3D modeling (Koch et al., 2019;Zhou et al., 2020;Li et al., 2022).However, high and steep landslides pose challenges due to their complex terrain, and without adaptive photography flight paths, the obtained images may suffer from missing information and poor photography angles.In the second task of multi-temporal 3D model comparison and analysis, direct digital elevation model (DEM) differential or dense point cloud registration can only identify the location of changes, but cannot accurately describe the changes in deformation points.While matching original images can provide accurate position relationships, redundant observations and the failure to identify the best pairs of images can result in imprecise descriptions of deformation points.
To address these issues, this paper presents a novel approach to measuring large-scale deformation of landslide bodies using nap-of-the-object photogrammetry.The method involves using UAVs to capture multi-temporal fine 3D models of high and steep landslide areas, and then matching the data to obtain accurate 3D deformation information of target points.This enables terrain deformation measurement.To obtain the most ideal images, a terrain-adaptive nap-of-the-object photography flight planning is first carried out for the high and steep landslides areas to design the optimal photography positions and directions.This ensures ultra-high-resolution images are obtained.Photogrammetry principles and technologies are then used to produce fine realistic 3D models.For terrain deformation measurement, an algorithm is proposed to obtain 3D correspondences by fusing DEM differential and 3D model texture matching.The 3D points variation vectors can be calculated, and the large-scale deformation measurement results of the landslide body can be formed after the vectors are aggregated.Overall, this approach provides a reliable and accurate method for measuring large-scale deformation of landslide bodies.This paper presents several significant contributions to the field of landslide deformation measurement: (1) A novel UAV flight planning method is proposed, which adapts to the terrain of high and steep slopes, allowing for the capture of fine 3D models of landslide bodies.
(2) An algorithm is proposed to obtain 3D correspondences by fusing DEM differential and 3D model texture matching, providing an accurate means of measuring deformation.
(3) The proposed method is applied to the Lijiebeishan landslide in Gansu Province, western China, and large-scale deformation measurement results of the landslide body are obtained, demonstrating the effectiveness of the approach in real-world scenarios.

FINE 3D MODELING
High and steep landslide bodies with complex terrains can be accurately modeled using UAV-based fine 3D modeling.This process involves the novel nap-of-the-object photogrammetry theory and technique proposed by Professor Zuxun Zhang of Wuhan University (Tao et al., 2019).This technique segments complex targets into multiple planar or curved surface elements with arbitrary slopes and orientations based on initial shape information.Each surface element is treated as an independent processing unit for close vertical photogrammetry at the best photography angle.As a result, the technique delivers high spatial resolution and accurate positioning.
To ensure the successful implementation of this technique, the acquisition of the most ideal high-resolution images is crucial for fine 3D modeling.The specific methodology involves utilizing initial terrain information to determine the optimal UAV photography positions and directions, plan the most suitable flight paths, guide the UAV to acquire high-resolution images of the high and steep landslide body, and employ photogrammetry principles to reconstruct the fine 3D model.The process and principles are depicted in Figure 1.

Segmentation of Terrain Units
The segmentation of terrain units refers to the process of dividing the study area into operational units using existing data (such as global SRTM data, DEM data, and rough terrain data obtained from UAV imagery processing).The process employs the region growing algorithm that considers various factors such as the similarity of adjacent point normals, the flatness of point sets, and the size of segmentation regions.In this study, the triangulation growing algorithm was employed for segmentation, which involves the following steps: (1) establishing a Delaunay triangulation network using existing data points; (2) calculating the normal vectors and local flatness (the average flatness of the three vertices) of all triangles in the network and sorting the triangles in ascending order of flatness; (3) starting with any triangle and growing the region based on the condition that the triangle normal vector angle ∆ ≤ , until the number of triangles in the region reaches the threshold , then fitting a spatial plane P using the vertices of the triangles contained in the region; (4) continuing to grow the region based on the conditions that the distance ∆ ≤ from the triangle center to plane P and the maximum length of the growth region ≤ , until the number of newly added triangles reaches , then refitting plane P and updating the surface equation and normal vector; (5) repeating steps (3) and ( 4) until no new triangles are found.
In the algorithm, the empirical values of , , , and are 15°, 50, 1000 m, and 50, respectively, and is 1/5 of the photography distance.

Trajectory Planning
Trajectory planning needs to be performed in a plane.For any segmented terrain unit, a spatial plane P can be obtained by fitting the terrain points, and a local coordinate system , , , for the terrain unit can be established as shown in Figure 2. Let the coordinate system of the object be , , , , the equation of the plane P can be fitted based on the initial terrain information using point-normal form, which can be expressed as: where , , = the normal vector of the plane P, and2 + 2 + 2 = 1 , , = the coordinate of the vertical foot point The axes of the local coordinate system are defined as the intersection of the spatial plane P and the horizontal plane ( = 0), whicn can be calculated as follows: (2) where = the normal vector of the horizontal plane is consistent with .and , follow the right-hand rule to constitute a three-dimensional Cartesian coordinate system, then The flight path is planned in plane P ' , which is parallel to plane P, with a distance of d.To minimize the time spent by the UAV on switching between flight lines during the photographing process, the minimum bounding rectangle of the target's outer contour is calculated.The long side of the rectangle is selected as the primary route direction, as it requires fewer flight lines to cover the top surface of the target.As shown in Figure 2, the rectangle in the plane P ' can be described by two sets of parameters: the midpoints 1 , 2 of the two short sides and the width w of the rectangle.Then the direction (i.e., the direction of the vector 1 , 2 ) of the long side of the rectangle serves as the primary direction of the flight path.
During photographing, principal axis of the camera is perpendicular to plane P. The distance between adjacent exposure points in the same flight line is ∆ , and the distance between adjacent flight lines perpendicular to the flight direction is ∆.Therefore, the number of flight lines required to photograph a terrain unit is: where = rounding upwards for The endpoints of these flight lines can be calculated according to the following formula: where = 1,2 ∈ [0, 1 ) = the unit vector in the direction after rotating the vector 1 , 2 by 90°clockwise The coordinate formula for each trajectory point, except the endpoints, is: where ∈ [1, 1 , 2 ∆ ) = the unit vector of the vector 1 , 2 In the next step, the rotation angle of the lens for each trajectory point is augmented by θ − 90°, where θ denotes the acute angle between the slope plane P and the horizontal plane.Finally, the calculated 1 flight lines are sorted, and the trajectory planning outcomes are obtained by combining them in a zigzag manner.

3D Modeling
After utilizing pre-planned flight paths to direct UAVs in capturing high-resolution images of high and steep landslide bodies, the acquired data can be processed to generate fine 3D models using photogrammetric principles and techniques.In this paper, the 3D reconstruction is processed with DPGrid, a software developed by Wuhan University, which incorporates mature theories and techniques (Xi and Duan, 2020;Gu, 2021).To provide concise and essential processing steps, redundant information is excluded in this section, and only the critical steps are presented as follows: (1) Aerial triangulation is employed to obtain the interior and exterior orientation elements of each image; (2) Dense matching is used to generate a dense point cloud; (3) An initial 3D mesh model is constructed; (4) The mesh is optimized to produce a fine 3D model.

DEFORMATION MEASUREMENT
After acquiring multi-temporal fine 3D models, deformation measurement can be accomplished through the comparison of two such models.A crucial aspect of 3D model comparison analysis is to identify correspondences between the two models.In this regard, this paper proposes an algorithm to obtain 3D correspondences by fusing DEM differential and 3D model texture matching.

The Texture Matching Algorithm Fused with DEM Differential
The fundamental concept of the texture matching algorithm fused with DEM differential is to first identify the deformation locations between two 3D models using DEM differential, and then seek the correspondences by texture matching between the deformed positions of the two 3D models based on the triangles.The algorithm mainly involves the following three steps: (1) DEM differential DEM differential is a technique that involves computing the disparity between the elevation values of corresponding coordinates in two DEMs from different time periods.This operation provides elevation change values, which can be used to identify deformation regions.To accurately identify the areas to be matched, a threshold is established based on the magnitude of these changes.The threshold can be determined by the magnitude of the expected deformation or the standard deviation, such as setting it to 5 cm or two times the standard deviation.
(2) Feature points detection In the previous 3D model, the Harris algorithm (C.Harris et al., 1998) is employed to extract features from regions with significant changes.Specifically, for each triangular element on the 3D model, the horizontal gradient and vertical gradient of the element are calculated using the first-order differential operator of the Gaussian kernel function to obtain the autocorrelation matrix 1 : The four elements in the autocorrelation matrix 1 undergo Gaussian smoothing to derive a new autocorrelation matrix 2 , as follows: where , = the two-dimensional Gaussian kernel function The corner response function CRF is then computed as follows: where r = a constant value between 0.04 and 0.06 To identify suitable feature points from the given 3D model, a criterion based on the CRF value is employed.Specifically, a pixel point is considered as a potential feature point if its CRF value exceeds a predefined threshold k, and it represents the maximum value within its defined neighborhood.To ensure the selection of only the most salient feature points, a nonmaximum suppression technique is applied to remove any extreme points that do not meet these criteria.This approach results in a set of highly relevant feature points for subsequent processing tasks.
(3) Similarity evaluation For each feature point extracted from the previous 3D model, the closest 3 to 5 points in the subsequent model are selected as candidate points for matching.To obtain the best candidate points, the probability relaxation matching method (Q.Zhou et al., 2022) is employed to screen the candidate points.Probability relaxation can resolve the matching using the local support provided by adjacent points and minimize the energy function formulated from the employed constraint (Zhang, 2005).In this study, an adaptive model is established using terrain as the constraint condition.We select the cubic surface to fit the disparity map with the support of several similar points around it and update the probability values by fitting residuals.
The disparity surface reflects the changes in terrain and is almost identical to the real terrain surface.Therefore, we optimize the local disparity surface by fitting the local terrain surface.The cubic surface equation used is given by: Where 0 = the initial probability of the j'th candidate point being the true matching point (, ) = the similarity between the current feature point and the j'th candidate point (, ) is defined as the weighted sum of the correlation coefficient and variance, and the calculation formula is as follows: The values of (, ) and • (, ) range from 0 to 1, and the value of (, ) ranges from 0 to 2.
The compatibility coefficient is an important parameter for the relaxation iteration, which expresses the reasonable degree of the j'th candidate point of the i'th point and the l'th candidate point of the k'th point being the true values at the same time.
The calculation formula is as follows: where = the residual error between the disparity value of the j'th candidate point of the i'th point and the fitted value λ = the adjustment parameter The smaller the value of is, the more suitable the candidate point fits the previously fitted local terrain surface, and more likely it is to be a correct matching point; .
In each iteration, the surface parameters get updated by refitting, and the compatibility coefficient is recalculated.The probability of the current round is updated from the probability of the previous round.The probability relaxation refers to the following equation: where n = the number of candidate points −1 = the matching probability of the j'th candidate of the i'th point in the r-1'th iteration = the relative probability of the r'th iteration = the true probability after normalization of = the set of similar neighbors of the current feature point 0 , 1 = the relaxation coefficient = the weight determined by the distance from the midpoint of the set to the feature point.is calculated as follows: The iteration ends when the probability change value does not exceed 1 × 10 −4 .Finally, the candidate point with the maximum probability and the residual between it and the fitted value not exceeding the threshold d is selected as the matching point.
By following these three steps, the matching of correspondences on the two-temporal 3D models can be accomplished.

Variation Vectors Aggregation
The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XLVIII-5/W1-2023 International Conference on Geomatics Education -Challenges and Prospects (ICGE22), 10-12 May 2023, Hong Kong SAR, China After obtaining 3D correspondences through matching, the 3D variation vectors ∆ between the correspondences can be calculated to derive the deformation of the landslide body in 3D space.The calculation formula is as follows: where 1 , 1 , 1 = the coordinates of the correspondences int the previous 3D model 2 , 2 , 2 = the coordinates of the correspondences in the subsequent 3D model Then the deformation magnitude dF between correspondences is calculated as: Finally, a cubic surface is used to fit the deformation of the 3D points.And after aggregation, the measurement results of largescale deformation of the landslide body can be obtained.The fitting formula is as follows:

EXPERIMENTS AND ANALYSIS
In this study, experiments were conducted on the Lijiebeishan landslide, which is situdtedd in Gansu Province, western China.The area is characterized by a high elevation, significant drop, complex terrain, and sparse vegetaion.Owing to the impact of rainfall and frost swelling, the landslide has exhibited continuous deformation over the years (X.Zhou et al., 2022).Notably, a large-scale mudslide occurred in August 2020, which blocked half of the Bailong River and caused the collapse of the road connecting Beishan Village, posing a threat to the safety of 751 people in 152 households (Di, 2022).

Trajectory Planning
To obtain the initial terrain of the landslide, a UAV was utilized to capture original images with a resolution of 0.2 meters using conventional photography.Subsequently, the images were processed using DPGrid software to generate a LAS point cloud.
With the assistance of the LAS point cloud, terrain segmentation and trajectory planning were conducted in accordance with the methodology outlined in section 2. The resulting partial flight planning outcomes are shown in Figure 3.

Production of 3D Models
The original images were processed using DPGrid software, resulting in the generation of two-temporal realistic 3D models of the Lijiebeishan landslide, as shown in Figure 5. Owing to the extensive scope of the landslide, only partial results are presented in this paper.Notably, the 3D model provides a detailed view of the surface cracks and fissures, highlighting fine information pertaining to the landslide morphology.

Deformation Measurement
The disparity between the two DEMs was derived via differential processing.Subsequently, semantic segmentation was excuted based on the standard deviation x, resulting in the outcomes depicted in Figure 6.The color scheme adopted for the figure illustrates the magnitude of the elevation modifications.Specifically, yellow-red indicates areas where the elevation has increased, green indicates areas where the elevation has remained relatively stable, and blue indicates areas where the elevation has decreased.The chromatic scale of the figure correlates with the extent of the elevation changes, with the hues at the two extremes signifying the most drastic changes.It is important to note that the drastic elevation changes are mainly caused by the landslide, and as such, the areas where the landslide has changed significantly can be identified through the approach utilized.

CONCLUSION
In this paper, we focus on the Lijiebeishan landslide in Gansu Province, western China, and propose a novel method for measuring large-scale deformation of landslide bodies based on nap-of-the-object photogrammetry.
The experiments demonstrate that the proposed method is effective and provides a solution to the limitations of traditional methods.It has several advantages over traditional methods.First, it enables the obtaining of detailed information about the target area, including the distribution of cracks, fissures, and other signs of instability.Second, it detects the three-dimensional deformation information of landslides.Third, it provides accurate and refined deformation data on a large scale.It has the potential to improve early warning systems and contribute to the safety of people and property.By accurately measuring the deformation of high and steep slopes, the risk of geological hazards can be reduced.The proposed method can also be applied to other areas where the terrain is complex and steep, such as cliffs, mountains, and canyons, to provide valuable information for geological research and engineering applications.

Figure 2 .
Figure 2. Local coordinate system establishment and nap-ofthe-object trajectory planning , s = non-negative integers Different from general surface fitting, this study treats all candidate disparity values of each feature point as weighted observations, where the weights are the probability values of the disparities.The initial probability value for each iteration is calculated as follows:

Figure 4 .
Figure 4. Partial images of the Lijiebeishan landslide

Figure 6 .
Figure 6.Differential results of two DEMsTo determine the deformation magnitude and trend of the landslide body in 3D space, the two-temporal 3D models are matched to obtain 3D correspondences.Subsequently, the 3D variation vectors between correspondences are calculated to obtain the 3D changes of the landslide body, as shown in Figure7.To verify the accuracy of the calculated results, the outcomes are compared with the GNSS-measured values at five designated points.The positions of these five GNSS receivers are shown in Figure8.

Table 1
compares the displacement values obtained from the 3D model calculation and the GNSS-measured values.The results indicate that while the displacement values obtained from the 3D model calculation are slightly less accurate, they exhibit good consistency overall with the GNSS-measured values.This affirms the effectiveness of the proposed method in measuring the deformation of the Lijiebeishan landslide.

Table 1 .
Comparison of displacement values obtained from 3D model calculation and GNSS measurement