Exploiting a Semi-Automatic Point Cloud Segmentation Method to Improve the Quality of Rock-Mass Characterization. The Cima Grappa Conservative Restoration Case Study

Abstract

This paper presents results from applying semi-automatic point cloud segmentation methods in the underground tunnels within the Military Shrine’s conservative restoration project in Cima Grappa (Italy). The studied area, which has a predominant underground development distributed in a network of tunnels, is characterized by diffuse rock collapsing. In such a context, carrying out surveys and other technical operations are dangerous activities. Considering safety restrictions and unreachable impervious tunnels, having approached the study area with the scan-line survey technique resulted in only partial rock mass characterization. Hence, the geo-mechanical dataset was integrated, applying a semi-automatic segmentation method to the point clouds acquired through terrestrial laser scanning (TLS). The combined approach allowed for remote performance of detailed rock mass characterization, even remotely, in a short time and with a limited operators presence on site. Moreover, it permitted extending assessing tunnels’ stability and state of conservation to the inaccessible areas.

1. Introduction

The Cima Grappa Military Shrine restoration project has been promoted and financed by the Italian Ministry of Defence of the Italian Presidency of the Council of the Ministers. It concerns a mountain area on the top of Monte Grappa (also called “Cima del Grappa”) in the Venetian Prealps context (Figure 1).

The Cima Grappa Military Shrine is the largest Italian military ossuary of the First World War. During the First World War, many Italian soldiers died. Most of them were buried in mass graves close to battlefields. Between 1920 and 1930, remains of the fallen were unearthed and reburied in a large war memorial. Politicizing the memory of the fallen was one of the primary purposes. The war memorials, designed by architects close to the regime, were constituted near former front lines and kept under the patronage of a special commission from the former Ministry of War [1].

The hereafter presented activities that insisted on the underground connection area between the Milan Barracks and the former underground infirmary, focus on securing the former infirmary rooms and tunnels. These areas were dug in a fractured limestone rock mass. The spaces are cramped, and the accessible areas are tunnels that span from 1.2 to 2.2 m in height and 0.9 to 2.4 m in width. The overall study area’s extension is 300 m² with a tunnel’s lengths of about 100 linear meters (Figure 2). During past years, several collapses have occurred in the underground infirmary and the access tunnels. Thanks to the restoration project, the site will become a new tourist area. A detailed geo-mechanical characterization of the former infirmary’s underground area has been necessary to identify and design the optimal stabilization and restoration interventions [2,3]. Discontinuities play a crucial role in determining rock mass stability [4,5], and we need to investigate extensive rock-mass portions [6,7,8,9,10] to perform sound statistical analyses. The manual surveying of fractures and planar rock mass surfaces is one of the most fundamental but time-consuming activities performed by surveyors.

The classical geo-mechanical survey [11] on its own would not have brought sound results [12]. Primarily because the tunnels’ roof was too high to be reached with classic instrumentation; and secondly, because the amount of data acquired through the classical scan-line survey would not have a good statistical significance, considering the geometrical complexity of the tunnels. Finally, the relatively poor conservative status of the underground area and some inaccessible spaces would not have allowed the technicians to access all the spaces. Furthermore, the restoration project has a strong conservative purpose and using non-destructive investigation methods is a compulsory request. The laser scanning technique was applied to collect high-density point clouds to identify the discontinuities main sets [13,14,15,16,17,18,19]. The geo-mechanics surveying community have been extracting fractures and planar surfaces orientations from 3D outcrop point clouds through semi-automatic and automatic methods developed in the last 15 years.

Using random sample consensus (RANSAC) algorithm-based methods to segment point clouds into subsets [20,21,22], deriving triangulated irregular networks from point cloud and grouping neighbouring polygons with a similar orientation to obtain planar features [23,24], and some other methods based on k-means clustering [25], moving sampling cube, point attributes, neighbouring points coplanarity testing, and principal component analysis (PCA), were developed and can be applied to gather rock mass geometrical settings instead of performing a classical geo-mechanic survey. Considering the narrow and complicated geometrical setting of the study area, a method designed with supervision capability was chosen and applied within the presented paper. Finally, kinematic analyses were carried out by applying 3D methodology [26,27,28,29].

2. Materials and Methods

