Multi-Scale Engineering Geological Zonation for Linear Projects in Mountainous Regions: A Case Study of National Highway 318 Chengdu-Shigatse Section

Abstract

In response to the challenges of long crossing distances and difficult site selection for linear engineering projects in mountainous areas, this article proposes a multi-scale engineering geological zoning (EGZ) method. This method is based on the linear engineering construction stage and transitions from regional EGZ to EGZ of key sections (areas with poor or worst engineering geological conditions). This method not only ensures the effect of EGZ but also reduces the workload. When carrying out the EGZ of key sections, the assessment ideas of geological disaster hazards were taken into consideration. An improved method for calculating the time probability and magnitude probability of disaster occurrence is proposed. Taking the National Highway 318 Chengdu-Shigatse section as an example, EGZ was carried out. Its results revealed that the Nyingchi section was the key section with poor and worst engineering geological conditions. EGZ of the key section showed that the areas with poor and worst engineering geological conditions were mainly distributed in the curved sections on the northern side of the linear project. The proposed method in this article provides guidance for EGZ for linear engineering projects in mountainous areas.

1. Introduction

With the implementation of national strategies such as the Belt and Road Initiative and the Western Development, many large-scale projects are gradually being launched in the southwestern region of China (including the Qinghai-Tibet Plateau, Yunnan-Guizhou Plateau and the Sichuan Basin). Among them, the construction of linear projects such as the Qinghai-Tibet Railway, Sichuan-Tibet Railway, and Sichuan-Tibet Highway [1], all these projects are of great importance for the development of the provinces in the southwestern region of China [2]. Affected by tectonic uplift and river cutting, the region is prone to frequent seismic activity, with rocks being highly fractured and intensely weathered and loose deposits being widely distributed [3,4,5,6]. The combination of unique and fragile geological conditions and sufficient precipitation made the southwestern region an area with frequent occurrence of geological disasters [7], which can seriously affect the safe construction and operation activities of the major existing and proposed linear projects [8]. Therefore, to ensure the long-term normal and safe operation of the proposed linear projects, it is necessary to select a suitable alignment for the proposed linear project by avoiding areas of potentially high risk [9,10].

Engineering geological zoning (EGZ) refers to the division of sections within the project outlines on the basis of similar engineering geological conditions and similar engineering geological problems as per the requirements of engineering planning [11]. The EGZ mainly integrates the stratum and tectonics of the area, lithotypes, hydrogeological conditions, geomorphological features, weathering, anthropogenic activities, as well as the distribution of geological hazards, etc. [10,12]. Reasonable EGZ cannot only provide the basis for selecting the engineering site but also guide the terrain assessment experts about the disaster prevention of the site once it is determined; thus, EGZ can effectively reduce the occurrence of disaster events during as well as after the project construction [13]. Lazaro et al. (2004) divided the Fortaleza Metropolitan Region, Brazil, into nine engineering geological zones according to the susceptibility to geohazards and foundation, excavation, and waste disposal conditions [14]. Shang et al. (2005) selected structural geology and ground stress, lithology, topography, hydrogeological conditions, and physical (natural) geological phenomena as evaluation factors and employed the interactive matrix method to perform EGZ along a 427 km section of the Sichuan-Tibet Railway spanning from Basu to Nyingchi [10]. Osipov et al. (2012) divided the Moscow territory into four categories based on the integration of data on the structural geodynamic, geomorphologic, and geologic structures, hydrogeologic conditions, and the occurrence of hazardous natural and anthropogenic induced processes and phenomena [15]. Considering geological conditions (i.e., lithology, weathering and rock mass structure) and historical disaster data, Qi et al. (2015) produced a zonation map for a 17 km highway in Beijing [9]. Xiao et al. (2018) used ten factors, such as elevation and slope, and employed various methodologies, such as decision tree analysis, support vector machines, back propagation neural network, and long- and short-term memory methods, to conduct a comprehensive susceptibility assessment of a 133 km segment of the China-Nepal Highway from Mengla to Friendship Bridge [16]. Ali et al. (2021) adopted a regional to site-specific approach, utilizing both quantitative and semi-quantitative methods, to evaluate a 200 km section of the Karakoram Highway, extending from Besham to Chilas based on the risk of rockfalls and debris flows [17]. Yang et al. (2021) carried out EGZ of the site after comprehensively analyzing the factors affecting the stability of the site, the factors determining the stability of the foundation, the degree of groundwater influence, and the factors of the drainage conditions of the site [18].

