Automatic Method for Detecting Deformation Cracks in Landslides Based on Multidimensional Information Fusion

Abstract

As cracks are a precursor landslide deformation feature, they can provide forecasting information that is useful for the early identification of landslides and determining motion instability characteristics. However, it is difficult to solve the size effect and noise-filtering problems associated with the currently available automatic crack detection methods under complex conditions using single remote sensing data sources. This article uses multidimensional target scene images obtained by UAV photogrammetry as the data source. Firstly, under the premise of fully considering the multidimensional image characteristics of different crack types, this article accomplishes the initial identification of landslide cracks by using six algorithm models with indicators including the roughness, slope, eigenvalue rate of the point cloud and pixel gradient, gray value, and RGB value of the images. Secondly, the initial extraction results are processed through a morphological repair task using three filtering algorithms (calculating the crack orientation, length, and frequency) to address background noise. Finally, this article proposes a multi-dimensional information fusion method, the Bayesian probability of minimum risk methods, to fuse the identification results derived from different models at the decision level. The results show that the six tested algorithm models can be used to effectively extract landslide cracks, providing Area Under the Curve (AUC) values between 0.6 and 0.85. After the repairing and filtering steps, the proposed method removes complex noise and minimizes the loss of real cracks, thus increasing the accuracy of each model by 7.5–55.3%. Multidimensional data fusion methods solve issues associated with the spatial scale effect during crack identification, and the F-score of the fusion model is 0.901.

1. Introduction

Rockfall, debris flows, and other landslide disasters are major threats to humans [1]. Landslides generally undergo long deformation development and evolution processes before a large-scale instability event occurs. As the ground surface can be affected by different types of stress, ground surfaces also exhibit different crack stages and matching characteristics [2]. Therefore, surface cracks, one of the main features of slope deformation, can provide effective forecasting information that is useful for the early identification of landslides, analyzing landslide instability models, and inferring the deformation direction [3,4]. However, as a result of the terrain conditions in mountainous areas, it is difficult to efficiently achieve accurate large-scale inspection results and refined surface crack measurements in these areas using manual, on-site inspection methods [5].

In recent years, the use of a variety of remote sensing data sources, such as Digital Orthophoto Maps (DOMs), point clouds, and Digital Surface Models (DSM) obtained by aerial remote sensing technology, has become a hot research topic for scholars involved in identifying cracks and landslides [6,7,8]. At present, the methods used to detect surface cracks based on remote sensing data mainly rely on manual visual interpretations and manual measurements [9]. However, these cracks are widely distributed, slender, and multibranched, and the background noise is often complex, causing these manual operations to exhibit certain errors and low efficiencies.

With the development of machine learning and information technologies, crack detection based on digital images has become the primary method for automatically extracting cracks. Effective results have been obtained from many studies aiming to identify cracks in concrete pavement, tunnel linings, building walls, and other single background types [10,11,12]. These studies use traditional image processing algorithms such as edge detection and threshold segmentation to obtain crack edge or texture features to achieve the purpose of crack identification [13]. However, when automatically identifying landslide deformation cracks under complex and changeable background conditions, the use of only image data containing 2D color information is particularly weak [14]. The main defect of these data involves the difficulty in solving the problem of distinguishing surface cracks from background noise (such as shrubs, rock accumulations, and light and shadow differences) [15], resulting in low model recognition accuracies. However, few studies on crack identification and extraction methods based on 3D morphological features have been conducted worldwide [16,17]. Among the existing studies, the main representative method is the use of point cloud data obtained by UAVs, as proposed by Al-Rawabdeh et al. [16], to detect landslide scarps. The study obtained results for bare surface data without vegetation, and no noise-filtering methods were performed with regard to shrubs, boulders, or tree stumps. In addition to traditional target detection methods, using deep learning methods such as the Convolutional Neural Network (CNN) to identify cracks is also an important means, and the recognition accuracy is higher in a certain type of specific scene [18]. However, due to the variability in landslide site environments and the expense of data collection, the model may suffer from underfitting due to insufficient sample size [6,19]. In summary, most current studies focus on crack feature detection in 2D digital images, ignoring the 3D morphological changes caused by different landslide crack scale effects. At the same time, there is no specific filtering algorithm for noise points such as shrubs and shadows in the scene of landslide cracks, which ultimately makes it difficult to comprehensively detect landslide cracks.

In this research, the main purpose is to use point clouds and digital images obtained by UAV-based photography and to achieve automatic identification of landslide cracks under small-scale sample training through image processing, filtering denoising, data fusion, etc. It is necessary to consider the crack image characteristics under different scale conditions and the coupling combination between different characterization indicators. The main innovation is to propose a set of automatic landslide crack detection workflows in which multidimensional remote sensing information is fused, taking into full consideration the recognition effect of landslide cracks under different dimensional data. In this process, a noise filtering method based on crack direction, frequency, and length is proposed to improve the accuracy of crack recognition. The rest of this article is organized as follows: Section 2 describes the basic situation of the study area, the data collection process, and the authenticity of the on-site investigation. Section 3 contains the technical methodological principles of this research, including initial crack extraction, noise filtering, and data fusion. Section 4 provides the crack identification results derived from the above methods in the demonstration area and discusses the different models. Section 5 presents the main conclusions of this research.

2. Description of the Datasets

2.1. Study Area

The mountainous areas in Southwest China feature well-developed valleys and multi-directional rivers, making them key regions for hydropower construction. However, the rise in the upstream water storage level has caused significant water pressure and the erosion of reservoir bank slopes, leading to the rapid development of crack deformation [9]. We selected WuLiPo as the study site, where deformation is occurring on the riverbank of the BaiHeTan Hydropower Station, the second-largest hydropower station in China (Figure 1). WuLiPo spans 300 m in length and 330 m in width, covering an area of approximately 7 × 10⁴ m². The slope is located on the east side of the Jinsha River in Qiaojia County, Shaotong City, Yunnan Province (Figure 1a). Since this reservoir area was impounded in mid-April 2021, the front-edge slopes closest to the riverbank have deformed and loosened (Figure 1b), gradually forming a chair-shaped topographical structure (Figure 1c). This has resulted in the formation of numerous scarps and tension cracks on the slope surface and highway pavement (now abandoned), which continue to expand over time.

