A morphological regularization-constrained fuzzy C-means clustering algorithm for tunnel rock discontinuity identification

1 Introduction

Rock discontinuities refer to various geological interfaces that separate rock masses into discrete blocks, including joints, fractures, bedding planes, and faults. Their development characteristics directly determine the mechanical properties and engineering response of rock masses (Zhang et al., 2016). The spatial distribution, orientation combinations, and physico-mechanical properties of discontinuities are key factors controlling geological hazards such as slope instability, tunnel collapse, and rockburst (Berisavljević et al., 2023; Li et al., 2023). As weak zones, potential sliding surfaces, and seepage pathways, discontinuities significantly affect the strength, deformation, and seepage behavior of rock masses, thereby controlling their macroscopic stability (Cook, 1992; Zhu et al., 2023; Zhang et al., 2026; Du et al., 2026). Therefore, accurate acquisition of the geometric information of rock discontinuities, especially orientation, spacing, and persistence, is fundamental for rock mass quality evaluation, stability analysis, and engineering design (Cao et al., 2025; Sun et al., 2026; Chen et al., 2024a; Zhang et al., 2025; Zhang et al., 2024). Traditional discontinuity mapping methods mainly rely on manual measurements using geological compasses and measuring tapes. These methods are labor-intensive, time-consuming, and difficult to apply to steep slopes or hazardous areas, where accessibility and operational safety are major concerns (Assali et al., 2014; Zhang et al., 2018). With the rapid development of surveying technologies, non-contact data acquisition techniques, such as three-dimensional laser scanning (LiDAR) and UAV photogrammetry, have been increasingly applied to rock discontinuity investigation (Abellán et al., 2006; Niethammer et al., 2012; Li et al., 2016; Deng et al., 2026; Ma et al., 2026). These techniques can rapidly acquire high-precision point cloud data of rock mass surfaces, providing rich data sources for the digital characterization of discontinuities. However, automatic and accurate identification and extraction of discontinuity information from massive and scattered point cloud data remain key challenges in current research.

Extensive studies have been conducted on the intelligent identification of rock discontinuities from point cloud data, resulting in multiple technical routes. Automatic identification methods based on clustering algorithms are currently among the mainstream approaches. Commonly used algorithms include K-means clustering (Chen et al., 2016; Zhou et al., 2024; Peng and Wang, 2025), density-based spatial clustering of applications with noise (DBSCAN) (Ren et al., 2025; Zhang et al., 2022), and fuzzy C-means (FCM) clustering. Among them, FCM has been widely applied to rock discontinuity identification because it introduces membership degrees and can describe the fuzzy characteristics of discontinuity boundaries. However, traditional FCM still has several significant limitations in rock discontinuity identification. It is sensitive to randomly selected initial cluster centers, requires the number of clusters to be predefined, easily converges to local optima, and is susceptible to outliers, resulting in unstable clustering results (Liu et al., 2025; Wang et al., 2025; Pu et al., 2025). In addition, when processing large-scale datasets, traditional FCM can also suffer from local optimum problems, which further affects clustering stability (Song et al., 2017; Cui and Yan, 2020). To overcome these limitations, various improved FCM strategies have been proposed. Li et al., (2025) combined an improved artificial bee colony (ABC) algorithm with FCM and used the global search capability of ABC to alleviate the local optimum problem, automatically determine initial cluster centers and membership degrees, and improve accuracy and stability. Wu et al. (2025) proposed an optimized algorithm combining density peak clustering and FCM (DPC-FCM), in which the DPC algorithm is used to determine global density peaks as the initial cluster centers of FCM, automatically determine the number of clusters, and remove abnormal points, thereby avoiding local optima and improving robustness. Song et al. (2017) improved FCM using quantum particle swarm optimization (QPSO), achieving global optimization and reducing sensitivity to initial cluster centers through particle swarm initialization, fitness calculation, and optimal-position updating. Gao et al., (2024) proposed an adaptive-neighborhood-information-based FCM algorithm, in which neighborhood information between samples and cluster centers is introduced to enhance the algorithm’s ability to perceive data structures and improve clustering stability.

Rock discontinuities not only exhibit similar orientation characteristics but also possess distinct spatial geometric characteristics. Point clouds belonging to the same dominant discontinuity set generally exhibit good local coplanarity, whereas the local geometric characteristics become more complex near discontinuity boundaries, intersection zones, and fractured regions (Battulwar et al., 2021; Chen et al., 2016; Kong et al., 2020). At discontinuity boundaries, point clouds from adjacent discontinuities may be mixed within the same local neighborhood, and clustering based solely on normal-vector orientation can cause boundary points to be absorbed into neighboring dominant sets, resulting in boundary expansion or the merging of fine strip-like discontinuities. In discontinuity-intersection zones, a single spatial neighborhood may simultaneously contain points from discontinuities with different orientations, causing the estimated local normal vectors to be jointly affected by multiple directional components and making it difficult for orientation-based clustering to maintain clear class boundaries. In fractured regions, increased point-cloud dispersion and surface undulation reduce the stability of local normal vectors and may also lead to scattered misclassification and cluster-center drift (Guo et al., 2017; Gu et al., 2024; Chen et al., 2024b; Zhou et al., 2026).

Existing improvements to FCM mainly enhance clustering stability through initial cluster-center selection, global search, cluster-number determination, and the use of neighborhood information. ABC-FCM and QPSO-FCM improve the initialization and local-optimum problems through swarm-intelligence optimization, DPC-FCM uses density peaks to determine initial cluster centers and assist in determining the number of clusters, and adaptive-neighborhood FCM enhances the perception of data-distribution structures through neighborhood information. These methods improve the initialization, search capability, or data-structure representation of FCM. However, in rock discontinuity identification, local planar geometric characteristics have not been directly incorporated into the membership calculation of the FCM objective function. For boundary-mixing zones, discontinuity intersections, and locally fractured regions, orientation, density, or neighborhood-distribution information alone is still insufficient to simultaneously describe the directional characteristics of a point and its local degree of coplanarity.

