**1. Introduction**

The identification of potential rockfall source areas is the first step in assessing rockfall susceptibility, hazard, risk, and determining rockfall disaster prevention and mitigation [1–7]. However, it is very difficult work to carry out in mountain areas, especially the steep and high-relief slopes, which are not accessible on site. Field investigation provides the most effective method to survey the distribution of potential rockfall source areas at a specific site [8]. Through field investigation, the engineering geology conditions controlling rockfall distribution, including the rock mass strength, orientation of structures, joint density, slope angle, relief, and the activity of tectonic faults, could be carefully studied on site [7,9,10]. Recent technologies, including unmanned aerial vehicles, terrestrial laser scanning, monitoring systems, photogrammetry, and point cloud analysis software tools (e.g., AgiSoft, Photoscan, and Coltop) [11–14], help researchers to acquire detailed information of the above conditions.

For the identification of potential rockfall source areas at a regional scale, the traditional field investigation is not as effective as that at a specific site because of its limited investigation scope and because it consumes much time and human resources. Hence,

**Citation:** Wang, X.; Liu, H.; Sun, J. A New Approach for Identification of Potential Rockfall Source Areas Controlled by Rock Mass Strength at a Regional Scale. *Remote Sens.* **2021**, *13*, 938. https://doi.org/10.3390/ rs13050938

Academic Editor: Paolo Mazzanti

Received: 7 February 2021 Accepted: 26 February 2021 Published: 3 March 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

some researchers have developed regional rockfall susceptibility models based on ArcGIS to estimate the spatial distribution of rockfall using the causal factors of rockfall including lithology, terrain, elevation, faults, rainfall, and earthquakes [15–17]. Based on the results of rockfall susceptibility assessments, the whole area could be zoned into different areas with different susceptible degrees, which provides useful guidance for regional land use and rockfall disaster prevention plans. Alternatively, some researchers have identified rockfall sources at a regional scale by remote sensing interpretation technologies using multi-temporal aerial photos, helicopter-based remote sensing imagery, and high-resolution digital elevation model (DEM) [5,6,18]. Subsequently, regional locations of rockfall source areas could be identified for further rockfall kinematic modeling and predicting regional rockfall hazards [19].

Rock mass strength is thought to be the basic controlling factor of slope stability [20,21]. A rock slope with a low rock mass strength fractured by different types of fractures is prone to rockfall [7,22]. Based on Culmann's two-dimensional slope stability model [23], a hillslope is susceptible to rockfall if its relief, slope angle, or both are larger than the threshold values [24–27]. This means that the relief and slope angles are a pair of important indicators that could be used to identify rockfall source areas on the slopes whose stability is dominantly controlled by the rock mass strength. Until now, previous studies rarely focused on the approaches combining the relief and slope angles to identify rockfall source areas controlled by the rock mass strength at a regional scale [5,28]. Focusing on this issue, comprehensive technologies, including helicopter-based remote sensing imagery, a DEM with 10 m resolution, images from Google Earth, and field work, were adopted in this study. Lastly, a new approach, the procedures, and the application criteria for identifying rockfall source areas at a regional scale are proposed and were applied in the study area.

#### **2. Materials and Methods**

The Wolong (WL) area of Tibet where ideal geological conditions exist to investigate the characteristics of rock mass outcrops and the distribution of rockfall source areas was selected for this study (Figure 1). To identify potential rockfall source areas at a regional scale, a new approach combining the relief and slope angles based on the rock mass strength is proposed.

#### *2.1. Study Area*

The lithology of the area is mainly diorite and granite with a small component of gneiss [29,30]. Because of the steep terrain and the widely distributed tectonic structures, the main type of slope failure in the WL area is rockfall. Based on the rockfall scars left on the steep slopes and the rockfall deposits distributed widely, rockfalls in the study area occur frequently, and the dominant kinematic failure modes are toppling and planer sliding (Figure 2).

In this study, based on a power grid project, the helicopter-based remote sensing imagery obtained in 2017 and the 10 m resolution DEM over the study area were acquired with the help of the research group State Grid Corporation of China. Using the data of the complete study area, images of Google Earth, and field work, a rockfall inventory including 235 rockfall scars on bedrocks and 109 rockfall deposits was prepared (Figure 1). The rockfall scars were identified based on the fresh bedrock color left on the scars (Figure 2). Because many of the rockfall deposits were removed by the Yarlung Tsangpo River, and some of the rockfalls are adjacent, it was not possible to track each rockfall's deposits. Rockfall deposits at the foot of the slopes were identified based on the shape of the deposit (e.g., pyramid) and the identifiable rockfall blocks (e.g., meters) left on the deposits.

