A two-stage hierarchical model for spatial location and evacuation allocation problem of urban earthquake shelters: a case study in Central urban area of Yangbi county, Yunnan province, China

Abstract The scientific planning of urban shelters can greatly reduce property losses and casualties, and can improve the efficiency and safety of post-disaster evacuation. Based on the hierarchical characteristics of shelters and the shortcomings of previous studies on evacuation demands and blocked roads, this study proposes a two-stage multi-hierarchy planning method for shelters. In the first stage, the construction cost of shelters is taken as the objective function to complete the spatial configuration of shelters. In the second stage, the evacuation allocation model for refugees is established, which minimizes the total evacuation distance. The planning method also establishes the calculation method for dynamic high-precision refuge demands, and considers the blockage of evacuation roads caused by damaged buildings. The traditional genetic algorithm is improved and hierarchical 0-1 coding is adopted to solve the spatial configuration model of hierarchical shelters. An “exchange” heuristic algorithm is proposed to solve the hierarchical evacuation allocation model. Finally, Yangbi County in Yunnan Province, China, is taken as an empirical case to verify the applicability and accuracy of the proposed method. When the algorithm iterates to 150 rounds, the fitness function value converges to 3040675, indicating that 61 new shelters should be constructed to satisfy the refuge demands, and the minimum construction cost of these shelters is 3,040,675 RMB. The minimum evacuation distances in the first, second, and third evacuation stages are determined to be 51645192, 17351731, and 16136090 m, respectively.