Based on the above considerations, this study proposes a morphological regularization-constrained FCM clustering algorithm. Unlike methods that mainly improve cluster initialization or the search process, the proposed method directly introduces the local planar residual into the FCM objective function as a morphological constraint, so that dominant discontinuity-set classification simultaneously considers normal-vector directional similarity and local spatial coplanarity. Points located within stable discontinuity surfaces generally have relatively small local planar residuals, whereas points with poor local geometric consistency in boundary, intersection, and fractured regions exhibit larger residuals and therefore incur higher assignment costs during membership calculation, thereby reducing their influence on cluster centers and class boundaries. On this basis, initial over-segmentation, geometric-quality screening, and similar-cluster merging are further employed to obtain the final dominant discontinuity sets, after which a plane-constrained improved region growing algorithm is used to extract individual discontinuities.

2 Data source and methodology

This section describes the proposed intelligent identification method for tunnel rock discontinuities based on the morphological regularization-constrained FCM clustering algorithm. The workflow consists of point cloud acquisition and preprocessing, classification of dominant discontinuity sets using MRC-FCM, cluster validity screening and merging, and individual discontinuity segmentation based on plane-constrained region growing. The overall strategy is to first identify dominant discontinuity sets and then extract individual discontinuities within each set.

2.1 Data acquisition and preprocessing

In this study, a tunnel face and its adjacent sidewalls are taken as the research object, and three-dimensional laser scanning is used to acquire point cloud data of the surrounding rock surface. Considering the large tunnel cross-section, pronounced local surface undulations, and potential safety risks associated with close-range manual measurements, three-dimensional laser scanning is more suitable for this study because it can rapidly acquire high-precision spatial geometric information under non-contact conditions.

During field data acquisition, multiple scanning stations were arranged to scan the tunnel face and adjacent sidewalls. Inter-station registration and point cloud stitching were then performed during post-processing to establish a 3D point cloud model in a unified coordinate system. The final raw point cloud contained approximately 5.40 million points, which can comprehensively reflect the distribution characteristics and spatial morphology of discontinuities on the surrounding rock surface of the tunnel face. The 3D point cloud model is shown in Figure 1. The region enclosed by the blue box is the clustering object selected in this study.

2.2 Morphological regularization-constrained FCM algorithm

Fuzzy C-means (FCM) is a classical clustering algorithm based on soft partitioning. Its basic idea is to introduce the concept of membership degree, allowing each sample to belong to multiple clusters with different membership values. Therefore, FCM is suitable for data with fuzzy boundaries. Let the point cloud normal-vector dataset be X = {X₁, X₂, …, X_n}, where Xi denotes the unit normal vector of the ith point. Let the cluster center set be V = {V₁, V₂, …, V_k}, and the membership matrix U = [u_ji]k × n satisfies the constraint given in Equation 1:

∑ j = 1 k u j i = 1 , u j i ∈ 0 , 1 (1)

The objective function is expressed in Equation 2:where m is the fuzzification exponent, and m = 2 is used in this study. Considering the directional equivalence of normal vectors of rock discontinuities, namely, that opposite normal vectors may correspond to the same discontinuity set, the sine-squared value of the angle between two unit normal vectors is adopted as the similarity distance instead of the conventional Euclidean distance, as given in Equation 3:

J U , V = ∑ i = 1 n ∑ j = 1 k u j i m d 2 X i , V j (2)

d 2 X i , V j = 1 − X i T · V j 2 (3)

The iterative update equations for the membership degree and cluster centers are given in Equations 4, 5, respectively:

u j i = ∑ l = 1 k d 2 X i , V j d 2 X i , V l 1 m − 1 − 1 (4)

V j = ∑ i = 1 n u j i m X i ∑ i = 1 n u j i m X i (5)

This model can describe the fuzzy boundaries involved in discontinuity classification; however, the clustering distance is still mainly determined by differences in normal-vector orientation and does not directly represent the local spatial coplanarity of the point cloud. Therefore, in discontinuity-boundary, intersection, and locally fractured regions, points with similar orientations but different spatial geometric relationships may still be misclassified. To further incorporate local geometric information of discontinuities, plane fitting is performed within the spatial neighborhood of each point, and the local planar residual is used to characterize the degree of local coplanarity of the point cloud.

To quantify the local geometric flatness of the point cloud, least-squares plane fitting is performed for the spatial neighborhood of each point. For an arbitrary 3D coordinate point, the neighboring point set within a radius of R = 0.05 m is selected, and the number of neighboring points is denoted as M. The plane is fitted using the least-squares method, as given in Equation 6:where p ¯ = 1 M ∑ q ∈ N i q, p i ¯ is the centroid of the neighborhood, q is a neighboring point and n i is the unit normal vector of the fitted plane. The distance from point p i to this plane is defined as the local planar residual in Equation 7:

min ∑ q ∈ N i n i T · q − p i ¯ 2 , n i = 1 (6)

r i = n i T · p i − p i ¯ (7)

The smaller the value of r i , the higher the likelihood that point i is located within the interior of a smooth discontinuity surface. Conversely, a larger r i indicates that the point is more likely to be located near a discontinuity boundary, fractured zone, cavity edge, or noisy region. Because this residual is determined solely by the 3D coordinates of the point, it can be precomputed before the clustering process.

On this basis, the local planar residual is introduced into the FCM objective function as a regularization term, and the total objective function of the morphological regularization-constrained FCM algorithm (MRC-FCM) is constructed as given in Equation 8:where λ is the regularization coefficient, which controls the strength of the spatial morphological constraint. When λ = 0, the model degenerates into the traditional FCM algorithm. When λ > 0, points with poor local geometric flatness are assigned a higher membership cost, thereby suppressing their contribution to cluster-center formation. In this study, r i denotes the normalized local planar residual before being introduced into the objective function. With fixed cluster centers, the membership update equation given in Equation 9 can be derived using the Lagrange multiplier method:

J M R C U , V = ∑ i = 1 n ∑ j = 1 k u j i m d 2 X i , V j + λ r i (8)

u j i = ∑ l = 1 k d 2 X i , V j + λ r i d 2 X i , V l + λ r i 1 m − 1 − 1 (9)

The cluster-center update equation is given in Equation 10, in which a weighted average is calculated and then normalized:

