Railway Alignment Optimization in Mountainous Regions Considering Spatial Geological Hazards: A Sustainable Safety Perspective

Sustainable railway construction and operation are threatened by densely occurring geological hazards in complex mountainous regions. Thus, during the alignment optimization process, it is vital to reduce the harmful impacts of geological hazards to a railway. However, current alignment-related studies solely consider such threats in existing geological hazard regions and, outside these regions, slight attention has been devoted to the assessment of potential hazardous impacts along the alignment. To this end, this paper proposes a novel railway alignment optimization model considering both existing and potential geological hazards based on quantitative geological hazard evaluation criteria from a sustainable safety perspective. More specifically, a geohazard zone classification method, within which an energy–slope model is integrated, is first developed. Three geohazard regions, namely the geohazard outbreak region, buffer region and fuzzy region, can then be obtained. Afterward, a spatial geological hazard assessment model is constructed considering the geological danger of three kinds of geohazards (debris flows, landslides and rockfalls) and railway construction vulnerability. This model is incorporated into a previous cost–hazard biobjective alignment optimization model. Finally, the effectiveness of the proposed model is verified by applying it to a real-life case of the Sichuan–Tibet railway. The results show that this method can effectively optimize mountain railway alignments by concurrently reducing geological hazards and costs, which is beneficial to railway safety and sustainable construction and operation.


Introduction
Railways are of crucial importance to national economies. For a railway project, the alignment design is the primary work that has profound impacts on the railway's sustainable construction and operation. However, alignment design is quite complex and laborious because of multiple design constraints and numerous influencing factors. In addition, human designers may often overlook many promising railway alternatives connecting the start and end points in the landscape since the time and resources of designers are limited.
To solve the above problems, numerous researchers have devoted considerable efforts to alignment optimization since the 1960s [1]. In the 20th century, some exact mathematical methods have been proposed, such as the calculus of variation [2] and enumeration [3]. In the 21st century, additional methods have been proposed for alignment optimization,   [23] mainly focused on the hazard evaluation of existing geological hazard regions. In fact, in complicated mountainous regions, the influencing ranges of a specific geological hazard are not constant but dynamically change, i.e., the existing geological hazard regions are very likely to further expand. In addition, there may be potential geological hazards outside existing hazard regions, which threaten the sustainable railway construction and operation. These factors limit the applicability of existing methods.
Given the gaps in previous methods noted above, we propose here a novel spatial geological hazard assessment model. The contributions of this paper include the following:

1.
A geological hazard zone classification method is developed to determine railway reachable regions based on an energy-slope model. More specifically, the railwayreachable region is divided into three categories, i.e., the geohazard outbreak region, buffer region and fuzzy region. Among them, the outbreak region includes three kinds of geohazards: rockfall, landslide and debris flow.

2.
Then, a spatial geological hazard quantitative assessment model considering geological danger in railway reachable regions and railway construction vulnerability is constructed for railway alignment optimization. An information value (IV) method is designed for the geological danger assessment process, and a simple score ranking approach is developed for evaluating construction vulnerability. 3.
Lastly, the above model is incorporated into a previous cost-hazard alignment optimization model   [23]). A real-world case of the Sichuan-Tibet railway is used to test the effectiveness of the present method. The detailed data analysis verifies that our method is promising for solving actual mountain railway design problems and is conducive for improving railway safety and sustainable construction and operation.
The rest of this article is organized as follows. Section 2 overviews the railway alignment optimization model and the search algorithm. Section 3 focuses on explaining the spatial geohazard assessment model. Section 4 presents the case study, and the last section summarizes this study.

Frameworks of Railway Alignment Optimization Model and Algorithm
An alignment refers to the 3D centerline of a railway, which can be determined by a set of points of intersection (PIs) between the two endpoints in the landscape. PIs include horizontal points of intersection (HPIs) (Figure 1a) and vertical points of intersection (VPIs) (Figure 1b) [12,24,25]. Thus, the alignment optimization is intended to find the best PIs for satisfying specific objectives, subject to multiple constraints. The basic model and search algorithm used in this work are briefly introduced below.

