Abstract
Our understanding of debris-flow initiation by slope failure is restricted by the challenge of acquiring accurate geomorphic features of debris flows and the structural setting of the rock mass in the remote mountainous terrain. Point cloud data of debris flows in Sabino Canyon, Tucson, Arizona, July 2006, with initiation by joint-controlled rock slope were obtained using multitemporal LiDAR scanning. Topographic changes were detected by comparing historical LiDAR scanning data of this area since 2005 by adopting open-source CloudCompare software. The results showed persistent scour and erosion in the debris flows after 2006. Point cloud data of joint-controlled rock in the initiation zone were generated by the means of photogrammetry using Pix4D software. The joint planes, the dip direction and the dip value of the joint plane, the joint spacing, and the joint roughness were therefore acquired by point cloud processing. Our study contributes a foundation for analyzing the relationship between the rock features, the generation of slope failure, and the initiation of debris flows.
1. Introduction
Debris flows that are often initiated in the mountainous terrain and triggered by heavy rainfall are significant geologic hazards that can lead to catastrophic damage across large regions [1–3]. Mitigation and early warning of the hazards posed by the debris flows require an understanding of the mechanisms leading to their initiation [4, 5].
Geologists have long recognized that debris flows were initiated by a slide, debris avalanche, or rockfall from a steep bank (>20°), synchronously or lag behind rainfall [1, 6, 7]. The instability of rock slope projected rock fragments directly into the steep upper channels, providing both the initial material and steep gradient for the debris flow, making debris flows heavy in terms of volume and weight and thereby leading to its rapid occurrence [7, 8].
The processes for the mobilization of debris flows from slope failure include the Coulomb failure within the sloping mass, liquefaction, and high pore-fluid pressures of the mass by the rapid undrained loading of rockslides and the conversion of the landslide translational energy to internal vibrational energy [1, 9, 10]. Gradient, angle of entry of failure into the channel, and the amount of in-channel stored sediment are regarded as the parameters that determine the initiation of debris flows by rock-slope failure [9, 10].
In order to figure out the initiation of debris flows, field observation, empirical models, experimental modeling, statistical and simulation models, and hydromechanical models have been adopted [11, 12]. However, these models are data-dependent. It always requires the estimation of a larger number of parameters, as well as knowledge of the initial conditions and boundary conditions [3]. Particularly for rockslide debris flow, the joints have profound effects on the rock-slope stability [13]. Joint parameters, that is, joint spacing, joint length, joint roughness, and joint orientation, influence the weathering velocity, the size of blocks, the rock strength, and the fissure-fluid pressure distribution during heavy rainfall events [14–16]. Subtle differences in joints can lead to significant differences in the erodibility and rock-slope stability, thus leading to the initiation of this type of debris flows [17].
However, challenges still exist in obtaining accurate and detailed features of the debris flows and rock mass in a safe way, since the debris flows usually originate in the remote mountainous terrain and/or steep rock slopes [17, 18]. Conventional survey methods present serious limitations for collecting the datasets required for the analysis and modeling of the debris flows. Therefore, technologies like the aerial photography, LiDAR (Light Detection and Ranging), and photogrammetry technologies are applied to extract the data in an unbiased and secure way [2, 19, 20]. Point cloud data generated in these ways could then be processed to extract the necessary parameters for rock-slope failure-initiated debris flows.
In the present work, the historical LiDAR datasets before and after the debris flow event in 2006 in the study area were used for detecting topographic changes. Point cloud data of the rock in initiation zones with joint sets generated by photogrammetry were processed for obtaining the detailed joint parameters. The aim of this research was to demonstrate the effectiveness of the LiDAR and point cloud processing in acquiring the geomorphic features of debris flows and the joint parameters.
2. Study Area
2.1. Setting
The study area is located on the Santa Catalina Mountains roughly 15 km northeast of Tucson, Arizona, the United States (Figure 1).
The Santa Catalina Mountains are a metamorphic core complex typical of the Basin and Range Province of North America, and the bedrock of the southern half of the range is almost entirely granitic [21]. Elevation ranges from 805 m at the mountain front at Sabino Canyon to 2800 m at Mount Lemmon [22]. The Santa Catalina Mountains are bounded on the south by the low-angle Catalina detachment fault and on the west by the high-angle Pirate fault. Offset along these faults occurred in two separate intervals of deformation and uplift. The first extension along an ∼240° azimuth was accompanied by tectonic tilting of an extension-parallel topographic ramp and by antiformal arching along a direction approximately orthogonal to extension. The second extension was oriented more nearly east-west compared to the earlier phase of extension [22].
The climate is semiarid across most of the Santa Catalinas, with mean annual precipitation ranging from 330 mm at the mouth of Sabino Canyon to 750 mm on Mount Lemmon. About 45% of rainfall falls during the summer monsoon season of July through September, typically in convective storms [23].
2.2. Debris-Flow Description
In July 2006, a week of extreme precipitation (recurrence interval > 1000 years for 4-day precipitation) in the study area caused 435 slope failures and spawned numerous damaging debris flows [22, 23]. Sabino Canyon, a heavily used recreation area, located at the southern flank of Santa Catalina Mountains (the drainage area is estimated to be 89 km2), had the most slope failures (Figure 1) [24]. At least 13 debris flows occurred in this area and resulted in the creation of the steep chutes down to the Sabino Creek. Numerous large boulders were entrained, adding mass to debris flows and compounding damage of the roads, bridges, and structures in Sabino Canyon including the destruction of structures and the roadway, and resulted in closed public access for months [24].
Digital orthophotography orthophotos produced by the Pima Association of Governments’ (PAG’s) Regional Remote Sensing Program show the debris-flow activities in Sabino Canyon (Figure 2).
There were more than fifty slope failures occurring in lower Sabino Canyon (Figure 2). Many of them aggregated debris-flow tracks and resulted in the creation of the steep chutes down to the Sabino Creek.
Debris flows were joint-controlled and oriented parallel and perpendicular to the extension-related joint sets and geomorphology in a range of settings (joint control on network geometry) [25].
Debris flows were typically initiated on a steep slope before developing into a rapid flow confined by a steep channel and ultimately deposited the material downstream on the Sabino Creek (Figure 3).
The topography of most slope failure surfaces was hollow chute. Colluvium and diluvium on these chutes, that is, debris, had already been removed by rain, flowing along the narrow (<4.5 m wide) and relatively long (about 457 m) chutes, all the way to the canyon bottom (the Sabino Creek). Bedrock was exposed after the debris migration process; therefore, the mean failure depth was consistently about the depth of colluvium and diluvium, 0.8 ± 0.4 m [24].
2.3. Rock Cliffs in Initiation Zone
The initiation zone of the debris flows was near vertical bedrock outcrops of the Wilderness Granite [21] and consisted of central core granitic intrusion flanked with a metamorphic core complex on the margins [25]. Steep talus slopes below these outcrops were covered with thin colluvium [24]. Each debris flow had one or multiple initiation points, finally gathered in Sabino Canyon (Figure 4).
The slope failures from July 2006 occurred between elevations of 1220 and 1830 m, and the failure surface in many cases was the bedrock.
Two nearly orthogonal, steeply dipping structural joint sets (WSW-ENE trending and SSE-trending) incised the bedrock into blocks (Figure 5) that mobilized the sediment and contributed to the initiation of the debris flows and also to the enormous potential and kinetic energy of the debris flows, leading to the scouring of material from the hillsides, channels, and downstream channels [23].
3. Change Detection of Debris Flows
After the unusual extreme rainfall event in late July 2006, the debris-flow activities in Sabino Canyon subsided, but it never meant being incapable of its recurrence, because of the unpredictable weather conditions. A little change of topography might attribute to the debris flows [26]. Monitoring the topographic changes in debris-flow areas was important to indicate the likelihood of slope failure and debris flows. LiDAR technology has the potential to precisely identify and quantify the topographic change [27], which has been applied in investigations to rapidly provide detailed topographic models [20, 28]. Therefore, we tried to precisely quantify the change in Sabino Canyon, by processing the point cloud generated from LiDAR scanning.
3.1. Point Cloud Data Preprocessing
Four high-resolution LiDAR datasets were collected by the Pima Association of Governments’ (PAG’s) Regional Remote Sensing Program before and after the debris-flow event, in 2005, 2007, 2011, and 2015, which imaged a 2.6 km2 swath of the ground covered study area.
The datasets were high-density collections (greater than or equal to 1 pulse per square meter (pls/m2)). Aggregate nominal pulse density (ANPD) was used to represent the point density of LiDAR datasets measured in the areas of Sabino Canyon Road, where the datasets were relatively flat and unchanged. The ANPDs of LiDAR dataset were 1.6301 pls/m2, 2.0517 pls/m2, 2.4804 pls/m2, and 3.0776 pls/m2, respectively. ANPDs of LiDAR data collected in 2007, 2011, and 2015 were up to the minimum acceptable quality level 2 (ANPD ≥ 2.0 pls/m2) [29].
The validity of comparing the topographic change of different LiDAR datasets was checked by calculating the error range of z values of a different LiDAR database in an area with a perfectly flat level surface, since dust could cause erroneous measurements of LiDAR data. We selected an overlapped ten meters long hard surface road in lower Sabino Canyon Road as the surrogates with a relatively subdued topographic change. The mean difference value ± 1 standard deviation was selected to be an error range (Table 1). A ±0.30 m overall error range was selected to represent the differenced value, which meant that when differentiating the LiDAR datasets, elevation changes less than ±0.3 m were not considered. Comparison and change detection were performed based on these datasets.
3.2. Change Detection
Following the quality testing of the data, we performed the topographic change analysis by means of the open-source CloudCompare software (version 2.11. Alpha, 2019). The CloudCompare software is a 3D point cloud (triangular mesh) editing and processing software, which supports various point cloud processing algorithms and provides a good comparison between the point cloud datasets.
The mesh to mesh distance computation in the CloudCompare is a unique way to compute the signed distances directly between the reference mesh and the compared mesh, which calculates the nearest neighbor Euclidean distance between the two meshes. This is regarded as a reflection of the topographic change.
The three steps that were taken before distance computation are as follows:(1)Importing and denoising the point cloud data: LiDAR datasets were imported, respectively, into CloudCompare software to clean the noise by using filtering scale value. We retained the vegetation instead of removing them because of the imprecise vegetation removal at present [30, 31]. Furthermore, the growth rate of the major vegetation in the study area, Saguaro, was less than 0.3 m in 10 years, which could be ignored compared with the accuracy of LiDAR datasets.(2)Meshing point cloud data and reconstructing slope surfaces: the extensive point data of each year were meshed by implementing the triangulated meshing tool to generate the triangulated irregular networks (TINs) by reconstructing the 3Dslope surface.(3)Registering and computing: two data to be compared were registered by applying the ICP algorithm, which aimed to minimize the Euclidean distance between the two 3D point clouds. The signed distance between every two meshes could be computed after registration to align with the real-life coordinates and for implementing a distance computation algorithm.
The topographic changed maps using the LiDAR data recorded the land surface before and after event is shown in Figure 6, where the previous year was set as reference data and the subsequent year represented a compared data. The positive distance meant that the compared mesh was outside the reference mesh and the negative distance meant that the compared mesh was inside the reference mesh. Color scale referred to the range of signed distance and the bulge on the right of the color scale represented the distribution of computation results. 700 represented the plotting scale, meter.
(a)
(b)
(c)
(d)
Compared with the large spatial scale, the change was quite insignificant. Most of the changes between the two meshes were concentrated on a small interval ranging from −1 m to 1 m, since the colluvium covered on the slopes was thin (<1 m). In order to clearly observe the topographic change relative to its area, we edited the color scale to concentrate it on an interval ranged from −1 m to 1 m to show the primary topographic change (Figure 6). The distance computation results exhibited evident changes before and after the debris-flow events that occurred in 2006, showing an extensive positive change in 2007 compared with 2005 Compared with 2007, a constant negative change occurred later, and a relatively bigger negative change happened during the period from 2007 to 2011 (Figures 6(b) and 6(c)), whereas a relatively inconspicuous change was exhibited from 2011 to 2015 (Figure 6(d)).
The topographic changes in debris-flow areas from 2006 to 2015 are illustrated in Table 2.
Debris-flow events that occurred in 2006 delivered the debris as well as the thin colluvium covered on the slope (around 1 m) from high elevation to low elevation, which exhibited a relatively flat gradient, thereby providing a deposit place for debris flows and leading to positive changes. Erosion of debris flows took place after 2007 due to the continuous weathering of the bedrock and the velocity of erosion was high during 2007 and 2011 and slowed down after 2011.
The largest and persistent positive or negative distance occurred along the Sabino Creek that was especially concentrated at the snout of the chute of debris flows, which was because the Sabino Creek provided a natural deposition place for all the debris flows as well as a scour and transportation for the debris.
At the initiation points, the rock cliffs of the debris flow in Sabino Canyon, especially the ones at upper and middle regions, had already reached the ridge, making a gentle gradient of the rock slopes on the initiation points. However, the ones located at lower Sabino Canyon had a steeper gradient, which was almost vertical, where most slope failures and debris flows occurred.
To show a more detailed change for all parts of debris flows, we segmented the point cloud in debris flows marked in Figure 2 and compared the change from 2007 to 2015 because of the fragmentary data of that area collected in 2005, as illustrated in Figure 7.
(a)
(b)
On comparing the changes in the detection results of the data collected in 2011 and those collected in 2007, a more visible change was seen and a stronger erosion occurred at the initiation area at the east side of the chute, about 1.2 m to 1.5 m, but a relatively positive change existed outside the debris-flow area due to probable deposition of the debris (Figure 7(a)). After 2011, the topographic changes seemed to be slowed down, with a relative change distance ± 1 m, and mainly concentrated in a range less than 0.3 m (Figure 7(b)).
A noticeable relative visible negative change occurred at the initiation points as well as the chutes of debris flows, which showed the continuously unstable condition of the rocks. A little change in this could cause the probability of the formation of debris flows.
From the field observations, it was found that the bedrock exposed on the steep rock cliffs was quite joint-controlled. Two joint sets incised the bedrock in the study area, which affected both the stability of the rock cliffs and the initiation of the debris flows. However, the historic LiDAR data were comprised of the topographic changes of the debris flows in an unbiased, rapid, and accurate manner from a remote distance. However, it lacked the detailed information and accuracy about the characteristics of the rocks in the initiation, which was necessary for analyzing the change of debris flows and the rock cliffs.
4. Joint Features of Rock at Initiation Zone
In the present work, we employed the photogrammetry technology to acquire a point cloud of rock, which has been widely used in geological disaster monitoring [31, 32]. The joint parameters were obtained using the CloudCompare software, which reduced the amount of manual labor involved and enabled us to use a mobile phone combined with a camera instead of the costly 3D laser scanner.
4.1. Point Cloud Data Generation by Photogrammetry
A rock block sample, closing to the initiation point of debris-flow chute, is shown in Figure 6(b) containing groups of joint sets selected in a talus rock slope (Figure 8(a)).
(a)
(b)
166 images were then taken using a mobile phone with the mobile phone built-in WGS 84 (EGM 96 Geoid) coordinate system. Subsequently, the images were uploaded to the photogrammetry Pix4D software (https://www.pix4d.com/), and the images were firstly calibrated with a 2.99% relative difference between the initial and optimized internal camera parameters. 62148 key points were then collected per image with an average 0.08 cm/1 cm ground sampling distance (GSD), and 26266.8 of them were matched per calibrated image. Then, the noise was filtered, and the surface smoothed point cloud of the block was densified and generated using the point cloud for processing automatically. The number of 3D densified points was 6203692 and the average density was 1028353 pls/m2, which was quite a high quality. Finally, the point cloud data were output with WGS 84/UTM zone 12 N (EGM 96 Geoid) coordinate system.
In order to acquire the joint parameters of the rock block, the point cloud file was imported to the CloudCompare software (7,080,405 points), the point cloud data were cleaned, and the vegetation was removed manually using a segmentation tool (4,559,945 points left) (Figure 8(b)).
4.2. Joint Parameters Measurement
4.2.1. Joint Identification and Joint Orientation
Joint surfaces can be detected by the means of the CloudCompare software using the “RANSAC Shape Detection” plugin, an interface to the automatic shape detection algorithm proposed by Schnabel et al. [33], which can extract shapes by randomly drawing minimal sets from the point data and constructing corresponding shape primitives.
The normal of the point cloud was computed first and 50,000 was selected as the number of samples per “plane” primitive to detect the joint face based on our cloud density and the size of the joint faces to be detected. 16 planes were detected, and 10 planes were picked up manually by removing those inaccurate distinctly. Accordingly, the joint parameters were measured or calculated. Dip value and dip direction were displayed after segmenting the rock blocks by cross section and converting the normal to dip and dip direction.
Three joint sets were identified (Figure 9), while it was noteworthy that the orientation of the two joint sets detected by RANSAC matched well with the field observations conducted in December 2018. However, the third bedding plane did not match. The reason for this could be the horizontal orientation of the bedding plane, which caused an error between the reality and the field measurement.
(a)
(b)
(c)
Detected results of the joints were compared with the field observation and displayed in Table 3.
4.2.2. Block Size
The joint spacing was measured by the distance detection between the joint planes after obtaining the section point cloud of each detection plane using the “Cross Section” tool. Joint density was approximately 1.5 times the SSE-NNW trending joint set compared to the WSW-ENE-trending joint set.
Block size was therefore calculated by computing the joint space and the results are shown in Table 4.
4.2.3. Joint Roughness
The normal vector is a three-dimensional parameter that represents a plane, which can be obtained in CloudCompare software. The roughness along a joint face in different azimuths may help characterize the roughness in a more realistic way, providing insight into the variation in shear strength along different azimuths [26, 34, 35].
In this study, we computed the dip direction and dip value of the orient normal vector of triangles formed by three adjacent points along a detected joint surface by using simple spherical statistics equations. The mean orient normal of the best-fitting plane was calculated by the calculation of the Fisher mean vector. The standard angular deviation of orientation normal between every polygon and the mean normal vector was assumed to represent the surface roughness.
The orientation of the normal was obtained and was plotted on a stereo net, wherein the scatter plot showed the point dense, and the contour plot displayed the distribution of orient normal (Figure 10).
(a)
(b)
(c)
(d)
The standard angular deviation of the orient normal along each detected joint surface could reflect the joint roughness, as shown in Table 2. Accordingly, the bedding planes exhibited the roughest surface, and the exposed joint surfaces were rougher than the fresh ones.
5. Discussion and Conclusions
The debris-flow activities in the Sabino Canyon were quite joint-controlled, and the initiation of these debris flows was related to the rock-slope failure. The point cloud data were recorded safely, were unbiased, and could be processed efficiently to exhibit the topographic condition of debris flows and detailed joint characteristics of rock slope in initiate zones. Furthermore, it presented opportunities for establishing a 3D model for further mechanic analysis and simulation. In this study, we showed that the point clouds generated by LiDAR scanning as well as photogrammetry could be processed and then be used not only in the analysis of the topographic changes in the large-scaled debris-flow activities but also in the measurement of the joint parameters of small-scaled rock block.
The LiDAR data of the study area were successfully processed, mapped, and precisely quantified for topographic changes, showing the location of sediment deposition and erosion after the debris-flow events.
Both the topographic change computation of large-scale LiDAR dataset and geometry characteristics obtaining small-scale LiDAR dataset demonstrated the capabilities and huge potential of LiDAR to precisely acquire information and assessment of hazards.
The spacing between the points and the mesh triangles is determined by scanning with different resolutions. High-resolution scanning is required to obtain more accurate data. In addition, the JRC (Joint Roughness Coefficient) is perhaps the most common empirical method for determining the joint roughness; however, visually assigning a JRC value to a joint profile is inherently subjective. The standard angular deviation between the orient normal of triangles surfaces generated by point and best-fitting plane of the joint faces is able to represent the joint roughness, and the bigger the angular deviation, the bigger the joint roughness.
The results are consistent with the actual situation, representing its feasibility for analyzing the debris flow initiated by rock-slope failure. However, the way to remove the vegetation accurately and to build 3D models of the debris flow initiated by rock slope should be done in future work.
Data Availability
The data of numerical results used to support the findings of this study can be obtained from the corresponding author upon request.
Conflicts of Interest
The authors declare that they have no known conflicts of interest that could have appeared to influence the work reported in this paper.
Acknowledgments
The authors would like to thank the Pima Association of Governments’ (PAG’s) Regional Remote Sensing Program for providing the laser scanner data for this research. This study was supported by the National Natural Science Foundation of China (51608083 and 41807272) and the Second Tibetan Plateau Scientific Expedition and Research (STEP) Program of China (2019QZKK0902).