#### *2.2. Framework of Identifying Potential Rockfall Source Areas*

By using ArcGIS and statistic tools, our new approach proposed in this paper includes the framework as follows:

**Figure 1.** Simplified geology map of the study area. Two zones are separated by the space partition for validation of the prediction result, with (**a**) the back calculation zone and (**b**) the verification zone. (**c**) Sample area presented in Figure 2.

(1) Calculate the rock mass strength parameters (*c*1, *ϕ*1) of the bedrock slope at the landscape scale and build the relief–slope (R–S) angle relationship curve in the study area (Figure 3a).

(2) Measure the present relief (*H*) and slope angle (*β*) of each specific area (*A*) (Figure 3b) with the potential rockfall over the slope areas.

(3) Calculate the limit relief (*Hc*) of the specific area (*A*) by the Culmann model [23] (Figure 3a). The Culmann model indicates that the relief of the slope is controlled by the rock mass strength, and the slope angle (*β*) has the following relationship with the limit relief (maximum height) (*Hc*):

$$H\_c = \frac{4C}{\gamma} \frac{\sin \beta \cos \varphi}{[1 - \cos(\beta - \varphi)]} \tag{1}$$

where *c* is cohesion, *γ* is the bulk density, and *ϕ* is the internal friction angle.

(4) Compare the present relief (H) and the limit relief (*Hc*) of the specific area to estimate its state of stability. When the present relief of the bedrock is larger than the limit relief (Figure 3b), the bedrock is prone to generate rockfalls. Hence, the area whose relief exceeds its limit relief is identified as the potential rockfall source area.

(5) Using the procedure from steps (1) to (4), all areas of the slopes over the study area are searched, and their states of stability are estimated. Eventually, all the potential rockfall source are identified.

**Figure 2.** Samples of the rockfall scars and deposits on the helicopter-based remote sensing imagery with the location c in Figure 1b and the historical rockfalls with kinematic failure modes of toppling and planer sliding and their deposits at the foot of the slopes.

**Figure 3.** Sketch map describing the approach. (**a**) Slope angle versus relief for hillslope. (**b**) Sketch map of estimating the state of stability for the specific area (A).

Different from the availability of rock mass strength tested in the laboratory, the rock mass strength of bedrocks at the landscape scale is very difficult to test on site because of the lack of suitable approaches [20,21]. Schmidt and Montgomery [24] proposed an approach to estimate the rock mass strength parameters (*c*, *ϕ*) using the relief and slope angles of historical rockfalls with the Culmann model. Previous studies have applied the approach to calculate the parameters of rock mass strength using the relief and slope angles of historical landsides or rockfall scars in some cases [25–27,31]. In reference to previous studies, data of the relief and slope angles of 235 historical rockfall scars were first extracted by ArcGIS. Then, the parameters of the rock mass strength (*c*, *ϕ*) at the landscape scale in our study area were calculated under the precondition that the bedrock relief is controlled by the rock mass strength.

#### **3. Results**

The geometrical characteristics of historical rockfalls were analyzed to define the specific area (A) used in the new approach for searching rockfall source areas. Applying our proposed approach, the rockfall source areas were identified in the study area and were zoned into three susceptibility classes.

#### *3.1. Geometrical Characteristics of Historical Rockfalls*

Before calculating the relief and slope angles of each specific area, a suitable value of the specific area (A) should be first determined. To define the value of A in this study, the geometrical characteristics of 235 historical rockfalls were first analyzed by ArcGIS. The values of the relief and slope angles and the areas of historical rockfalls were measured separately. According to the statistical results, the relief of historical rockfalls is mainly distributed between 40 and 130 m, the slope angle is generally larger than 45◦, and the area of each historical rockfall scar is generally less than 9000 m2. All three groups of data show Gaussian distribution characteristics (Figure 4). The mean area of historical rockfalls is 5217 m2, which was adopted as the specific area (A) for searching rockfall source areas of the slopes over the study area.

(**a**). Relief frequency statistics

(**b**). Slope angle frequency statistics

**Figure 4.** *Cont*.

