Mitigating the Vulnerability of a High-Speed Railway–Air Network by Optimizing the Location of Integrated Transportation Hubs

High-speed trains and flights are the backbone of interregional transportation. A coupled high-speed railway–air transportation network (HRATN) allows these two modes of transportation to substitute for each other when the HRATN is threatened by extreme events. During train and flight operations, the area of influence and occurrence of extreme events are difficult to predict. Reasonably selected integrated transportation hubs can mitigate the vulnerability of HRATNs to extreme events. Optimizing the locations of integrated transportation hubs makes full use of the transport resources of both modes and facilitates the transfer of passengers, thereby decreasing the losses caused by extreme events. In this research, a tri-level model is established to reduce the potential worst loss by optimizing the location of the integrated transportation hub. The column-and-constraint generation algorithm based on duality is used to solve the model. For an illustrative example case, the enumeration method and sensitivity analysis show that building more integrated hubs can mitigate HRATN vulnerability but decreases unit investment income. This trivial conclusion verifies the validity of the model constructed in this paper. For a real-world case, comparing the corresponding models that do not consider extreme events reveals that the proposed strategy does not conflict with the decision-making in the normalized situation, thus verifying the usability of the present research.


I. INTRODUCTION
High-speed rail and air transportation have become the backbone of interregional passenger transportation systems in China. In terms of both passenger volume and infrastructure, high-speed rail and air transport have achieved considerable scale. High-speed rail and air transport carried 2.16 billion and 420 million passengers, respectively, in China in 2020. The number of civil-certificated schedule airports in mainland China reached 260 in 2020, and the operation length of high-speed rail reached 38,000 km. Given their important The associate editor coordinating the review of this manuscript and approving it for publication was Yu Wang . role in interregional passenger transportation, failures can cause serious economic and social losses. Considering the high-speed rail and air transportation system as a coupled transportation system can mitigate vulnerability when an unfavorable event occurs. This research aims to reduce the impact of extreme events on the high-speed railway-air transportation network (HRATN) by optimizing the locations of integrated transportation hubs.
In this research, a HRATN is defined as a holistic transportation system that provides travel services to passengers by one or two transportation modes. Passengers choose the transportation mode depending on their individual preferences. In the HRATN, reasonably selecting the location of the integrated transport hub provides convenience for passengers transferring between different modes and facilitates the full use of different modes of transport in the intermodal network. Shanghai Hongqiao Hub in China and the Frankfurt Hub in Germany are examples of integrated transport hubs that couple two transportation networks and enable passengers to easily complete the transfer between the two modes. However, the construction of this type of transportation hub requires the full cooperation of railway and airline companies and large government investment. Therefore, there are few integrated transportation hubs. Mainland China has 260 civil aviation airports and more than 1,000 high-speed rail stations, but only 13 high-speed rail stations interoperate with an airport. Given the existing high-speed railway infrastructure network and aviation network, the reasonable selection of new nodes as coupling hubs of the two types of transportation is a practical problem that planners need to consider.
The selection of an integrated transport hub in a multimode transportation network should consider not only the economics of transportation but also the vulnerability of the network to extreme conditions [1]. For example, during the smog weather in North China in December 2017, a large number of passengers who needed to fly to Beijing had to continue their journey by train after landing in other cities. If the high-speed rail track is interrupted due to geological disasters, passengers can take a flight at the airport to complete their travel. Therefore, the rational selection of integrated transportation hubs and new airlines can reduce the impact of extreme events on the transportation network.
This article defines a local extreme event as a disaster event that affects the operation of transportation infrastructure in a spatial area and for which the probability distribution of occurrence is difficult to accurately estimate. Such a local extreme event can disturb the operation of transportation components and cause serious HRATN losses (the specific influence of local extreme events will be discussed in detail in the section ''Problem description and network representation'').
One feature of a local extreme event is the spatial influence area. Extreme events have a certain spatial impact area, such as the heavy rain in London on December 24, 2019, the smog weather in North China in December 2017, and the eruption of the Eyjafjalla volcano in Iceland on April 14, 2010, and affect the transportation facilities of a city, region, country, or even multiple countries. An extreme event that interferes with a certain station or a certain railway section will also affect adjacent transportation facilities. This type of extreme event can be viewed as a special local extreme event with a small radius. The geological disaster of the tunnel on China's Lanzhou-Xinjiang high-speed rail in December 2018 caused the downward and upward three stations of the tunnel to be closed. In the Brussels railway station bombing in June 2017, the station was attacked and closed, and trains could not enter or leave the station. This outcome is equivalent to the destruction of all tracks connected to the station. In some studies, scholars have analyzed the impact of a large number of single-component failures on transportation networks [2], [3], [4], [5]. However, these studies fail to consider the spatial proximity of components. To compensate for this defect, this article focuses on the impact of local extreme events with a spatial area of influence on the transportation network.
Another feature of local extreme events such as earthquakes is the inability to predict their occurrence probability and spatial impact area. Some local extreme events can be anticipated based on past data or experience, such as flooding that often occurs in summer, but such predictions are only inferences based on past knowledge, and it is difficult to obtain the probability distributions of location and time. Consequently, it is difficult for planners to prepare before extreme events occur. Many papers have analyzed the impact of different unfavorable events on the transportation network under a given probability distribution [6], [7], [8], [9], [10], [11]. However, as Li pointed out [12], such analyses can yield misleading results when the potential complexity leads to incorrect assumptions. Even for analyses that do not rely on an accurate probability distribution, such as scenario-based reliability research [13], [14], [15], [16], it is necessary to quantify the cognition of existing experience as a weight that describes the possibility of occurrence of different scenarios. The results of this method also depend on the understanding of the possibility of different scenarios.
As the Chinese government continues to invest in high-speed rail and civil aviation infrastructure, more cities will be selected for the construction of integrated hubs. One of the goals of our research is to reasonably select the locations of these hubs and make HRATNs more closely coupled. In addition, since the consequences of extreme events are severe and the probability of occurrence is unpredictable, this research takes as its objective minimizing the potential worst losses rather than the expectation of losses when optimizing the locations of integrated hubs. Additionally, it is necessary to consider how to avoid severe losses due to extreme events in the planning stage of the HRATN, that is, in the hub location selection stage. Therefore, another goal of this paper is to mitigate the vulnerability of HRATNs by optimizing the locations of the integrated hubs.
In this paper, a tri-level model is established that takes the spatial area affected by extreme events into consideration. To solve this model, the column-and-constraint generation (C&CG) algorithm is used. Through this model, the locations of integrated transportation hubs are optimized to reduce the impact of extreme events on passenger travel.
The potential contributions of this research are as follows: 1. This research considers the range of influence of extreme events as a more realistic spatial region rather than several nonadjacent facilities.
2. To overcome the problem that the probability of extreme events is difficult to accurately predict, this research aims to minimize the worst loss rather than the expectation of loss to develop a strategy to mitigate the vulnerability of HRATNs.
3. This research mitigates the vulnerability of HRATNs by optimizing the location of integrated hubs. This strategy can mitigate the vulnerability of a HRATN while strengthening the coupling of HSR and air transport networks at the planning level.
The remainder of this paper is organized as follows. The ''Problem description and network representation'' section describes the details of the problem in this study, the ''Model presentation'' section presents our model, and the ''Algorithm'' section introduces the C&CG algorithm used in this study. Verification and sensitivity analysis of our model are performed in the ''Illustrative example case study'' section. In the ''Real-world case study'' section, we provide advice on the selection of integrated transportation hubs in mainland China based on our model. In the ''Results and Discussion'' section, we summarize the results of the calculations in the ''Illustrative example case study'' and ''Realworld case study'' sections. The last section provides the conclusions and future research directions.