In correspondence to inaccessible areas and areas prone to collapsing events, the classic geo-mechanical survey (Figure 3) scan-line detection method [30] of the rock mass was substituted with measures obtained from laser scanning datasets [31,32]. The geo-mechanical analysis results were combined with the analysis of the point clouds acquired by the laser scanner. In particular, a semi-automatic extraction of rock mass structural data from high-resolution laser point clouds was adopted [33,34,35]. This section first introduces the classical geo-mechanical scan-line method (Figure 3) and then uses segmentation procedures to extract structural data from point clouds.

2.1. Classic Manually Performed Geo-Mechanical Scan-Line Survey

A geo-mechanical survey of the calcareous rock mass surfacing in the infirmary rooms was performed. It consisted of the following steps, (a) geo-mechanical characterization of intact rock, (b) geo-mechanical characterization of discontinuities, and (c) geo-mechanical characterization of the rock mass. Skipping a detailed description of a geo-mechanical survey [36] beyond the present paper’s scope, the essential parameters to describe joints, cracks, and discontinuities are reported. In particular, cracks orientation, number of families of discontinuity, spacing, persistence, opening, roughness, surface resistance, degradation, filling, and filtration [37].

2.2. Remote Geo-Mechanical Survey

Using terrestrial laser scanning (TLS) for rock mass discontinuity characterization [38,39,40,41] was reported in several scientific works [42,43]. The DiAna3D algorithm is applied to extend the rock-mass geo-mechanical assessment to the less accessible areas. The method is more rapid than the standard scan line survey, and it can minimize the operator’s exposure to rock collapse risk. Thanks to the laser scanning data high spatial resolution, it is possible to use point clouds to integrate the traditional survey data and obtain information on the rock mass geometry remotely if compared to the manual magnetic compass method.

In the specific case, a calculation algorithm was applied to the point clouds to extract accurate geometric information, valid for the rock mass reconstruction and characterization. In particular, through DiAna, the discontinuity plans were identified semi-automatically, and their spatial orientations were calculated.

The analysis was carried out in correspondence with the limestone rock mass’s main outcrop areas and extracted 2761 discontinuous surfaces. However, a surface of discontinuity that is not perfectly planar and has a high persistence may be fragmented into a certain number of "sub-surfaces". The number of discontinuities extracted is overestimated.

The laser scanning technique was used to generate point clouds and measure geometrical parameters as persistence, spacing, roughness, and discontinuities orientation.

2.3. Acquisition Campaign

The laser scanning acquisition campaign involved the underground rooms of the infirmary behind the Milan barracks. The used laser scanner model is a Riegl LMS-Z420i ground-based (

Table 1. Laser scanner hardware parameters.

Parameters	Value
Measurement Range	0.2–1000 m
Accuracy	10 mm
Laser class	1
Minimum angle step width	0.004°

).

A series of high-resolution laser scans were carried out from 13 different scanning positions starting from the square in front of the Milan barracks to mapping the underground hospital area. The obtained point clouds were then georeferenced by aligning them to that acquired externally. The scans were carried out mainly with the vertical instrument setup; a horizontal configuration was also used in correspondence with the major collapses to allow the complete acquisition of the tunnel roof. For each scanning position, low-resolution framing scans and detail scans were acquired. A total of over 40 million points were thus acquired, with the average angular resolution of 0.056°, a minimum acquisition distance of 1.5 m, and a maximum acquisition distance of 20 meters. One point/cm² is the average point cloud density of the processed point clouds after optimization and unification. Before proceeding with the laser scan, some targets/reflectors were positioned within the scene (Figure 4).

Up to 15 cylindrical reflectors (retroreflector targets) with a pre-determined size of 50 × 50 mm were deployed, performing each scanning acquisitions. Reflectors’ use is linked to the possibility of the acquisition software recognizing them in the various scans and combining them into a single point cloud, minimizing the problems related to shaded areas. Thanks to the reflector’s high reflectivity, their positions were measured with the highest TLS resolution.

The TLS scans registration accuracy was improved by applying the multi-station adjustment algorithm (MSA). It was developed by the laser manufacturer and based on the recognition of homologous plans materialized for each scan (which therefore require a scan overlap between them) and their iterative alignment, up to minimizing the position differences (Figure 5).

The overlapping value varied from 35 to 60%. Further, 60% could be considered a relatively high value. It is the partial consequence of the very articulated tunnels’ morphology, characterized by many bumps, depressions, and irregularities of the rock mass surface. Indeed, the consequent massive number of shaded areas led to complicated scanning geometries and made designing scan positions effectively a problematic task.