Railway Alignment Optimization Model
Based on our previous studies on railway alignment design, the coordinates and horizontal curve radius of HPIs and station points, as well as design elevations of VPIs are used as design variables [13]. The alignment cost, which consists of construction and operation costs [12], and the geological hazards of debris flows, landslides and rockfalls, are combined as the objective, while the geometric (i.e., code requirements), construction (mainly for bridges and tunnels), location (e.g., forbidden zones) and geological hazard constraints must be satisfied [15].
The framework of this railway alignment optimization model is shown below (Table 1). Detailed formulations of this model can be found in our previous studies [12,13,15,23] and are not duplicated here. In this study, we mainly focus on improving the geological hazard evaluation part of this model, i.e., the geological hazard objective function.
Overall objective function: Geological hazard function: Minimize the total geological hazard value F Hazard (M A )of a railway.
Railway cost function: Minimize the comprehensive costF Cost (M A ) of a railway.
Where X i and Y i denote an HPI's coordinates; R i denotes horizontal curve radius; K j and H j denote the station and design elevation of each VPI; m and n are the numbers of HPIs and VPIs, respectively. S and E denote the coordinates of the start and end points, respectively; (x k , y k , e k ) denotes the coordinates of k th station. C B , C T , C E , C L and C R denote costs spent on bridge, tunnel, earthwork, length-dependent construction and right-of-way expenses, respectively; C A , C M and C G denote annual operating costs related to the deflection angles at HPIs, alignment length and gradients, respectively; ∆ (cost/year) is a capital recovery factor used to convert total construction cost to annual construction cost, which is usually set as 0.065 in China.

Cost-Hazard Bi-Objective Optimization Algorithm
In our previous studies [14,15,23,26,27], two-stage search algorithms were devised to solve the above model. At the first stage, a 3D-DT search algorithm was developed to search for railway corridors in the broad study area (area-corridor stage). At the second stage, we proposed a hybrid PSO-GA to obtain final alignment solutions in a specific railway corridor (corridor-alignment stage). To find tradeoff solutions considering alignment costs and geological hazards simultaneously, a multi-criteria tournament decision (MTD) method is incorporated into the 3D-DT and a crowding distance computation (CDC) is integrated into the PSO-GA. In the above studies, geological hazards regions are assumed to be static. Only when a railway crosses existing geohazards regions does the procedure start to evaluate the geological hazard. The dynamic development and potential hazards to a railway are overlooked. Thus, in this study, a novel spatial geohazard assessment model is designed to quantitatively evaluate the overall geological hazards of a railway during the alignment optimization process.

The Spatial Geohazard Assessment Model for Railway Design
The alignment geological hazard value used in this study contains two components: geological danger and construction vulnerability. The geological danger value estimates the degree of threats in different kinds of geohazard regions based on a geological hazard zone classification method. The spatial vulnerability value indicates the resistance abilities of an alignment to various geological hazards by evaluating the type and elevation difference from the ground of its construction components [23]. The main flow of the evaluation process is shown in Figure 2.

Geological Hazard Zone Classification
Digitization of the design information is the first step in railway alignment optimization. In our work, all kinds of geographical information required for railway alignment optimization are stored discretely in a regular lattice, which consists of a set of grids, to form a comprehensive geographic information model (CGIM, Figure 3) [13].

Determining the Railway Reachable Region
For railway alignment optimization, the two endpoints should usually be pre-determined in the landscape. Then, due to the constraints of the maximum allowable circuity coefficient (γ max ) and maximum gradient (G max ), the reachable region (R-region) of a railway alignment is therefore limited. More specifically, the R-region of a railway can be determined by analyzing the horizontal and vertical attainability. Detailed explanations are provided below.
Below, S (r S , c S , H S ) and E (r E , c E , H E ), respectively, denote the grids that contain the start and end points, where r S and r E are, respectively, the row of the start and end points; c S and c E are, respectively, the column of the start and end points; H S and H E are the start and end points' design elevations, respectively; L S and L E are the linear distances from a random grid n (r,c) to the start and end grid, respectively.

Horizontal Attainability Determination
In order to meet the requirement of γ max , n (r,c) must satisfy the following conditions (since otherwise n (r,c) is unreachable): where ω is the grid width. The horizontal reachable region is shown in Figure 4a.