II. PROBLEM DESCRIPTION AND NETWORK REPRESENTATION
In this paper, the vulnerability of the network is mitigated by optimizing the locations of the integrated hubs and newly opened airlines. As mentioned above, the local extreme events considered in this paper are disaster events whose probability and location cannot be accurately estimated. For local extreme events whose occurrence can already be accurately judged, corresponding accident prevention can be carried out based on the prediction results. How to deal with events that can be accurately predicted is outside the scope of this article.

A. SPATIAL INFLUENCE RANGE OF LOCAL EXTREME EVENTS
A common method for analyzing the impact of the spatial area is to divide the study area using a regular or irregular grid [18], [19], [20]. However, the results are affected by the location of the regular grid, as moving the grid will change the results. Therefore, the accuracy of this method depends on the uniqueness and completeness of the segmentation results [17]. The maximum coverage set (MCS) proposed by Ouyang can meet the above requirement and appropriately approximate the area affected by local extreme events [21], [22], [23], [24]. The MCS is defined as the set of all components that can be covered by a circle with radius r unless a new component is added [17]. The advantage of this method is that it uses the spatial position information of nodes and edges to obtain the unique division of the study area for a given influence range. In addition, the number of subareas obtained by this method is limited, and no transportation facilities will be missed. This division method is based on the following geometric features of a circle: (1) For a given radius, the circular area centered on the node in the network is unique. Checking all nodes in the network, all the circular areas that satisfy this relationship can be obtained as shown in Fig. 1 (a). (2) For a given radius, the circular areas that can pass through two nodes are certain. Checking all node pairs in the network, all the circular areas that satisfy this relationship can be obtained as shown in Fig. 1 (b). (3) For a given radius, the circular area that can pass through one node and tangent to one edge is certain. By checking all pairs of nodes and edges in the network, all the circular areas that satisfy this relationship can be obtained as shown in Fig. 1 (c). (4) For a given radius, the circular area that can be tangent to two edges is certain. Checking all pairs of edges in the network, all the circular areas that satisfy this relationship can be obtained as shown in Fig. 1 (d).
Refer to [17] for the detailed method of determining the position of the center of the circular area. By further filtering and removing the circular areas covering the same transportation facilities, we can obtain all the MCSs that can cover all transportation facilities in the network within the given influence range, as shown in Fig. 1 (e).
In addition to the circular MCS used in this article, other shapes can be used to divide the study area if needed to accurately predict the influence area of local extreme events. However, as mentioned above, it is very difficult to accurately predict or approximate the scope of local extreme events before they occur. As a preliminary study, this article does not focus on other shapes.

B. THE INFLUENCE OF LOCAL EXTREME EVENTS ON HRATN VULNERABILITY
Local extreme events have diverse causes but all increase the travel times of passengers in the affected area and their wait times at airports, stations or trains. For example, after a flood breaks through a railway track section, passengers must wait in the train or station for line maintenance; when a railway tunnel is threatened by geological disasters, speed will be reduced to ensure safety; airports will wait for heavy rain and windy weather to subside before allowing flights to take off and land normally. When a local extreme event has a very severe impact on the transportation system and flights and trains cannot be resumed in a short time, airlines and railway companies will cancel trains and flights.
To maintain generalizability, the specific causes of extreme events are not discussed in this study. In this research, local extreme events are defined as emergencies for which the probability of occurrence cannot be predicted and that will cause flights or trains to fail to operate normally within a period of time and have a certain spatial influence range. Furthermore, according to the different modes of transportation affected by extreme events, this paper classifies extreme events into three categories: affecting only aviation transport, affecting only high-speed rail transport, or affecting both transport modes.
There is no standard definition of vulnerability in research on the vulnerability of transportation systems, but the general understanding is that the vulnerability of the system reflects the ability of the network to maintain a certain level of performance in the face of perturbation [25], [26], [27]. The flow-based cost is a common metric used to measure the performance of transportation systems. The flow-based cost not only identifies the path changes caused by the disturbance but also reflects the direct impact of the disturbance on the transportation facilities by setting different use costs of the transportation facilities [28], [29], [30], [31]. Flowbased costs more fully reflect the changes in the ability of the transportation network to meet traffic demand when there is a disturbance. On this basis, this research uses the relative change in flow-based general cost when the HRATN will or will not be threatened by local extreme events to measure the vulnerability of the network. We set the total generalized travel cost of passengers as P 0 when local extreme events will not affect the HRATN and P E to consider the impact of local extreme events. Then, the vulnerability of the transportation network when the transportation network is threatened can be expressed as follows: In this study, we consider travelers' different transportation mode preferences. We define passengers into two categories: one category has a higher time value and is more inclined to choose to travel by plane. The other category has a lower time value and prefers to travel by train. This method of studying travelers' mode choice preferences based on time value has been widely used in similar studies [32].
Passengers are rational and choose a travel path that minimizes the sum of generalized costs. The passenger's generalized costs are calculated based on the time cost, ticket price, and discomfort cost. A passenger's choice of transportation mode reflects not only the travel time and ticket price but also factors such as subjective travel comfort, which causes passengers to choose different transportation modes. To reflect this choice tendency, this paper adopts a method similar to that used in a previous study [32] and introduces a discomfort coefficient γ g m to describe the discomfort cost of passenger g on transportation mode m. The unit of this parameter is yuan/min. In addition, this article introduces the parameter α g (unit: yuan/min) and the fare β m (unit: yuan/km), which describe the time values of different passengers and the ticket price, respectively. This article uses the above three parameters to calculate the travel costs of passengers under different travel paths.