2.2. Data Acquisition and Description

UAV photogrammetry is important for landslide investigation due to its high resolution, large-scale imaging ability, rapid mapping ability, and inexpensive cost [5,20,21]. Compared to airborne light detection and ranging (LiDAR), the data acquisition cost is greatly reduced when UAV photogrammetry is employed. Additionally, it offers greater flexibility in alpine and canyon areas with significant relief. However, its main limitation is the difficulty in obtaining terrain information in densely vegetated areas. With the introduction of SFM (Structure From Motion), multiple overlapping aerial images can capture the 3D features of the target through feature matching and 3D reconstruction [22], which provides an efficient means for obtaining landslide images and terrain information in low-vegetation mountainous areas. There are only a small number of shrubs or branches in the study area, which do not cover the cracks. Therefore, we adopted UAV aerial photography and SFM to complete the image data acquisition and processing. In this research, we used the PIX4D mapper (from URL (accessed on 28 October 2024): https://www.pix4d.com/) software to implement the specific work of aerial image processing.

The FeiMa D200 multirotor drone was selected for the field flight (Figure 2b). The UAV hardware parameters and flight information are listed in

Table 1. UAV hardware parameters and flight information.

Hardware Parameters of the Drone		Field Flight Configuration
Drone model	FeiMa D200 Drone	Ground spatial resolution	1.5 cm
Load model	SONY ILCE-6000	Overlap rate between flight bands	80%
Sensor size	23.5 × 15.6 mm	Course overlap rate	85%
Effective Pixels	6000 × 4000	Aerial photography altitude	80 m
Lens parameters	20 mm (focus)	Number of routes	13
Differential mode	PPK/RTK and its fusion mode	Route length	3.9 km
Control radius	5 km	Aerial photography area	0.07 km2

. The terrain-following flight route was planned through preliminary terrain detection (Figure 2c). The 3D point cloud (with a density of approximately 205/m²) was generated in the WuLiPo deformation area through dense image matching (Figure 2e). Next, a DOM with a ground resolution of 3.3 cm was generated through Triangular Irregular Network (TIN) grid construction and automatic texture mapping (Figure 2d). In the collection process, a total of six control points were selected to spatially correct the image (Figure 2a). The plane error of the image before correction was ±0.02059 m and the elevation error was ±0.01697 m. The modeling accuracy could satisfy the subsequent crack identification and evaluation requirements.

2.3. Description of Ground Truth

As shown in Figure 3a, the landslide occurred on a reverse-inclined rock slope with structural characteristics of upper hard rock and lower soft rock. The upper part consists of hard and thick quartz sandstone and dolomite limestone. According to the on-site investigation, the layered rock mass exposed on the surface was broken due to the cutting of two other groups of fractures affected by valley downcutting and unloading [9]. The landslide body formed a potential sliding surface along the interface between the highly weathered rock mass and the medium-weathered rock mass, and a landslide scarp of nearly 200 m was formed on the rear edge. There are a large number of cracks in the middle of the landslide. To evaluate the accuracy and precision of the method, we conducted an interpretation of crack images and field investigation verification in the Wulipo deformation zone from May to July 2021. This investigation identified a total of 133 large-scale tension cracks, accounting for 6.68% of the study area, and a total of 256 local shear cracks on slope surfaces and highways, accounting for 1.31% of the study area. The survey results are shown in Figure 3b. Most of the slope surface cracks are primarily formed in the Quaternary coarse-grained soil layer, exhibiting localized stepped displacements (Figure 3a). The maximum dislocation width is 8.2 m (Figure 3c-④). A large number of shear cracks on the road surface are mainly distributed in the slip areas on both sides of the landslide, resulting in a steep slope of nearly 4.37 m (Figure 3c-⑥). Finally, we numbered each crack and measured their geometric parameter information.

3. Methodology

The DOM and point cloud obtained through UAV photogrammetry are used as data sources. The method can be mainly divided into the following 5 steps. First, we completed the pre-processing of data to achieve data standardization. Second, we analyzed the multidimensional image features and spatial scale effects of surface cracks and defined the scope and type of crack identification. Third, we adopted 6 automatic landslide crack identification models, including roughness, eigenvalue ratio, slope, gray threshold, edge gradient, and supervised classification models, representing remote sensing data sources with different dimensions. Next, the preliminarily extracted cracks were repaired using morphological methods, and 3 filtering methods (the orientation, frequency, and length methods) were proposed to remove errors caused by vegetation and mountain shadows. Finally, the minimum-risk Bayesian probability decision-making method is proposed to perform decision-level fusion on the results of the 6 models listed above. The technical process of this research is shown in Figure 4.

3.1. Data Pre-Processing

Most researchers tackling image-based crack extraction tasks rely on RGB values as the primary features to distinguish cracks from the background surface. In complex mountainous areas, obtaining optical image data necessitates a large amount of information due to low contrast and the challenging working environment [6]. These conditions result in significant noise, small target images, and substantial influence from light and weather conditions. This research employs image contrast enhancement, filtering denoising, and grayscale image pre-processing methods. The purpose is to minimize the noise caused by ground objects or aerial photography angles and normalize the image data.

The Contrast-Limited Adaptive Histogram Equalization processing method (CLAHE) [23] is selected to enhance the image contrast. The principle of this method involves using a fixed-size matrix to collect samples covering the entire image range, setting a certain pixel frequency threshold, and then evenly distributing the regions in which the threshold is exceeded to each gray level. The block-suppressed contrast histogram equalization method increases the differences between cracks and background surfaces under different light intensities and prevents excessive noise amplification, improving the harmoniousness of the entire image. For the filtering and denoising tasks, we adopt bilateral filtering [24]. This algorithm uses 2 Gaussian filters, space and edge intensity difference, to smooth local noise while ensuring the pixel edge features are unchanged; in addition, further crack feature mining often requires only single-variable or binarized image information, so image grayscale processing is necessary. The purpose of this task is to change the crack pixel values from three channels to one channel. This process not only reduces the feature dimensions but also removes irrelevant variables. The weighted conversion grayscale value method is used to convert the three-channel RGB values into a single-channel grayscale value [25]. The above operation process is shown in Figure 5.

During the collection of point cloud data, noise is often generated, which can affect the experimental results. For example, smooth surfaces of soil or rock layers can cause light reflections, while flying birds or dense branches in the air can impact crack extraction accuracy. Therefore, it is necessary to remove these noises using a local distance mean filtering algorithm [26], a common point cloud filtering method that will not be elaborated upon here.

3.2. Crack Identification Model and Principle

3.2.1. Types of Surface Cracks Observed on Landslides

In the early deformation stage, cracks are mostly caused by slope sliding. The width of the cracks in the slope varies under different accumulated deformations. For example, in the WeiHe Basin of China, the initial width of deformation is about 2–10 mm. As the slip surface advances, this deformation width can generally reach 1000–3000 mm [27]. The manifestations of cracks of different sizes in both 2D and 3D spaces are also quite variable.

In the 2D image, the fresh surfaces of cracks formed in slopes composed of different materials exhibit varying colors, but they all produce the same black shadow under sunlight. The shadow is an important feature for identifying cracks [6]. Cracks with different openings produce varying black shadow ranges under different lighting angles. Specifically, cracks with larger openings exhibit almost no black shadows, even if the sunlight is squint, this feature makes the subsidence areas of these cracks indistinguishable from the surrounding flat surfaces (as shown for cracks No. 1 to 3 in Figure 6a). The above reasons make it difficult to separate these cracks from the background surface.

In the 3D scene, the point cloud density generated by the UAV photogrammetry generally ranges from 20–200/m², and the ground resolution of the DSM ranges from 10–50 cm. A point cloud containing 3D information can effectively reflect surface crack fluctuations and is not affected by the angle of illumination. However, for small deformation cracks, especially those with widths smaller than 10 cm, the identification effectiveness is poor (as shown for cracks No. 4 to 5 in Figure 6b).

Different opening sizes will cause different size effects on the cracks in terms of shape and texture. Landslide cracks can be divided into mode I (opening), mode II (sliding), and mode III (tearing) based on the formation mechanism [3]. In addition, there are also some collapsible or closed cracks on the surface near some landslides [13]. This article mainly considers the recognition effect of different data types on different crack size effects. Therefore, the landslide crack deformation size type is agreed as follows: cracks caused by tensile stress, which include tension cracks, grooves produced by the trailing edges of landslides, and scarps caused by sliding crack deformation, are called “large deformation cracks” (as shown in Figure 6c). The cracks caused by shear stress, expansion, and contraction, which include shear cracks on both sides of the landslide, front-edge bulging and tearing, and collapsible ground cracks, are called “small deformation cracks”.

Due to the scale effects of cracks, it is difficult to comprehensively identify complex slope deformation cracks from any single data source. Considering that the resolution of aerial images captured by drone cameras is approximately 2–5 cm, small deformation cracks can produce noticeable black shadows under most lighting angles, making their identification from ortho-images advantageous. Conversely, point cloud data reflect the 3D steepness information of large deformation cracks. Therefore, we leverage the strengths of both data sources to identify cracks of different sizes.

3.2.2. Crack Extraction Based on Point Cloud

Based on the research foundation provided by Al-Rawabdeh et al. [16], three identification indicators in a point cloud, the roughness, eigenvalue rate, and slope, were obtained to identify large deformation cracks. In terms of data indexing, the K-D tree data structure is widely used for nearest-neighbor searches and approximate nearest-neighbor searches of key data in multidimensional spaces due to its space partitioning ability [28]. The advantage of the K-D tree over the octree is that it is more adaptable, can handle data of any dimension, has higher memory efficiency, especially in low-dimensional cases, usually has better search performance, and performs well when data are unevenly distributed. The K-D tree is used to collect local point cloud (X, Y, and Z) coordinates in a 3D space in this research. This process begins from the first point, collects point cloud data at the specified search radius, and calculates the identification indicators (as shown in Figure 7). The definition of each identification indicator is described below.

The roughness indicator can be defined as the irregularity of the terrain surface. Clear surface roughness differences can be found between undisturbed areas and cracks. Equation (1) can be used to calculate the standard deviation of the Z values of the point cloud under the local search radius to represent the roughness ($r_i$) value of the crack region:

$$r_i=\sqrt{{\displaystyle \frac{\sum {\left(Z-\overline{Z}\right)}^2}{n-1}}}$$

(1)

where $r_i$ is the roughness at point $i$, $\overline{Z}$ is the average height at all points within the search range of point $i$, and $n$ is the number of all point clouds within the local search range.

The eigenvalue ratio can describe whether a point cloud set has directionality. If the spatial distribution of surface cracks shows a local downward dislocation trend, the corresponding point cloud will show subsidence and concavity on the flat slope. Principal Component Analysis (PCA) is used to calculate the eigenvalue ratios of the 3 principal vectors of all points in the local neighborhood. If a crack has strong ductility in a fixed direction, the eigenvalue ${\lambda }_1$ (the largest eigenvalue) will be much larger than ${\lambda }_2 \mathrm{o}\mathrm{r} {\lambda }_3$. At the same time, the ${\lambda }_2$ or ${\lambda }_3$ values will be close to 0. The calculation steps of the local point cloud eigenvalue ratios ($d_i$) are described as follows. First, we calculate the covariance matrix ($CO{V}_{3\times 3}$) using Equation (2). The eigenvalues $\left({\lambda }_1,{\lambda }_{2,}, \mathrm{a}\mathrm{n}\mathrm{d} {\lambda }_3\right)$ and eigenvectors $\left(\begin{array}{ccc}\overrightarrow{e_1},& \overrightarrow{e_2}, \mathrm{a}\mathrm{n}\mathrm{d}& \overrightarrow{e_3}\end{array}\right)$ can be obtained by matrix transformation using Equation (3) and are arranged in descending order. Equation (4) can then be used to calculate the eigenvalue ratio ($d_i$).

$$CO{V}_{3\times 3}={\displaystyle \frac{1}{n}}{\sum _{\mathrm{i}=1}^n\left(\left[\begin{array}{c}{P}_{\mathrm{iX}}\\ P_{iY}\\ P_{iZ}\end{array}\right]-\left[\begin{array}{c}\overline{P_{\mathrm{cX}}}\\ \overline{P_{cY}}\\ \overline{P_{cZ}}\end{array}\right]\right)\left(\left[\begin{array}{c}{P}_{\mathrm{iX}}\\ P_{iY}\\ P_{iZ}\end{array}\right]-\left[\begin{array}{c}\overline{P_{\mathrm{cX}}}\\ \overline{P_{cY}}\\ \overline{P_{cZ}}\end{array}\right]\right)}^T$$

(2)

$$CO{V}_{3\times 3}=W\Lambda W^T=\left[\begin{array}{ccc}\overrightarrow{e_1}& \overrightarrow{e_2}& \overrightarrow{e_3}\end{array}\right]\left[\begin{array}{ccc}{\lambda }_1& 0& 0\\ 0& {\lambda }_2& 0\\ 0& 0& {\lambda }_3\end{array}\right]\left[\begin{array}{c}{\overrightarrow{e_1}}^T\\ {\overrightarrow{e_2}}^T\\ {\overrightarrow{e_3}}^T\end{array}\right]$$

(3)

$$d_i={\displaystyle \frac{{\lambda }_3}{{\lambda }_2}}$$

(4)

where $P_{\mathrm{iX}}$, $P_{iY}$, and $P_{iZ}$ are the $X$, $Y$, and $Z$ coordinate values of point $i$, respectively, $\overline{P_{\mathrm{cX}}}$, $\overline{P_{cY}}$, and $\overline{P_{cZ}}$ are the X, Y, and Z coordinates of the centroid of the point cloud set, respectively, and $CO{V}_{3\times 3}$ is the covariance matrix.

The slope indicator, which refers to the undulation size of the fitting plane in the point cloud neighborhood set, is an important parameter in slope stability analyses. A sudden change in the gradient between two consecutive steep slopes means that a scarp is present on the slope [29]. The slope angle at each point can be estimated using Equation (5). The fitting plane utilized herein is the coplanar components of $\overrightarrow{e_1}$ and $\overrightarrow{e_2}$, as processed by PCA in the previous step. The 3 eigenvectors are perpendicular to each other, and the slope (${\theta }_i$) can be calculated using Equation (5):

$${\theta }_i=\left|{\mathit{tan}}^{-1}{\displaystyle \frac{\sqrt{N_x^2+N_y^2}}{N_Z}}\ast {\displaystyle \frac{180}{\pi }}\right|$$

(5)

where ${\theta }_i$ is the inclination angle of cloud characterizing point i and $Nx$, $Ny$, and $Nz$ are the components of the smallest eigenvector ($\overrightarrow{e_3}$) about the X-axis, Y-axis, and Z-axis, respectively, in the local point cloud set.

3.2.3. Crack Extraction Based on Digital Images

According to the expression of pixels, digital images can be divided into color images and binary images, and the crack identification methods applied to different pixel units can be divided into three categories: (1) grayscale image threshold segmentation, (2) edge detection methods, and (3) crack supervised classification based on machine learning [30].

The gray threshold segmentation is the simplest crack identification method. Cracks often exist as black pixels in images. Following the grayscale processing task, a binary image can be segmented using a certain gray value threshold to achieve crack extraction (as shown in Figure 8a,b).

Edge detection is the process of judging whether a pixel belongs to the edge of an object by calculating the size of the gray gradient. Edge detection can be used to effectively detect the outline of an object. Therefore, in the field of computer image processing, various edge detection operators have been developed and are often used to identify cracks or abnormal objects [21]. In this research, the Sobel edge detection operator is used to extract surface cracks. The Sobel operator is a first-order derivative function, and its main principle is an edge-enhancement method that causes the extreme value to appear at the edge by weighting the adjacent gray levels in the local window [31]. In this process, two convolution kernels with a size of 3 × 3 in the horizontal and vertical directions are used to perform convolution calculations on each binary pixel; Then, Equations (6)–(8) are used to calculate the pixel gradient value ($\nabla f$) to represent the size of the edged feature. When the $\nabla f$ value is relatively large, the pixel value changes relatively fast and is biased toward the edge of the contour (as shown in Figure 8c,d):

$$g_x={\displaystyle \frac{\partial f}{\partial x}}=(P7+2P8+P9)-\left(P1+2P2+P3\right)$$

(6)

$$g_y={\displaystyle \frac{\partial f}{\partial x}}=\left(P3+2P6+P9\right)-\left(P1+2P2+P7\right)$$

(7)

$$\nabla f=\sqrt{{g_x}^2+{g_y}^2}$$

(8)

where P1–P9 represent the local-window pixel values, $g_x$ is the horizontal change gradient, $g_y$ is the vertical change gradient, and ∇ f is the first-order pixel gradient.

The maximum likelihood supervised classification based on RGB values is used to complete the image distinction task according to the image features (cracks and surrounding objects). This method is a statistical method based on Bayesian probability. The principle involves obtaining parameters such as the mean and variance in each category by counting the area of interest to determine the classification function. Then, each pixel in the image is substituted into the classification function of each category, and the category exhibiting the largest return value of the function is used as the attribution category of the scanned pixel [32]. Figure 9 provides the spatial distributions of different ground object pixels following the application of the maximum likelihood estimation supervised classification.

3.3. Crack Pixel Repair and Filtering

Following the preliminary extraction process, crack images may still have some noise. Thus, a small portion of the cracks will inevitably be misidentified. We propose 3 filtering algorithms based on the orientation, length, and frequency of cracks. These algorithms can filter out the noises outside the crack area and retain the real pixels unchanged. Notably, to facilitate matrix operations, the data utilized for the crack pixel repair are the binarized raster data following the initial extraction process from the various models described above (a crack is given a value of 1, while the background is given a value of 0). In particular, the extracted crack point cloud is also interpolated into raster data with the same dimensions as the image pixels.

• Morphological repair

In the centers of the initially extracted crack pixel regions, a small portion of the real crack values may be misidentified, especially when using the pixel gradient model. This model can only identify crack edge features, leading to the overlook of a large number of positive samples. Therefore, the morphological closure operation method is used to repair the pixel [33]. The principle of this method involves applying an expansion morphological operation followed by an erosion operation on the extracted crack binary grid. This process fills gaps in the foreground object while maintaining the overall position and shape of the object. The main steps of this method are described as follows: (1) A fixed-size rectangular sampling frame is used to traverse each binarized pixel (as shown in Figure 10a,b). (2) Expanding pixels are identified by calculating the “|” (OR) operation of all the pixels in the sampling frame with each pixel as the center (that is, when there is at least one true value in the sample, the true value is assigned to the center point) (as shown in Figure 10c). (3) Corroding pixels are identified by performing the “&” (AND) operation on the pixels in the expanded sampling frame (that is, when there is at least one false value in the sample, a false value is assigned to the center point) (as shown in Figure 10d,e).

2.

Orientation filtering

Cracks have a certain directional ductility, meaning they extend in thin branches along one or more orientations. This characteristic is a key feature that distinguishes cracks from background noise (hill shadows, bushes, etc.). The ability of PCA to evaluate crack-directionality ductility has been discussed above. Using this feature, each real crack pixel is examined with a fixed-size rectangular frame. The Eigenvalue ratios (${\mathsf{\lambda }}_1/{\mathsf{\lambda }}_2\left({\mathsf{\lambda }}_1>{\mathsf{\lambda }}_2\right)$) of the truth pixels can be calculated as the orientation index of each center pixel, and a threshold can be set according to the image crack feature to conduct refiltering to remove pixel points that do not exhibit strong directionality (as shown in Figure 11a,b). Additionally, the real crack pixels are partially filtered if a reduction in the orientation index at the branch junction occurs. If this happens, the morphological repair can be used to perform the repair process once again.

3.

Length and frequency filtering

The excessive use of the local sampling frame filtering method causes distortion in small cracks, especially at the ends of small cracks. Therefore, refined crack-filtering methods should be based on analyzing each crack as a whole to distinguish the crack from background noise. The filtering steps applied in this work are described as follows:

First, we use the Density-Based Spatial Clustering of Applications with Noise (DBSCAN) algorithm [34]. We cluster pixels according to the following parameter requirements: Eps =$\sqrt{2}$ (represents the neighborhood radius when making density accessibility judgments, which is fixed to ensure that the diagonal distance between two adjacent pixels in a crack is not exceeded),and min_samples = 2 (represents the minimum number of samples in a cluster, which is fixed to ensure that each pixel of a crack has at least two adjacent crack pixels within the cluster search radius). These parameter values are set to ensure that adjacent pixels with at least one common endpoint are clustered; then, all densely concentrated binary pixels are classified into a single crack (as shown in Figure 11f). Second, taking each clustered crack as an index unit, the pixel frequency and the minimum circumscribed rectangle diagonal length of each crack are counted and used as the filter indicators of the frequency and length, respectively. Finally, we set a certain refiltering threshold according to the derived distribution characteristics (as shown in Figure 11g,h). The overall analysis of the fitted crack pixels is performed using the DBSCAN clustering index method to prevent the risk of small-crack end regions being filtered out. This method addresses the disadvantages associated with filtering methods such as the mean filtering method, in which a fixed convolution kernel size is used as a window (as shown in Figure 11b–d).

3.4. Multidimensional Information Fusion

As mentioned above, the types and shapes of cracks identified from different dimensions can vary significantly. Therefore, it is necessary to use the mutual compensation abilities of various algorithms to achieve the optimal recognition effect. This process is known as multidimensional information fusion.

At present, image information fusion levels can be mainly divided into three levels: the pixel level, feature level, and decision level [35]. Decision-level information fusion is the most reliable information fusion method. The decision-level information fusion model can effectively overcome the influence of the errors caused by the limitations of the individual information sources on the true and false target evaluation results. Meanwhile, this method permits the absence of some information sources, thus reducing the possibility of global misjudgments arising due to local errors. In this article, we propose a decision-level Bayesian probability fusion method based on minimum risk.

Due to the differences in the data structures of point clouds and images, the recognition results from different dimensions need to be integrated into the same format. Raster data (the file format is “.tif”) are a data model of geographic spatial information that can be used to describe continuous geographic phenomena, such as terrain, land cover, etc. [36]. Therefore, the point cloud feature indicators (Eigenvalue ratios, Roughness, and Slope) of the same range can be converted into raster format, and then the same number of rows, columns, and corner point coordinates as the digital image (grayscale, edge gradient, and RGB) can be set to complete the data registration. For 3D point cloud data, in this work, we use the TerraScan module in Terrasolid software (from URL (accessed on 28 October 2024): https://terrasolid.com/) to perform point cloud Kriging interpolation and generate grids with equal resolution.

Bayesian decision theory is a basic statistical pattern recognition method in which classifiers are generated based on probabilistic data analyses [37]. This fusion method considers the probability of each feature combination relationship being assigned to each decision category. Therefore, it is not necessary to consider the number of features or the affiliation of the fusion information. In fact, each feature combination only needs to be independent of each other combination, thus allowing the full-feature information fusion of the 6 crack identification models used to be realized.

The Bayesian decision principle involves judging the posterior probability that the input feature pattern belongs to each decision category. However, this method does not take into account the risks introduced by each decision. Therefore, we added the risk coefficient generated by each decision in the final selection process to obtain the optimal crack information fusion results under different risk conditions. The applied implementation steps are described as follows:

(1) The number of input features ($n$) and permutation and combination types are determined to obtain the number of multidimensional information fusion sample features ($x_i,\mathrm{w}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{e} i=1,2,3\dots ,n$).

(2) The output category type is determined by obtaining the final Bayesian decision classification result (${\omega }_j,\mathrm{w}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{e} j=1,2,3\dots ,m$).

(3) According to the statistics obtained for different decision results in various scenarios through existing historical cases, the prior probability distribution ($P\left({\omega }_j\right)$) of the decision result is obtained.

(4) The statistical probability ($P\left(x_i|{\omega }_j\right)$) values of the different feature distributions derived under each decision in the sample function are obtained.

(5) The Bayesian posterior probabilities ($P\left({\omega }_j|x_i\right)$) of different classification results (${\omega }_j$) are obtained under the condition of each feature combination ($x_i$) using Equation (9).

$$P\left({\omega }_j|x_i\right)={\displaystyle \frac{P\left(x_i|{\omega }_j\right)·P\left({\omega }_j\right)}{\sum _{c=1}^{m}P\left(x_i|{\omega }_c\right)·P\left({\omega }_c\right)}}$$

(9)

(6) The risk coefficient (${\lambda }_{\alpha ,\beta }$, where ${\lambda }_{\alpha ,\beta }$) is determined; this coefficient represents the loss that may be caused when the decision classifies the $x$ feature as the ${\omega }_{\alpha }$ class but $x$ actually belongs to the ${\omega }_{\beta }$ class (where $\alpha ,\beta \in \left(1,2,3\dots m\right)$).

(7) Finally, it is necessary only to calculate the risk of $x$ falling into each class according to Equation (10) and select the classification $R\left({\omega }_{\alpha }|x\right)$ with the smallest risk as the $x_i$ fusion decision.

$$R\left({\omega }_{\alpha }|x\right)={\sum }_{\beta =1}^m{\lambda }_{\alpha ,\beta }·P\left({\omega }_{\beta }|x\right)$$

(10)

4. Results and Discussion

4.1. Remote Sensing Data Pre-Processing

To obtain the final 2D image data, we use the OpenCV-python computer vision image processing library (Version: 4.9.0.80) to complete the DOM pre-processing task. This open-source platform provides image-processing algorithm functions, including pixel editing, geometric transformation, image segmentation, etc., which can complete the image data preprocessing in this work. For 3D point cloud data, we use the TerraScan module in Terrasolid software (Version: 10.09.01.01) to denoise the point cloud. The results are shown in Figure 12.

4.2. Initial Extraction

Based on the six crack indicator identification methods proposed in Section 3.2, the point cloud and orthophoto image collected in the demonstration area are processed as follows:

(1) To better highlight terrain changes, and considering that a smaller radius increases sensitivity to noise such as isolated stones, deadwood, and branches, we choose a neighborhood search radius of 1 m. The search radius is adjustable and generally roughly the same as the average width of the cracks to be identified. To obtain the roughness, eigenvalue ratio, and slope, Equations (1), (4) and (5) are used (the results are shown in Figure 13a–c). To facilitate the subsequent matrix operations, the point cloud is converted into a raster with the same unit pixel size as the DOM (the results are shown in Figure 13d–f).

(2) The pre-processed grayscale image is shown in Figure 13h. The OpenCV-python image processing library is used to complete the calculations described by Equations (6)–(8), and the calculated gradient value is scaled to an 8-bit image depth (the results are shown in Figure 13g). For the RGB images, the maximum likelihood method toolbox of the ArcGIS software (Version: 10.3) is used. By combining the RGB images with the characteristics of the WuLiPo image, the pixels are divided into 8 categories and 2000 pixel points are collected manually for each sample category. The sampling positions and prediction results by the maximum likelihood method are shown in Figure 14.

In the six crack identification models listed above, except for the supervised learning classification, all methods are pixel-oriented crack data binarization processing methods. Therefore, the results only need to pass a threshold value to complete the two-classification task to obtain crack pixels (i.e., crack and non-crack pixels). The Receiver Operating Curve (ROC) test can be used to analyze the predictive ability of a classifier under different thresholds [38]. The vertical coordinate of the curve is the True Positive Ratio (TPR), that is, the proportion of positive samples correctly predicted by the model as positive cases under different thresholds: TPR = TP/(TP + FN). The abscissa is obtained by calculating the False Positive Ratio (FPR), that is, the proportion of false positives in all inverse samples derived under different thresholds: FPR = FP/(FP + TN). The meanings of the abbreviations TN, FP, FN, and TP are shown in Figure 15.

In this research, 100,000 pixels representing both positive and negative samples are randomly selected from the test set of large and small deformation cracks, and the 5 binary classification models described above are applied to draw ROC curves (the result is shown in Figure 16).

Figure 16 shows that when identifying large deformation cracks, the AUCs of the slope and roughness classification models are both above 0.8, indicating a good identification effect. The classification effect of the eigenvalue ratio model is poor. The main reason for this result is that this model is highly sensitive to uneven unidirectional changes. Consequently, it has a strong recognition effect on unclassified boulder blocks and low vegetation, resulting in its relatively low recognition accuracy. In the recognition of small deformation cracks, the grayscale segmentation model exhibited the best recognition effect, with an AUC value reaching 0.83. The pixel gradient model has a poor predictive ability, mainly because the edge gradient has an optimal recognition effect only on the boundary crack pixels, while the recognition effect on the black pixels in the middle of the crack and the surrounding regions is not optimal.

On a ROC curve, the closer a point is to (0, 1), the better the segmentation threshold corresponding to that point is, meaning that while ensuring a high positive rate, the false positive rate should be minimized to the greatest possible extent. Therefore, the point with the smallest distance ($\sqrt{x^2+y^2}$) from (0, 1) on the curve is the optimal segmentation threshold of the corresponding model. Figure 16 shows that when the best segmentation threshold is obtained for various models, the TPR is lower than 80%. To better meet actual engineering needs, it is necessary to reduce the certain threshold standard for the segmentation task. Finally, the corresponding threshold of each model under the true case rate of 0.8 is used as the segmentation threshold, which is set following the adjustment process to complete the semantic segmentation of cracks. The extraction results are shown in Figure 17.

4.3. Repair and Filtering

Figure 17 shows that, with an increase in the TPR, the noise of the segmented false cracks gradually increases, and the overall recognition precision of the model begins to decline. Therefore, using the pixel-repairing and filtering method proposed in Section 3.3, the morphological repair is conducted using a 5 × 5 Sampling range for the operation. Filter processing utilizes three indicators, including the eigenvalue ratio (${\mathsf{\lambda }}_1/{\mathsf{\lambda }}_2$), crack frequency, and length, where the filtering threshold is 20% of the maximum value of each indicator. The results are shown in Figure 18.

To evaluate the accuracy of each model after the repair filtering step, the results extracted from each model are compared with the test sets derived from manual field checking. The precision, TPR, and FPR values of the crack identification results are counted based on 500,000 pixels in a typical scene area in WuLiPo (the number of positive and negative samples is set to be the same). The F-score is applied to comprehensively characterize the precision and TPR values. This metric can be used to evaluate the overall characteristics of each model before and after repair filtering. The above precision and F-score calculation methods are shown in Equations (11) and (12). The evaluation results are shown in

Table 2. Two-class evaluation results of crack identification model before and after repair filtering.

Category	Method		TPR	FPR	Precision	F-Score
Small deformation cracks	Grayscale value segmentation model	Initial	81.32%	32.54%	71.42%	0.760
		After filtering	71.23%	21.5%	76.81%	0.739
	Supervised classification model	Initial	79.67%	29.64%	72.88%	0.761
		After filtering	72.31%	17.34%	80.66%	0.763
	Gradient value segmentation model	Initial	79.84%	59.63%	57.24%	0.667
		After filtering	82.56%	10.32%	88.89%	0.856
Large deformation crack	Eigenvalue Ratios segmentation model	Initial	77.34%	48.32%	61.54%	0.685
		After filtering	73.54%	34.86%	67.84%	0.706
	Roughness segmentation model	Initial	82.45%	31.98%	72.05%	0.769
		After filtering	76.35%	18.86%	80.19%	0.782
	Slope segmentation model	Initial	78.46%	27.64%	73.95%	0.761
		After filtering	77.36%	13.56%	85.09%	0.810

and Figure 19.

$$Precision={\displaystyle \frac{TP}{TP+FP}}$$

(11)

$$\mathrm{F}\text{-}\mathrm{Score}=2·{\displaystyle \frac{Precision·TPR}{Precision+TPR}}$$

(12)

The filtering algorithm proposed in this article reduces the FPR of each model by 27.9–82.7%, thus greatly reducing the crack identification error rate. Inevitably, due to the complexity of the background, a small number of real cracks are still filtered out. A small decrease in the TPR can be observed, of which the value ranges between 1.4% and 12.4%, but the accuracy rate is increased by 7.5% to 55.3%. Especially for the edge detection algorithm model, the morphological repair can not only “fill” the positive samples inside the edge contour but the filtering algorithm based on the crack feature can also reduce 82.7% of the noise caused by the boundary shadows of trees and rock blocks. Moreover, the precision metric has increased by 55.3%. Table 2 shows that the F-score of the model exhibited an improvement following the filtering process.

Among the identification methods of small and large deformation cracks, gradient value segmentation and slope segmentation models have the best effect. The crack results extracted by the two models are superimposed and compared with the manual investigation results (as shown in Figure 20).

4.4. Multidimensional Information Fusion Results

4.4.1. The Result of Fusion

Six different models are selected for deformation cracks, and these models are all independent of each other. Therefore, the six studied models can be used together as feature samples. According to the fusion method described in Section 3.4, the detailed decision-level fusion steps are described below.

(1) The decision-level fusion process considers the output of each model, so the judgment results of the six models for a single pixel can have values of only 0 (non-crack) or 1 (crack). The permutation and combination relationships of the six models can be used to determine the total number (N) of all feature combinations: $N=2\times 2\times 2\times 2\times 2\times 2=64$.

The number “2” in the calculation of $N$ above represents the possible value of each model result, that is, the sample feature number $x=\left\{x_i\left|x_i=\left[A,B,C,D,E,F\right],A~F\in \left(\mathrm{0,1}\right)\right.\right\}$, where A~F refers to the roughness, slope, eigenvalue ratio, grayscale segmentation, supervised classification, and edge gradient values output by each model. The results of each feature displayed in the WuLiPo image are shown in Figure 21a.

(2) The prior probability is a preexisting known quantity in Bayesian probability theory. Therefore, according to the investigation of crack development characteristics and historical data in the demonstration area, $P\left({\omega }_1\right)$ = 8.27% and $P\left({\omega }_2\right)$ = 91.72%, where ${\omega }_1$ and ${\omega }_2$ refer to the judgment results of cracks and noncracks, respectively.

(3) Different feature distribution probabilities ($P\left(x_i|{\omega }_j\right)$) are then obtained. When the results of ${\omega }_1$ and ${\omega }_2$ are known, the distribution probabilities $P\left(x_i|{\omega }_1\right)$ and $P\left(x_i|{\omega }_2\right)$ of $x_i$ can be calculated from the training set.

(4) According to Equation (9), the Bayesian posterior probabilities $P\left({\omega }_1|x_i\right)$ and $P\left({\omega }_1|x_i\right)$ of different classification results (${\omega }_j$) under the condition of each feature combination $x_i$ are calculated; the results are shown in Figure 21b.

(5) The decision risk coefficient (${\lambda }_{\alpha ,\beta }$) is then substituted. There are four types of decision risk coefficient combinations considered, that is,${\lambda }_{1,1},{\lambda }_{\mathrm{1,2}},{\lambda }_{2,1},\mathrm{a}\mathrm{n}\mathrm{d} {\lambda }_{2,2}$, among which ${\lambda }_{1,1}$, ${\lambda }_{2,2}$ = 0. To analyze the impacts of different risk factors on the decision-making results, a total of six risk conditions are considered:${\lambda }_{1,1}:{\lambda }_{\mathrm{1,2}}:{\lambda }_{2,1}:{\lambda }_{2,2}$ = $\left(0:1.5:1:0\right)$, $\left(0:1.2:1:0\right)$, $\left(0:1:1:0\right)$, $\left(0:1:2:0\right)$, $\left(0:1:4:0\right)$, and $(0:1:7:0)$.

(6) Finally, through Equation (10), we compare the risk values of $R\left({\omega }_1|x_i\right)$ and $R\left({\omega }_2|x_i\right)$ and perform the decision fusion task (as shown in Figure 22).

4.4.2. Evaluation of the Fusion Results

The same test set as that used in the previous section was used here to calculate the F-scores of the fusion methods. The results are shown in

Table 3. Statistics of the multidimensional data fusion results.

Fusion Method	Type	TPR	FPR	Precision	F-Score
The Layer stacking fusion	Slope + Edge Gradient	80.81%	12.74%	86.38%	0.835
The Bayesian probability fusion method based on minimum risk	${\lambda }_{12}:{\lambda }_{21}=1.5:1$	1.24%	0.13%	90.51%	0.024
	${\lambda }_{12}:{\lambda }_{21}=1.2:1$	61.59%	3.24%	95.03%	0.750
	${\lambda }_{12}:{\lambda }_{21}=1:1$	76.67%	8.97%	89.53%	0.826
	${\lambda }_{12}:{\lambda }_{21}=1:2$	83.73%	10.64%	87.15%	0.864
	${\lambda }_{12}:{\lambda }_{21}=1:4$	92.36%	12.36%	85.37%	0.901
	${\lambda }_{12}:{\lambda }_{21}=1:7$	93.52%	34.32%	71.96%	0.822

. The fused model was more robust than any individual recognition model. In addition, it can be seen from the table that the Bayesian decision-level fusion results are not necessarily better than the layer-stacking fusion results. The various risk coefficients exhibit great differences in the fusion results. It can be seen from Figure 23 that as the value of ${\lambda }_{12}:{\lambda }_{21}$ decreases, TPR begins to increase, indicating that as the risk ratio decreases, the true-positive rate increases. When the risk coefficient decreases to 0.5, the true-positive rate exceeds the optimal effect of the optimal layer superposition fusion (as shown in Figure 23a). As the true-positive rate increases, the false-positive rate begins to increase due to the decrease in the risk coefficient. When the risk coefficient exceeds 0.25, the probability of judging non-cracks as true cracks suddenly increases, causing the accuracy rate to begin to decline (as shown in Figure 23b,c). At this time, the F-score is maximized with a value of 0.901 (as shown in Figure 23d). Therefore, the crack identification effect using Bayesian probability fusion changes with the change in the risk coefficient. When the crack identification accuracy is given priority, the risk coefficient should be appropriately increased; when the crack identification recall rate is given priority, the risk coefficient should be minimized. If both are considered comprehensively, the risk coefficient should be maintained at 0.25 according to the results of this research.

5. Conclusions

In the research, we investigated the key technologies related to landslide cracks based on multidimensional image data, including their automatic identification, extraction, repair, filtering, and fusion. The major conclusions are as follows:

(1) A methodological system for automatically identifying cracks and collecting crack information based on fusing multidimensional information was constructed. The proposed system addressed the shortcomings of the currently available individual remote sensing data identification technologies and could be applied to identify and investigate large-scale, multidimensional complex landslide deformation cracks. In the pixel unit-oriented automatic identification algorithm, starting from the crack image features of different scales, six different crack identification models were used and tested using the calculated ROC values. The derived AUC values were all between 0.6 and 0.85.

(2) Due to the influence of light angles, shrubs, boulders, and branches, there are often noise points in the identification of cracks, which can mislead people’s judgment of the instability characteristics of the landslide. The morphological repairing and filtering algorithm proposed in this research could effectively reduce the background noise of cracks and improve the overall accuracy of the applied models. Based on morphological repair and orientation, frequency, and length filtering tasks, which could sufficiently remove complex noise points (thus reducing the FPR value) and minimize the loss of real cracks (thus reducing the TPR as much as possible), the precision of each model increased by 7.5–55.3%, making the crack identification results more realistic and reliable.

(3) We proposed Bayesian probability decision-making with minimum risk methods to significantly improve the crack identification ability of any individual model. The highest F-score value obtained for the studied model following this fusion process was 0.901, and the risk coefficient could be adjusted according to the actual needs of a project to obtain improved and more practical results.

The investigation and statistical analysis of landslide cracks conducted in past work were often time-consuming and laborious, especially the crack information measurement work conducted in obviously deformed and damaged highly dangerous areas. The use of remote sensing technology to generate multidimensional data that can be used to automatically identify cracks can achieve an efficient operation method. These technologies can be used to solve the difficulties and bottlenecks encountered in traditional geological disaster investigations and are critical for further promoting the early identification of geological disasters and for ensuring the intelligence and automation of monitoring and early warning systems. As mentioned in the introduction, the environmental conditions of the site are important factors that determine whether the model can accurately identify landslide cracks. The crack identification process proposed in this research has multiple threshold and filtering methods to choose from. These choices are not absolute. Therefore, we systematically summarize the principles and ranges of threshold selection in the method section. However, there may still be problems such as noise filtering transitions and difficulty in processing large datasets. Therefore, determining how to perform adaptive parameter selection and large data cluster processing algorithms may be the focus of research in the next step. In addition, considering the data costs and sparse vegetation in the field, in this research, we adopted the UAV oblique photography method to obtain the target point cloud. However, in densely vegetated fields, the use of airborne LiDAR would effectively “penetrate” the vegetation layer to obtain real surface data, and this could improve the crack identification accuracy [39]. Therefore, if conditions permit, the use of LiDAR to obtain point cloud data is also a good choice.