Vertical Attainability Determination
For each grid within the horizontal reachable region, we first generate a so-called spatial voxel set, which is obtained by extending the grid's ground elevation upward and downward to determine several feasible ω × ω × ε voxels, where ε is the height of each voxel.
The upper boundary (H max ) is determined by considering the constraints of i max , γ max and the maximum allowable height of bridges (H maxb ), as shown in Equation (4). In addition, given that the geological hazards studied here are surface geological hazards, if the railway is buried in a tunnel and the overlying soil layer reaches a certain thickness (H maxt , assumed to be 25 m in this study), the influences of geohazards are deemed to be negligible. Therefore, the lower boundary (H min ) is determined according to Equation (5). The vertical reachable regions are shown in Figure 4b.
H max > H n > H min (6) where H min , H max denote the minimum reachable depth and height at grid n (r,c), respectively; H n is the ground elevation at the grid.

Geological Hazard Zone Classification
After obtaining the railway R-region, we divide it into three kinds of subregions, i.e., geological outbreak regions, buffer regions and fuzzy regions.
In this paper, the regions where existing rockfalls (U R ), landslides (U L ) and debris flows (U D ) are located are called the geohazard outbreak regions (O-regions). Rockfalls, landslides and debris flows are different types of geological hazards, which are obtained by intersecting type of movements and type of materials. According to [23], a rockfall refers to the geological phenomenon of the sudden detachment of rocks and soil masses on a steep slope under gravitational forces. A landslide refers to the natural phenomenon in which a soil mass slides down the slope of the mountain along a certain soft surface under influences of river erosions, rainwater immersions or earthquakes. Debris flow refers to a special sudden flood with heavy sediments and stones that is triggered by torrential rain, snow or glacier melting.
The grids beside the O-region are threatened by the existing geohazards. The rock and soil masses may expand outward and destroy the railway constructions when the geological hazard breaks out. These grids affected by the O-region constitute the buffer region (B-region), whose main hazard source is the moving masses. The remaining grids outside the B-region are basically not affected by the O-region. However, these grids may be impacted by potential geohazards because the information on potential geohazards is generally difficult to obtain, especially for mountainous regions. Therefore, we cannot directly determine such geological hazards, and the corresponding regions are defined as fuzzy regions (F-regions). A zone classification diagram of a study area is shown in Figure 5. At the alignment design stage, O-regions are explicit and can be obtained from geological survey and historical data. Thus, the main aim of the geological hazard zone classification is to determine the B-region, which reflects the transition between the Oregion and F-region. To achieve this, a geological hazard energy-slope model is applied to generate the B-region. Finally, except for the O-regions and B-regions, other areas in the R-region are F-regions.

Energy-Slope Model
Generally, the soil movement process is regarded as a process of energy conversion [28]. During this process, part of the gravitational potential energy is converted into kinetic energy while the rest is dissipated. The energy dissipation of the mass is mainly caused by the internal and external frictions and the damping motion. This process is expressed in Equation (7): where h, β and m are the elevation, sliding direction and weight of the soil mass, respectively; W 0 and W ∆H are the initial energy and gravity potential energy conversion, respectively; W f and A k are the energy spent on overcoming the total friction and damping motion, respectively. The computing method for Equation (7) is provided below.

Geological Hazard Energy-Slope Model
The run-out distance and impact force of mountain geological hazards can be determined by using an energy-slope model, which is a generic model explaining the movement in specific geomorphologic processes [28]. The soil and rock masses of rockfalls, landslides and debris flows slide from the top to the bottom along multiple single-slope sections. In addition, in the alignment optimization process, the computation accuracy is determined by the CGIM grid width [29]. Thus, based on the energy-slope model, each single-slope section is evaluated by analyzing the energy conversion of all the grids traversed by the slope section.
The single-slope motion process can be abstracted as the horizontal motion (Figure 6a,b) and vertical slide (Figure 6c). In the horizontal motion, the mass moves from the nth grid to the (n + 1) grid along the horizontal direction (β n ) and the maximum vertical slope (α n ). The vertical slope (α n ) between the nth and (n + 1) grids can be computed using Equations (8) and (9). After determining the horizontal motion direction, the vertical sliding process of the soil mass can be simulated (Figure 6c). Taking Figure 6a as an example, this process can be specified as in Equations (10)- (13) according to the law of energy conservation, which is expressed in Equation (7). The case shown in Figure 6b can be similarly obtained.
From the n th grid to the (n + 1) grid: Among them, where F f is the friction, L n is the single slope sliding distance, µ i is the internal friction coefficient computed as µ i = tan ϕ [28] and µ e is the external friction coefficient computed as µ e = e −7.13+9.33 NDVI [30], where ϕ is the angle of internal friction and the NDVI (Normalized Difference Vegetation Index) reflects the vegetation growth status and vegetation coverage of a region, which can be obtained from remote-sensing images, as discussed in Vicente-Serrano et al. (2016) [31]. When the soil mass moves from the starting position and then reaches the nth grid, the energy conversion can be computed with Equation (14), as shown in Figure 7. Denoting that: The kinetic energy of the mass can be determined when it reaches the nth grid: The process of energy conversion is slightly different for different types of geological hazards: (1) For rockfalls and landslides, the very slow starting velocity of the soil mass can be neglected. (2) For rockfalls, the air friction of the falling rock is negligible. Therefore, the boundary conditions for different geological hazards are set as follows: The above procedure will be used in processing all the geological hazard regions in the study area to generate their corresponding B-regions. Therefore, the F-regions can also be determined and the geological hazard zone classification of the study area can be completed.