C. NETWORK REPRESENTATION
In Fig. 2, we show the network structure of the HRATN defined in this paper. In the HRATN, there are three types of nodes. (1) City nodes. As shown by node a in Fig. 2, city nodes are the origin and destination nodes of passenger trips. The set of city nodes is N city . (2) Transport terminal nodes. As shown by nodes a r and a a in Fig. 2, transport terminal nodes correspond to high-speed railway stations and airports. The airports can be divided into hub and spoke airports. The set of high-speed railway stations is N r , and the set of airports is N a . (3) Logic transport terminal nodes. To facilitate modeling, this paper defines logic transport terminal nodes, as shown by a r and a a in Fig. 2. The logic transport terminal nodes are established in the candidate cities of integrated transport hubs. The set of logic transport terminal nodes is N l .
Based on whether a city is an alternative integrated transport hub city and the type of airport in the city, the cities in this article are divided into three categories. (1) Nonintegrated transport hub candidate cities. Such cities cannot be selected as new integrated transport hub cities. In reality, some cities, such as Beijing and Shanghai, may already have integrated transportation hubs. (2) Spoke cities. The airports in such cities are spoke airports. (3) Hub cities. The airports in such cities are hub airports. In this paper, there are logical transport terminals in the latter two types of cities, which are integrated transportation hub candidate cities, and the set is N P .
The edges in this paper are composed of two categories: transportation edges and logic edges. The costs of these edges differ.
1. Transportation edges. Transportation edges connect high-speed railway stations and airports in different cities, which are shown as (a a , d a ) and (a r , d r ) in Fig. 2, respectively. For type−g passengers, if the transport mode edge m with line length d is selected, the cost of this edge is cost The capacity of transportation VOLUME 10, 2022 edges is limited. The air transportation network in this paper is a hub-spoke network. The hub airports are connected by trunk airlines, and connections are not limited by distance; an example trunk airline is denoted by (d a , e a ) in Fig. 2. The set of trunk airlines is A a trunk . If the distance between a hub airport and a spoke airport is shorter than D branch , the two airports are connected by branch airlines, as (b a , d a ) shows in Fig. 2. The set of branch lines is A a branch ; there are no airlines between spoke airports. Due to the advantages of scale, a trunk airline has a higher transportation capacity and lower ticket price than a branch airline, which means that Cap trunk > Cap branch and β air trunk < β air branch . If an airport is selected as an integrated transport hub, a new branch airline can be opened to connect adjacent spoke airports within D branch , as shown by (a a , b a ) in Fig. 2. The set of all potential new branch airlines is A a P . 2. Logic edges. The logic edges in this paper facilitate the modeling of passenger travel in the integrated transportation network. Therefore, the capacity of logic edges is unlimited.
(1) Access edges. Access edges are used to connect city nodes, transport terminal nodes and logic transport terminal nodes. The access edges are divided into two types. The access edges connected to the transport terminal nodes and logic transport terminal nodes are the hub selection edges, as shown by (a a , a a ) and a r , a r in Fig. 2. The set of this type of access edge is A hub acc . Another type of access edge is the logic edge that connects the city and transport terminal nodes, such as (a a , a) and (a, a r ) in Fig. 2. The set of this type of access edge is A mode acc . This type of access edge describes how the passengers depart from and arrive at a city. The leaving access edge, such as (a, a r/a ), describes passengers choosing a transport mode and preparing to leave a city. The set of leaving access edges is A mode acc,out . The arriving access edge, such as (a r/a , a), describes passengers arriving at a city by a certain transport mode. The set of arriving access edges is A mode acc,in . (2) Logic transfer edges. Edges such as (a a , a r ) in Fig. 2 are used to connect terminal nodes in the same city. Passengers use this type of edge to transfer between different modes of transport. The set of logic transfer edges is A E trans . (3) Logic transfer edges of the potential integrated transport hub. Edges such as (a a , a r ) in Fig. 2 are used by passengers to complete transfers between different modes of transportation if city a is an integrated transportation hub candidate city. The set of logic transfer edges of the potential integrated transport hub is A u,P trans . The cost of logic edges is set to describe the security check time and the transfer time. Taking a plane requires a stricter safety check. Therefore, passengers need to arrive at an airport much earlier before flight take-off. We denote the times of the safety checks for taking a train and a plane as t train check and t plane check , respectively, where t train check < t plane check . When passengers need to transfer to another transport mode in a city, their transfer time will depend on whether the city has an integrated transportation hub. The travel time between a station and airport will be shorter when the station and airport are part of an integrated transportation hub. We denote the travel time between the station and airport in an integrated transportation hub as t hub and denote this travel time in a city without an integrated transportation hub as t nonhub , where t hub < t nonhub .
Notably, once a passenger arrives at the transport terminal of the destination city, the trip is finished. Therefore, the cost of arriving access edge is set as 0. In addition, the advantage of an integrated transport hub is set as a shorter transfer time in this research. Hence, the cost of the hub selection edge is also set as 0. In summary, the costs of different types of logic edges are shown in Table 1. 3. Super edge. In this research, a super edge is defined as a virtual edge that will not be affected by any local extreme events and has unlimited capacity. The super edge directly connects two city nodes, as shown in (a, d) in Fig. 2. The super edge reflects the unmeet traffic demand that may arise after the HRATN is affected by a local extreme event. The cost of the super edge is set to a value much larger than that of the other edges, so passengers will only choose the super edge when other routes cannot be used. This choice means that other routes cannot meet the requirements of the constraints and that passengers cannot complete the trip. In this study, the cost of the super edge between two cities is set to 50 times the cost of the lowest-cost path.