Introduction
With the increasing complexity of natural and social systems, the frequency and types of disasters, especially earthquakes, which are highly uncertain, are increasing. The scientific planning and construction of urban shelters are significant to urban disaster prevention (Hu et al. 2012); they can greatly reduce property losses and casualties and improve the efficiency and safety of post-disaster evacuation. Two types of methods that solve the spatial configuration of shelters have been proposed. The first type of method includes the use of the analytic hierarchy process (Nappi and Souza, 2015), the weighted linear regression method (Kar and Hodgson 2008), the TOPSIS method (Chu and Su 2012), the fuzzy aggregation method (Ahmadi et al. 2021), and the grey evaluation method (Ma et al. 2011), which are integrated into a comprehensive evaluation system to complete spatial allocation. Candidate shelters are evaluated by their geological characteristics (e.g., distance from active faults) and socio-economic characteristics, which have been determined as quantifiable indicators of risk during evacuation (Chu and Su, 2010). The second type of method refers to the spatial allocation of public facilities such as hospitals, schools, and fire stations, based on which the locations of shelters are determined from candidate facilities by maximizing or minimizing the indicators. Researchers have proposed the classical set coverage model (Aly and White 1978), the maximum coverage model (Church and Velle 1974), the P-center model (Tsai and Yeh 2016), and the P-median model (Hakimi 1964). Based on the classical models, the objective functions and constraint conditions have been improved. The single-objective model (Ye et al. 2012) has been widely used to establish the location model, the main optimization objectives of which are the minimization of the construction cost, the number of shelters, or the total evacuation distance. Due to the complexity of the location models of shelters, researchers have further established multi-objective optimization models by integrating multiple optimization objectives, such as the construction cost (Praneetpholkrang et al. 2021), the evacuation distance (Hu et al. 2012), time (Sherali et al. 1991), the safety of evacuation roads (Zhang et al. 2015), and disaster risks (Barzinpour and Esmaeili 2014). After development and evolution, the competitive location problem (Fern andez et al. 2007), the coordination of the rescue resources (Yi and Ozdamar 2007;Ghasemi et al. 2022), the disaster evolution process (Gama et al. 2016), the rationality of the behavior of refugees (Ng et al. 2010), and traffic congestion (Sherali et al. 1991) have also been considered in the spatial allocation models of shelters.
To improve the computational efficiency, a series of optimization algorithms, such as the genetic algorithm (Hu et al. 2014), the simulated annealing algorithm (Ng 2010), the particle swarm optimization algorithm , the ant colony algorithm (Arnaout 2013), the tabu search algorithm (Salman and Y€ ucel 2015), and the epsilon-constraint method (Khalili-Damghani et al. 2022), have been used to solve multi-objective and multi-constraint shelter location models.
The spatial allocation model of shelters has experienced development from simple to complex, from single factors to complex urban environmental factors, and from the macro level to the micro level. However, the temporal and spatial variations of refuge demands, the hierarchy of shelters, and different evacuation stages have seldom been considered. Hierarchical facility location models (HFLMs) (Şahin and S€ ural 2007) have been widely used in education systems (Teixeira and Antunes 2008), medical systems (Şahin et al. 2007), and the distribution of relief goods (Widener and Horner 2011). However, few studies have applied HFLMs to the hierarchical modeling of emergency shelters. Chen et al. (2013) and Li et al. (2017) assumed that evacuees are assigned to different types of shelters and minimized the total evacuation distance. Zhao et al. (2018Zhao et al. ( , 2019 considered the minimization of the total evaluation time and the total shelter area as the objective functions, and proposed the MPSO-GA combined algorithm to solve the location and evacuation allocation model. The existing hierarchical models of shelters are characterized by the following disadvantages. (1) Due to the separation of workplaces and residences in the central area of the city, the population distributions are different during daytime and nighttime. For example, the number of people in the commercial area reaches its peak value during daytime and decreases at nighttime. Thus, refuge demands assessed by the permanent population (Yin et al. 1990;Fan et al. 2014) cannot reflect the temporal and spatial changes of the population, and cannot satisfy the refuge demands during the peak period. Second, the scale of calculation is not sufficiently refined. At present, the practice of planning earthquake shelters is generally based on experience and statistical data (Yin et al. 1990), and does not take into account the actual damage state of each building after an earthquake; thus, the refuge demands are not precise.
(2) The reduction of the effective width and blockade of evacuation roads after an earthquake has not been considered. Previous models ignored the possibility of obstacles on evacuation roads and assumed that all evacuation roads are passable, which is inconsistent with the authentic evacuation process. The evacuation distance is always the objective function for a shelter location model, and the total evacuation distance will be inaccurate if blocked roads are not considered, as pedestrians will be considered to use blocked roads. (3) When formulating the evacuation allocation scheme, current algorithms cannot realize the equilibrium between the service distance and capacity of shelters. Current algorithms follow the principle of the "minimum distance priority" or "maximum capacity priority" (Daskin and Jones 2008), which are prone to getting stuck in the local optimal solution.
Based on the hierarchical characteristics of shelters and the shortcomings of previous studies on evacuation demands and blocked roads, this study proposes a hierarchical spatial allocation model of shelters with the minimization of the construction cost as the objective function, as well as a hierarchical evacuation allocation model with the minimization of the evacuation distance as the objective function. Based on the Code for Design of Disasters Mitigation Emergency Congregate Shelter (Ministry of Housing and Urban-Rural Construction of the People's Republic of China 2015), shelters are divided into three levels: emergency evacuation and embarkation shelters (ESs), resident emergency congregate shelters (RSs), and central emergency congregate shelters (CSs). The first-stage model considers the coverage and capacity constraints of shelters, and takes the total construction cost as the objective function to complete the spatial allocation of shelters. Once the shelter locations and levels are determined, the corresponding allocation for refugees is then addressed. The secondstage model establishes the evacuation allocation model for refugees, which minimizes the total evacuation distance. Based on the hierarchy of shelters, the evacuation process is divided into three stages. The proposed model also considers the temporal and spatial variations of refuge demands and establishes a calculation method to determine the dynamic high-precision refuge demands. Moreover, this study also considers the blockage effect of evacuation roads caused by buildings damaged due to earthquakes. To solve the spatial allocation model of hierarchical shelters, hierarchical 0-1 coding is adopted to improve the traditional genetic algorithm. To solve the hierarchical evacuation allocation model, an "exchange" heuristic algorithm is proposed. The central urban area of Yangbi County, Yunnan Province, China, is taken as an empirical case. The detailed planning of hierarchical earthquake shelters is carried out with ArcGIS software, which involves the spatial distribution of shelters and the evacuation allocation scheme.

Study area
Cangshanxi Town, Yangbi County, Yunnan Province, China, was selected as the study area ( Figure 1). Yangbi County, which is located in the west of Yunnan Province, is squeezed by the Indian Plate and the Eurasian Plate. On May 21, 2021, a magnitude 6.4 earthquake occurred in Yangbi County. The epicenter (118 25 0 67 00 E and 24 25 0 02 00 N) was located in Cangshanxi Town, Yangbi County, and the focal depth was 8.0 km (Zhang et al. 2022). Three people were killed and 34 people were injured by the earthquake. The study area was the epicenter of the earthquake, and is an earthquake region with intensity 8. There are 2218 buildings in the study area, including residential buildings, schools, hospitals, office buildings, large shopping malls, etc. The study area has a high composite function, and reflects the assembly of several industries such as trade, finance, business, entertainment, and services. Thus, the population distributions during daytime and nighttime are different.

Data source
The hierarchical model requires the acquisition of spatial data, population data, and information about buildings in the study area. The basic information of the buildings, including the spatial location, building area, number of floors, floor height, structure type, and construction time, were derived from the Urban Rural Development Bureau of Yangbi County and the China Seismological Bureau. The spatial data were sourced from OpenStreetMap (Haklay and Weber 2008) and Google Maps, and remote sensing images and vector data were provided by the Urban Rural Development Bureau ( Figure  2). The detailed criteria for the adaptability evaluation of land resources for the construction of earthquake shelters include the slope and geological conditions, as well as the distance to a river or lake, oil or gas station, and high-tension line (Table 1). After the evaluation of the safety and adaptability of land resources, 113 alternative shelters were selected ( Figure 3). According to the value of the effective area, the candidate shelters were further divided into 10 CSs, which provide refuge services from level 1 to level 3, 39 RSs, which provide refuge services for levels 1 and 2, and 64 ESs, which provide refuge services for level 1. The shortest road network distances between the demand point and shelter, and those between different hierarchies of shelters, were determined based on the Dijkstra algorithm (Dijkstra 1959). The Dijkstra algorithm is a classical algorithm for searching the shortest path from one node to other nodes. The main characteristic of the Dijkstra algorithm is that it takes the starting point as the center and expands outward layer by layer (breadth-first search, BFS) until it extends to the end point. Table 2 reports the parameters for different hierarchies of candidate shelters.
To determine more accurate population data for which plots are considered as units, the population over 24 hours was determined based on the statistics of Chinese mobile signaling data from 15 June 2021, to 28 June 2021 ( Figure 4). However, the number of people in the plot based on Chinese mobile signaling data does not represent the total population, as the proportion of the population using Chinese mobile communication equipment in the study area is unknown. The nighttime population determined by the mobile signaling data was corrected with reference to the permanent population, after which the proportion of the population that accounts for China's mobile phone users in the study area was determined. Based on the permanent population data provided by the Statistical Yearbook of Yangbi City (2020), the  population determined by the mobile signaling data was corrected based on proportion. It was found that the number of Chinese mobile phone users in Cangshanxi Town accounts for 17.52% of the permanent population. The population at 4 a.m. was taken as the refuge demands at nighttime, and the population at 12 am was taken as the refuge demands at daytime (Table 2). Figure 5 demonstrates the population of each plot during daytime and nighttime. The refuge demands in different evacuation stages were then determined based on the calculation method established in Figure 5.

Framework construction and indicator selection
Based on the Code for Design of Disasters Mitigation Emergency Congregate Shelter, shelters were divided into three categories from a low grade to a higher grade, namely emergency evaluation and evacuation shelters (ESs), resident emergency aggregate shelters (RSs), and central emergency aggregate shelters (CSs). Based on Chinese codes, a three-level hierarchical location model and evacuation allocation model are proposed in this study.  Urban emergency shelter planning is divided into two stages: the first stage decides which shelters should receive refugees, and the second stage decides how to allocate refugees to the selected shelter. The spatial allocation of shelters is related to the location of emergency facilities, and the evacuation allocation of refugees is related to the emergency distribution. Most emergency facility location models evaluate the evacuation efficiency based on the evacuation time or distance (Toregas et al. 1971). The longer the evacuation distance, the greater the risk of disasters. Uncertain events, such as aftershocks, traffic congestion, and falling objects, will not only delay the evacuation progress, but will also expose the evacuees to unexpected risks, especially during an earthquake. Therefore, when considering the hierarchical structure of shelters, the most challenging work is determining how to minimize the total evacuation distance and satisfy the refuge demands in different evacuation periods. In addition, the construction cost of shelters is an important factor considered by government managers in planning practice. The budget for the construction of shelters is limited, and the maintenance cost of shelters is relatively high. Therefore, based on the principles of fairness (satisfying the refuge demands of all people), efficiency (completing the evacuation safely and quickly), and economy (minimizing the construction cost), the total evacuation distance and the construction cost of shelters are taken as the objective functions to establish a two-stage, three-level spatial allocation and evacuation allocation model for earthquake shelters. In the first stage, under the limitations of the capacity and service distance, the spatial configuration scheme of shelters with the minimum construction cost is determined. Once the shelter locations and levels are determined, the corresponding allocation for refugees is addressed. In the second stage, the optimal three-level evacuation allocation scheme is determined by taking the total evacuation distance as the objective function. Figure 6 illustrates the location model of the three hierarchies of emergency shelters, the coverage radius at different hierarchies, and the evacuation allocation process for refugees at different stages. As shown in Figure 7, all residents are allocated and transferred to the shelters immediately after an earthquake, which is the first stage of emergency evacuation. In the first stage, candidate shelters can be ESs, RSs, and CSs. With the continuation of time, refugees will require a higher level of refuge service. If the existing shelters cannot provide the required services, they will be transferred to a higher level of shelter. During this period, refugees will experience two re-locations (from ESs to RSs and CSs, and then from RSs to CSs). In the second evacuation stage, when ESs are no longer able to satisfy the refugees' demands for relief resources, refugees will be transferred from ESs to higher-level shelters (RSs or CSs). In the third evacuation stage, similarly, refugees in RSs will be transferred to CSs. If refugees can initially be assigned to a shelter that can satisfy higher refuge demands in the future, there is no need for resettlement, which indicates that higher-level shelters can provide the same refuge services as lower-level shelters. For example, CSs can be regarded as RSs and ESs, and RSs can be regarded as ESs, but the opposite is not true.
The model consists of the following four steps. 1. Basic data acquisition. Basic data, including population data, spatial data, and building information, are obtained by data processing. 2. Determining blocked roads. Based on the elastic-plastic analysis (Xiong et al. 2016;Lin 2017) of building groups, the inter-story drift ratios and velocities are determined, and the blocking effects on evacuation roads caused by collapsed buildings and falling non-structural components are quantified. Then, the blocked evacuation roads are identified. In addition, the reduction of the effective area of shelters caused by collapsed buildings and falling non-structural components is also considered. 3. Establishing a two-stage, three-level spatial allocation and evacuation allocation model of shelters. In the first stage, under the limitations of capacity and the service distance, the spatial allocation scheme of shelters with the minimum construction cost is determined. In the second stage, the optimal three-stage evacuation allocation scheme with the shortest evacuation distance is determined. 4. All candidate shelters are open spaces. 5. Solution algorithm design. An improved genetic algorithm with hierarchical 0-1 coding and an "exchange" allocation heuristic algorithm are proposed to solve the three-level spatial allocation and evacuation allocation model. The framework for the research method is illustrated in Figure 8.

The model assumptions
The three-hierarchy location model and evacuation allocation model are subject to the following four assumptions.
1. The location of refuge demands in each demand point is concentrated on the centroid. 2. When refugees are allocated to the assigned shelters, they follow the distribution routes recommended by the government, rather than randomly choosing potentially dangerous routes during the evacuation process. If the residents in demand point i are assigned to shelter j, it is assumed that the residents of demand point i evacuate to shelter j along the shortest accessible route. 3. In the first and second stages, people complete the evacuation on foot. In the third stage, vehicles are used to transfer people from RSs to CSs. 4. To efficiently and orderly evacuate refugees to the shelters, and to reduce the confusion in the emergency evacuation process, the decision-maker regards all the refugees in one demand point as an entity, and allocates them to the same shelter. In the first stage, refugees in the same demand point are assigned to the same ES, RS, or CS. In the second stage, refugees in the same ES are assigned to the same RS or CS. In the third stage, refugees in the same RS are assigned to the same CS.
To provide adequate refuge services, the spatial location and evacuation scheme should satisfy the following two constraints.
1. The limitation of the coverage radius. To rapidly complete the evacuation, it is necessary to evacuate refugees to neighboring shelters. Therefore, the distance between each demand point and its designated shelter should be within the service radius. 2. The limitation of the capacity of shelters. To ensure the basic needs of daily life, the refuge area required by refugees should not exceed the effective area of shelters. Table 3 illustrates the capacity and distance limitations of the three hierarchies of shelters (ESs, RSs, and CSs). Note: The table layout displayed in 'Edit' view is not how it will appear in the printed/pdf version. This html display is to enable content corrections to the  First stage: Establishing the hierarchical spatial allocation model of shelters In the first stage, the spatial allocation model of shelters with a hierarchical structure is established as follows.
The objective function is as follows.
The constraint conditions are as follows.
In Eqs.
(1)-(17), I ¼ fiji ¼ 1, 2, 3 … , mg is the set of demand points; for each demand point i, i 2 I. Moreover, J ¼ fjj j ¼ 1, 2, 3 … , ng is the set of candidate shelters; for each candidate shelter j, j 2 J. The candidate shelters are divided into three categories: E is the set of candidate ESs, R is the set of candidate RSs, C is the set of CSs, e represents ESs, r represents RSs, and c represents CSs (E, R, C 2 J, 8e 0 2 E 0 , 8r 0 2 R 0 , 8c 0 2 C 0 ). Moreover, k represents the hierarchy of shelters and the evacuation stage (k ¼ 1 refers to ESs and the first evacuation stage, k ¼ 2 refers to RSs and the second evacuation stage, k ¼ 3 refers to CSs and the third evacuation stage). K is the set of the hierarchy of shelters, w j represents the effective area of candidate shelter j, and A k represents the refuge area per capita in evacuation stage k (in the first stage, A 1 ¼ 1 m 2 ; in the second stage, A 2 ¼ 2 m 2 ; in the third stage, A 3 ¼ 4.5 m 2 ). Additionally, D k represents the coverage radius of shelters in evacuation stage k (D 1 ¼ 500 m, D 2 ¼ 1500 m, D 3 ¼ 5000 m), and U k represents the construction cost per unit area of shelters. The construction costs per unit area of ESs, RSs, and CSs are respectively 5, 15, and 40 RMB (Chen et al. 2013). Moreover, d k ij represents the road network distance from demand point i to shelter j, and d 1 i.e., d 1 ir, and d 1 ic, respectively represent the road network distances from demand point i to an ES, RS, and CS in the first stage. d 2 er and d 2 ec, respectively represent the road network distances from an ES to an RS or CS in the second evacuation stage, and d 3 rc represents the road network distance from an RS to a CS in the third evacuation stage. C ij is the matrix for the distance constraint. If the road network distance from demand point i to shelter j does not exceed the distance limitation, the value of c ij is 1; otherwise, it is 0. Additionally, P k j,daytime and P k j,nighttime, respectively represent the numbers of refugees received by shelter j at evacuation stage k during daytime and nighttime. Yk j is the binary location decision variable. If j is selected to construct a new shelter at level k, the value of Yk j is determined as 1; otherwise, it is 0. Finally, X k ij is the binary assignment decision variable. If refugees from demand point i are assigned to shelter j at evacuation stage k, the value of X k ij is determined as 1; otherwise, it is 0.
Equation (1) is the objective function of the spatial allocation model, which minimizes the total construction cost of shelters. Equations (2)-(17) are constraint conditions. Equations (2)-(4) represent the calculation to ensure that at each stage, refugees from each demand point are assigned to the same candidate shelters. Equation (2) ensures that refugees from the same demand point are assigned to the same ES, RS, or CS in the first evacuation stage, Eq. (3) ensures that refugees in the same ES are assigned to the same RS or CS in the second stage, and Eq. (4) ensures that refugees in the same RS are assigned to the same CS in the third stage. Equations (5)-(7) respectively determine the numbers of people accepted by candidate shelters during daytime in the first, second, and third evacuation stages. Equations (8)-(10) respectively determine the numbers of people accepted by candidate shelters during nighttime in the first, second, and third evacuation stages. Equations (11)-(13) represent the capacity limitation of the candidate shelters. Equations (14)-(16) are the calculation to ensure that refugees can only select shelters within the distance limitation. Finally, Eq. (17) represents the constraints for binary decision variables.
Second stage: Establishing the hierarchical evacuation allocation model The solution of the hierarchical model is divided into two types: location decision variables (Yk j) and allocation decision variables (X k ij). In fact, when the value of Yk j is determined, the value of X k ij is also determined. Because the spatial allocation model does not consider the optimization of the evacuation distance, the shelters are prone to evacuate to farther shelters. To ensure that refugees can quickly evacuate to the designated shelter at each stage, this section optimizes the total evacuation distance and formulates the optimal evacuation allocation scheme under the service distance and capacity limitations.
The objective function is as follows.
The total evacuation distance in the first stage during daytime and nighttime is as follows.
The total evacuation distance in the second stage during daytime and nighttime is as follows.
The evacuation distance in the third stage during daytime and nighttime is as follows.
The objective function for the minimization of the overall evacuation distance is as follows.
The constraint conditions are as follows. X In these equations, I ¼ fi j i ¼ 1, 2, 3 … , mg is the set of demand points; for each demand point i, i 2 I. Additionally, J 0 ¼ fj 0 jj 0 ¼ 1, 2, 3 … , n 0 g is the set of selected shelters; for each selected shelter j 0 , j 0 2 J 0 . The selected shelters are divided into three categories: E 0 is the set of selected ESs, R 0 is the set of selected RSs, C 0 is the set of selected CSs, e 0 represents an ES, r 0 represents an RS, and c 0 represents a CS (E 0 , R 0 , C 0 2 J 0 , 8e 0 2 E 0 , 8r 0 2 R 0 , 8c 0 2 C 0 ). Moreover, d k ij' represents the road network distance from demand point i to selected shelter j 0 , and d 1 i.e.,', d 1 ir', and d 1 ic' respectively represent the road network distances from demand point i to the selected ES, RS, or CS. Additionally, d 2 e'r' and d 2 e'c' respectively represent the road network distances from the selected ES to RS or CS in the second evacuation stage, and d 3 r'c' represents the road network distance from the selected RS to CS in the third evacuation stage. The meanings of the other parameters are the same as those in the spatial allocation model, but j, e, r, and c are respectively transformed into j 0 , e 0 , r 0 , and c 0 .
Equation (18) is the objective function of minimizing the total evacuation distance during three evacuation stages. Equations (19)-(21) are the calculation to ensure that at each stage, refugees from each demand point are assigned to the same selected shelters. Equations (22)-(24) determine the numbers of people accepted by the selected shelters in each evacuation stage during daytime. Equations (25)-(27) determine the numbers of people accepted by the selected shelters in each evacuation stage during nighttime. Equations (28)-(30) represent the capacity limitations of the selected shelters, and Equations (31)-(33) are the calculation to ensure that refugees can only select shelters within the distance limit. Finally, Eq. (34) represents the constraints for the binary decision variables.
It should be noted that refugees complete evacuation on foot in the first and second evacuation stages. In the third evacuation stage, refugees enter the CS by vehicles. To quantify the overall evacuation efficiency, the difference between the evacuation speeds of pedestrians and vehicles is considered. In Eq. (18), b represents the ratio of the evacuation speed of pedestrians to the vehicle speed, which is considered 0.2 in this study.

Identifying the blockage of evacuation roads
The nonlinear time-history analysis of building groups provides basic data to determine the distribution of falling objects. However, there are many buildings in the city, and it is difficult to obtain detailed information such as architectural drawings. Moreover, traditional finite element modeling would lead to a huge workload, which is not suitable for seismic analysis in urban areas. Based on the nonlinear multidegree-of-freedom (MDOF) shear layer model (Xiong et al. 2016) and the bending shear coupling model (Xiong et al. 2017), You Simulator v2.0 (Lin et al. 2020;Wang 2020) is adopted in this study to complete the nonlinear and elastic-plastic analysis of building groups. The basic information about the buildings, including the number of stories, height, building areas, structure type, and construction time, is inputted into You Simulator v2.0 to complete the nonlinear time-history analysis. The seismic responses of the buildings, including the absolute displacement, velocity, acceleration, and inter-story drift ratios, are determined. Then, the damage state of each building, namely basically intact, slightly damaged, moderately damaged, severely damaged, or collapsed, is determined. The next step is to determine the influence ranges of the collapsed buildings and falling non-structural components. According to a historical survey of building structures damaged during earthquakes, the impact range of debris after a building collapses does not exceed half of the building height (Du 2007). Therefore, when the building structure collapses, the impact range is determined as 1/2 of the building height. On the other hand, even if the building structure does not collapse, non-structural components such as infilled walls will fall off during an earthquake. In this study, the failure criterion of masonry infilled walls is determined by the standard proposed by Xu et al. (2016): where D wall,i is the inter-story drift ratio of the ith floor, and D a is determined by the Minimum design loads for buildings and other structures (ASCE 2010). After the failure of non-structural components, the fragments will be thrown horizontally. The influence ranges of the collapsed building structures and non-structural components are determined by Eq. (36): The building is damaged The building is not damaged i ¼ 1, 2, 3, … .n; j ¼ 1, 2, 3, … .m,where H represents the total height of the building. v i,j represents the velocity of the jth floor at time step i, h j represents the height of the jth floor, g represents the acceleration of gravity, n represents the total time step after the failure of non-structural components, and m is the number of floors. This method (Figure 9) reveals the influence ranges of collapsed buildings and falling non-structural components on evacuation routes. By drawing the buffer range in ArcGIS, the effective width is determined and the blocked evacuation roads are identified.

Estimation method of the high-precision and time-varying evacuation demands
The number of refugees during an earthquake depends on two factors, namely (1) the population quantity and distribution and (2) the damage state of buildings. The current practice of earthquake shelter planning determines the refuge demands via an empirical method (Xie et al. 1996). The actual damage state of each building during an earthquake is not taken into account, so the number of people who actually need refuge service is not determined.
To overcome the shortcomings of traditional methods, this study proposes a calculation method for dynamic high-precise refuge demands. "Dynamic" refers to the variation of the population over 24 h based on the statistics of mobile signaling data, which can satisfy the full-time refuge demands. Additionally, "high-precision" refers to the refined number of refugees in each building. First, the seismic damage state of each building is determined based on seismic elastic-plastic analysis. According to the seismic damage state, the buildings in which people need refuge service are determined. In the first evacuation stage, all people need refuge services and evacuate to shelters. In the second and third stages, the residents inside the buildings in an intact state do not need refuge service, whereas the residents inside buildings in other damage states require refuge service. Second, the dynamic population data of each building are determined. The FEMA P-58 report provides the curve of the variation of the population density with time for different building functions (FEMA 2012). Although the population densities of China and the U.S. are different, the proportion of the population in buildings with different functions can be calculated. The population in each building in the plot is then determined based on the proportional distribution of the total population of the plot. This method realizes the spatio-temporal coupling between the seismic damage state of buildings and the population distribution during an earthquake, and the accuracy of refuge demands is improved from the plot scale to the building scale.
The calculation method of refuge demands includes the following steps.
1. The population in demand area i (P i (t)) at moment t is determined based on mobile signaling data. 2. The population density (d b (t)) of building b at moment t is determined based on the population density of buildings with different functions provided by FEMA P-58. 3. The damage state of each building in the demand area j is determined based on seismic elastic-plastic analysis. The parameter k i is introduced. If building b does not suffer destruction, the value of k b is determined as 0; otherwise, it is 1. 4. The number of refugees in building b at moment t is determined by Eq. (37): where the number of buildings in demand area i is n (b ¼ 1, 2,., n), p b,i (t) represents the population of building b at moment t, m b represents the number of floors of building b, and A b represents the architectural area of building b. 5. The evacuation demands in demand area i at moment t are calculated as follows.
6. The refuge demands of the three-stage evacuation are determined based on the following assumptions. In the first stage, all people need refuge services. In the second and third stages, the residents inside non-damaged buildings do not need refuge services.
The detailed calculation process of the dynamic high-precision refuge demands is illustrated in Figure 10. Therefore, Figure 5 which illustrates the population distribution at daytime and nighttime also represent the number of refugees in the first stage. The number of refugees in the second and third stages is determined through the methodology in Figure 10.

Design of the solution algorithm
Researchers have implemented various methods to solve location and optimization problems. These methods can usually be divided into two types, the first of which uses an exhaustive method or mathematical solvers (Teixeira and Antunes 2008) to solve small-scale problems, but large-scale optimization problems cannot be solved. The second type of method uses imprecise heuristic algorithms, such as the genetic algorithm, simulated annealing algorithm, particle swarm optimization algorithm, ant colony algorithm, or tabu search algorithm, to determine the near-optimal solutions in large-scale scenarios. However, the application of the traditional heuristic search rules of optimization algorithms to the proposed model will lead to a high proportion of infeasible solutions and the requirement of additional computational work. The traditional genetic algorithm is improved in this section. According to the characteristics of the hierarchical location, hierarchical 0-1 coding is adopted to determine the location scheme of shelters with the minimum construction cost. Moreover, an "exchange" heuristic method is developed to optimize the total evacuation distance of the evacuation allocation scheme.

Solution algorithm of the hierarchical location model
The specific solution process of the improved genetic algorithm is explained as follows.
1. Chromosome coding. Real coding is a method widely used to represent the feasible solution, for which k represents the hierarchy of the candidate shelters (Boffey et al. 2003). However, under real coding, the selection of hierarchies of candidate shelters is the same, and many infeasible solutions will appear. To overcome these disadvantages and improve the calculation speed, the hierarchical 0-1 coding method is adopted in this study ( Figure 11). The mapping rules for the location scheme and gene coding are as follows. Each chromosome represents a location scheme of a shelter, one chromosome has n gene spots, and one gene spot represents a candidate shelter. The value of gene codes is 0 or 1; 1 means that the candidate shelter is selected, and 0 means that it is not selected. For the hierarchical 0-1 coding method, one chromosome is divided into three segments based on the three hierarchies of candidate shelters ( Figure 6). The difference between hierarchical 0-1 coding and ordinary 0-1 coding is that if the ith gene in the third segment is coded as 1, the ith gene in the second segment and the first segment is coded as 1. If the ith gene in the second segment is coded as 1, the ith gene in the first segment is also coded as 1; conversely, it is not. 2. Initial population generation. First, this procedure determines the hierarchy of shelters (ESs, RSs, or CSs) to divide the segments. The selection probability is set as 50%, and a chromosome (i.e., a set of shelter location schemes) is then generated. This random process is repeated N times to produce an initial population with a population size of N. After trial calculation, the population size in this study was determined as 200. 3. Determining the fitness function value of the location scheme. The fitness function is the criteria for screening chromosomes (i.e., the location scheme of shelters). In the spatial allocation model, the fitness function is the construction cost of shelters. When determining the fitness function value, a penalty function should be introduced to consider the situation in which the shelter location scheme does not satisfy the constraint conditions. If the location scheme does not satisfy the constraint conditions, 1,000,000 RMB is added to the fitness function value. 4. Selection operation. Roulette wheel selection and the elitism strategy are adopted in the selection operation. First, the fitness function value of each chromosome is determined by Eq. (1). The selected probability of each chromosome is proportional to its fitness function value. Second, the selection operation follows the elite retention strategy. The chromosome with the best fitness function value is Figure 11. Schematic diagram of hierarchical 0-1 coding.
directly copied into the next generation without mutation and crossover, which ensures that the optimal fitness function value of the generation is superior to that of the parent generation. 5. Crossover operation. Multi-point crossover can improve the diversity of solutions, while single-point crossover can accelerate the convergence of the algorithm and improve the computation speed. Therefore, the hybrid crossover strategy of multi-point and single-point crossover is adopted in this study. First, 20 crossover points in the parent chromosomes are randomly determined to exchange gene blocks. Then, single-point crossover is adopted to determine the crossover position of two matching chromosomes generated by multi-point crossover, and the progeny chromosomes are obtained. Figure 12 illustrates the process of the hybrid crossover operation. 6. Mutation operation. The progeny chromosomes not only inherit a part of the parent chromosomes, but also mutate with probability p m . First, a random number r d 2 (0,1) is generated for a gene point. If the value of r d is less than the probability p m , the value of the gene code is changed. The mutation operation is conducted on each gene point until a new chromosome is generated. Figure 13 illustrates the process of the mutation operation. 7. Algorithm termination judgment. Considering the calculation time and accuracy, the iterative termination conditions are determined as follows: (1) the maximum number of iterations, and (2) the optimal fitness function solution does not improve in future generations. When the algorithm terminates, the chromosome with the optimal fitness function value is determined as the solution for the spatial allocation model.

Solution algorithm of the hierarchical evacuation allocation model
The evacuation allocation model is also solved by the genetic algorithm, which is different from the solution algorithm of the location model in the following two aspects: (1) the use of hierarchical real number coding, and (2) the "exchange" heuristic search method.
1. Hierarchical real number coding. The spatial allocation model selects a group of shelters. The numbers of ESs, RSs, and CSs are respectively denoted as N es , N rs , and N cs . The selected shelters are renumbered. The serial number of the selected ESs is set from 1 to N es , the serial number of the selected RSs is set from (N es þ 1) to (N es þ N rs ), and the serial number of the selected CSs is set from (N es þ N rs þ 1) to (N es þ N rs þ N cs ). The chromosome has three segments (Figure 14), and the length of the chromosome is determined as the sum of m (the number of demand points), N es (the number of ESs), and N rs (the number of RSs). The value of the gene is the serial number of the selected shelters. The range of gene encoding in the first level is from 1 to (N es þ N rs þ N cs ), the range of gene encoding in the second level is from (N es þ 1) to (N es þ N rs þ N cs ), and the range of gene encoding in the third level is from (N es þ N rs þ 1) to (N es þ N rs þ N cs ). The value range of genes illustrates that high-level shelters can provide the same refuge services as lower-level shelters.
2. The "exchange" heuristic search method. The allocation variables (X k ij) generated by the spatial allocation model are able to guide the evacuation allocation of refugees. The allocation algorithms used in previous studies include "capacity priority" algorithms and "distance priority" algorithms. The allocation schemes based on "capacity priority" allocate refugees to the shelters with the largest remaining capacity. However, this scheme may assign refugees to a distant shelter rather than a neighboring shelter. The allocation schemes based on "distance priority" allocate refugees to the nearest shelter; this only satisfies the local optimization of demand points with high priority, and ignores the overall optimization. This section proposes a new "exchange" heuristic search after the "capacity priority" and "distance priority" heuristic search.   Figure 15 illustrates the process of the "exchange" allocation heuristic algorithm. First, the demand points are ranked in ascending order based on the number of covered shelters. If the number of shelters that cover two demand points is the same, they are arranged in descending order according to the refuge demands. This method gives priority to the demand points with the minimum number of covered shelters and the maximum refuge demands. Second, the first heuristic search is conducted following the "capacity priority" and "distance priority" rules. The "capacity priority" rules mean that the refugees from each demand point will be allocated to the shelter with the maximum remaining capacity. The "distance priority" rules mean that the refugees from each demand point will be allocated to the nearest shelter. By comparing the total evacuation distances determined by the two heuristic search rules, the evacuation allocation schemes with the relative shortest evacuation distance are determined. Third, an orderly local exchange operation is implemented on the allocation scheme solved by a "capacity priority" or "distance priority" heuristic search. The nearest shelters of each demand point constitute a set L. The exchange operation is the exchange of the shelters allocated to the two demand points. For example, demand point i 1 is initially assigned to shelter j 1 and demand area i 2 is assigned to shelter j 2 . After the exchange operation, demand point i 1 is assigned to shelter j 2 , and demand point i 2 is assigned to shelter j 1 . The prerequisite for performing the "exchange" operation is that if shelter j 1 belongs to set L, then demand point i 1 is assigned to shelter j 1 . If shelter j 1 does not belong to set L and inequality (39) is true, the shelter allocation of i 1 and i 2 will be exchanged; i 1 will be allocated to j 2 , and i 2 will be allocated to j 1 . Otherwise, if inequality (39) is not satisfied, the original allocation scheme is maintained.
where P i1 and P i2 , respectively represent the refuge demands from demand points i 1 and i 2 . Moreover, d i1,j1 , d i1,j2 , d i2,j1 , and d i2,j2 , respectively represent the road network distances from demand point i 1 to shelter j 1 , from demand point i 1 to shelter j 2 , from demand point i 2 to shelter j 1 , and from demand point i 2 to shelter j 2 .

Results and analysis
Based on the estimation method of the high-precision and time-varying evacuation demands established in Figure 10, the numbers of people who need the refugee services in the three-hierarchy location model are determined. Figure 16 demonstrates the numbers of people who need the refugee services in three levels. In the first stage, the numbers of people who need the refugee services is the total population of each demand block. In the second stage, the demand points are selected shelters in the first stage. In the third stage, the demand points are selected shelters in the second stage. Based on seismic elastic-plastic analysis, the seismic damage state of each building in the study area, namely basically intact, slightly damaged, moderately damaged, severely damaged, or collapsed, was determined ( Figure 17). After the occurrence of an earthquake with a intensity of 8, 126 buildings collapsed completely, and 360 buildings were seriously damaged. The influence range of damaged buildings after the earthquake was determined by the method presented in Figure 9, after which the blocked evacuation roads were identified (Figure 18).
In the spatial allocation model, the minimum number of shelters under the limitations of the capacity and service distance was determined. In the genetic algorithm, 200 groups of location schemes were randomly established as the initial population, and the crossover and mutation probabilities were respectively determined as 0.9 and 0.3. Figure 19 illustrates the relationship between the fitness function value and the number of iterations. After trial calculation, the fitness function value was found to reach the minimum value when the algorithm iterated to 150 rounds. The fitness function value converged to 3040675, indicating that 61 new shelters should be constructed to satisfy the refuge demands, and the minimum construction cost of the shelters is 3,040,675 RMB ( Figure 19). Moreover, 8 CSs, 22 RCs, and 31 ESs were found to jointly provide different levels of refuge services in the three evacuation  stages after the earthquake (Figure 20). Based on the "exchange" heuristic algorithm designed in Figure 15, a three-stage evacuation allocation scheme (Figure 21) was determined. The evacuation distances in the first, second, and third evacuation stages were respectively determined to be 51645192, 17351731, and 16136090 m.
As for the number of refugees received by each shelter, Figure 22 demonstrates the number of refugees at each shelter in the first, second and third stages respectively.

Comparison of two situations considering road blockages or not
Based on the seismic elastic-plastic analysis, the seismic damage state of each building was determined, and the blockage of evacuation roads caused by damaged buildings was identified. The location model and evacuation allocation model of multi-hierarchical shelters consider the influence of blocked roads during an earthquake. The results of the traditional location and evacuation allocation model without considering the effect of blockages were compared with the results of the proposed method. Table 4 and Figure 23 compare the construction cost and total evacuation distance with or without considering the blockade of evacuation roads. Without considering the blockade of evacuation roads, the construction cost and total evacuation distance were respectively found to be 85.4 and 87.0% of the values determined by the proposed method. Because refugees will avoid blocked roads even if the initial route is the shortest, it is necessary to consider the blockade effect of evacuation roads in the shelter location model and evacuation allocation model.

Comparing the results of different heuristic search methods
Based on the traditional "capacity priority" and "distance priority" heuristic search methods, this study proposed an "exchange" heuristic search method. Table 5 illustrates the results of the evacuation distance determined by the three different allocation algorithms. The evacuation distance based on the "capacity priority" rules was longer than that based on the "distance priority" rules. Moreover, the evacuation distance based on the "distance priority" rules was longer than that based on the "exchange" heuristic search. Thus, the proposed "exchange" heuristic search method effectively reduced the total evacuation distance.

Sensitivity analysis
In this study, the seismic intensity was determined as the key indicator in the sensitivity analysis. With the change of the seismic intensity, the damage state of building groups and the total length of blocked roads change. Figure 24 illustrates the damage state of buildings under different seismic intensities. Eight buildings were seriously damaged during an earthquake with intensity 7. During an earthquake with a magnitude of 8, 126 buildings completely collapsed and 360 buildings were seriously damaged. During an earthquake with a magnitude of 9, 772 buildings in the study area collapsed and 584 buildings were seriously damaged.
Based on the "Identifying the blockage of evacuation roads" section, the influence range of the damaged buildings and the blocked evacuation routes under different seismic intensities were determined. As is illustrated in Figure 25, during earthquakes with magnitudes of 7, 8, and 9, the length of blocked roads respectively accounted for 4.62, 13.73, and 51.26% of the total length of the evacuation roads. During an earthquake with a magnitude of 9, some pedestrians could not complete the evacuation because the demand points were surrounded by blocked roads. Therefore, the spatial allocation of shelters was carried out only for earthquakes with magnitudes of 7 and 8.  Table 6 compares the construction costs and evacuation distances under different seismic intensities. With the increase of the seismic intensity, the construction cost and evacuation distance were found to increase.

Discussion
Considering the different hierarchies of shelters, different evacuation stages, and the dynamic variation of refuge demands, a two-stage multi-hierarchy shelter location and evacuation allocation method was established in this study. In the first stage, the spatial allocation of shelters is completed with the minimum construction cost as the objective function. In the second stage, a three-stage evacuation allocation scheme is formulated with the total evacuation distance as the objective function. Compared with the traditional shelter location and evacuation allocation methods, the proposed method makes the following four contributions.
1. A calculation model for dynamic high-precision refuge demands was provided.
Based on mobile signaling data, the temporal and spatial variations of the population model are established to satisfy the full-time refuge demands. Moreover, based on seismic elastic-plastic analysis, the damage state of each building during an earthquake is determined. The damage state of a building and its population are coupled in the temporal and spatial dimensions. The accuracy of refuge demand is improved from the plot scale to the building scale.

2.
A two-stage multi-hierarchy refuge location and evacuation allocation model was presented. This study proposed a hierarchical spatial allocation model with the minimization of the construction cost as the objective function, as well as a hierarchical evacuation allocation model with the minimization of the evacuation distance as the objective function. Based on the Code for Design of Disasters Mitigation Emergency Congregate Shelter, shelters were categorized into three levels with different service distances and capacity limits. Based on the hierarchy of the shelters, the evacuation process was divided into three stages, which satisfies different refuge demands in different evacuation stages.  3. The blockage effect of roads after an earthquake was considered. Based on seismic elastic-plastic analysis, the seismic damage state of each building during an earthquake is determined. The influence range of collapsed buildings or nonstructural components is then determined, and the blocked roads are identified.
If the blockage effect of roads is not considered, refugees will enter the shelter along the shortest routes, which is inconsistent with the situation of disaster occurrence. After an earthquake, some roads are blocked and no pedestrians can get through, so refugees choose a relatively longer route to complete the evacuation. The consideration of the blockage effect of roads realizes the interaction between the disaster environment and shelters, and improves the authenticity and accuracy of the model. 4. A reasonable heuristic algorithm was used. The hierarchies of candidate shelters are the same under traditional real coding, which cannot satisfy the setting that high-level shelters can provide low-level refugee services. In this study, 0-1 hierarchical coding was adopted to solve the hierarchies of the location model. To solve the three-hierarchy evacuation allocation model, the "exchange" heuristic allocation method was proposed for use after the "distance priority" and "capacity priority" rules. The proposed "exchange" heuristic search method effectively reduces the overall evacuation distance.  Although a test was completed based on a relatively small-scale case to verify the applicability and accuracy of the proposed model, the proposed model is also suitable for large-scale problems. However, the implementation of the model also has limitations. Because this study is applicable to earthquake disasters and focuses on shelter planning, it cannot be extended to other types of natural disasters. For example, earthquake shelters are always built in open outdoor areas, while during other  disasters, such as floods, refugees enter indoor areas. Therefore, the criteria for choosing shelters should conform to the specific characteristics of disasters and different evacuation space environments. Another point to consider in future research is the uncertainty of evacuation. As the psychological status of refugees will affect their behavior, the effects of the evacuation behavior (herd behavior, dynamic selection of evacuation routes, preference of shelters, etc.) on refuge demands and the evacuation process should be further studied.

Conclusions
Based on the hierarchy of shelters, the evacuation process was divided into three stages, which satisfies different refuge demands in different evacuation stages. This study proposed a two-stage and multi-hierarchy shelter planning method. The firststage model considers the limitations of the service distance and capacity, and takes  the total construction cost as the objective function to complete the spatial allocation of shelters. Once the shelter locations and levels are determined, the corresponding allocation for refugees is addressed. The second-stage model establishes the evacuation allocation model for refugees, for which the minimization of the evacuation distance is taken as the objective function. The proposed planning method also considers the dynamically changing refuge demands and improves the accuracy and scale of refuge demands. The blocking effect of evacuation roads caused by damaged buildings during an earthquake is simulated, so the accuracy and scientificity of the proposed model are improved. Cangshanxi Town, Yangbi County, Yunnan Province, China, was selected as the study area to verify the applicability and accuracy of the proposed two-stage multihierarchy shelter planning model. This study improved the traditional genetic algorithm and adopted the hierarchical 0-1 coding method to solve the spatial allocation model of hierarchical shelters. The crossover probability and mutation probability were respectively determined as 0.9 and 0.3. After trial calculation, when the algorithm iterated to 150 rounds, the fitness function value reached the minimum value at 3,040,675 RMB. The findings indicate that 61 new shelters, including 8 CSs, 22 RCs, and 31 ESs, should be constructed to satisfy the refuge demands. The three-stage evacuation allocation scheme was determined based on the "exchange" heuristic algorithm. The evacuation distances in the first, second, and third evacuation stages were determined to be 51645192, 17351731, and 16136090 m, respectively.
Without considering the road blockage effect, the construction cost and total evacuation distance were respectively found to be 85.4 and 87.0% of the values determined by the proposed method. It is necessary to consider the blockage effect of evacuation roads in the shelter location model and evacuation allocation model. The evacuation distances determined based on "capacity first" and "distance priority" rules were longer than those based on the proposed "exchange" heuristic search. The proposed "exchange" heuristic search method effectively reduces the total evacuation distance.
Although this study made contributions in terms of evacuation demands, blocked roads, the hierarchical shelter location model, and an advanced heuristic algorithm, the development of this model remains characterized by some deficiencies. The proposed model does not consider the shelter location under multiple hazards, the evacuation behaviors of different pedestrians, or the change of the degree of congestion on evacuation roads during the evacuation process. The proposed shelter location model also does not concentrate on the emergency evacuation process. It is therefore suggested that further research be carried out in these respects.