(**c**). Area frequency statistics

**Figure 4.** The geometrical statistics of historical rockfalls.

#### *3.2. Locations of Identified Rockfall Sources*

A clear inverse relationship between the relief and slope angles of 235 historical rockfalls (Figure 5) enabled us to calculate the rock mass strength at the landscape scale in this study. According to previous studies [25–27,31], the minimum and maximum parameters of rock mass strength can be estimated using the lower envelope and the upper envelope of the R–S curves obtained by data of the relief and slope angles of failed slopes. The upper envelope of the R–S curve (Figure 5) represents the maximum strength of the rock mass, and the lower envelope represents the minimum.

**Figure 5.** Rock mass strength fitting result and the relief–slope (R–S) angle relationship in the study area.

Applying the results of 235 historical rockfalls to Equation (1) by the optimization algorithm (Figure 5), the minimum and maximum rock mass cohesions (*c*) in the study area are 28 Pa and 270 kPa, respectively, and both internal friction angles (*ϕ*) are 23◦. Using the upper envelope and the mean value of the upper envelope and the lower envelope, with the lower envelope corresponding to each slope angle, three R–S relationship curves were built as the threshold to determine if each specific area is stable or unstable (Figure 5). The R–S curves could be regarded as three different estimates: the aggressive estimate, the moderate estimate, and the conservative estimate, corresponding to the upper envelope, the mean value, and the lower envelope curves, which could be selected by different aims or rockfall disaster prevention and mitigation strategies.

Based on our approach, to identify the rockfall source area is to identify the area whose present relief exceeds the limit relief corresponding to its slope angle (Figure 3). Hence, we calculated the limit relief corresponding to each slope angle (Table 1) using the three R–S curves presented in Figure 5. Considering the actual slope distribution of the historical rockfalls (Figure 4b), we mainly focused on the slopes with angles larger than 45◦. Comparing the present relief and the limit relief in each specific area of the slopes in the study area, we obtained the rockfall source areas in ArcGIS (Figure 6). More rockfall source areas are distributed at the lower parts of the slopes whose slope angles are relatively bigger, which is probably affected by the intense incision of the Yarlung Tsangpo River [32].

**Table 1.** Limit relief of the specific areas corresponding to the slope angle.


Note: for the slope angle, the front part is open interval, and latter part is closed interval. For example, 46 represents the range of (45, 46]. C-U: R–S curve by upper envelope; C-M: R–S curve by mean value; C-L: R–S curve by lower envelope.

**Figure 6.** *Cont*.

**Figure 6.** The identification results of the rockfall source areas in the study area. The results of (**a**–**c**) correspond to high (I), medium (II) and low (III) rockfall susceptibility classes, respectively.

#### *3.3. Zoning Map of Rockfall Susceptibility*

In rockfall risk analysis, the rockfall susceptibility assessment is the first step to carry out [16,28,33]. The rockfall susceptibility map helps to highlight the spatial distribution of potentially unstable slopes [34], which is usually zoned into different susceptibility classes, e.g., high, medium, and low susceptibility, to represent different rockfall susceptible degrees of slopes. Based on the sketch in Figure 3, if *H* > *Hc*, the area A on the slope is unstable, which means it is prone to rockfall in the future. The larger H is than Hc, the higher the possibility of rockfall in area A, and hence the higher rockfall susceptibility of area A. To provide a reference for the susceptibility study following the identification of the rockfall source areas by our approach, the rockfall source areas were zoned into three susceptibility classes, high (I), medium (II), and low susceptibility (III) areas, using the upper envelope, the mean envelope, and the lower envelope in Figure 5. In this way, the regional rockfall susceptibility maps were produced in the study area (Figure 6).

The slope angle and the elevation of all the susceptibility classes were analyzed, and their frequency distributions were obtained. According to the results (Figure 7), the rockfall source areas with different susceptibility classes have different ranges of slope angle. The rockfall source areas within the high susceptibility class are mainly distributed on the slopes with the angles of 60–66◦, those of medium susceptibility are distributed on the slopes with the angles of 54–61◦, and those of low susceptibility are distributed on the slopes with the angles of 46–55◦. However, based on the distribution statistics of elevation in Figure 7, no obvious relationship between elevation and susceptibility classes was observed in this study.

**Figure 7.** The distribution of slope and elevation of the identified rockfall source areas within different susceptibility classes.