V j = ∑ i = 1 n u j i m X i ∑ i = 1 n u j i m X i (10)

It should be noted that MRC-FCM explicitly suppresses the influence of non-coplanar points during soft partitioning, making the clustering within discontinuity surfaces more compact and the boundary transitions more concentrated. Traditional FCM requires the number of clusters k to be specified in advance, which is often difficult in practical engineering applications. Considering that the number of dominant discontinuity sets in tunnel surrounding rock is usually limited, the clustering procedure begins with an initial over-segmentation strategy. Specifically, Kmax is set to be larger than the possible true number of discontinuity sets. After one MRC-FCM iterative process reaches convergence, the clusters are screened according to their geometric quality. The weighted average local planar residual of the jth cluster is defined in Equation 11:where r j ¯ is the weighted average local planar residual of the jth cluster. When r j ¯ > γ r , with γ r being an empirical threshold, the cluster is regarded as being dominated by boundary points, transition points, or noise points and is removed. When the angle between the normal vectors of two cluster centers, θ jl = cos − 1 V j T · V l, is smaller than the prescribed threshold, which is set to approximately 10° in this study, the two clusters are treated as the same dominant discontinuity set and are subsequently merged.

r j ¯ = ∑ i = 1 n u j i m r i ∑ i = 1 n u j i m (11)

Through the three-step procedure of over-segmentation, pruning, and merging, invalid clusters with poor geometric quality can be removed while preserving the complete identification capability of discontinuity sets. Consequently, more stable dominant discontinuity-set classification results that are more consistent with engineering geological significance can be obtained. The overall workflow of the algorithm is shown in Figure 2.

To evaluate the effectiveness of dominant discontinuity-set classification, commonly used fuzzy clustering validity indices are adopted, including the partition coefficient (PC), partition entropy (PCE), Xie-Beni index (XB), and Fukuyama-Sugeno index (FS). These four validity indices are calculated using Equations 12–15: where V ¯ is the mean of all cluster centers, and d 2 V j , V ¯ denotes the distance between the jth cluster center V j and the mean cluster center V ¯ . Generally, a larger PC value indicates better separation among clusters, whereas smaller PCE, XB, and FS values indicate superior clustering performance. Comprehensive comparison of these four indices enables quantitative evaluation of clustering results under different numbers of discontinuity sets, thus helping determine the optimal number of dominant discontinuity sets.

P C = 1 n ∑ i = 1 n ∑ j = 1 k u j i 2 (12)

P C E = − 1 n ∑ i = 1 n ∑ j = 1 k u j i ⁡ ln ⁡ u j i (13)

X B = ∑ i = 1 n ∑ j = 1 k u j i m d 2 X i , V j n ⁡ min ⁡ d 2 V l , V k , l ≠ k (14)

F S = 1 n ∑ i = 1 n ∑ j = 1 k u j i m d 2 X i , V j − d 2 V j , V ¯ (15)

2.3 Plane-constrained improved region growing algorithm

After the identification of dominant discontinuity sets, each set generally still contains multiple spatially separated individual discontinuities, and the clustering results alone cannot adequately represent the differences among individual discontinuities in spatial position and extent. Meanwhile, rock discontinuities generally exhibit approximately planar characteristics at a local scale, but different dominant discontinuity sets may be spatially close in intersection zones. If region growing is performed directly on the entire point cloud, candidate points may cross the actual boundaries between different dominant discontinuity sets, allowing discontinuities with different orientations to participate in the growth of the same region. Therefore, point-cloud subsets corresponding to the dominant discontinuity sets are first extracted according to the MRC-FCM clustering results, and neighborhood relationships are constructed independently within each set so that the search range of candidate points for region growing is restricted to the corresponding dominant discontinuity set. On this basis, a plane-constrained improved region growing algorithm is introduced. The fitted plane of the current growing region is used as an overall reference, and candidate-point inclusion is controlled in terms of directional consistency, spatial coplanarity, and overall planar stability, thereby extending the procedure from dominant discontinuity-set identification to spatial segmentation of individual discontinuities.

For the plane fitted to the current growing region, let n denote its normal vector and x ¯ denote its center. The coordinate of the candidate growing point is denoted by x j, and its normal vector is denoted by n j.

• Directional consistency

Directional consistency is defined in Equation 16 to describe the degree of agreement between the candidate point and the overall orientation of the current fitted plane:where θ j denotes the angle between the normal vector of the candidate growing point and the normal vector of the current fitted plane. A smaller θ j indicates that the point is more consistent with the current discontinuity in orientation.

θ j = cos − 1 n T · n j (16)

• 2. Spatial coplanarity:

The distance from the candidate point to the current fitted plane is defined in Equation 17 and is used to determine whether the candidate point satisfies the coplanarity requirement in spatial position:where ε j denotes the orthogonal distance from the candidate growing point to the current fitted plane. A smaller ε j indicates that the point is closer to the plane of the current region.

ε j = n T · x j − x ¯ (17)

• 3. Overall planar stability:

To ensure the overall planarity of the growing region, the change in fitting error before and after the inclusion of the candidate growing point should also be examined. Assuming that the current region contains k points, the root mean square error of plane fitting for the current region can be expressed as given in Equation 18:

R M S = 1 k ∑ i = 1 k n · x i − x ¯ 2 (18)

After the candidate growing point x j is included in the current region, the updated fitting error is given in Equation 19:

R M S j = k · R M S 2 + ε j 2 k + 1 (19)

Accordingly, the fitting error increment caused by the candidate point is defined in Equation 20:

Δ R M S j = R M S j − R M S (20)

A smaller Δ R M S j indicates that the inclusion of the point causes less disturbance to the overall planarity of the current region. Conversely, a larger Δ R M S j indicates that the point is not suitable for inclusion in the current discontinuity.

The candidate point x j is included in the current planar region only when it simultaneously satisfies the threshold requirements for θ j, ε j and Δ R M S j . Through the above three geometric criteria, candidate-point inclusion is simultaneously constrained by directional differences, local spatial coplanarity, and the overall plane-fitting quality of the growing region. This allows region growing to further account for the local planar geometric relationships of individual discontinuities within the same dominant discontinuity set. The overall workflow of the plane-constrained improved region growing algorithm is shown in Figure 3.