III. MODEL PRESENTATION
This paper establishes a tri-level model to improve the reliability of high-speed railway-aviation networks (TMIRHAN) to reduce the loss caused by local extreme events. The mathematical notations used in this research are listed in the Table 2.
The outer-level model optimizes the locations of the integrated transportation hubs, the middle-level model finds the worst potential local extreme events, and inner-level model allocates passengers. For ease of description, we introduce the variables in the model here. For the defender, we define the variables w u , w p (u a ,j) and w T ,u (i,j) . They denote whether city u will be selected as an integrated transportation hub, whether potential airline (u a , j) that takes off from the airport in city u can be opened and whether the logic transfer edges of the potential integrated transport hub in city u can be used. We denote the variables of the defender as a vector w. (i,j) denote whether the transportation edges and the potential airline will be affected by the type c extreme event s. x c,e (i,j) and x c,p (i,j) denote whether the transportation edges and the potential airline will be affected by the type c extreme event. We denote the variables of the attacker as the vector x. For the passengers, we define the variable y as the passenger flow of type g between (s, t), which is assigned on edge (i, j). We denote the variable y s,t,g i,j as vector y.
As mentioned earlier, local extreme events are divided into three types in this research, namely, events that can affect the operation of high-speed railways, aviation and both modes. The set of types of extreme events is Event, where Event = {Air, Rail, All} . In this study, all logical edges are within the city and are not affected by extreme events. A potential airline is not allowed to open if it is within the scope of an extreme event that affects air transport. Therefore, depending on the type of extreme event c, the composition of the set of edges affected by or not affected by the extreme event varies, as shown in Table 3.

A. OBJECTIVE FUNCTION
For a certain type of extreme event c, the goals of the defender, the attacker, and the user can be quantified by the generalized travel costs of the passengers in the HRATN. This total cost includes the costs of the edges that are affected or not affected by the local extreme event and the loss caused by the unmet traffic demand, as shown in equation (2): In equation (2), cost g (i,j) is the generalized travel cost of edge (i, j) for the type g passenger. pen g (i,j) is the penalty value of edge (i, j) for the type g passenger if this edge is affected by the local extreme event. In this research, this value reflects that passengers need a longer time to complete the trip: In this equation, d (i,j) is the length of edge (i, j), and v m is the speed of transportation mode m. We use XPen to reflect the severity of the local extreme events. As mentioned above, in this study, the impact of local extreme events on a HRATN is reflected as delays in the affected area, that is, a longer travel time. In equation (3), α g + γ g m is the time cost per unit time, and d/v m is the travel time of transportation mode m passing through edge (i, j) under normal conditions. We use the parameter XPen to reflect the impact of local extreme events on the transportation network; that is, when affected by local extreme events, it will take XPen times longer than normal to complete the trip. VOLUME 10, 2022  To facilitate the following description, we use the vectors y and x to represent the variables in equation (2), and the objective function can be expressed as f (x, y).

B. PASSENGER FLOW ASSIGNMENT MODEL WITH CAPACITY CONSTRAINTS IN THE HRATN
The proposed passenger flow assignment model with capacity constraints for a high-speed railway-air network (PFAMHAN) simulates passenger choices regarding transport modes and routing in the HRATN. For a given type of extreme weather c, the PFAMHAN is listed as follows: In the above constraints, optimization goal (4) means that the passenger flow is allocated with the goal of minimizing the total travel generalized cost of passengers. Equation (5) is the passenger flow constraint. Equations (6) and (7) ensure that passengers who depart from city s and arrive at city t will not pass through other city nodes.
Equation (8) is the constraint on the total travel time of passengers, which means that for g−type passengers between any OD pair, the total travel time should be less than a value corresponding to the travel time requirements in the transportation network. In the formula, TimeCon st,g is the total travel time of g−type passengers between OD pair st without the influence of local extreme events. This parameter is obtained in advance by solving the PFAMHAN that fixes w and x to 0 and removes constraints (8) and (11). The parameter Pre reflects the planner's tolerance for the impact of local extreme events; that is, when the total travel time caused by the local extreme events will not exceed Pre times that of the normal condition, the impact is acceptable.
Equation (9) reflects the impact of different local extreme events on different transportation modes. In the formula, F Total is the total travel demand in the HRATN. This formula specifies that if the potential airline (u a , j) will be affected by local extreme events (x c,p (u a ,j) = 1), then even if this potential airline is selected as an openable route (w p (u a ,j) = 1), passenger flow will not be allocated on this edge. For local extreme events that only affect the operation of highspeed railways, no potential airline will be affected. However, for local extreme events that will affect air transportation or both modes of transportation, the potential airlines may be affected. In the next section, we will further discuss how to determine whether a potential airline will be affected by a local extreme event.
Equation (10) shows that if city u is selected as a new integrated transport hub, passengers can transfer through the logic transfer edge of the potential integrated transport hub. Equation (11) is the capacity of edges.
This research defines the super edge to reflect the traffic demand that cannot be met. Without the super edge, local extreme events may cause the travel time of passengers to increase significantly, and constraint (8) cannot be satisfied, which makes the model infeasible. The super edge can avoid this situation. The super edge is set as a virtual edge that has no capacity limitation and has much higher usage. The super edge is only directly connected to city nodes. In the objective function, the optimization aim is to minimize the total generalized cost. Therefore, the passenger flow will be allocated to the super edge only when other edges cannot be selected, which corresponds to the scenario in which local extreme events prevent part of the travel demand from being met.
As mentioned above, the anticipated impact of local extreme events on the transportation network is a longer travel time, which will increase the generalized travel cost. Although this method does not directly reflect the impact of local extreme events on the capacity of the transportation edges, it can still reflect the reduction in transportation capacity. This model allocates passenger flow with the goal of minimizing the total travel generalized cost. If the transportation edges are affected by local extreme events, then the passenger flow will be allocated as little as possible to these edges to avoid the impact of local extreme events. The difference in flow allocated to this edge before and after the impact of the local extreme event reflects the capacity drop caused by VOLUME 10, 2022 pen g (i,j) in equation (3). If the local extreme event has a very large impact on the transportation edges, then XPen will be so large that once edge (i, j) is selected, constraint (8) cannot be satisfied. Consequently, the actual ability of the edge will be completely reduced to 0.
In contrast to directly changing the capacity value of transportation facilities, we use the increasing costs to reflect the external disturbance. This method can describe transportation delays caused by external disturbances. In the study of the shortest path interdiction problem, this method is called partial interdiction of the transportation network. Combined with virtual traffic flow or super edges, this method can also reflect the complete interdiction of transportation edges caused by larger disturbances. Therefore, this method has been widely used in related research [15], [33], [34].
We use vector Y to represent the feasible region of variable y. From equations (8), (9) and (10), we can see that the feasible region Y is influenced by w and x. Therefore, we denote equations (5)- (12), which give the feasible region of variable y, as Y(w), where y ∈ Y(w, x).