Geological Danger Evaluation
For different geological hazard subregions, i.e., O-regions, B-regions and F-regions, different methods are used to estimate the degree of threats by using geological danger values.

Geological Danger of O-Regions
Referring to our previous study [23], in an O-region, a susceptibility value S can be obtained by analyzing the environmental factors such as slope angle, slope aspect, NDVI value, area and Gravelius index [32]. Thus, in this work, we use S to quantify the geological danger in an O-Region. Then, S is normalized as shown in Equation (21): where D O (r, c) is the geological danger value of the grid (r,c) which is in O-regions; S(r, c) is the susceptibility value of the grid (r,c); S max and S min are the maximum and minimum susceptibility values in O-regions.

Geological Danger of B-Regions
When a railway passes through a B-region, it may be threatened by geohazards due to the potential expansion of the O-region. As explained in Section 3.1.2, we can quantify such danger with an energy-slope. As the energy of the rock and soil masses dissipates, their kinetic energy also decreases, and hence, in a B-region, the destructive power of a specific kind of geohazard is reflected by the kinetic energy [28]. Therefore, the kinetic energy obtained in Section 3.1.2 can be normalized with Equation (22) and used as the danger value in B-regions.
where D B (r, c) is the geological danger value of the grid (r,c), which is in B-regions; E n (r, c) is the kinetic energy value of the grid (r,c); E max and E min are the maximum and minimum kinetic energy values in B-regions. In addition, if a grid is located in the intersection regions of several B-regions, its danger value can be determined by adding the kinetic energy of every intersected B-region.

Geological Danger of F-Regions
The grids in an F-region may be far from O-regions and are basically not affected by existing geohazards. However, there may be potential and unknown geohazards in F-regions. In such conditions, we propose using a susceptibility value to measure the geological danger in the F-region. In this analysis, rockfall, landslide and debris flows have similar hazard-generating factors and, hence, for simplicity, can be treated as generalized landslides.
In this field, many researchers have used different models to assess the susceptibility of geohazards on different scales according to different factors, such as an artificial neural network model [33], a bivariate statistical approach [34], a Bayesian model [35], a machine learning method [36] and an information value model (IVM) [37,38]. For large-scale railway design problems, we adopt the IVM because it can assess the geo-susceptibility of a broad area very quickly. The assumptions of the IVM are adopted according to [39,40]: first, the areas that have experienced geological hazards in the past are likely to experience them in the future; second, the areas with a similar set of geo-environmental conditions to that of the failed areas (i.e., the regions which have experienced geo-hazards in the past) are also likely to fail in the future.
The occurrence of landslides is influenced by various factors (x i ) [41]. In this study, the environmental parameters include altitude, slope angle, slope aspect, NDVI value, distance to rivers, distance to faults and lithology [42,43]. The basic information of these environmental parameters is shown in Table 2. Moreover, the contributions of different factors to landslide occurrences are different and can be measured using the Information Value (IV): The greater the IV, the greater the contribution of this factor to the occurrence of specific geohazards. In this study, we use the area of the landslide as a parameter in computing the information value. The information value of x ij (the jth subclass of the factor x i ) is determined (as shown in Figure 8) as follows: where N ij is the area of a landslide in subclass x ij (km 2 ); N is the total area of a landslide in the R-region (km 2 ); S i is the total area of subclass i (km 2 ); S is the total area of the R-region (km 2 ). The area data can be obtained from a geohazard region map. As shown in Figure 8, according to the value of the x i factor in the grid n(r,c), the IV of the x i factor can be determined with Equation (23). Then, the total IV for each grid (r,c) can be computed by summing all the IVs of each factor layer. The susceptibility value can thus be obtained by normalizing the IV of each grid in the F-region using ArcGIS 10.2. As shown in Figure 9, after the above process, different geological danger values are set for different geohazard regions. These geological danger values D (r,c) are then recorded in the CGIM, as shown in Figure 10.