The three criteria constrain candidate-point inclusion at different levels. Normal consistency evaluates the directional compatibility between a candidate point and the current region and restricts points with large directional differences from entering the current discontinuity. The point-to-plane residual further evaluates the local coplanarity of the candidate point with respect to the currently fitted plane and is used to distinguish points with similar orientations but spatial positions that do not belong to the same discontinuity. The RMS residual increment reflects the change in the overall plane-fitting quality after a candidate point is incorporated and controls the stability of the fitted plane during region expansion. Only candidate points satisfying all three criteria are incorporated into the current region, thereby forming progressive constraints from directional consistency and local spatial coplanarity to overall planar stability.

The geometric information represented by the three criteria has different emphases. Normal-vector orientation alone cannot further distinguish discontinuities with similar orientations but spatially separated locations, whereas the point-to-plane distance alone cannot fully represent the directional characteristics of a candidate point or its influence on the overall plane-fitting quality after incorporation. Therefore, the three criteria constrain candidate points in terms of orientation, local spatial position, and the overall planar characteristics of the growing region. Dominant discontinuity-set classification and the three geometric criteria operate at different levels: the former restricts the source of candidate points for region growing and limits direct cross-set expansion between different dominant discontinuity sets, whereas the latter further controls candidate-point inclusion within the same dominant discontinuity set to extract spatially separated individual discontinuities.

3 Results and analysis

3.1 Clustering results of rock discontinuities

To verify the applicability of the proposed morphological regularization-constrained fuzzy C-means clustering algorithm to tunnel rock discontinuity identification, a 3D laser scanning point cloud of a tunnel face was selected as the research object. According to the data acquisition results described above, the original point cloud contained approximately 5.40 million points. Considering the large data volume of the entire tunnel face point cloud, as well as the presence of local noise and boundary transition zones, direct clustering of all points would not only involve high computational costs but also make it difficult to highlight representative discontinuity combinations. Therefore, while preserving the spatial distribution characteristics of discontinuities on the tunnel face as much as possible, a representative region with concentrated discontinuity development and clear geometric morphology was selected as the clustering object. The cropped region contained 581,816 points.

The cropped original point cloud and the corresponding ground-truth labels are shown in Figure 4. Figure 4a presents the original point cloud, whereas Figure 4b shows the ground-truth labels manually established according to the spatial morphology, local orientation characteristics, and spatial continuity of the discontinuities. The rock mass surface in the study region exhibits distinct strip-like and planar structural features. Discontinuities with different orientations intersect with each other, and local surface undulations and boundary-transition characteristics are also present, indicating relatively complex spatial geometric characteristics. The detailed procedure used to establish the ground-truth labels is further described in Section 3.2.

On this basis, the proposed MRC-FCM algorithm was applied to classify the point cloud in the study region into dominant discontinuity sets, and the input parameters are listed in Table 1.

TABLE 1. Parameter settings of the morphological regularization-constrained FCM algorithm.

The discontinuity-set classification result obtained by MRC-FCM is shown in Figure 5. The point cloud in the study region was ultimately divided into five dominant discontinuity sets, which exhibit relatively good spatial continuity and banded distribution characteristics. Plane fitting was subsequently performed on the point cloud of each discontinuity set identified by MRC-FCM, and the corresponding representative dip directions and dip angles were calculated, as listed in Table 2.

TABLE 2. Orientation parameters of discontinuities identified by clustering.

3.2 Evaluation of classification performance

To quantitatively evaluate the discontinuity-set classification results obtained by MRC-FCM, the manually established ground-truth labels shown in Figure 4b were used as the reference. The labels were generated by considering the spatial morphology of the original point cloud, local orientation characteristics, and the spatial continuity of discontinuities. Dominant discontinuity sets were first identified according to the normal-vector orientations and spatial extents of continuous planar regions, with particular attention paid to discontinuity-intersection zones, edges of narrow strip-like structures, and locally fractured regions. For points with ambiguous boundaries, the labels were adjusted according to the overall orientation and spatial continuity of neighboring discontinuities. Five point-wise ground-truth classes were finally established and used for the subsequent classification-performance evaluation.

To further evaluate the reliability of the classification results, Figure 6 presents the ground-truth labels, the classification results of the MRC-FCM model, and the misclassification map for the cropped region. In the misclassification map, green points represent correctly classified points, whereas red points represent incorrectly classified points. As shown in the misclassification map, the errors are mainly concentrated in discontinuity intersection zones, along the edges of narrow strip-like discontinuities, and in local boundary transition zones, whereas large continuous discontinuity surfaces are dominated by correctly classified points. Overall, the model can effectively preserve the spatial continuity of discontinuities, and the residual errors mainly arise from the geometric mixing of different discontinuities within local neighborhoods.

Figure 7 shows the confusion matrix of the classification results. The classification accuracies of the five discontinuity sets are 95.05%, 94.89%, 93.88%, 94.03%, and 90.72%, respectively. The accuracy of each discontinuity set exceeds 90%, indicating that the model has high reliability in discontinuity classification.

In addition, global accuracy (GA) and mean Intersection over Union (mIoU) were further calculated to evaluate the classification performance, and the results are listed in Table 3. The GA value is 93.75%, and the corresponding mIoU is 87.99%. These two indicators further confirm the good performance and robustness of the MRC-FCM algorithm in discontinuity classification.

TABLE 3. Classification evaluation indicators of discontinuity sets.

To further evaluate the rationality of MRC-FCM-based discontinuity-set classification, the PC, PCE, XB, and FS indices were used to comprehensively analyze the clustering results under different numbers of discontinuity sets. The validity evaluation results of the discontinuity-set classification based on the four indices, namely, PC, PCE, XB, and FS, are shown in Figure 8. It can be seen that PC gradually increases and reaches its peak when the number of discontinuity sets is five, after which it begins to decrease. PCE gradually decreases and reaches its minimum value at five discontinuity sets, followed by an increase. XB obtains its minimum value at five discontinuity sets and then slightly rebounds. FS shows a certain fluctuating trend but also reaches its minimum value at five discontinuity sets. These four curves indicate that the clustering performance is optimal when the data are divided into five discontinuity sets. Therefore, the optimal number of discontinuity sets for this dataset is determined to be five.