C. SPATIALLY LOCALIZED INTERDICTION MODEL OF HRATN WITH RESOURCE CONSTRAINTS
The spatially localized interdiction model of a HRATN (SLIMHAN) with resource constraints indicates that for a limited number of potential spatial areas affected by local extreme events, a virtual attacker develops an attack strategy that can maximize the loss of the HRATN. For a given type of extreme event c, the SLIMHAN is given by the following: x s,c,e x c,p Equation (14) shows that the number of extreme events that can occur at the same time is B A . In equation (15), d s,c,E (i,j) is the relation parameter between the transportation edges and the spatial area affected by the extreme event. For a type c extreme event s with an impact range of d, if edge (i, j) can be affected by this extreme weather and is within the corresponding range, d s,c,E (i,j) is 1; otherwise, it is 0. Equation (15) builds the relationship between the extreme event s and edge (i, j). Equation (16) indicates that for type c extreme events, even if only one of the extreme events affects edge (i, j), edge (i, j) will be affected; otherwise, it will not be affected. Equations (17) and (18) show a similar relation between the potential airlines and the local extreme event. The parameter d s,c,P (i,j) is the relation parameter between the potential airline and the spatial area affected by the extreme event. The values of d s,c,E (i,j) and d s,c,P (i,j) are precalculated according to the location information of the center of the MCS and the edges. For a type c extreme event s with an impact range of d, if edge (i, j) can be affected by this extreme event and is within the corresponding range, d s,c,P (i,j) is 1. For a local extreme event that only affects the operation of trains, the value of d s,c,P (i,j) for all potential airlines is 0. For a local extreme event that can affect air travel or both modes, the value of d s,c,P (i,j) is calculated according to the distance between the edges and the center of the MCS. This parameter indicates that different types of extreme events can affect different types of transport modes. For the convenience of description, we denote the feasible region described by constraints (14)- (19) as X, where x ∈ X.

D. HRATN DESIGN MODEL WITH RESOURCE CONSTRAINTS
In the HRATN design model with resource constraints (HANDM), the planner optimizes the locations of the integrated transportation hubs and selects the newly opened alternative branch airlines. For a given type of extreme event c, the HANDM equation set is as follows: In equation (21), w u is a binary variable. If city u is selected as an integrated transport hub, w u is 1; otherwise, it is 0. Equation (22) indicates that B D cities can be selected as integrated transport hubs among the integrated transportation hub candidate cities. Equations (22) and (23) indicate that if city u is selected as an integrated transport hub, a potential new branch airline that departs from city u can be opened, and the passengers can transfer through the logic transfer edges of the potential integrated transport hub in city u.
We denote the feasible region of w, which is given in equations (21)- (24), as W, where w ∈ W.

E. THE OVERALL FORM OF TMIRHAN
In summary, for a given type of local extreme event c, the overall mathematical form of TMIRHAN can be expressed From constraints (5)-(10), we can see that the feasible region of the inner variable y varies with the outer and middle variables w and x. Although w and x do not directly affect the feasible region of each other, the value of the variable y affected by the two can evaluate the solution of w and x through the objective function. For different combinations of w and x, different passenger flow assignments y can be obtained by solving PFAMHAN. After obtaining the best defense strategy w * and the worst attack strategy x * , the passenger flow assignment y * that can minimize the loss of the HRATN can be obtained by solving PFAMHAN when the worst local extreme event occurs.

IV. ALGORITHM
The model proposed in this research has a tri-level nonlinear structure and cannot be solved by commercial solvers. For this type of model, the decomposition algorithm proposed by Alderson et al. [7], [35] and Zeng and Zhao [34] is a common solution method, and Zeng et al. named this algorithm the column-and-constraint generation algorithm (C&CG algorithm). This algorithm has been widely used to solve such tri-level models [15], [33], [34], [36], [37]. The convergence of this algorithm has been proven in previous studies [33], [34], [36]. In this section, we introduce this algorithm.
As described above, TMIRHAN can be expressed as follows: For a givenx k , the max operator of the middle level is temporarily ignored, and equation (25) can be expressed as follows: Optimizing equation (26) gives w and y k , which can minimize the total general cost of passengers f (x k , y k ). This is an ideal condition for the givenx k and does not consider other possible x. Other unconsidered local extreme events may cause greater losses to the system. Although this situation is not practical, it provides a lower bound (LB) for the total generalized cost of the system. Furthermore, we suppose we obtain h local extreme events and the corresponding setX H ,X H = x 1 ,x 2 , . . . ,x h . Then, TMIRHAN can be expressed as follows: The min operators of the innermost and outermost models have the same objective function and optimization direction. For everyx h inX H , we can solve min independently in advance. Therefore, the innermost and outermost models can be combined as follows: WhenX H includes all possible x, the w obtained from equation (28) is the best w * that can minimize the loss caused by x. If we enumerate all combinations of w and x, w * can be obtained just by comparing the objective function value f (x, y). When the number of combinations of w and x is small, the enumeration is still acceptable, but in reality, the explosion of combinations makes this method impractical. To solve this problem, we can slack equation (28) as a relaxation main problem: Solving this relaxation main problem, the defense strategŷ w that considersX H is obtained. As mentioned above, it is an ideal situation to formulate a defense strategyŵ for the known partialX H . Other x that are not considered may cause greater losses. Therefore, V M P is an LB of the generalized cost.
Similarly, for a givenŵ, TMIRHAN can be decomposed into a subproblem expressed by the following formula to obtain x: The above formula expresses that with a known defense strategyŵ, the attacker expects to find an x that can maximize the generalized cost. The above formula only considers the givenŵ, which is also an idealized condition reflecting the maximum loss that may be caused by local extreme events when the protection strategy has not yet reached the optimum. Other protection strategies may reduce the generalized cost caused by this x, but for the givenŵ, V S P is still an upper bound UB).
Solving equation (32) yields a knownx. Adding the newly generatedx intoX H and adding the corresponding Y(w,x) to the constraints (30)- (31) in the main relaxation problem as a new set of cuts allows the main relaxation problem to VOLUME 10, 2022 be solved iteratively. As the generated cut increases, the UB and LB gradually approach each other. When B U = B L , the obtained w * is the optimal defense strategy for a given set of local extreme events X. In summary, this algorithm solves the original problem min w∈W max x∈X min y∈Y(w,x) f (x, y) by decomposing the original problem into a relaxation main problem and subproblems. In the following, we will show the details of the subproblems and the relaxation main problem in this research.