Spatial Construction Vulnerability
The construction vulnerability reflects the different resistance of different railway constructions with respect to specific geohazards [44]. To estimate the construction vulnerability, a simple but effective score ranking method is proposed. More specifically, the vulnerability is determined by considering the type and the elevation difference from the ground of different alignment constructions [23], including bridges, tunnels, fills and cut sections. In this process, the type of a construction is determined by checking the fill height or cut depth (∆h), which equals the difference between the railway design elevation (e d ) and the ground elevation (e g ), as specified below: Tunnel, where ∆H B denotes the fill-bridge threshold height; ∆H T denotes the cut-tunnel threshold depth; ∆H B-max is the maximum allowable bridge height. This process is illustrated in Figure 11. Generally, bridges are safer than cut-and-fill sections, and hence the preliminary vulnerability weight (V P ) for bridges is set at 2.0, while that for cut-and-fill sections is set at 3.0. In addition, tunnel sections can avoid the adverse effects of ground hazards well, and hence their vulnerability weights are 1.0. Then, each type of construction can be further classified according to its elevation difference from the ground surface (∆h).
As shown in Table 3, bridges can be classified as general bridges, high bridges and super-high bridges based on two threshold heights of ∆H B−1 and ∆H B−2 . The buried depth of shallow tunnels is less than ∆H min . When the buried depth of a tunnel section exceeds ∆H min , it is deemed a deep tunnel section, which is not affected by surface geological hazards. Lastly, the earthwork sections whose fill height or cut depth are within the range of (∆H C , ∆H F ) are deemed general earthwork sections. Otherwise, a cut section whose excavation height exceeds ∆H F is called a deep cut, and similarly, a fill section whose height exceeds ∆H C m is deemed a high embankment. Constructions of the same type but different elevation differences from the ground also have different vulnerabilities. Each construction with a different elevation difference is assigned a correction factor c (1, 2 and 3) according to its resistance abilities to geological hazards. After determining V P and c, the final vulnerability is determined as follows: Afterwards, the ground elevation of each grid cell in the CGIM can be extended upward and downward to generate the longitudinal space, which can be divided into a series of stacked voxels (ω × ω × ε), thereby extending the geological hazard attributes of the ground surface to 3D spaces. The geological hazard value of each grid can be determined with Equation (27): where r and c are the row and column of each spatial voxel, respectively. The spatial geological hazard map is shown in Figure 12.

Railway Geohazard Determination
To determine the geohazard value of a railway, we first convert the alignment into a set of station points, as shown in Figure 13, based on the CGIM grids that the railway traverses [45]. Then, the total hazard value of a railway can be obtained by accumulating the geological hazard values of all station points, as expressed in Equations (28) and (29): where (x k , y k , e k ) T denotes the coordinates of the kth station point on the alignment. r k ,c k denote the row and column of the kth grid. The steps in computing the total geological hazard value are as follows: STEP 1: Convert the railway alignment into a set of station points according to the ω of the CGIM grids that the alignment traverses. STEP 2: Interpolate the ground elevation of the stations based on a digital elevation model. STEP 3: Compute ∆h for all stations. STEP 4: Determine the type and locations of railway constructions according to ∆h of each grid, i.e., Formula (24). STEP 5: Accumulate the geological hazard values of all the station points using Equation (28) to obtain the total hazard value of an alignment.

Case Profile
The Sichuan-Tibet Railway is one of the most complex railway projects in the world [12] and the most important railway project in China. Here, we use as a case study the Lulang-Linzhi section, which is one of the most difficult sections of the Sichuan-Tibet Railway.
High mountains and deep valleys fill this region, and 62 geological hazard regions are designated in this area. The topographic map, the geological hazard map, the lithology map and the NDVI map of the study area are shown in Figures 14-17, respectively. These data are supported by the China Railway First Survey and Design Institute Group Co. Ltd. The linear distance and elevation difference between the start and end points are 39 km and 414 m, respectively. The constraints are shown in Table 4. The unit costs for alignment optimization are shown in Table 5. The construction threshold heights for vulnerability evaluation are shown in Table 6.