3.3 Discontinuity segmentation results

After the dominant discontinuity sets were classified, the plane-constrained improved region growing algorithm was separately applied to the five dominant discontinuity sets. Point-cloud subsets corresponding to each set were first extracted according to the MRC-FCM clustering results, and neighborhood relationships were independently established within each set so that points belonging to different dominant sets did not participate in the candidate-point search of the same growing region. Normal consistency, point-to-plane residual, and RMS residual increment were then used as within-set growing criteria, and candidate points satisfying these conditions were progressively incorporated into the growing regions. The specific thresholds used in the plane-constrained improved region growing algorithm are listed in Table 4.

TABLE 4. Detailed parameters of the plane-constrained improved region growing algorithm.

After segmentation of each dominant discontinuity set, 477, 365, 380, 348, and 376 individual discontinuities were extracted from the five dominant discontinuity sets, respectively, giving a total of 1,946 individual discontinuities. Figure 9a presents the dominant discontinuity-set classification results, Figure 9b presents the individual-discontinuity segmentation results obtained using the plane-constrained improved region growing algorithm, and Figure 9c shows the pole-density plot of discontinuity orientations. Based on the dominant discontinuity-set classification, the large-scale discontinuity regions were further divided into multiple spatially separated individual discontinuities, while local small-scale and fractured discontinuities were also correspondingly segmented.

To further evaluate the accuracy of orientation calculation for the segmented individual discontinuities, 15 representative individual discontinuities were selected from each of the five dominant discontinuity sets, giving a total of 75 samples for comparative analysis. For each selected discontinuity, the discontinuity extent was manually identified and delineated in the original point cloud, and plane fitting was performed on the manually selected point-cloud data to obtain the corresponding dip direction and dip angle. The absolute differences in dip direction and dip angle between the automatically extracted results and the manually interpreted results were then calculated to evaluate the preservation of orientation parameters after region growing. The statistical results are presented in Table 5.

TABLE 5. Statistical errors in orientation parameters of individual discontinuities.

As shown in Table 5, the orientation parameters obtained for the different discontinuity sets are generally close to those obtained through manual interpretation. The mean absolute dip-direction errors of Sets 1–5 range from 1.21° to 2.83°, whereas the mean absolute dip-angle errors range from 0.91° to 2.35°, with some differences among the discontinuity sets. Set 1 shows relatively small overall errors, Sets 3 and 4 exhibit relatively larger dip-direction errors, and Set 5 mainly shows an increase in dip-angle error. For all 75 individual discontinuities, the overall mean absolute errors in dip direction and dip angle are 1.91° and 1.62°, respectively, indicating a relatively low overall error level.

The differences in errors among the discontinuity sets are mainly related to the spatial scale of individual discontinuities, boundary completeness, and local point-cloud undulation. For discontinuities with relatively complete boundaries and distinct planar characteristics, the fitted plane remains relatively stable during region growing, resulting in small orientation differences between the automatically extracted and manually interpreted results. In locally fractured regions, areas with indistinct boundary transitions, or regions where point clouds from adjacent discontinuities are mixed, the distribution of points involved in plane fitting changes, which may cause a certain deviation in the fitted-plane normal vector. Small-scale discontinuities and surfaces with pronounced local undulations are also more sensitive to boundary and anomalous points. Overall, the automatically extracted results are in good agreement with the manual interpretation, indicating that the plane-constrained improved region growing method can spatially separate individual discontinuities while reasonably preserving their overall planar characteristics and basic orientations. In the preceding segmentation procedure, the dominant discontinuity-set constraint limits the inter-set search range of region growing, whereas the three within-set geometric criteria further control the spatial expansion of individual discontinuities. Together, these two levels form a continuous procedure from dominant discontinuity-set identification to individual-discontinuity extraction.

4 Discussion

4.1 Advantages and robustness of the proposed method

For large-section tunnels, the continuity, boundary locations, and spatial combinations of discontinuities on the tunnel face directly affect the block cutting pattern and potential instability modes of the surrounding rock. Therefore, the clustering results should not only distinguish different discontinuity sets but also preserve the authenticity of their spatial distribution as much as possible. The comparison results between different clustering methods and the ground-truth labels are shown in Figure 10. As shown in Figure 10a, the discontinuities generally exhibit interlaced strip-like and sheet-like distributions. Different discontinuity sets show good spatial continuity, and relatively clear boundaries are still maintained in the central intersection zone and near the dominant strip-like structures on the right side. Compared with the ground-truth labels, the DBSCAN clustering results show obvious erroneous merging and fragmentation in local areas. Some strip-like structures that should belong to different discontinuity sets are connected into the same discontinuity set, whereas some continuous strips are divided into scattered patches. This indicates that DBSCAN has difficulty in stably recovering the true spatial combination of discontinuities under the current data conditions. Figures 10c,d show the clustering results of the traditional FCM method when the neighborhood point numbers are set to 25 and 30, respectively. When the neighborhood point number is 25, the FCM algorithm can divide the discontinuities into five discontinuity sets. However, obvious local bias occurs in the complex intersection zone in the lower-left part and near several slender strip-like structures. This is manifested as abnormal absorption of some local regions by the same discontinuity set, indicating that the algorithm falls into a suboptimal solution in local areas. When the neighborhood point number increases to 30, the clustering results instead show obvious over-smoothing characteristics. The originally fine strip-like structures become significantly blurred, obvious clustering errors occur, boundary ranges expand, and local small-scale discontinuities are absorbed into adjacent dominant discontinuity sets, making it difficult to preserve the true spatial morphology of the discontinuities.