All mega linear projects, i.e., National Highway 318 Chengdu-Shigatse section, always traversed through diverse geological and geomorphological units, offering complex geological engineering conditions all along its length, ranging from several hundred to several thousand kilometers [10]. Principally, EGZ is assumed to be carried out all along the entire project corridor, but in such cases, the accuracy of the EGZ evaluation results is often limited due to the restrictions of the regional scope, the precision of available data, and the workload. Therefore, scholars often only choose a specific section of mega linear engineering projects [9,10,16,17]. Based on the various existing problems in the current large linear projects EGZ and referring to its regional stability evaluation and geological disaster hazard assessment, this paper proposes a multi-scale EGZ method for linear projects in mountainous areas. For the planning and feasibility study stages of large-scale linear projects, the EGZ of regions and key sections is carried out, respectively. The key section is the area with poor or worst engineering geological conditions traversed by the proposed highway alignment. When evaluating key sections, refer to the idea of geological disaster hazard assessment and carry out EGZ based on the spatial probability, time probability and magnitude probability of disaster occurrence. This paper innovatively proposes to calculate the time probability and magnitude probability of the geological disaster based on the number and the average area of historical disasters in the slope unit and superimpose them with the spatial probability of disaster occurrence to obtain the EGZ results of key sections.

Based on the proposed method, this paper illustrates the multi-scale EGZ of large-scale linear projects in mountainous areas using the example of the Chengdu-Shigatse section of National Highway 318. Although the highway has been completed, the results of the EGZ can still guide the highway maintenance crew about the prevention of future disasters along the route and even for the realignment of certain sections of the mentioned highway.

2. Study Area

The Chengdu to Shigatse section of National Highway 318 starts from Chengdu, Sichuan Province in the east, passes through Ya’an, Litang, Mangkang, Bomi, Nyingchi, Lhasa, and finally reaches Shigatse in Tibet, with a total length of more than 2000 km (as shown in Figure 1). This highway is the most convenient land transportation channel for economic and cultural exchanges among surrounding provinces and regions in Southwest China. The area’s terrain is characterized by a gradual decrease in elevation from west to east, with a maximum vertical difference of over 5000 m. The section from Chengdu to Ya’an is located in the Sichuan Basin, with a hilly landscape. The section from Ya’an to Shigatse is located in the Qinghai-Tibet Plateau, with a rugged terrain of high peaks and deep valleys. The highway passes through various mountain ranges, including the Erlang Mountains, Zheduo Mountains, Hengduan Mountains, Nyainqentanglha Mountains, the Himalayas, and some major rivers such as the Yangtze River, Lancang River, Nujiang River, and Brahmaputra River [19]. The landforms are characterized by a high alpine relief with the uplifting of hills and trenching of streams [10]. In addition, due to the abundant precipitation, seasonal movement of glaciers, and the freeze-thaw action, the rocks in the study area disintegrated into fragments and blocks. Anthropogenic activities such as slope cutting and mining also worsen the natural surface processes. The superposition of the above unfavorable factors has led to the widespread development of natural disasters such as rock falls, landslides, and debris flows in the proposed project area [12,20], as shown in Figure 2.

The Chengdu-Shigatse section of National Highway 318 is situated in a tectonically active area. The alignment of the highway traverses many active mountain fold belts and some famous regional and local scale faults and shear zones, i.e., the N-E trending Longmenshan fault zone in the East, the N-S trending Jinshajiang fault zone and Lancangjiang fault zone nearly in the middle and the E-W trending Yarlung Zangbo River fault in the West [21], as shown in Figure 3. The available instrumental seismic data provides several pieces of evidence of past earthquake events that occurred in the area along the project corridor. On the national seismic intensity zoning map, the basic intensity of earthquakes along National Highway 318 is level VI–VII [22]. The lithology and strata encountered along the 318 highway corridor are quite complex and range in age from the Upper Proterozoic Sinian to the Cenozoic Quaternary.

Due to the influence of terrain and monsoon circulation, there are significant climate differences in the study area. The eastward regions of the Erlang Mountains have a mid-subtropical monsoon humid climate with high rainfall, while the westward regions of the Erlang Mountains have a plateau cold type continental climate zone with large temperature differences and concentrated rainfall from May to September while the regions near the Yarlung Zangbo Grand Canyon have a warm and humid climate with relatively high rainfall [21].

3. Methods and Data

3.1. Methods

Based on the quantum of work and different tasks, large-scale linear projects can be divided into different stages, such as the planning stage and the feasibility stage. This paper proposes a multi-scale EGZ approach, i.e., EGZ of the regional sections and EGZ of the key sections, corresponding to the mentioned different stages of linear projects, Figure 4 shows the generalized flowchart of the multi-scale EGZ method. The EGZ of the regional sections corresponds to the initial planning stage of linear projects. At this initial planning stage, a regional-level EGZ (less precise EGZ) can assist in the process of selection of relatively suitable alignment for the linear projects. The latter, i.e., the EGZ of the key sections, corresponds to the feasibility engineering stage of linear engineering projects. Here, key sections are the areas along the linear projects having poor and worst engineering geological conditions. Therefore, for all the key sections along the linear project, a more precise EGZ should be developed on a local scale. The precise and correct development of EGZ always helps project managers in selecting an area of relatively stable engineering geological conditions, which then ultimately reduces the possibility of geological disasters like landslides and rock falls along the project corridor.