Results Analysis
Using the proposed method and the software "Vizrail 2016", which was specifically developed by our research team for railway alignment design and has been widely used by many Chinese railway design companies, such as China Railway First Survey and Design Institute Group Co. Ltd. and China Railway Eryuan Engineering Group Co. Ltd., an alignment solution is generated automatically (denoted as A C ) and compared with the best railway alignment designed manually by very experienced human designers (denoted as A M ) in the China Railway Eryuan Engineering Group Co. Ltd. The horizontal and vertical alignment comparisons are shown in Figures 18 and 19, respectively. The detailed data analysis is shown in Table 7. It can be observed from these results that:    Figure 18b is located downstream from a landslide. Thus, although the tunnel is shortened, it may be threatened by geohazard expansions. In addition, A C also reduces the hazard value in the F-regions by a total of (1255 − 1127)/1255 = 10.2% compared with A M . Therefore, overall, A C is (1431−1264)/1431 = 11.7% safer than A M . Overall, the model-produced solution is less expensive than the manually designed one by 13.9% while reducing the geological hazard value by 11.7%, which confirms that the proposed method is effective in solving this realistic problem.

Conclusions
Frequent geological hazards in complex mountainous regions greatly threaten the safety and sustainable construction and operation of a railway. It is very important to minimize the geological hazard of a whole railway during a railway alignment optimization process. Existing studies simply treat geological hazards as forbidden areas or only focus on evaluating existing geological hazards that are assumed to only occur within static regions but overlook potential geohazards. That is, the range of hazards may expand outward from existing geohazard regions. Meanwhile, the complex geological conditions of mountainous regions may create potential geological hazards. Moreover, a railway combines different types of construction with varied elevation differences from the ground surface and, hence, its vulnerabilities to specific hazards are difficult to determine. Thus, mountain railway design that considers geological hazard impacts poses a challenging problem.
To solve the above problem, this paper proposes a quantitative spatial geohazard assessment model for railway alignment optimization. Then, this model is incorporated into a previous cost-hazard alignment optimization model. The following are the main contributions of this study:

1.
A geological hazard zone classification method is developed based on an energy-slope model. More specifically, the railway-reachable region (R-region) is divided into three categories, i.e., geohazard outbreak region (O-region), buffer region (B-region) and fuzzy region (F-region). Among them, the O-region includes three kinds of typical geological hazards, i.e., rockfalls, landslides and debris flows.

2.
A spatial geological hazard quantitative assessment model considering geological danger in the railway reachable region and railway construction vulnerability is developed for railway alignment optimization. An information value model (IVM) is designed for the geological danger assessment process, and a simple score ranking approach is developed for evaluating construction vulnerability. 3.
The above model is incorporated into a previous cost-hazard alignment optimization model. A real-world case of the Sichuan-Tibet railway is used to test the effectiveness of the present model. The experimental results show that the model-produced alignment reduces the cost by 13.9% and the hazard value by 11.7% compared with the best manually designed alignment, which verifies that the proposed method is promising for optimizing complex and realistic mountain railway alignments and is beneficial to railway safety as well as sustainable construction and operation.
In the future, the impacts of underground geological hazards, such as caverns, on railway designs should be investigated. Moreover, while this paper mainly focuses on developing the relations for a railway alignment optimization model, improved search algorithms may be found through future research.
Author Contributions: Conceptualization, H.P., J.X., T.S. and W.L.; data curation, J.W. and J.H.; formal analysis, J.X., P.S. and T.S.; funding acquisition, H.P.; investigation, J.X.; Methodology, H.P., J.X., P.S. and T.S.; resources, W.L., J.W. and J.H.; software, H.P., J.X., T.S., W.L., J.W. and J.H.; supervision, H.P. and P.S.; validation, J.X.; writing-original draft, H.P., J.X., P.S. and T.S.; writing-review and editing, H.P., J.X., P.S. and T.S. All authors have read and agreed to the published version of the manuscript. Data Availability Statement: Data available on request due to restrictions, e.g., privacy or ethical. The data presented in this study are available on request from the corresponding author. The data are not publicly available because all design data for the case study presented here, such as topography map, geological hazards map and lithology map, are provided by the China Railway First Survey and Design Institute Group Co. Ltd. These data were provided to support this research project, but detailed information cannot be publicly available because the proprietary rights belong to that company.