Overall, Figure 10 shows that both DBSCAN and traditional FCM have certain limitations in complex intersection zones, fine strip-like structures, and boundary transition regions. DBSCAN mainly suffers from local erroneous merging and strip fragmentation, whereas traditional FCM tends to fall into local optima when the neighborhood parameter is small and tends to produce over-smoothed results when the neighborhood parameter is large. In contrast, the proposed MRC-FCM algorithm better preserves the continuous distribution and boundary morphology of discontinuity sets. This indicates that, after introducing morphological constraints, the clustering results not only satisfy normal-vector directional similarity but also better reflect the spatial coplanarity and geometric continuity of discontinuities, thereby showing better stability and applicability in discontinuity-set classification.

This result also indicates that discontinuity-identification errors are mainly concentrated in boundary, intersection, and fractured regions where both directional information and local spatial relationships vary simultaneously, while the introduction of the local planar residual provides an additional geometric constraint for these regions beyond a single orientation feature.

Following the comparison among different clustering methods, the robustness of the proposed method was further evaluated with respect to changes in point-cloud quality in terms of point-cloud density, data completeness, and coordinate perturbation. First, the original point cloud was randomly downsampled to retain 75% and 50% of the original points, respectively, thereby simulating different point-cloud-density conditions. Second, a continuous region was removed from the original point cloud to simulate local data loss caused by occlusion or limited viewing angles during actual three-dimensional scanning. In addition, zero-mean Gaussian random noise with standard deviations of 0.0025 m and 0.005 m was separately added to the three-dimensional coordinates of the original point cloud to simulate different levels of coordinate perturbation. The same discontinuity-identification procedure was applied to all perturbed datasets, and the results are shown in Figure 11.

As shown in Figure 11, the main spatial distributions of the dominant discontinuity sets are generally preserved after the point-cloud density is reduced. When 75% of the original points are retained, the identification result remains generally close to that obtained from the original point cloud. The spatial positions, extension directions, and overall morphologies of the major strip-like discontinuities remain essentially stable, with only limited dispersion occurring in some fine discontinuities and boundary regions. When the point cloud is further reduced to 50% of the original data, the local point density decreases substantially, the continuity of some fine discontinuities and intersection regions declines, and local changes along the boundaries become more evident. Nevertheless, the main distributions of the dominant discontinuity sets remain clear, and the overall spatial pattern of the major discontinuities does not change substantially. These results indicate that reduced point-cloud density primarily affects local boundaries and fine discontinuities, whereas larger discontinuities with distinct planar characteristics retain relatively stable identification results.

Under the local point-cloud-loss condition, the removed region forms an obvious spatial gap, and the continuous distribution of discontinuities within and adjacent to the missing region is affected to some extent. However, the main distributions of the dominant discontinuity sets in the remaining regions are generally consistent with those obtained from the complete point cloud, without obvious large-scale class reorganization. This indicates that the influence of local data loss is mainly concentrated near the missing region, whereas discontinuities located farther away can still be relatively stably classified using the remaining local geometric information.

The results under different Gaussian-noise conditions show that when the noise standard deviation is 0.0025 m, the dominant discontinuity sets generally retain spatial distributions similar to those obtained under the noise-free condition. The spatial positions and overall morphologies of the large central discontinuities and major strip-like regions are well preserved, and the changes are mainly concentrated near local boundaries and fine discontinuities. When the noise standard deviation increases to 0.005 m, the spatial positions, extension directions, and large-scale morphologies of the major discontinuity sets can still be identified, although class mixing increases at some intersection and boundary regions and the local continuity of fine discontinuities decreases to some extent. As the coordinate perturbation increases, the spatial distribution of points within local neighborhoods changes, affecting both normal-vector estimation and planar-feature calculation. Consequently, boundary-transition zones and regions with more complex local geometries exhibit greater sensitivity to noise.

Overall, the results under different data perturbations show that point-cloud downsampling, local data loss, and random coordinate noise mainly affect fine discontinuities, intersection zones, and boundary-transition regions, whereas their influence on larger dominant discontinuities with distinct planar characteristics is relatively limited. When the point-cloud retention ratio is reduced to 50%, a continuous local region is missing, or Gaussian coordinate noise with a standard deviation of 0.0025–0.005 m is introduced, the main spatial positions, extension directions, and overall distributions of the dominant discontinuity sets are still reasonably preserved. Overall, the proposed method shows good adaptability to a certain degree of point-cloud-density reduction, local information loss, and random coordinate perturbation, although the identification stability of fine discontinuities and boundary regions decreases as the data quality deteriorates.

4.2 Parameter sensitivity analysis

In MRC-FCM, the local neighborhood radius R, regularization coefficient λ, and upper limit of the initial cluster number Kmax control the spatial scale of local geometric-feature calculation, the strength of the morphological constraint, and the degree of initial over-segmentation, respectively. The values R = 0.05 m, λ = 0.5, and Kmax = 8 were adopted in the preceding dominant discontinuity-set classification. In this section, the effects of continuous variations in these parameters on the clustering results are further analyzed to evaluate the adopted parameter settings.

The local neighborhood radius R determines the spatial range included in local plane fitting. When the neighborhood is relatively small, fewer effective neighboring points participate in plane fitting, and the estimated local normal vectors and planar residuals become more susceptible to point-cloud undulations and isolated points. As R increases, the local geometric information becomes more stable; however, an excessively large neighborhood may cross the actual boundaries between adjacent discontinuities and cause geometric information from different discontinuities to be jointly included in plane fitting. A radius of R = 0.05 m allows the local neighborhood to contain sufficient geometric information while maintaining the representation of local discontinuity characteristics. Therefore, different neighborhood scales ranging from 0.02 to 0.10 m are examined in this section, with 0.05 m retained as the value adopted in the preceding analysis.

The regularization coefficient λ controls the relative contribution of the normal-vector directional feature and the local planar residual. Specifically, λ controls the strength of the morphological constraint relative to the directional-distance term. When λ = 0, the local morphological constraint is inactive. As λ increases, the influence of local planar characteristics on membership assignment gradually increases. The local planar residual is mainly introduced to reduce the influence of boundary points, fractured points, and points with poor local geometric quality on cluster centers. The value λ = 0.5 allows the morphological constraint to participate in clustering at a moderate strength. In this section, λ is further varied from 0 to 1 at an interval of 0.25 to investigate the classification performance as the morphological constraint gradually increases.