3.1.1. Regional Engineering Geological Zoning

Regional EGZ is the process of carrying out EGZ at a regional scale, i.e., dividing large areas into zones based on similar engineering geological conditions and issues and then assessing the engineering geological condition grade, similar to regional stability evaluation. Therefore, the work content and methods of regional stability evaluation can be used as a reference for regional EGZ [23]. Regional stability evaluation refers to the stability of the current crust and its surface within a certain regional scope under the comprehensive action of internal and external dynamic forces and anthropogenic activities, as well as the interaction between this stability and engineering construction [24]. It can be considered that areas with poor engineering geological conditions are more prone to geological disasters.

During the development of regional EGZ of linear engineering projects, all alternative routes should also be included in the outlines of the project corridor, and then during evaluation, significant factors such as geologic structures, topography, geomorphology, and seismic activity that can affect the engineering geological conditions in the study area should be considered. Grid units with a side length of 30 m can be used as evaluation units. Regional EGZ will be carried out through the qualitative and semi-quantitative Analytic Hierarchy Process (AHP) method [18,25,26]. This method was first proposed by Saaty (1980) [27] and is considered a commonly used approach that combines qualitative and quantitative analysis. This particular method decomposes the objective or criteria from top to bottom and determines the weight of each factor by solving the comparison matrices. It has been widely applied in geological disaster hazard assessment [28,29,30]. Due to its diverse applications, this paper will not go into further detail.

Based on the AHP method, the weight values of each evaluation factor can be calculated and then combined with the level values of each factor so that the regional EGZ index for each grid unit can be obtained. The formula is as follows:

$$RZI=\sum _i^n{W}_i\ast {\mathrm{C}}_{ij}$$

(1)

where $RZI$ is the value of the regional EGZ index, $W_i$ is the weight of the $i\mathrm{th}$ factor, ${\mathrm{C}}_{ij}$ is the level value of the $j\mathrm{t}\mathrm{h}$ class of the $i\mathrm{th}$ factor.

On the basis of the regional EGZ index, the entire project area under investigation can be divided into four regions: (1) good engineering geological conditions, (2) medium engineering geological conditions, (3) poor engineering geological conditions, and (4) worst engineering geological conditions. With the results of this evaluation, the process of selection of relatively suitable alternative route schemes starts.

3.1.2. Key Section Engineering Geological Zoning

EGZ of key sections is carried out based on regional EGZ findings. On the basis of regional EGZ results, a tentatively suitable route is selected for further precise investigations from the list of initial alternative route options. As mentioned earlier, the key sections are the areas with poor or worst engineering geological conditions along the proposed route. As compared to regional EGZ, the key section EGZ has the characteristics of a small evaluation scope and high accuracy. Furthermore, the evaluation parameters of the key sections of EGZ are also quite different from those of regional EGZ evaluation parameters. The nearby ridgelines of the steep mountains on both sides of the project corridor are set here as the evaluation area because slope-related disasters outside these ridgelines shall not affect the proposed 318 National Highway project. Geological disasters have a great negative impact on engineering construction. The EGZ cannot be completely equated with disaster hazard assessment because the EGZ focuses on the engineering geological conditions while the disaster hazard assessment focuses on the spatial probability, time probability and magnitude probability of disaster occurrences. However, there is a certain relation between the two. Some scholars have carried out engineering geological zoning based on the susceptibility level of disasters and other factors [14,15,18]. Areas with poor engineering geological conditions are always prone to geological disasters. So, based on this assumption, using the method of geological disaster hazard assessment, the EGZ of key sections along the project corridor has been carried out. Here, the geological disaster hazard refers to the occurring probability of a disaster of a certain intensity or magnitude within certain time period, which is related to the spatial probability, time probability, and magnitude probability of the disaster [31,32]. The key section EGZ should be conducted based on the above three aspects of a disaster. The EGZ index of key sections can be calculated using the below Equation (2).

$$KZI=P\left(S\right)\times P\left(T\right)\times P\left(M\right)$$

(2)

where KZI is the value of the key section EGZ index, P(S) is the probability index of a disaster occurrence, P(T) is the time probability index, and P(M) is the magnitude probability index.

