Rocks naturally begin as lateral bedding of sedimentary layers (Oberender and Plan, 2015). However, under the action of tectonic stress, rock folds and fractures, resulting in discontinuous rock mass structure separated by weakness planes, i.e. joint, fault and bedding surfaces (Assali et al., 2014, Oberender and Plan, 2015). Geological structures influence the shape of the landscape and determine the degree of gravitational hazards along steep mountainside and rock slopes. However, the essential geotechnical information required to determine potential slope instability is not always visible on the Earth’s surface, except for the region where the outcrop is exposed through the processes of erosion, landslide, or at rock cuts along transportation corridor (Meisina et al., 2006, Alejano et al., 2010, Sacchini et al., 2016). Planar characteristics of the rock convey a wealth of information about the state of the rock. Measurements of strike and dip are used to represent the two-dimensional (2D) maps, and rocks are oriented in three-dimensional (3D) space in order to describe the directional aspects of tilted bedding planes and faults (Roncella and Forlani, 2005). Usually, geotechnical engineers determine how deformations and numbers of folds in rock occur by exploring the landscape to gather and record the orientation information of strike and dip of rock mass with compass-clinometer, scanline survey and measuring tape (Turanboy and Ulker, 2008, Lee et al., 2013, Assali et al., 2014). The traditional approach is still commonly used, however it has limitations as the process of data gathering relies on visual examination and manual measurement which depends on the experience and judgement of the engineers (Roncella and Forlani, 2005, Lee et al., 2013). This often leads to biased results and consequently poor assessments of hazard susceptibility rating and potential hazard of unstable rocks (Fekete et al., 2009, Sturzenegger and Stead, 2009).
Evolution of remote sensing imaging technologies (close-range photogrammetry and terrestrial laser scanning (TLS)) has made it possible to acquire 3D data of the terrain surface with high accuracy in less time and over wider area coverage (Gratchev et al., 2013, Kim et al., 2013, Assali et al., 2014, Matasci et al., 2015). They are also capable of surveying inaccessible steep topography and rock slopes safely. This technology development has also changed the way that geotechnical engineers collect and interpret information about rock geostructure for slope stability assessment (Sturzenegger and Stead, 2009, Wang et al., 2014, Matasci et al., 2015). In addition, there has been growing efforts to develop efficient computer programs to map rock discontinuity from 3D models, such as Split-FX (Slob et al., 2005), COLTOP (Jaboyedoff et al., 2007) and PlaneDetect (Lato and Voge, 2012), and also from 3D point cloud (e.g. FACET (Dewez et al., 2016)). These parallel developments in hardware and software capabilities have accelerated solution-focused geotechnical applications of close-range photogrammetry and TLS on the Earth’s surface over the years.
Technological advancement in TLS presents a unique opportunity that has widely attracted multidisciplinary scientific studies underground, where photogrammetric technique is rarely possible due to darkness (Gikas, 2012, Canevese et al., 2013). Since early 2001, laser scanning devices have been extensively used in the subsurface environment to directly observe and map underground features (Idrees and Pradhan, 2016). The subterraneous space is unique because it displays vertical and overhang slopes, and thus provides good overview of the quality of rock outcrop, condition and characteristic of the geology. This is currently being used especially in tunnels where laser scanning has gained popularity for operational, management and geostructural applications (Lemy et al., 2006, Gosliga et al., 2006, Fekete et al., 2009, Fekete et al., 2010, Gikas, 2012, Fekete and Diederichs, 2013, Pejić, 2013, Zhang et al., 2016).
The first study to demonstrate the potential of laser scanning technique in tunnel was reported in 2006 (Lemy et al., 2006). They used Zoller + Fröhlich Imager 5003 laser scanner in tunnel at different times and determined the displacements of excavation surface by comparing the two point clouds. They highlighted the advantages of the technique with respect to access restriction, efficient time management and data completeness. Also in 2006, a similar study was conducted along RandstadRail, north of Rotterdam, where Gosliga et al. (2006) detected deformation in the bored tunnel using a statistical testing procedure. These two studies described the earliest investigations in tunnel using laser scanning for as-built 3D documentation (Gikas, 2012, Pejić, 2013) and deformation assessment (Gosliga et al., 2006).
Geotechnical applications in tunnel attracted attention in 2009–2010 following a number of laser scanning campaigns in drill-and-blast tunnels between a group of researchers from Queen’s University (Kingston, Ontario, Canada) and the Norwegian Geotechnical Institute (NGI) (Oslo, Norway) (Fekete et al., 2009, Fekete et al., 2010). These studies emphasized on geostructure and stability of the tunnels by analysing the rock mass and excavation information obtained from 3D point clouds in order to improve excavation operations. Furthermore, Lato and Voge (2012) developed automatic 3D model-based rock discontinuity extraction software, PlaneDetect, and demonstrated its capability using data collected from some tunnels in Oslo. The study continued in 2013 in two active drill-and-blast tunnels as part of the railway expansion project near Oslo. In the study, Fekete and Diederichs (2013) presented a workflow, from data collection to design outputs, for integrating Lidar-derived point cloud data and rock mass modelling in geotechnical tunnel analyses for stability evaluation. The authors summarized the Oslo tunnels experience into two application domains: (a) operational applications composed of support evaluation, scaling assessment, potential leakage mapping and evaluation of design specification, and (b) geotechnical applications such as analysis of structurally controlled overbreak, structural discontinuity evaluation, discontinuity spacing and 3D models, surface characterization, and geological feature identification.
Similarly, Gikas (2012) utilized laser scanning technique to track excavation and construction activities in two tunnel construction sites in Greece. The author accurately determined the geometry, profile, and volume of the tunnels from 3D scanned data. Again in Niksic – Podgorica (Montenegro), Pejić (2013) proposed an optimal solution for scanning along tunnels with survey control network design for inspection and as-built documentations. The author demonstrated the effectiveness of the “arbitrary georeferencing” approach using 12 tunnel sections, totalling 3.5 km along the single railway track. He concluded that the proposed approach yields optimal registration error size within allowable limits.
Unlike in tunnel where human intrusion of the rock often results in rock surface slope and discontinuity that may not always reflect the natural state of the outcrop (Fekete and Diederichs, 2013), caves display a natural state of geological rock underground. Several dozens of caves have been documented with TLS for different applications over the last decade (Dublyansky et al., 2016, Idrees and Pradhan, 2016, Fabbri et al., 2017, Hoffmeister, 2017). However, apart from the work of Lyons-Baral (2012), studies have rarely been reported on cave stability assessment. For this, this study tries to investigate the influence of rock geostructure on cave channel development and to assess stability of the cave by determining areas susceptible to different failure types using surface discontinuity data extracted from laser scanning point cloud.
Materials and methods
Gomantong hill is a tower-like limestone outcrop that sits within Sabah Forest Reserve (118°04′E and 5°32′N), some 31.4 km south of Sandakan, Sabah (see Fig. 1). The hill is a typical synclinal feature approximately 300 m in thickness above the floodplain (Lundberg and McFarlane, 2012). From geological point of view, the limestone contains dense lithic fragments of grey sedimentary rock rich in organic matters. The bedding exists as laterally thin and well-defined layers that alternate with bed of thin grey-green fossil and sometimes interpenetrated lithologically by sandstones and shales of the Oligocene Labang Formation (Noad, 2001, Lundberg and McFarlane, 2012). Evidences obtained from nanofossil dating of samples taken from marls, mudstones, and benthic foraminifera, in the region, show that the formation was formed about 23 million years ago, i.e. between the Late Oligocene and Early Miocene. The detritus particles are believed to have recrystalized during the Oligocene sediment formation, resulting in varying sizes of crystal structures in the uplifted rock (Hutchison, 2005, Lundberg and McFarlane, 2012).
From spelaeological point of view, the hill is penetrated by complex cave systems that consist of two major passages, one above the other. Entrance to the lower cave (Simud Hitam) is along one major fault which opens to the base of the hill that coincides with level of the bank of a small stream nearby, as opposed to the upper level cave (Simud Putih) which has its access situated some 85 m above the floor of the opening to the lower cave (Lundberg and McFarlane, 2012).
Laser scanning survey and data processing
In July 2014, the upper section of Gomantong cave was scanned with FARO Focus3D laser scanner (McFarlane et al., 2015). The scanner employs the phase-shift technology with 360° horizontal field of view and 305° vertical field of view. It is a high accuracy system with theoretical range noise of ±2 mm for a distance of 0.6–120 m indoor or outdoor furthered with low ambient light and normal incidence to a 90% reflective surface. Like many other common phase-based scanners, Focus3D emits a laser light in the short infrared field (905 nm) with laser power of 20 mW at measurement speed of up to 976,000 points per second (Idrees and Pradhan, 2017). During the fieldwork, scanning was achieved using the 1/4 scan resolution mode (244,000 points/s) to optimize data collection time. Scanning was done from different scan stations to ensure sufficient overlap of adjacent scans to minimize occlusions and vertical orientation sampling bias (Sturzenegger and Stead, 2009, Assali et al., 2014).
After the field work, the scans were aligned in FARO SCENE-5.5 software (www.faro.com/downloads/training/Software) using the iterative closest point (ICP) algorithm (Besl and McKay, 1992). The registration quality yields root mean square error (RMSE) of 0.0044 m (4.4 mm), which is acceptable because the value is less than the average point spacing of 6 mm. The aligned scans were sampled at every 7th point and exported as individual scan file in ASTM E57 File Format (Huber, 2011). In CloudCompare software (CloudCompare, 2016), the decimated scans were combined into a single discontinuous point cloud, followed by filtering and removal of unwanted points. Moreover, the point cloud was spatially homogenized by sub-sampling points at 5 cm interval in both axes. Before assessment of slope stability, the point cloud was divided into six sections according to the cave passage alignment (see Fig. 2).
Automatic extraction of planar surfaces
FACET, a structural geology plug-in embedded in CloudCompare, was applied to automatically extracting planar surfaces from 3D point clouds. The tool can strategically partition points into sub-cells, and then compute the basic planar objects which are progressively combined into larger polygonal surfaces based on user-defined geometric criteria of co-planarity (Dewez et al., 2016). In the FACET, point segmentation is possible with either Kd-Tree or Fast Marching algorithm; the only difference between the two is that the former utilizes irregular space division while the latter employs regular approach. But, considering the irregular shape of the conduit, we used the Kd-Tree algorithm for point partitioning.
Kd-Tree requires defining three key geometric enforcement constraints (maximum distance, maximum angle and minimum number of points). During segmentation, the tool recursively subdivides the point cloud into quarter cells down until all the points in the cell match the best fitting plane given in the root-mean-square threshold (maximum distance). When the cells have no more point left or the number of points is less than six, the subdivision stops. Thereafter, the minimum point criterion is imposed; any sub-cell with points below the threshold will be considered as a planar surface. Then the program is executed up the tree in the opposite direction during which adjacent cells are merged if they share a common dip and dip direction using the maximum angle parameter. Additionally, the normal distance between them should be less than the maximum distance. Details of the principle, function and capabilities of FACETS can be found in Dewez et al. (2016). With this tool, several hundreds to thousands of facets were extracted from each cave section, providing the spatial orientation (dip and dip direction) required for stability analysis.
Kinematic analysis for cave slope stability assessment
The first step in rock stability assessment is to perform a kinematic analysis. Kinematic analysis involves interpreting the spatial interactions of the geological discontinuities with respect to the slope to determine the likelihood of potential for the various modes of rock slope failures – planar sliding, wedge or toppling (Sturzenegger and Stead, 2009, Babiker et al., 2014). Following the method described by Goodman (1989), kinematic analyses were performed to identify potential unstable slope using RocScience Dips v6.0 (Rocscience, 2010). The mean orientation of joint set is compared with slope angle and friction angle to evaluate the potential of planar or toppling failure. In addition, it also computes all potential plunge directions of intersection lines to estimate the potential for wedge failure (Admassu, 2012). The orientation data were plotted on stereonet, and the pattern of orientations and spatial interactions between discontinuity sets (pole clusters) were analysed. Assali et al. (2014) reported that discontinuity density influences the rate of slope erosion. Therefore, we grouped the surface discontinuities into joint sets by observing pole concentration using contour plot. According to Rocscience (2010), pole density less than 4% should be regarded with suspicion unless the quantity of data involved is relatively large. In this work, pole density less than 1% was not considered because cave section contains several hundreds to thousands of poles.
Based on Markland’s test (Admassu, 2012), the probability of any failure mode occurring depends on the relationship between the geometry of slope and friction angle (a geomechanical characteristic). According to Rocscience (2010), planar failure is likely to occur when a discontinuity dips in the same direction (within 20°) with the slope face, and at an angle inclined less steeply than the slope angle but greater than the friction angle. Also, toppling failure may take place when a steeply dipping discontinuity is parallel to the slope face (within 30°) and dips into it (Hoek and Bray, 1981). A wedge failure may be induced when the line of intersection of two discontinuities, forming the wedge-shaped block, plunges in the same direction as the slope face and the plunge angle is less than the slope angle but greater than the friction angle along the planes of failure. The kinematic analysis was performed to analyse the potential for plane, toppling and wedge failures using friction angle of 35° assumed (Babiker et al., 2014).
Results and discussion
Planar facets and orientation
The outcome of segmentation is 2D polygonal planar facets that carry important information about each plane such as orientation (dip/dip direction), quality (e.g. RMSE) and geometric properties (horizontal extent, vertical extent, surface extent and surface area). The result can be exported as Shapefile and/or CSV delimited text file into any geospatial information system (GIS) or geological analysis software. In the present case, the point cloud was segmented according to the passage structure (see Fig. 2) and each section was analysed independently. Fig. 3 shows one of the planar facets extracted from point cloud of the Main hall.
A large number of facets were extracted from each segment. A general overview of each section of the cave is presented in the stereogram (see Fig. 4). It can be observed that the pole plot appears patchily due to discretization of points into planar facets. The solid red line is the strike while red dotted lines indicate slope aspect. Important information useful for examining the relationship between the direction of cave passages and rock structure is the mean surface orientation and best–fit plane of each segment.
Rock structure and fold analysis
Orientation-based geological data (mean dip/dip direction) obtained from the extracted facets provide information for general assessment of rock folding. Thousands of planes were extracted from the point cloud of the respective cave section. The mean orientation of the six cave sections shows that the rock dips between 26° and 80° and the dip direction is from 75° to 322° (see Table 1). This variation in dip angles also reflects in the channel slope of each passage section which varies from gentle slope to cliff landform.
Average orientation of the respective cave section was plotted on stereonet. This allows visualizing the dip/dip directions as a network of planes on stereographic projection for interpreting the overall rock geostructure (general orientation of the folded rock and specific direction of the individual cave channel). Combination of the global best–fit planes of the entire cave (and the six cave sections) and the great circles of the mean orientation of the respective sections provides useful information to describe folded bedding planes (Fig. 5). According to Matasci et al. (2015), orientation of the fold axis normally coincides with the trend and plunge of the normal vector of the best–fit plane. The global best–fit plane trends east-west containing all beddings parallel to the pole vectors of the folded planes and perpendicular to the main east-west, west-east and south oriented great circles of the average orientations (Fig. 5a). From geological point of view, this pattern of slope facing planes indicates that the rock possibly folds inward from the east and west (and possibly south), producing radial network of rock fractures (and faults) typical of synclinal feature (Lundberg and McFarlane, 2012).
Another interesting finding is the link between the alignment of the passages and rock structure. For example, the entry passage and chamber trend north-south geographically (Fig. 5b); however, they have east and west slope facing mean orientations, respectively (Fig. 5a). Similar dip-controlled passage development is visible in Shaft 1 which runs almost parallel to the bedding and perpendicular to the slope face. It also can be noticed that the passage of the Main hall is neither parallel to the slope face nor perpendicular to it; instead, it intercepts the slope face at an acute angle, forming a wedge, east of the planes. Perhaps, the development of the Main hall conduit may have been influenced by geological bedding plane.
Shafts 2 and 3 are completely different. Both passages are nearly parallel (Fig. 5b) and slope towards the east (Fig. 5a) but trend in northeast-southwest and southeast-northwest, respectively. The interaction between the slope and rock folding structure indicates a probable cave passage initiated by lateral slip faults. Evidence of this is seen in the severely fractured roofs of the two passages and the rock piles on the floor (Idrees and Pradhan, 2017). Structural fold analysis reveals the interconnectivity of the cave passages, obviously linking the east-west trends of the Main hall that forms the largest cave (Fig. 5b). From geostructural point of view, it is evident that the development and pattern of the cave network are motivated by the rock structural characteristics (fold, bedding and faults).
Potential slope failures
The direction and pattern of clustering of discontinuities indicate where lineaments are concentrated. Such area increases water inflow and thereby accelerates weathering process; rock failure may occur if lineation is unfavourably oriented (Kohlstedt, 2006, Nassir et al., 2010). What determines if it is favourably oriented or not and which type of slope failure may occur are a function of the relationship between the dip directions of discontinuities and the slope face.
Clustering the obtained orientation sets (Table 2) allows estimating the amount of planar sliding, toppling and wedge sliding failures in each cave section. It was impossible to identify if a discontinuity set is bedding, faults or joints in the cave; however, examination of the mean dip angle of each group provides a distinction between two types of discontinuity sets: bedding and joint (fault) sets. The general characteristic of the joint structure is homogeneous across all the sections. Except for Shaft 3, each cave section has at least one bedding set (see JS1 (Entry hall), JS2 (Chamber), JS5 (Main hall), JS1 (Shaft 1) and JS4 (Shaft 2)). It can be understood that they have slope angles varying from 7° to 18° and gently slope to the west and north-west between azimuths 273°–326°. Besides, the bedding sets form a somewhat persistent roof structure as revealed in their variability statistics, indicating that the orientation data share a similar spatial bonding. The second category of joint sets has steep slope angles (48°–87°), typically of steep slope and cliff topography, without defined dip direction like the bedding sets.
The spatial distributions and significance of potential risk of planar sliding, toppling and wedge sliding in each section of the cave are presented in Fig. 6, Fig. 7, Fig. 8. With pole friction cone of 30° and range of slope dip direction of ±30° as recommended by Goodman (1989), only the cave chamber and Shaft 3 (Fig. 6b and f) are entirely free from any form of slope failure mode. The potential of planar sliding occurrence is present in Fig. 6a, c, d and e (Entry hall, Main hall, Shaft 1 and Shaft 2), although at varying degrees of vulnerability. Compared to Shaft 2 which has significantly high-risk level, particularly for the discontinuity set JS2 with 32.49% risk (see Table 3), the probability of planar sliding for the Entry hall, Main hall, and Shaft 1 are relatively low. The Main hall and Shaft 2 are primarily favourable for this type of failure because they have kinematically unstable poles with discontinuity orientations that plunge in the east-west direction of the cave channel development (Fig. 5b), providing rock mass with free space to slide if frictionally unstable (Rocscience, 2010).
There is significantly high risk of toppling failure in the cave, specifically within the Entry hall, Main hall, Shaft 1 and Shaft 2 (see Fig. 7a, c, d and e) with very steep slopes (>70°). In Fig. 7, toppling zones are defined by the region between the slip limit plane, stereonet perimeter and the lateral limits. The main joint sets that can result in high-risk toppling failure are essentially the east-facing slopes (JS3 of the entry and Main halls with 75.26% and 43.07%, respectively) and west-facing slopes (JS3 in Shaft 2 with 24.82%). The other joint sets have very low risk of toppling occurrence. There is apparent low-risk toppling in the north- and south-facing slopes of the cave sections.
The number of wedge sliding failure in the cave is considerably less than the number of other failure types. In Fig. 8, the percentage of intersecting planes within the crescent shaped area in the plane friction cone (measured 30° from the perimeter of the stereonet) to the total intersection represents the degree of risk (Rocscience, 2010). Obviously, only the Entry hall and Shaft 1 passages (see Fig. 8a and e) are at risk of wedge failure. Based on our analysis, the most important wedges are formed by intersection of JS2 – JS4 which makes one-tenth (10%) of all the intersections (Fig. 8a). Also, the intersection of planes JS1 – JS2, JS3 – JS2, and JS4 – JS2 make 50% wedge failure in Fig. 8e. Other cave sections show no risk at all; however, over time, the near rectangular intersections in Fig. 8c and f may be accounting for cave instability, particularly at the roof sections.
Analysis of the surface discontinuity orientation derived from laser scanning data provides important information for understanding the rock geostructure and its influence on the development of the cave passages. The limestone dissolved in all directions, creating a typical phreatic tube of circular or oval shape with scalloped walls and ceilings. This is obviously steered by strike and dip within the plane of the geological strata and accelerated by compression-induced forces (tensile cracking, shattering and shearing). Under the gravity of rock mass on the top of weak soluble carbonate base, the exerted forces develop contraction of the rock blocks, thus creating cave channels of different sizes. Generally, this may vary according to the intensity of the localized weathering activities since the impact of exterior factors such as running water, high relief, gravity, and the lithotype may not be constant (Sacchini et al., 2016).
The geometry of discontinuity sets with respect to the direction of the cave channels has significantly influenced the spatial distribution and orientation of potential slope failures. Kinematic tests indicated the relative severity of each section of the cave to sliding, toppling and wedge failures, highlighting the dominant failures that can undermine the stability of the cave. The location with the highest risk level is the opening section of the Entry hall with toppling probability of 75.26%, while for the Main hall, the toppling risk level is 43.07%, and for the Shaft 2, the planar sliding, toppling and wedge sliding probabilities are 32.49%, 24.82% and 50%, respectively. We noticed that all the high-risk discontinuity sets fall within steep-facing slopes, which is greater than 70° (see Table 2). It is emphasized that the results obtained are mainly based on analysis of the discontinuity orientation sets, however, other properties of the rock such as lateral joints, tension cracks, geological strength index, persistence, trace length, and shear strength (Barton, 2013) are not considered. Extensive field investigations of the structural characteristics of the rock by experienced structural geologists will improve the interpretation of the kinematic test results with respect to the degree of stability and/or instability of the cave.
There has been increasing use of TLS survey technology in caves in recent years, which has motivated broader scientific enquiries in the underground environments. Analysis of extracted rock surface discontinuities from point cloud collected in Gomantong cave using this technology allowed us to measure the rock deformation, evaluate its effect on cave network development, and assess cave stability. The point-based geostructural rock characterization is much optimized in terms of processing time and computational efficiency. Moreover, it overcomes the issues of interpolation errors and concave/convex surfaces by digital elevation model-based technique. Stereonet display of the mean surface orientation and global best–fit planes allowed us to understand the outcrop having a bedding plane that mostly dips 290° (west direction) with a major north-south strike and northeast-southwest fault, all of which influenced the development of the passage network. Also, kinematic tests pinpointed the spatial distributions, directions and probability of instability arising from different failure mechanisms.
The contribution of rock mass and hydraulic properties in the study area is not modelled in this work; therefore, the severity of the various failure modes detected in the cave can be determined through field investigations by experienced structural geologists. Combination of the cave model and surface topographic model will allow the thickness of the rock mass between the terrain surface and the passage conduits to be measured. This will allow identification of passage sections that are close to the surface and at risk of collapse through either surface depression that collects water or root of trees that penetrates rock cracks that develop into seepage causing rock decay. Stability assessment is critical to the survival of the cave to sustain resource exploitation (of edible bird-nest) and to ensure the safety of bird-nest harvester as well. Again, since the cave is becoming more and more opened to tourists, their safety must be guaranteed by identifying risk zones and taking necessary measures to protect the cave.
Source: Geostructural stability assessment of cave using rock surface discontinuity extracted from terrestrial laser scanning point cloud
Authors: Mohammed Oludare Idrees, Biswajeet Pradhan