Cross-sections of the infirmary point cloud are reported in Figure 6 to better understanding the former infirmary’s spaces and shapes.

2.4. Point Clouds Segmentation for Remote Geo-Mechanical Survey

As mentioned in the section above, a semi-automatic segmentation method, DiAna3D, was used to extract geometric characteristics of discontinuities and relative planes. By processing point clouds, some of the rock-mass parameters [44,45] are block size, the number of sets, persistence, orientation (Figure 7), spacing/frequency (and derived RQD), and scale-dependent roughness can be determined.

The three-dimensional approach of the segmentation method is hereafter shortly reported. The first step is aimed at identifying the principal planes on the cloud. A searching cube contains a selection of a point cloud subset. The best-fitting plane is identified for the searching cube volume by applying the least square method. The plane is described with the general equation (Equation (1)) end its direction cosines characterize it $\left(l,m,n\right)$:

$$ax+by+cz+d=0$$

(1)

The algorithm computes the selected points’ standard deviation for the best fitting plane. Simultaneously, the operator can designate a threshold value (s_t) depending on the point cloud resolution, roughness of the rock slope, and searching cube dimension. A cluster is formed with the points extracted from the point cloud when s < s_t. A box analysis can be performed; in this case, the searching cube is moved along the geographic axes according to a regular pattern to process the entire point cloud. A discrete analysis can also be processed; in this case, the searching cube is centred at every point of the point cloud. For each valid determined cluster, the correlate plane orientation is calculated by applying Equations (2) and (3):

$$\alpha =arctan\left(\frac{m}{l}\right)+Q$$

(2)

$$\beta =arctan\left(\frac{n}{\sqrt{l^2+m^2}}\right)$$

(3)

where α and β are the dip and dip/direction of the plan; $\left(l,m,n\right)$ are the direction cosines of the plane; $Q$ is a constant which depend on $l,m$.

Although the method is quite far from moving with a compass and clinometer on the rock face, taking measurements on discontinuity planes allows local face orientation measurements.

The operator can choose a variation on cube size, and it returns details that can be observed on the point cloud.