Kmax represents the upper bound of the candidate-cluster number allowed during the initial over-segmentation stage. In the preceding analysis, PC, PCE, XB, and FS were used to comprehensively evaluate different numbers of clusters, and the results showed that five discontinuity sets provided the best overall clustering performance. Because MRC-FCM follows an initial over-segmentation, geometric-quality screening, and similar-cluster merging procedure, the initial number of clusters should be moderately larger than the effective number of discontinuity sets. This prevents potential discontinuities from being prematurely merged during the initial stage and provides sufficient representation space for boundary regions, transition zones, and small candidate discontinuities. Therefore, Kmax = 8 was adopted in the preceding analysis. In this section, Kmax is further varied to investigate the influence of the degree of initial over-segmentation on classification accuracy and the final number of effective discontinuity sets.

To quantitatively analyze the effects of the three parameters, a one-factor-at-a-time strategy was adopted, in which only one parameter was varied while the other two were kept at the values used in the preceding analysis. Specifically, R was set to 0.02, 0.04, 0.05, 0.06, 0.08, and 0.10 m; λ was set to 0.00, 0.25, 0.50, 0.75, and 1.00; and Kmax was set to 4, 6, 8, 10, and 12. Global accuracy (GA), mean intersection over union (mIoU), and the number of effective discontinuity sets after screening and merging were used to comprehensively evaluate the clustering results under different parameter settings. The specific parameter settings are listed in Table 6. The clustering results obtained under different values of R, λ, and Kmax are summarized in Tables 7–9, respectively.

TABLE 6. Parameter settings for the MRC-FCM sensitivity analysis.

TABLE 7. Clustering results under different R values.

TABLE 8. Clustering results under different λ values.

The results obtained under different neighborhood radii show that GA and mIoU generally increase first and then decrease as R increases. At relatively small neighborhood scales, the number of effective points participating in local plane fitting is limited, making the estimated normal vectors and planar residuals more sensitive to local point-cloud undulations and causing local dispersion within the same discontinuity. As R increases to 0.04–0.06 m, the calculation of local geometric characteristics becomes more stable, and both the classification accuracy and the final number of effective discontinuity sets remain relatively stable. When R is further increased, the probability that the neighborhood crosses the boundary between adjacent discontinuities increases, leading to mixing of local geometric information from different discontinuities and a corresponding decrease in classification performance.

Variation in the regularization coefficient λ reflects the strength of the local morphological constraint during clustering. When λ = 0, the local planar residual does not participate in membership adjustment, and clustering is mainly controlled by normal-vector orientation, resulting in relatively weak constraints in boundary-transition regions and areas with poor local geometric quality. As λ increases, local coplanarity information gradually contributes to clustering, and both GA and mIoU increase accordingly. When λ is further increased, the local morphological term becomes stronger, and some points that exhibit local surface undulations but still belong to the same discontinuity are more strongly constrained, causing the classification performance to decrease. Therefore, an appropriate local planar residual provides complementary information to a single orientation feature, whereas an excessively strong morphological constraint may reduce the local continuity of discontinuities.

Kmax mainly affects the ability of the initial over-segmentation stage to represent potential discontinuity sets. When Kmax is smaller than the number of effective dominant discontinuity sets, some discontinuities are merged during the initial stage, resulting in an insufficient final number of discontinuity sets and decreases in both GA and mIoU. When Kmax is increased above the effective set number, geometric-quality screening and similar-cluster merging allow the final number of effective discontinuity sets to remain relatively stable. Further increasing Kmax increases the number of initial candidate clusters but also produces more redundant and boundary-transition clusters, providing only limited improvement in the final classification performance.

The above single-factor analyses characterize the clustering responses when each parameter is varied independently. To further examine the spatial clustering characteristics when multiple parameters simultaneously deviate from their suitable ranges, two representative parameter combinations were selected for comparison of under-segmentation and over-segmentation.

When the parameter values are set as shown in Table 10, a large R combined with small λ and Kmax values tends to cause under-segmentation. When R is too large, the local neighborhood can easily cross the true boundaries of discontinuities, causing the geometric features of different discontinuities to be averaged during plane fitting. This results in blurred boundaries and the merging of fine strip-like discontinuities into adjacent dominant discontinuities. Meanwhile, when λ is too small, the morphological constraint becomes insufficient, and the model becomes closer to traditional FCM. In this case, clustering mainly depends on normal-vector similarity, and mixed points near discontinuity boundaries are more likely to be absorbed into adjacent large clusters. When Kmax is insufficient, the model lacks adequate representation capability at the initial stage, and discontinuities with different orientations or adjacent spatial positions can easily be merged into the same discontinuity set. Therefore, the results usually exhibit under-segmentation, characterized by an insufficient number of discontinuity sets, boundary diffusion, and missed identification of small-scale discontinuities, leading to over-smoothed clustering results, as shown in Figure 12.

TABLE 9. Clustering results under different Kmax​ values.

TABLE 10. Under-segmentation parameters of the morphological regularization-constrained FCM algorithm.

Conversely, when the parameter values are set as shown in Table 11, a small R combined with large λ and Kmax values tends to cause over-segmentation. When R is too small, the number of neighboring points used for local plane fitting is insufficient, making the normal vectors and local planar residuals more sensitive to noise. As a result, scattered patch-like distributions may appear within the same discontinuity, and intra-set continuity decreases. When λ is too large, the model overemphasizes local planar consistency, causing some points with slight surface undulations but still belonging to the same discontinuity to be excluded from the dominant discontinuity set. This leads to intra-cluster shrinkage and local fragmentation. At the same time, an excessively large Kmax introduces more invalid small clusters and boundary transition clusters, making the subsequent cluster pruning and merging process more complicated. Therefore, although the model can identify more fine local features under this condition, the overall continuity of dominant discontinuities is weakened, as shown in Figure 13.

TABLE 11. Over-segmentation parameters of the morphological regularization-constrained FCM algorithm.