A. FORMULATION OF THE SUBPROBLEM
As mentioned before, for a givenŵ, the subproblem can be expressed as max x∈X min y∈Yŵ,x) f (x, y), which is combined with the middle-level model SLIMHAN and the inner-level model PFAMHAN. This subproblem is a bilevel model and cannot be solved by a commercial solver.
To solve this problem, we apply the strong duality theorem to combine the middle and inner models of the subproblem into a single-level model. This method is a common way to solve the bilevel subproblems in such tri-level models [13], [15], [28], [33], [34]. For a givenŵ, the super edges with no capacity limitation can ensure the feasibility of the inner model regardless of how the local extreme event affects the HRATN. Therefore, the strong duality theorem can be applied to obtain the dual model max π f (ŵ,x, π ), which has the same optimal value. π is the dual variable of the constraints of PFAMHAN. We can see that the inner dual model and the middle model have the same optimization direction. Thus, the inner dual model and the middle model can be combined into a single-level model max which is the dual subproblem of the original subproblem. Specifically, for a given type of local extreme event type c andŵ, we set the dual variables in equations (5)-(11) as α (s,t),g i , β (s,t),out i , β (s,t),in i , θ s g t, γ ( i, j), δ ( i, j) u and π ( i, j). Then, we obtain the dual subproblem: There are four nonlinear terms in the dual subproblem. The first two nonlinear terms are the max operator in equations (16) and (18). These two equations can be linearized as the following linear equations: Another two nonlinear terms are the product of a binary variable and a continuous variable in the objective function and equation (33), which are x c,p (i,j) ·γ (i,j) and x c,e (i,j) ·θ (s,t) g . We can use the big-M method to linearize those two nonlinear terms.
We introduce two new continuous variables, h (i,j) = x c,p (i,j) · γ (i,j) and e s,t,g , which are also added to the dual subproblem: The values of M 1 and M 2 can be set as follows: for each edge in A a P , we set the value of M 1 as the maximal total generalized cost if all passengers select this edge, which means M 1 (i,j) = max{cost a (i,j) , cost r (i,j) } · F Total . The value of M 2 is also set as the maximal total generalized penalty: M 2 (i,j) = max{pen a (i,j) , pen r (i,j) } · F Total . Therefore, the overall linearized dual subproblem is as follows: By linearizing this objective function, we can obtain the main problem (58)-(67), as shown at the bottom of the next page.

C. THE OVERALL ALGORITHM STEPS
By combining the above subproblems and the main problem, the flowchart of the C&CG algorithm used in this article can be obtained as follows: Step 1: Initialize parameters such as the UB and LB, establish the initial defense strategy, and generate a set of local extreme events of a given type and radius of influence.
Step 2: Solve the dual subproblems, and add the obtained attack strategyx to obtain the attack strategy set X h .
Step 3: If the UB value is greater than the objective function value V SP of the dual subproblem, the UB is updated with V SP , and the currentx is recorded as the worst local extreme event x * .
Step 4: If the gap between the UB and LB is sufficiently small, that is, UB − LB ≤ ε · LB, or the number of iterations exceeds the set value T max , then the currentŵ is the optimal defense strategy w * , and the algorithm terminates; otherwise, the algorithm proceeds to Step 5.
Step 5: Solve the main problem, and obtain the current defense strategyŵ and the objective function value V M P.
Step 6: If LB < V MP , update LB with V MP .
Step 7: Determine whether the algorithm is terminated according to Step 4. If the termination condition is not met, enterŵ into Step 2, and continue the calculation.
Step 8: Fix w and x as w * and x * , and solve PFAMHAN to obtain y * .

V. ILLUSTRATIVE EXAMPLE CASE STUDY A. CASE PRESENTATION
The geographical layout of the case explored in this study is shown in Fig. 3. The coordinates of the cities in Fig. 3 are listed in Table 4.
In this case, each city has high-speed railway stations and airports. Among them, cities 4, 14, and 21 are integrated transportation hubs, and the airports in those cities are hub airports. The other cities are integrated transportation hub candidate cities with the logic transfer edges of potential integrated transport hubs and transport hub selection access edges. The airports in cities 1, 7, 11, 13, and 23 are hub airports. The airports in the other cities are spoke airports.
There are 76 directed railway edges, 100 existing air edges and 44 potential airlines in the network. The different transport edges are straight lines, and the distance is the straight-line distance between the corresponding cities. Table 5 summarizes the attribute calculation standards and capacity values of different edges in the network as follows. A summary of the attribute parameter values used in this article is given in Table 6.
The experiments were performed on a laptop with an Intel Core i5-7200U 2.5-GHz CPU and 16 GB memory. The algorithm was programmed with Python 3.7, and the solver was Gurobi 9.0.
In Table 7, we show the calculation times of local extreme events with a radius of 200 km that affect both air and high-speed rail transportation under different combinations of B A and B D . Table 7 shows that as B A and B D increase, the feasible region of the solution gradually increases, and the required calculation time also increases. Because this study considers the optimization of the location of integrated transport hubs at the planning level, the efficiency of the algorithm used in this paper can meet the relevant requirements.