The geological disaster susceptibility index can replace the spatial probability index of disaster occurrence. Susceptibility refers to the probability of disasters controlled by geological conditions in a certain area [31,33,34]. It is the degree to which an area can be affected by future slope movements, i.e., an estimate of “where” disasters are likely to occur [31,32]. In mathematical language, disaster susceptibility is the probability of spatial (geographical) occurrence of slope failures in a given set of geo-environmental conditions [31,32]. The susceptibility assessment results can be obtained by analyzing disaster data and geological background information of the disaster-prone areas, i.e., topography, lithology, geologic faults, folds and shear zones, hydrology and drainage pattern, vegetation, etc. Slope units can be used as the evaluation unit. This unit division method divides the study area into hydrological regions bounded by drainage and divide lines [35] and has proven to be a reliable unit division method for landslide susceptibility assessment [36]. Slope units can be generated through the Digital Elevation Model (DEM) and the hydrological analysis function of ArcGIS 10.2 (ESRI, Redlands, CA, USA) software [37].

Common susceptibility assessment methods include information quantity, deterministic factor, vector machine, logistic regression (LR) methods [33,34,38], etc. This specific article employs the LR method, which can be used to analyze the relationship between a binary dependent variable and a series of independent variables and can also be used to predict the occurring probability of a disaster. Therefore, it is often used for susceptibility assessment of geological disasters [39,40,41]. Mathematically, the LR model can be expressed as follows:

$$L\left(S\right)=\frac{1}{1+e^{-\left({\alpha }_0+{\alpha }_1{x}_1+{\alpha }_2{x}_2\cdots +{\alpha }_n{x}_n\right)}}$$

(3)

where ${\alpha }_0$ is a constant, n is the number of independent variables, ${\alpha }_i$ (i = 1, 2, …, n) are the slope coefficients of the model and $x_i$ (i = 1, 2, …, n) represent the independent variables. L(S) is the dependent variable, and the value range is 0~1, which can be interpreted as the probability of event occurrence.

The time probability of disasters is estimated by assuming that slope failures are independent random point events in the time domain [32,42]. Some scholars use the Poisson model to obtain the time probability of disasters in a given period [42,43]. Therefore, when there are multiple given periods, multiple time probability disaster maps can be produced. As a result, the workload will be relatively huge when performing hazard calculations. Therefore, this paper proposes to determine the time probability of disaster occurrence based on the number of disasters that have occurred along the slope units in the past. For the disasters obtained through investigation in the study area, it can be assumed that they all occurred within the same X years, assuming that slope unit A has experienced “m” disasters and slope unit B also has experienced “n” disasters, if m > n, then the disaster recurrence period of slope unit A is m/X, which is greater than the disaster recurrence period of slope unit B, n/X. Therefore, it can be assumed that the probability of the occurrence time of disasters in slope unit A is greater than that in slope unit B.

Some scholars use the relationship between disaster volume and cumulative frequency to calculate disaster magnitude probability [43]. Guzzetti et al. (2005) used the probability density function of the landslide area to predict the probability of a specific landslide area within each slope unit [31]. This paper simplifies this method and determines the magnitude probability of disaster occurrence based on the average area of disasters in the slope unit in the past. Assuming that slope unit A has experienced “m” disasters with an average area of “x”, and slope unit B has experienced “n” disasters with an average area of “y”, if x < y, it can be considered that the probability of the occurrence magnitude of the disaster in slope unit A is smaller than that in slope unit B.

After determining the spatial probability index, time probability index, and magnitude probability index of each slope unit, the EGZ of the key section can be carried out, and the KZI of each slope unit in the key section can be obtained. The key section should also be divided into four regions: (1) good engineering geological conditions, (2) medium engineering geological conditions, (3) poor engineering geological conditions, and (4) worst engineering geological conditions.

3.1.3. Remote Sensing Interpretation of Disasters

During performing the EGZ of key sections of linear engineering projects, it is necessary to use the disaster inventories of the study area. Sometimes, due to various factors, it becomes difficult to obtain geological disaster data, such as the slope aspect and elevation of each disaster, through physical field surveys. Therefore, this paper employs the advanced applications of remote sensing for acquiring disasters in this study area all along the project corridor for further interpretations and analysis. Later on, the findings of this study were validated through a field survey.

The remote sensing interpretation process of disasters follows the principle of gradual transition from an area with a high research level and detailed disaster inventories to areas with a low research level and lack of detailed disaster inventories. The interpretation method used in this article mainly refers to the disaster interpretation methods proposed by Nichol and Wong (2005), Singhroy (2009), and Yao (2021) [44,45,46]. First, perform operations such as band synthesis, image enhancement, color calibration, spatial correction, and georeferencing on remote sensing images. The remote sensing interpretation signs of disasters are determined based on the distribution locations of existing disasters and their image characteristics. Subsequently, the location of the suspected disaster was analyzed through color, morphology, etc., to draw the boundary of the disaster.

It is worth mentioning that disasters in the study area include landslides, rock falls, debris flows, etc. As the occurrence of debris flow is closely related to the volume of unconsolidated materials in the channel, the susceptibility assessment method of debris flow is quite different from that of landslides or rock falls. Therefore, debris flow has not been considered in this study.