Based on the two parameter conditions, the effects of the MRC-FCM parameters on the clustering results have clear geometric meanings. A too-small R weakens the stability of local plane fitting, whereas a too-large R may cross the true boundaries of discontinuities. A too-small λ leads to insufficient morphological constraint, whereas a too-large λ causes excessive constraint and results in fragmented discontinuity sets. A too-small Kmax leads to under-segmentation, whereas a too-large Kmax increases the number of invalid clusters and boundary clusters.

For the tunnel face point cloud selected in this study, the recommended parameter combination is shown in Table 12. When R = 0.05 m, λ = 0.5, and Kmax = 8, GA and mIoU reach 93.75% and 87.99%, respectively, and five effective dominant discontinuity sets are obtained after screening and merging. The classification performance and final effective set number also remain generally stable within the neighboring parameter ranges. These results indicate that this parameter combination provides a suitable balance among local geometric-feature representation, morphological constraint, and initial discontinuity-set representation.

TABLE 12. Recommended parameters of the morphological regularization-constrained FCM algorithm.

For point clouds with different geometric characteristics, the parameters can be determined progressively according to the local geometric scale, morphological-constraint strength, and effective number of discontinuity sets. A candidate range of the neighborhood radius R can first be established according to the local scale of the point cloud. The regularization coefficient λ can then be adjusted after normalization of the local planar residual. Finally, PC, PCE, XB, and FS, can be used to estimate the candidate effective number of discontinuity sets, on the basis of which Kmax is progressively increased until the number of effective discontinuity sets obtained after screening and merging becomes stable. This procedure provides a parameter-selection strategy for point clouds with different characteristics.

Under the combined effect of these parameters, the MRC-FCM algorithm can identify discontinuity sets in large-section tunnel faces more stably than traditional FCM and DBSCAN, especially under conditions involving multi-directional intersections, blurred boundaries, and obvious local fragmentation.

4.3 Comparative validation

To further validate the reliability of the proposed method in discontinuity orientation identification, the open-source point cloud dataset of the Ouray roadcut slope in Colorado, USA, was selected for comparative analysis. This dataset has been used by Riquelme et al., 2014 and Kong et al., 2020 for automatic rock discontinuity identification, providing a good basis for algorithm comparison. The method proposed by Riquelme et al. mainly uses PCA-based coplanarity testing, KDE-based kernel density estimation, and DBSCAN clustering for discontinuity identification. In contrast, the method proposed by Kong et al., 2020 employs IRPF normal-vector calculation, CFSFDP density peak clustering, Fisher’s K-based noise removal, density-ratio clustering, and RANSAC plane fitting to achieve automatic discontinuity identification and geometric parameter extraction. In this study, the MRC-FCM algorithm was applied to the same dataset, and the orientations of the five identified dominant discontinuity sets were compared with the results of the above two methods, as listed in Table 13. The original point cloud of the open-source dataset and the clustering results obtained by the proposed method are shown in Figure 14.

TABLE 13. Comparison between the proposed algorithm and other algorithms on the open-source dataset.

As shown in Table 13, the orientations of the five discontinuity sets identified by the proposed method are generally close to those obtained by Riquelme et al., 2014 and Kong et al., 2020. The dip-direction and dip-angle errors of all discontinuity sets are less than 4°, which falls within the acceptable error range for automatic discontinuity identification from point clouds. Overall, the identification results of the proposed method on the open-source Ouray roadcut slope dataset are in good agreement with those of existing representative algorithms, indicating that MRC-FCM has good reliability and stability in dominant discontinuity-set classification and orientation extraction.

5 Conclusion

To improve the automatic identification accuracy of rock discontinuities from 3D point clouds of tunnel faces, this study proposed a morphological regularization-constrained FCM clustering algorithm (MRC-FCM). The proposed method was combined with a plane-constrained improved region growing algorithm to achieve classification of dominant discontinuity sets, individual discontinuity segmentation, and orientation extraction. The main conclusions are as follows:

• MRC-FCM introduces a local planar residual constraint into traditional FCM so that the clustering process simultaneously considers normal-vector similarity and spatial coplanarity. The proposed method reduces the influence of boundary mixed points, fractured points, and noise points on cluster centers, thereby improving the continuity and boundary-preserving capability of discontinuity-set classification.
• A representative point-cloud region from a tunnel face was selected for validation. MRC-FCM divided the study region into five dominant discontinuity sets, with classification accuracies of 95.05%, 94.89%, 93.88%, 94.03%, and 90.72%, respectively, while the overall GA and mIoU reached 93.75% and 87.99%. Under point-cloud downsampling, local data loss, and Gaussian coordinate-noise perturbations, the main spatial positions and overall distributions of the dominant discontinuity sets were still reasonably preserved, indicating that the proposed method has good adaptability to a certain degree of point-cloud-quality variation.
• Based on the dominant discontinuity-set classification, the plane-constrained improved region growing algorithm was further used to extract 1,946 individual discontinuities. For 75 representative individual discontinuities, the mean absolute errors in dip direction and dip angle were 1.91° and 1.62°, respectively. The pole-density plot showed five distinct density concentration zones that were generally consistent with the dominant discontinuity-set classification, indicating that the extracted individual discontinuities reasonably preserved their overall planar characteristics and basic orientations.
• Compared with DBSCAN and traditional FCM, DBSCAN is prone to erroneous merging and strip fragmentation in complex intersection zones, whereas traditional FCM tends to fall into local optima when the neighborhood parameter is relatively small and to produce over-smoothed results when the neighborhood parameter is relatively large. In contrast, MRC-FCM better preserves strip-like discontinuity boundaries and local fine structures. Further comparison using the open-source Ouray roadcut slope point-cloud dataset shows that the orientations of the five dominant discontinuity sets identified by the proposed method are generally close to those reported by Riquelme et al. (2014) and Kong et al. (2020) with the dip-direction and dip-angle differences of all sets being less than 4°. These results further demonstrate the reliability of MRC-FCM in discontinuity-set classification and orientation extraction.

Generative AI statement

The author(s) declared that generative AI was not used in the creation of this manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Footnotes

Edited by: Zhen Zhang, Taiyuan University of Science and Technology, China

Reviewed by: Zhengzheng Cao, Henan Polytechnic University, China

Meilu Yu, Anhui University of Science and Technology, China