B. COMPARISON OF THE TOTAL GENERALIZED COST UNDER DIFFERENT w u
As mentioned in the algorithm introduction, since the numbers of w and x are limited, enumerating all possible w and x can also find the optimal w * in theory. In this section, we consider the case of B D = 1, B A = 1 and a radius of influence of 75 km and use enumeration to verify the  = 0 ∀ (s, t) ∈ OD, i ∈ N city |i ={s, t}, h = 1, 2, . . . , t    correctness of the model and algorithm under different types of local extreme events. When B D = 1 and B A = 1, the number of combinations of the integrated transport hub candidates and the potential local extremes is |N P |·|S c |. In this case, this number of combinations is still acceptable in terms of computational cost.
We separately consider the local extreme events that can affect the operation of flights, trains and both modes. Fixing w u and O s as the city in N P and the local extreme event in S c and solving PFAMHAN, we can obtain the total generalized travel cost for different combinations of w u and O s . For the three types of local extreme events, we denote the total generalized travel cost under the different combinations of w u and O s as cost u O c s . In addition, we apply the models and algorithms used in this research to obtain the optimal w u : the integrated transportation hub should be built in cities 12, 17 and 12 for the three types of local extreme events. Next, for the different types of local extreme events, we fix w u as 12, 17, and 12, fix O c s as the local extreme event in S c , and solve PFAMHAN to obtain the total generalized travel cost cost * O c s . In Fig. 4, we show the relative difference in the generalized cost between when other cities are selected as the integrated transportation hub and when the optimal w u is selected for different local extreme events O c s . This value is calculated as , and we map this value as a heatmap in Fig. 4. As shown in Fig. 4, regardless of what types or which local extreme events occur, the defense strategy obtained by our model will always make the generalized travel cost of the system no larger than in other cases. This result confirms that the model and algorithm in this research can achieve the expectation of minimizing the maximum loss.

C. THE WORST LOCAL EXTREME EVENT AND CORRESPONDING DEFENSE STRATEGY
Solving TMIRHAN when B D = 0 and B A = 1 can yield the local extreme event that has the greatest impact on the HRATN without establishing a new integrated transport hub. In Fig. 5, we identified three types of local extreme events with different ranges (radii of 200 km, 175 km, and 150 km) that can influence the operation of high-speed railways, airplanes and both modes. As shown in Fig. 5, the worst local extreme events can be distributed in different areas according to the influence radius and type of event.
Considering the influence of new integrated transportation hubs, the corresponding defense strategy can be obtained by solving TMIRHAN. In Fig. 6, we identified defense strategies for three types of extreme events when B D = 3 and B A = 1. The influence radius of the local extreme events is set as 200 km. Fig. 6 shows that the defense strategies for different types of extreme events are also different. For example, although cities 4, 14, and 21 were selected as the new integrated transportation hubs to protect against extreme events, the new branch airlines that are opened are different. VOLUME 10, 2022

D. THE EFFECTS OF VARIATIONS IN THE INFLUENCE RADIUS OF LOCAL EXTREME EVENTS, B D AND B A
Next, we analyze the effects of B D and B A on the vulnerability of the HRATN system for local extreme events of different types and radii by metric Normalized vulnerability. It should be noted that for local extreme events that only affect the operation of high-speed railways, the MCSs for radii of less than 150 km coincide with the MCSs with larger radii, so the calculation results are the same and will not be shown.
In Fig. 7, we analyze the influence of different B D values on system vulnerability when the radius of influence changes. The curve of B D = 0 in the figure shows that the influence of local extreme events on the system gradually decreases as the radius of influence decreases. By comparing the curves with B D = 0, 3, 6, and 9 in the figure, we find that as B D increases, the loss of the system decreases regardless of the radius of the local extreme event. In addition, we find that the value of the curve of B D = 0 is much higher than the      8 shows that at the same radius, regardless of whether a defense strategy is adopted, the loss of the HRATN system will increase with increasing B A ; that is, the greater the number of local extreme events that occur at the same time, the greater the impact on the network system. Regardless of the value of B A , the curve of B D = 0 is always higher than the curve corresponding to B D = 3; thus, the defense strategy can reduce the impact of local extreme events on the network system and mitigate the vulnerability of the HRATN system. In addition, a comparison of Fig. 8(a), (b) and (c) indicates that the value of the curve in (c) is higher than the values of the corresponding curves in (a) and (b); therefore, compared with the local extreme event that affects a single mode of transportation, the local extreme event that affects both modes of transportation will cause greater losses to the network system.

VI. REAL-WORLD CASE STUDIES A. CASE DESCRIPTION
In this section, we consider a larger provincial level case that includes 35 cities. Among these cities, Beijing, Shanghai, Shenzhen, Guangzhou, Wuhan, Tianjin, Zhengzhou, Changchun, Shijiazhuang, Guiyang and Lanzhou have built integrated transportation hubs. The diagram of this case is shown in Fig. 9.
In this case, if a city has multiple airports or stations, they are combined into one airport or station. The latitude and longitude of the city are used as the location information of the airports and railway stations. For cities without integrated transportation hubs, the travel time between stations and airports can be obtained by searching the Baidu Map, and this value is set to 0.5 hours in the integrated hub. The relevant capacity of high-speed railway and air transportation in the case is calculated according to the number of trains and flights operating in the area for a week. Both the train and flight schedules are from January 6 to January 12, 2020. The timetable data for the trains with the prefix G, D, or C are obtained from the ticket book website of the China Railway Corporation. The flight data used here were obtained from Ctrip, a popular tourism website in China. We assume that each train has a capacity of 600 people (the capacity of CRH2 8-carriage trains) and each flight has a capacity of 300 people (the capacity of Airbus A-330). The capacity of potential airlines is set to seven flights per week. The total number of travel requests is the total number of passengers carried by trains and planes in Guangdong, Fujian, Hubei, Hunan, Anhui, Jiangsu and Jiangxi Provinces, which were covered by this case during the first week of the Chinese Spring Festival in 2019. These data were obtained by consulting Chinese provincial government websites. The travel demand between cities in the region is calculated according to the GDP ratio of each city. The values of other parameters in this case are the same as those in the illustrative example case study. It should be noted that the detailed geographic information of the high-speed railway track and the flight path data are difficult to obtain due to the safety requirements of transportation companies. Therefore, in this case, a straight line is used instead. In the actual application of government planning departments, real geographic data, security checks, transfer times and other parameters can be used to obtain more accurate results.