The extraction process of the discontinuity planes consisted of 3 main steps: In the first one, organized point clouds were analyzed (Figure 8), looking for the points around which the local geometry is sufficiently planar (standard deviation of the points concerning the interpolation plane, below a user-defined threshold (blue areas in Figure 9a; starting from these points, step 2 consisted of a progressive germination algorithm that allowed to gradually add all the neighbouring points, co-planar with the initial nucleus (Figure 9b), and to define its external 3D polygon (Figure 9c,d). The last step allowed to combine the results obtained from the individual scan positions by eliminating redundancies in the overlapping areas. For each identified plane, its persistence was calculated by measuring the maximum linear distance between two points of the cluster of points that compose it. The 3D polygons that delimit the discontinuity planes with persistence greater than 10 cm; thus, identified are shown in Figure 10a,b.

Once cluster orientations are determined (Figure 11), they are plotted on stereographic projections (Figure 12).

In this way, defining a discontinuity set as a series of discontinuities, for which orientation is broadly the same, discontinuity sets can be easily identified and extracted from the point cloud. A common ID is assigned to the clusters belonging to the same set. This last-mentioned process is manually driven to have more substantial supervision of the set recognition procedure and guarantees a real-time integration with possible validated in-field data.

By family of discontinuities, we mean the set of discontinuities with a similar orientation (sub-parallels), generally of the same type and, presumably, originating with the same genetic mechanism.

Four families of main discontinuities were recognized (Figure 12), whose respective modal planes’ orientations (

Table 2. Values (dip and dip direction) of joint families recorded by operators using classic scan-line survey and values of joint families extracted by the remote geo-mechanical survey.

Scan-Line Survey	Remote
JN1 059/78°	JN1: 65°/76°
JN2 002/87°	JN2: 001°/82°
JN3 115/75°	JN3: 121°/75°
BG 250/09°	BG: 245°/13°

).

2.5. 3D Kinematic Stability Analysis of the Rock Mass

After doing the geo-mechanical characterization of the intact rock, the discontinuity system, and the rock mass, the data obtained were used to carry out a 3D kinematic analysis of the rock mass [46]. The stability of the blocks composing a rock mass is essentially controlled by the discontinuities’ geometric and geo-mechanical properties. The term “kinematic analysis” refers to analyzing blocks of the rock’s behaviour under their weight and with purely attractive sliding resistance along the discontinuity surfaces [47]. The analysis allows investigating the following instability mechanisms (Figure 13): plane failure (PF), wedge failure (WF), block toppling (BT), flexural toppling (FT), and free fall (FF).

The kinematic analysis of slope stability is a valuable tool for understanding rock mass behaviour. It represents a valid connection point between the technical geological survey and the design phase. The analysis is generally carried out by graphically verifying certain geometric conditions on the lower hemispheric projection. Equiangular projection is usually used as it allows to maintain the relations between the angles.

The basic assumptions for kinematic analysis are:

• the rock mass is divided into blocks by flat and infinitely persistent discontinuities;
• the shear resistance along the discontinuity planes is purely due to friction;
• the block system is subject to weight force only.

The term “kinematic analysis” refers to analyzing individual blocks of the rock’s behaviour under its weight and with purely attractive sliding resistance on discontinuous surfaces.

This type of analysis applied to the slopes constitutes a valuable tool for understanding rock mass behaviour. They allow for a rapid assessment of the propensity to failure for the various breaking mechanisms in the sectors of the slopes analyzed.

The kinematic analysis of the different failure mechanisms can be performed using geometric relationships amongst discontinuity’s planes, intersection lines between different discontinuity planes, and the excavation face plane.

These can be expressed as vector relations between the vector units that identify these planes, expressed in polar coordinates using the azimuth $\propto$ and the inclination on the horizontal $\beta$.

In particular, the planes of discontinuity are identified by the versors of the lines of maximum slope $p\equiv \left({\alpha }_p;{\beta }_p\right)$, or their normal $n\equiv \left({\alpha }_n;{\beta }_n\right)$ or poles;

The intersection lines between pairs of plans are identified by versors $i\equiv \left({\alpha }_i;{\beta }_i\right)$.

While the excavation line is identified by the versor of the maximum slope line $f\equiv \left({\alpha }_f;{\beta }_f\right)$.

By convention, all versors units are considered with the extreme in the lower half-space of the horizontal plane.

A quantitative approach can be pursued by defining a kinematic hazard index for each kinematic mechanism, which expresses the relative probability that a given mechanism failure may occur based on the structural and geo-mechanical conditions of the rock mass. Such indexes can be calculated by counting the following quantities:

• Npf = number of poles of the discontinuities that satisfy the conditions for planar sliding;
• I_wf = number of the intersections of the discontinuities that satisfy the conditions for the sliding of wedges;
• N_bt = number of the poles of the discontinuities that satisfy the conditions for direct overturning;
• I_bt = number of intersection lines that satisfy the conditions for direct overturning;
• N_ft = number of the poles of the discontinuities that satisfy the conditions for bending overturning.

Indicating with N the total number of discontinuities sampled in the domain and with I = 0.5 (N2-N), the number of all possible intersections, it is possible to calculate the hazard kinematics indexes reported in

Table 3. Hazard kinematic indexes [49], N = number of poles satisfying a specific failure condition; I = number of intersections satisfying a specific failure condition.

Instability Mechanism	Hazard Index
Planar failure	$C_{pf}=\frac{N_{pf}}{N} \left(\%\right)$
Wedge failure	$C_{wf}=\frac{I_{wf}}{I} \left(\%\right)$
Flexural toppling	$C_{ft}=\frac{N_{ft}}{N} \left(\%\right)$
Block toppling	$C_{bt}=\frac{N_{bt}}{N}·\frac{I_{bt}}{I}\left(\%\right)$

.

2.6. Extraction of Detailed Meshes

The overall point cloud, obtained by combining the information acquired from each scanning position, was triangulated, getting a continuous surface from discrete data of a high point density.

The resulting 3D model is shown in Figure 14.

From the detailed mesh obtained by applying the Poisson surface reconstruction method [49], it was then possible to extract contour lines with an equidistance of 10 cm, as shown in Figure 15.

3. Results

Thanks to the performed survey campaigns, a considerable amount of the 3D data of the underground spaces was obtained. A series of elaborations were processed, capitalizing on the high-resolution point clouds acquired by laser scanning.

Works were explicitly aimed at:

(a)

Extracting a digital 3D model of the limestone rock-mass;

(b)

Integrating the traditional geo-mechanical survey by identifying the main discontinuities affecting the rock mass and dividing them into families;

(c)

Providing detailed geometric reconstruction of the collapse events that occurred;

(d)

Identifying main failure mechanisms of the rock mass in correspondence of the tunnels and underground spaces (kinematic analysis).

The high-resolution point clouds acquired with the laser scanner were also used to reconstruct the position and volumes of the main collapse phenomena that have affected the underground rooms under investigation over the years. In particular, six main events were identified, the location of which is shown in Figure 16.

The corresponding volumes were calculated using laser scanner data for each collapse, considering the detachment niches and accumulations (Figure 17).

In particular, collapse 1 (0.4 m³), collapse 2 (1.2 m³), collapse 3 (1.7 m³), collapse 4 (6.0 m³), collapse 5 (1.5 m³), collapse 6 (7.0 m³).

Figure 18 shows the results of the graphic processing relating to the quantitative kinematic analysis of the rock mass understudy, as described in Section 2.3. The most probable failure mechanisms are wedge failure and, subordinately, planar sliding by observing above mentioned figures. The maximum values of the hazard indices reach 35% for the individual mechanisms. These data testify to a high propensity for instability of the investigated tunnels and allow us to identify which sectors of the rock mass are most likely to generate phenomena of kinematic instability.

4. Discussions

By applying point clouds segmentation, remotely extracted orientations of planes were measured. Orientations of joint families modal plans were identified (BG: 13°/245°, JN1: 76°/069°, JN2: 76°/001°, JN3: 75°/125°) and a stereographic projection was populated as presented in Figure 12.

From the comparison between the stereographic projections (Figure 14 and Table 2), a good correspondence between the joint families measured by scan line survey and the ones extracted applying the semi-automatic method is noticeable. The discrepancy stays within 10% for the studied areas. Furthermore, a good correspondence between the main collapse and main joint families is evident, as Figure 16 shows.

The definition of pre-assigned size value to the searching cube of the DiAna3D algorithm to efficiently perform the semi-automatic family extraction at a different scale and with a different desired detail could be an alternative solution to implement.

There are many remote sensing applications to support restoration projects of cultural heritage sites, including laser scanning and photogrammetry [50,51,52,53,54,55]. Performing considerable quantities of measurements, including inaccessible areas with fewer efforts than manual measurements, is one of the main advantages. It allows for drastically reducing the time it takes to perform the survey and enables a broader and more detailed view of the context. Moreover, the reduction of operators’ time exposure to risk-prone areas and the capability to acquire valuable data from inaccessible spaces completes the benefits list.

However, to survey the whole area, a group of 3 operators spent more than three days. Considering the relatively small extension, the long time spent appears to have resulted from the highly complex surface’s morphology and the cramped spaces. Although static TLS is an appropriate technique for this kind of environment, a mobile and presumably remoted operating system, e.g., a TLS mounted on a robotized platform, would have performed the same in less time. Investigating results obtained by applying other segmentation methods [56] to the same point cloud dataset could be further researched to compare different performances and identify the most suitable method.

5. Conclusions

A laser scanning survey was effectively exploited to remotely perform rock mass characterization in the presented context, a conservative restoration project for a high-value cultural heritage underground area.

As the studied area has inaccessible portions, two different types of survey were carried out, the traditional manual geo-mechanical survey and the remote geo-mechanical survey.

Main joint families were extracted by applying the semi-automatic classification algorithm DiAna3D. A dataset of the whole investigation area’s joints was defined by integrating results from the traditional geo-mechanical survey and semi-automatic classification.

The acquired high-resolution point clouds have made it possible to extract the discontinuities that pervade the rock mass dividing it into small blocks. This process, combining the results of traditional geo-mechanical surveys, has allowed us to characterize the rock mass from a quantitative perspective.

The traditional geo-mechanical survey effectiveness (scan-line survey) is minimal since the investigated area is inaccessible or unreachable. As in the proposed case, it is very articulated and intensely fractured.

Performing a remote geo-mechanical survey was allowed to carry out a high-quality geo-mechanical rock-mass characterization quickly compared to the traditional method and maintain the survey activities under an acceptable level of risk exposure.

Comparing results from traditional scan line surveys and the remote one, in terms of joint family orientations, points out the remote method’s good reliability, even in narrow and unconventional environments as relatively small underground tunnels.