The main remote sensing interpretation characteristics of disasters are: the shape of landslides usually has elliptical and irregular polygonal boundaries, armchair-shaped back walls and other characteristics; the image color and brightness of landslides or rock falls are greatly different from those of the surrounding environment. Images of landslides or rock falls that have just occurred often show mottled light gray tones or bright white elongated strips (Figure 5a). For old landslides or rock falls, the color and brightness of these disasters in the image are basically consistent with the surrounding environment. However, due to the migration of materials, the landform has been transformed, such as cavities (negative topography) on the hillside, bulging of the front edge of the slope, changes in river channels, etc. (Figure 5b), disasters can be interpreted based on the above landform characteristics.

After desk study and remote sensing-based disaster interpretation, it is mandatory to carry out field investigations for validation and verification of the disaster inventories. After the field survey, corrections will be made to the disaster inventories if required.

3.2. Data

The 30 m DEM data used during the regional EGZ process in this study were downloaded from NASA’s Land Processes Distributed Active Archive Center (https://lpdaac.usgs.gov/ (accessed on 10 February 2022)) while the 12.5 m DEM data used during key section EGZ were downloaded from ALOS PALSAR DEM products (https://search.asf.alaska.edu (accessed on 5 June 2022)). The lithological information has been gathered from the 1:200,000 geological map (https://geocloud.cgs.gov.cn (accessed on 2 August 2021)) while information regarding geologic faults in the project area has been gathered from China Earthquake and Fault Information System (https://data.activetectonics.cn (accessed on 8 October 2022)). The Bouguer gravity anomaly gradient data comes from the International Center for Global Earth Models (ICGEM) (http://icgem.gfz-potsdam.de/home (accessed on 16 November 2022)) [47]. Geothermal heat flow data comes from the China Heat Flow Database (https://chfdb.xyz/ (accessed on 22 November 2022)) [48]. Surface deformation is derived from data collected by Wang and Shen (2020), which can be downloaded from the Harvard Dataverse website (https://doi.org/10.7910/DVN/C1WE3N (accessed on 26 November 2022) for free public access) [49]. The seismic data came from China Earthquake Networks Center (https://news.ceic.ac.cn/ (accessed on 12 December 2022)) and was processed through ArcGIS 10.2 software. The rainfall data has been downloaded from the National Meteorological Science Data Center (http://data.cma.cn/ (accessed on 18 August 2021)).

The remote sensing data used in this article include Google Earth images, Gaofen-2 (GF-2) satellite data and Landsat 8 satellite data from Geospatial Data Cloud (https://www.gscloud.cn (accessed on 9 April 2021)). Google Earth is an interactive geographic data browser that allows high-precision, long-distance navigation in a three-dimensional virtual environment to observe the shape of disasters from different angles. GF-2 satellite images have extremely high resolution (0.8 m) and can identify the boundaries of disasters in detail. Google Earth and GF-2 data are jointly used during disaster interpretation in the study area. In some areas, the Google Earth and GF-2 imageries were not found helpful in identifying disasters; therefore, Landsat 8 imagery and other graphics generated by the DEM were used for the disaster interpretation along the highway corridor. Figure 6 shows the interpretation results of a disaster located in the study area (latitude 30.001°N, longitude 95.039°E) using different remote sensing images.

4. Results

4.1. Regional Engineering Geological Zoning

4.1.1. Evaluation Factors

Based on previous research [50,51,52,53,54] and collected data, the evaluation factors selected for regional EGZ in this paper include fracture density, relief amplitude, Bouguer gravity anomaly gradient, geothermal heat flow, surface deformation, and seismic impact. The faults referred to in this paper are all active, have inherited activity and may be active in the future [55]. The evaluation factors are defined as follows:

• Fracture density (m/km²) (C1) refers to the length of fractures per unit area. The activities of faults disrupt the crustal rock layers, resulting in differential movement of adjacent blocks, thus affecting the stability of the crust [55,56]. The larger the value of C1, the more serious and adverse effects it will have on the region’s stability.
• Relief (C2) (m) is the difference in altitude between the highest and lowest points in a certain area. It can represent the depth of tectonic cutting and the degree of surface erosion and can characterize the intensity of tectonic activities in a region [57].
• Bouguer gravity anomaly gradient (Eotvos) (C3) is the derivative of Bouguer gravity anomaly. Differences in the density of underground rocks and changes in the nature and morphology of geological formations cause Bouguer gravity anomalies. Regions with high Bouguer gravity anomaly gradient values are mostly located in deep fault zones with poor crustal stability [58].
• Geothermal heat flow (mW/m²) (C4) is a phenomenon in which heat energy is transmitted from the earth’s interior to the surface. Its distribution has good correspondence with the distribution and magnitude of earthquakes. The stability of the crust is better in areas with low heat flow values and poorer in areas with high heat flow values [58].
• Surface deformation (mm/yr) (C5) reflects the movement of the crust. The larger the value of the surface deformation, the more active the crustal movement and the worst the stability of the crust.
• Degree of seismic impact (M) (C6) indicates the degree of earthquake influence in the area. The closer the distance to the earthquake epicenter, the more obvious the crustal movement and deformation, and the more unstable the region. The degree of the seismic impact on the project corridor has been obtained by overlying the historical earthquake epicenters data. Each evaluation factor is shown in Figure 7.

4.1.2. Construction of Evaluation Model

Through consultation with other experts in the field of engineering geology, the judgment matrix required for the AHP method was constructed, as shown in

Table 1. The importance scale of the two evaluation factors.

	C1	C2	C3	C4	C5	C6	Weights
C1	1	2	3	2	3/2	5/4	0.260
C2	1/2	1	3/2	1	3/4	1/2	0.125
C3	1/3	2/3	1	3/4	1/2	1/3	0.085
C4	1/2	1	4/3	1	2/3	1/2	0.120
C5	2/3	4/3	2	3/2	1	5/6	0.177
C6	4/5	2	3	2	6/5	1	0.233

, which was solved to obtain the maximum eigenvalue of the matrix as 6.01, the eigenvectors as (1.562, 0.752, 0.511, 0.723, 1.062, 1.399), the weights of the factors as (0.260, 0.125, 0.085, 0.120, 0.177, 0.233), CI is 0.0018, CR is 0.0015, CR is less than 0.1, and the matrix passes the consistency test, so the weights of evaluation factors obtained by AHP method are valid.

Using the natural break method, each factor was divided into five levels, corresponding to values from 1 to 5, as shown in

Table 2. Classification of the different levels of evaluation factors.

Level	C1	C2	C3	C4	C5	C6	Value
First	[0, 3.59)	[0, 166.07)	[0.23, 12.01)	[46.92, 78.42)	[6.26, 12.97)	[5.2, 5.50)	1
Second	[3.59, 9.91)	[166.07, 332.15)	[12.01, 27.79)	[78.42, 115.66)	[12.97, 19.07)	[5.50, 5.63)	2
Third	[9.91, 17.25)	[332.15, 498.22)	[27.79, 50.69)	[115.66, 147.16)	[19.07, 23.46)	[5.63, 5.80)	3
Fourth	[17.25, 26.48)	[498.23, 691.98)	[50.69, 82.53)	[147.16, 170.07)	[23.46, 27.61)	[5.80, 6.10)	4
Fifth	[26.48, 43.74]	[691.98, 3542.92]	[82.53, 164.71]	[170.07, 253.13]	[27.61, 37.48]	[6.10, 6.65]	5

. The higher the level, the worse the geological engineering conditions.

4.1.3. Evaluation Results

Using Equation (1), the regional EGZ index (RZI) of each unit in the study area was obtained. The engineering geological conditions in the study area have been divided by employing the natural break method into four levels based on the RZI, that four levels are good (G), medium (M), poor (P), and worst (W), as shown in Figure 8.

During this study, it has been observed that regions with good engineering geological conditions are mainly distributed in the eastern part of the evaluation area, regions with medium engineering geological conditions are mainly distributed in the western part of the evaluation area, and regions with the poor and worst engineering geological conditions are mainly distributed in the central part of the study area (Figure 8). After capturing the regional engineering geological conditions of the study area, all the alternative route options will be screened, and ultimately, in the end, a single route option will be selected. The National Highway 318 is the route selected after studying the prevailing regional engineering geological conditions and other factors comprehensively along the proposed project corridor. It can be seen that the areas with poor and worst engineering geological conditions crossed by the National Highway 318 route are mainly concentrated near Nyingchi City, as shown in the blue rectangle in Figure 8.

4.2. Engineering Geological Zoning of Key Sections

Based on the aforementioned regional EGZ results, the proposed route is selected. As mentioned earlier, key sections in this paper refer to all the areas along the proposed highway route with poor or worst engineering geological conditions. The EGZ of the key sections along the proposed highway route was carried out between the ridgelines of the mountains on both sides of the 318 highway, as shown in Figure 9 (the blue rectangular area in Figure 8).

4.2.1. Disaster Inventory

A total of 780 disasters with minimum and maximum area of 0.11 km² and 1.37 km² have been identified during this study along the 318 National Highway highway, both with the help of remote sensing interpretation and field survey (Figure 9). Furthermore, with the help of the hydrological analysis module of ArcGIS 10.2 version, a total of 2511 slope units were identified in the key sections along the proposed highway route. Three hundred ninety slope units out of 2511 were found to have geological disasters. Therefore, these 390 slope units were declared as disaster units in this study (Figure 10a).

4.2.2. Evaluation Factors

In this particular disaster susceptibility evaluation method, seven factors were selected: (a) distance from the fault, (b) elevation, (c) slope, (d) annual average rainfall, (e) distance from the rivers, (f) engineering geological rock group, and (g) slope aspect. Figure 10b–h shows all these factors.

4.2.3. Construction of Evaluation Model

When using the LR method for developing an evaluation model, both positive and negative samples are required, i.e., both disaster units and non-disaster units. Three hundred ninety slope units located more than 500 m away from the disaster units were randomly selected as non-disaster units. Figure 10a shows the distribution of disaster and non-disaster units used here in this paper for establishing the disaster susceptibility evaluation model along the proposed highway corridor. In order to reduce the effect of sampling while establishing the mentioned evaluation model, both the disaster and non-disaster units were divided equally into five parts. Among them, four parts were selected as training data for establishing the evaluation model each time, the remaining 1 part was used as prediction data for testing the accuracy of the evaluation model, and a total of five evaluation models were established.

4.2.4. Evaluation Results

Spatial Probability

The spatial probability of disaster occurrence can be obtained from the results of the disaster susceptibility evaluations. The training and prediction results of the evaluation model established by the LR method are shown in

Table 3. The accuracy rate of each evaluation model established by the LR method.

Model	Category	Unit Type	Evaluation Result		Accuracy	Overall Accuracy
			Non-Disaster Unit	Disaster Unit
Model # 1	Training	Non-Disaster Unit	280	32	89.74%	88.62%
		Disaster Unit	39	273	87.50%
	Prediction	Non-Disaster Unit	74	4	94.87%	89.10%
		Disaster Unit	13	65	83.33%
Model # 2	Training	Non-Disaster Unit	281	31	90.06%	88.78%
		Disaster Unit	39	273	87.50%
	Prediction	Non-Disaster Unit	69	9	88.46%	87.82%
		Disaster Unit	10	68	87.18%
Model # 3	Training	Non-Disaster Unit	285	27	91.35%	89.90%
		Disaster Unit	36	276	88.46%
	Prediction	Non-Disaster Unit	69	9	88.46%	85.26%
		Disaster Unit	14	64	82.05%
Model # 4	Training	Non-Disaster Unit	286	26	91.67%	89.90%
		Disaster Unit	37	275	88.14%
	Prediction	Non-Disaster Unit	62	16	79.49%	84.62%
		Disaster Unit	8	70	89.74%
Model # 5	Training	Non-Disaster Unit	275	37	88.14%	88.62%
		Disaster Unit	34	278	89.10%
	Prediction	Non-Disaster Unit	69	9	88.46%	87.18%
		Disaster Unit	11	67	85.90%

. The average accuracy of the training data of the five models is 89.17%, while the average accuracy of the prediction data is 86.79%. As the average accuracies of both the training and prediction data are more than 85%, it indicates that the effect of these models is good.

Using model # 3 (highest overall accuracy) to evaluate the remaining slope units, the susceptibility index for each slope unit was obtained, and the slope units were classified into four categories based on their susceptibility index. The evaluation results show that there are 633 slope units with extremely high disaster susceptibility, 405 with high disaster susceptibility, 404 with medium disaster susceptibility, and 1069 with low disaster susceptibility in the study area (Figure 11). High and extremely high disaster susceptibility slope units are mainly distributed in the middle of the key sections along the proposed highway route.

Time and Magnitude Probability

The time probability and magnitude probability of a disaster occurrence along slope units are divided into five grades based on the number and average area of disasters that occurred along each slope unit. The number of disasters of each grade and the corresponding time probability, the average disaster area, and the corresponding magnitude probability are shown in

Table 4. Time probability and magnitude probability of disaster occurrence in slope units.

Level	Number	Time Probability	Average Area (km2)	Magnitude Probability
First	0	0.2	0	0.2
Second	1	0.4	0~0.05	0.4
Third	2~3	0.6	0.05~0.15	0.6
Fourth	4~5	0.8	0.15~0.35	0.8
Fifth	≥6	1	≥0.35	1

. Figure 12 shows the slope unit’s time and magnitude probability of disaster occurrence.

Zoning Results

KZI of each slope unit can be calculated by employing Equation (2) and then using the natural break method for division of the slope units into four categories. As a result, there are 2158, 186, 140, and 27 slope units with good, medium, poor, and worst engineering geological conditions in the key sections along the 318 National Highway route (Figure 13).

5. Discussion

5.1. Zoning Ideas

The multi-scale EGZ method presented in this paper first evaluates the entire region containing all the alternative route options to serve the route selection process accordingly. After screening down and selecting a comparatively better route option among all the alternative route options through regional EGZ, areas having poor or worst engineering geological conditions within the selected comparatively better route option were further selected for key section EGZ.

Regional EGZ corresponds to the planning stage of linear engineering projects. During the regional EGZ stage, many alternative route options were considered and studied based on the regional engineering geological conditions of those areas. Due to the large scope, the evaluation accuracy of regional EGZ is always relatively low. The main purpose of regional EGZ is to provide a qualitative evaluation for the selection of comparatively better route options. This is similar to regional stability evaluation analysis [59], so the relevant research on regional stability evaluation is referred to here in the regional EGZ, too.

The EGZ of key sections corresponds to the feasibility study stage of linear engineering projects. At this stage, the route has been basically determined, and the evaluation scope is further narrowed down. The area within the ridgeline on both sides of the proposed route is taken as the study area because even if a disaster occurs in the area outside the ridgeline on both sides, it will generally not affect the engineering work inside the ridgeline. Referring to geological disaster hazard assessment, the EGZ of key sections is carried out with respect to the spatial, time, and magnitude probability of the disaster occurrence. The worse the geological engineering conditions, the greater the threat of potential disasters [10]; therefore, it becomes necessary for such areas to strengthen the monitoring or preventive measures.

5.2. Zoning Units

Grid units are selected for conducting regional EGZ, while slope units are selected for conducting key section EGZs. The selection of zoning units varies depending on the purpose of the zoning at different scales [60].

The assessment unit division is rapid and convenient when using grid units for calculations, making it the most commonly used evaluation unit [60,61,62,63]. The division of slope units is closely related to topography and geomorphology, and it is necessary to use the hydrological analysis module of ArcGIS 10.2 software to divide slope units according to the terrain in the study area. Moreover, it is necessary to manually adjust the unreasonable areas in the divided slope units based on actual terrain, where the workload would be relatively large [33]. The occurrence of disasters is closely related to topography and geomorphology; slope units can reflect the characteristics of the terrain and are, therefore, used more frequently in small-scale, high-precision assessments [64].

5.3. Zoning Indicators

The selection of evaluation factors is closely related to the ease of acquisition and its evaluation purpose [10]. In regional EGZ, due to the large evaluation scope, factors related to crustal stability under internal and external geological processes such as fault activity, topography and geomorphology, geophysical fields, surface deformation, and seismic activity can be screened. High accuracy is required to evaluate disaster hazards in key sections, so factors with significant spatial variations are selected, such as distance to fault, elevation, slope, average annual rainfall, distance to rivers, engineering geological rock group, slope aspects, etc.

5.4. Limitation Analysis

This paper employs the geological disaster hazard assessment method to perform EGZ in key sections. While EGZ and geological disaster hazards exhibit differences in evaluation factors and details, it is justifiable to utilize the method of disaster hazard assessment for EGZ in light of the significant adverse impact that geological disasters can have on linear engineering projects, particularly in disaster-prone regions characterized by unfavorable engineering geological conditions. Furthermore, to distinguish this approach from regional EGZ, this paper adopts the method of disaster hazard assessment as a means to carry out the EGZ of key sections. This constitutes a novel endeavor and needs to be verified in practice.

Only landslides and rock falls are considered here during the geological disaster hazard assessment and interpretation in this paper. Debris flows, and land subsidence will surely have an impact on the geological disaster hazard assessment results of the proposed highway route. In the future, disasters such as debris flows, and land subsidence should also be taken into consideration as much as possible.

5.5. Future Prospects

Underground engineering projects are gradually increasing with the development of engineering construction, such as many sections of the Sichuan-Tibet Railway comprising tunnels [13]. For such linear engineering projects that include underground construction, in addition to conducting surface EGZ, underground profile EGZ should also be carried out. The geological engineering conditions on the underground profile are evaluated by analyzing the changes in geological structures, engineering geological rock groups, seismic activity, hydrogeological conditions, geostress, geothermal, groundwater, and other factors. Through new technologies such as machine learning, surface engineering, underground profile engineering and geological zoning results combine to construct a three-dimensional EGZ model, which will provide great safety assurance for the construction of linear engineering projects.

6. Conclusions

This paper proposes a multi-scale EGZ method suitable for linear engineering projects in mountainous areas and introduces this method with the example of the Chengdu- Shigatse section of National Highway 318. The regional EGZ result with large-scale and less accuracy shows that the Nyingchi section is an area with the worst or poor engineering geological conditions for National Highway 318, and it is taken as a key section for small-scale and high-precision EGZ. By determining the spatial probability, time probability, and magnitude probability of disaster occurrence, it is found that areas with the worst and poor engineering geological conditions in the key sections along the highway route are mainly distributed in the curved sections in the north, posing a greater threat to the highway. The method proposed in this paper can effectively reduce the workload while performing the EGZ of linear projects, and the evaluation results can provide a reference for the construction of linear engineering projects and the adjustment after local damage of the route.