B. THE LOCATION OPTIMIZATION OF THE INTEGRATED TRANSPORTATION HUB CONSIDERING EXTREME EVENTS
In Fig. 10, we show the selected integrated transportation hub when B A = 1 and B D = 1 and the influence radius of the local extreme events is set as 75 km. The calculation times for these three cases are 1383.58 s, 3733.12 s and 68754.46 s, respectively.
As shown in Fig. 10, Jinan is selected as the integrated transportation hub when local extreme events can affect the operation of air travel and both modes, whereas when local extreme events can affect train operations, Huangshan is selected as the integrated transportation hub. Jinan is the capital city of Shandong Province in China. Currently, Jinan's integrated transportation hub construction plan is in progress, and the connection between Jinan Airport and the high-speed railway station is expected to be completed in 2025. Huangshan is an important city in Anhui Province. Huangshanbei Station is the hub station of the Hefei-Fuzhou High-speed Railway and Hangzhou-Huangshan High-speed Railway. Huangshan Tunxi National Airport is the second largest civil airport in Anhui Province.
Because of the huge investment and the indispensable coordination between the rail and aviation enterprise, the selection of integrated transportation hubs is government led in China. When a decision-maker needs to make a decision in those cities, a more comprehensive consideration should be performed. In this section, we attempt to provide advice on the selection of integrated transportation hubs according to the model proposed in this research. We further consider the situation when 2-3 integrated hubs can be selected, and the calculation results are listed in Table 8.
We can see that the selected integrated transportation hubs may differ according to the types of local extreme events. As shown in Table 8, Jinan, Shenyang, Nanjing, and Huangshan are selected many times. Among these cities, Nanjing and Shenyang are among China's twelve trunk airports and have not yet undergone integrated transformation. Shenyang is a high-speed rail hub connecting the railway networks in northeastern and northern China. Shenyang is FIGURE 10. The newly selected integrated transport hub and potential airlines in the provincial level case. The spatial local extreme events in (a), (b) and (c) are extreme events that can affect aviation, high-speed railways and both modes, respectively. the provincial capital of Liaoning and is a central city in northeastern China. Nanjing is the provincial capital of Jiangsu, which is one of the most prosperous provinces in China. The Nanjingnan railway station is the largest high-speed railway station in Asia and is a critical station in several high-speed railway lines. Therefore, Nanjing and Shenyang have solid foundations for high-speed rail and aviation transportation. If the high-speed rail stations and airports in these two cities can be integrated, more passengers can be served when local extreme events occur, and the vulnerability of the HRATN can be mitigated. According to the calculation results in this paper, these cities should be prioritized when planning new national-level integrated hubs.

C. THE LOCATION OPTIMIZATION OF THE INTEGRATED TRANSPORTATION HUB WITHOUT CONSIDERING EXTREME EVENTS
Setting the value of B A in constraint (14) to 0 and fixing the value of x to 0, TMIRHAN will degenerate into an integrated hub location model for minimizing the total generalized cost of passengers when extreme events are not considered (IHLM). IHLM is given by the following: The min-min operator in equation (58) has the same objective function and optimization direction and can be combined as a single-level model: This model is a classic location-routing problem (LRP). Although the various factors considered in LRP models describe different scenarios [38], [39], [40], some common basic constraints reflect the basic scenarios of the LRP model, such as the flow balance constraints and routing-location coupling constraints in equations (9) and (10).
In this section, we compare the calculation results of TMIRHAN and IHLM to verify the practicality of the vulnerability mitigation strategy proposed in this research. In the validation, we also use the example in Fig. 9 and consider extreme events that only affect air transport. The comparison results are listed in Table 9.
The comparison of the calculation results of the two models shows that when two cities need to be selected for integrated hub construction, both models choose Huangshan and Shenyang. When three cities are chosen for integrated hub construction, TMIRHAN and IHLM choose to build integrated hubs in Huangshan and Nanjing, respectively, and both select Shenyang and Jinan as the remaining integration hubs. According to the above comparison, considering the worst possible losses caused by extreme events will have a certain impact on the selection of an integrated hub. However, considering the consequences of extreme events will not unduly interfere with the decision on the location of the integration hub compared to normal conditions. Thus, the method proposed in this paper for mitigating the vulnerability of the HRATN by optimizing integrated hub location selection does not conflict with the hub location selection under normal situations and can be applied to actual hub location selection optimization.

VII. RESULTS AND DISCUSSION
To verify the strategy proposed in this research to mitigate the vulnerability of the HRATN, an illustrative example case and a real-world case were analyzed. In the illustrative example case, the method of enumeration and comparison were first used to verify that the location plan obtained by the model constructed in this paper can realize the goal of minimizing the potential worst loss. On this basis, a sensitivity analysis of the influence range of extreme events, the number of new integrated hubs that can be selected, and the number of extreme events that may occur simultaneously was performed. The numerical experiments showed that when one value of B D or B A is fixed, increasing the other value will lead to a significant increase or decrease in the loss of the HRATN (for example, if B A is fixed, increasing B D will reduce the losses of the HRATN, and if B D is fixed, increasing B A will cause a larger loss of the HRATN). However, as this value increases further, the difference in loss decreases. This law of change is in line with the law of diminishing returns on investment in economics. If B D and B A are regarded as resources available to HRATN managers and virtual attackers, the benefits of using these resources will be more significant when they are used in small amounts. However, as the use of these resources increases, the rate of return per unit resource will gradually decrease. This trivial result is not counterintuitive and justifies the results of this study. In addition, this paper validated the model on a real-world case with real background and compared it with a model for the classic LRP model that does not consider the influence caused by extreme events. The numerical calculation results showed that the strategy proposed in this research does not deviate too much from normal decision-making, thus also verifying the viability of the strategy proposed in this paper.

VIII. CONCLUSION
HRATNs have become the backbone of China's passenger network. A multimodal transportation network composed of these two modes of transportation can provide passengers with more secure and convenient transportation services. This research strengthens the connection between high-speed rail and aviation transport networks and mitigates the vulnerability of HRATNs by optimizing the locations of integrated hubs.
In this research, the impact range of extreme events is considered as a spatial range rather than as several nonadjacent facilities. On this basis, the impact of extreme events on the HRATN is reflected in the fact that the travel time of passengers increases. Then, with the goal of minimizing the worst loss, a tri-level model is constructed to optimize the location of the integrated hub, and the C&CG algorithm is applied to solve this model. To verify the correctness and practicality of the strategy proposed in this research, numerical verifications are carried out in an illustrative example case and a real-world case. The calculation results show that the strategy proposed in this paper has a certain practicality.
However, as a preliminary study, there are still some defects in this paper that need to be further improved upon in future work. First, we assume that the local extreme event has the same influence on all relevant edges; however, in reality, the impacts on different edges will differ. In addition, this paper does not consider the construction cost of the integrated hub. This construction cost includes not only the fixed cost of building the hub but also the necessity of building integrated comprehensive hubs in different cities. In future research, we will incorporate this factor into the outer layer model. The consideration of the above factors will be reflected in the constraints and objective functions of the model. Finally, we consider the possibility of different extreme events and describe the real scene with a set of extreme event scenarios with uncertain probability, which makes the strategy in this paper more practical.