Emergency shelters location-allocation problem concerning uncertainty and limited resources: a multi-objective optimization with a case study in the Central area of Beijing, China

Abstract

The special planning for shelters is an integral part of overall city planning and is a key issue in solving the urban security problem. During the planning of earthquake emergency shelters, it is usually necessary to prioritize the designated shelters to avoid wasting resources. However, in previous research initiatives, special attention has been given to the general location-allocation model without prior consideration of the designated shelters. This article presents a multi-objective supplemental location-allocation optimization model for earthquake emergency shelters that addresses the concerns related to the designated shelters. This model is believed to remain more consistent with the knowledge of shelter planners and policy makers by incorporating the uncertainties of the temporal dynamics of population distribution and the spatially heterogeneous distribution of evacuees under various earthquake damage scenarios by using the elliptical attenuation model for seismic intensity into the location-allocation problem. A case study of the central area of Beijing, China is presented. A modified particle swarm optimization algorithm is applied to solve the supplemental location-allocation model and the general location-allocation model, and the solutions obtained from the two models are then compared. The results show that the solutions of the supplemental location-allocation model significantly minimized the shelter area by prioritizing the designated shelters. By contrast, without prior consideration of the designated shelters, the shelter areas significantly increased, which slightly decreased the evacuation distance in the general location-allocation model. The policy makers could select solutions from the obtained minimal-area scheme using the supplemental location-allocation model that ensures acceptable evacuation distances and minimal budgets by prioritizing the designated shelters during the phased planning of shelters in Beijing.

Keywords:

• Earthquake emergency shelter planning
• spatial-temporal dynamic variations of evacuees
• multi-objective optimization
• supplemental location-allocation model
• Beijing

1. Introduction

Most cities in China have dramatically expanded over the last several decades, and the immeasurable value losses in light of frequent disasters may be enormous (Ding et al. 2005). Construction of emergency shelters has proven to be one of the most efficient ways to reduce human casualties and property losses during a disaster. For example, after the Great East Japan Earthquake on March 11, 2011, emergency shelters played an extremely important role, with approximately 250,000 people taking temporary refuge (BBC 2011). Typhoon Haiyan struck the Philippines in 2013, devastating communities and claiming more than 6,300 lives. These local communities are in need of emergency shelters that can withstand severe storm conditions. Therefore, GlobalMedic initiated a project rooted in disaster preparedness in 2015, and four evacuation shelters were constructed in highly vulnerable communities (GlobalMedic 2016). Powerful Typhoon Mangkhut slammed into the Philippines in 2018; however, there were no major damages or casualties in the region (CBS 2018), where residents poured into emergency shelters some of which were already filled to capacity (RNZ 2018). As an integral part of the overall urban planning, urban emergency shelters as safe spaces in a city, their detailed and special planning is necessary to make evacuation as efficient as possible, which has been a key issue to solve the urban security problem. Hence, the planning of emergency shelters has to be established, implemented, monitored, maintained and constantly updated by local city decision makers (Scheer et al. 2012).

The planning of disaster emergency shelters in cities must consider the interests of both planners and residents. Consequently, the multi-objective location-allocation model has received increasing attention for the planning of disaster shelters. In addition, several objectives are generally selected to minimize evacuation distance or time, facility area, construction costs and disaster risks in the shelter location-allocation problem. For instance, Alçada Almeida et al. (2009) selected four objectives to minimize the total distance from all of the population to shelters, the total risk of the primary path being impassable, the fire risks at shelters and the total time from the shelters to the university hospital; then, the authors proposed a multi-objective p-median model for evacuation during major fires and solved it by using a web-based decision support system. Tzeng et al. (2007) developed a deterministic multi-objective model for earthquake relief distribution with regard to minimizing the total costs and total travel time and maximizing the satisfaction of demand points in a planning period and solved it using a fuzzy multi-objective programming method. Doerner et al. (2009) presented a multi-objective optimization model for hurricane disaster shelters with the objectives of minimizing the weighted average of the minimum and maximum distances, the risk of tsunami events and costs and proposed a heuristic approach based on the non-dominated sorting genetic algorithm-II (NSGA-II) to provide optimal solutions. Yushimito et al. (2012) looked to maximize coverage of affected regions while minimizing human suffering through the use of a social cost function to create an uncapacitated nonlinear location optimization model that was solved by a Voronoi-based heuristic algorithm, and they presented an application to a real-world problem with data from Hurricane Katrina. Bozorgi-Amiri et al. (2012) proposed a robust multi-objective location-allocation optimization model in a stochastic environment by using discrete scenarios, which yielded more flexibility and practical aspects in real-world cases. The model was formulated as a mixed-integer nonlinear programming to minimize the sum of the expected total cost (which includes costs of location, procurement, transportation, holding, and shortage) and the variance of the total cost. Furthermore, a modified particle swarm optimization (PSO) algorithm was introduced to solve the model that exhibited successful performance in any problem size. In a subsequent study, the model added an objective of maximizing the satisfaction levels of the affected areas by minimizing the sum of the maximum shortages in the affected areas using a case study for disaster planning for earthquake scenarios in the region of Iran (Bozorgi-Amiri et al. 2013). Zhang et al. (2013) proposed a bottleneck Steiner tree based multi-objective optimization model to find the location of emergency facilities by minimizing the generalized cost of the system and minimizing the maximal distance between facilities and demand points. The model is solved using a population-based pattern-matching algorithm called cellular stochastic diffusion search. Barzinpour and Esmaeili (2014) developed a multi-objective mixed-integer linear programming model, with the objectives of maximizing the population coverage and minimizing the construction costs and traffic costs by using a virtual zoning approach to achieve humanitarian and financial goals. Hu et al. (2014) developed a bi-objective model for earthquake emergency shelter location and allocation to minimize the total evacuation distance and the total cost, subject to capacity and contiguity constraints, and the model was solved with a non-dominated sorting genetic algorithm based on NSGA-II was developed to solve the bi-objective model. Based on the scenarios, Rodríguez-Espíndola and Gaytán (2015) checked the flood areas by using GIS and proposed a model with the objective of minimizing distance and costs, and solved it using CPLEX 11.0 with the weighted-sum method. Zhao et al. (2015) introduced a multi-objective allocation model for earthquake emergency shelters, with the objectives of minimizing the total weighted evacuation time and the total area of shelters, and solved it using a modified PSO algorithm combined with the simulated annealing (SA) algorithm. In the follow-up work, the relationships between evacuation population size, earthquake emergency shelter capacity, and evacuation time were explored (Zhao et al. 2017), and the solutions of the model with four objective groups were given and compared using safety, capacity and investment evaluation index (Xu et al. 2018). Trivedi and Singh (2017) developed a hybrid multi-objective location-allocation model for emergency shelters based on analytic hierarchy process (AHP), fuzzy set theory and goal programming approach. It was solved by converting all objectives into a single objective function, and the results were validated by considering two real case studies of Nepal earthquakes in 2015. Zhao et al. (2017) proposed a location optimization model for emergency rescue facilities that consisted of three objectives representing service capacity, global efficiency and equity in the case of large-scale urban environmental disasters. Appropriate spatial representation and encoding strategy were designed and coupled with NSGA-II to solve the model. Based on a hypothetical disaster scenario, a case study of the Chaoyang district in Beijing was presented to demonstrate the methodology and explore the conflicting objectives. Ma et al. (2018) developed a multi-objective location-allocation model for earthquake emergency shelters using Wenchang of Hainan province, China as a case study. The two objectives were to minimize the total area of all shelters and to minimize the total evacuation distance, with the constraints of capacity and distance, and the model was solved by the modified PSO algorithm.

In recent years, there has been a growing number of studies that have focused on the multi-objective location-allocation problem for disaster emergency shelters. Nevertheless, uncertainty is an inseparable part of the location-allocation problem, which should be fully considered. Therefore, a multi-objective, multi-mode, multi-commodity, and multi-period stochastic model was presented by Najafi et al. (2013) that was capable of handling real-life emergency relief conditions under capacity and demand uncertainties. Jalali et al. (2016) developed a reliable bi-objective facility location model with provider-side uncertainty of facilities that was employed by a set of scenarios. In addition, a meta-heuristic algorithm called multi-objective biogeography-based optimization was employed to find a near-optimal Pareto solution, which was used to verify the solution that was compared with a non-dominated ranking genetic algorithm (NRGA) and a multi-objective SA algorithm. To incorporate uncertainty into a multi-objective model, Haghi et al. (2017) developed a robust multi-objective optimization model for logistics planning by defining a set of scenarios that accounted for the uncertainty parameters. The ε-constraint method was utilized, and the MOGASA and NSGA-II algorithms were proposed for the large-size instances to solve the model.

On the one hand, research has focused on considering the multi-objective location-allocation problem as purely static processes, whereas there has been mild interest in formulating and solving dynamic multi-objective optimization problems, especially for disaster emergency shelters. Because, most real-world optimization problems involve objectives, constraints, and parameters that constantly change with time (Deb and Karthik 2007). For example, the dynamic changes, quick responses and unpredictable events are the main features of the emergency location-allocation problem, which make them quite different from the traditional location-allocation systems (Zhang et al. 2013). In addition, the standard for urban planning on earthquake resistance and hazardous prevention (GB50413-2007) has clearly stated that it is advisable to consider the dynamic diurnal variations of citizens when estimating the number and distribution of evacuees (PRC Ministry of Construction 2007). Therefore, the uncertainties in this article incorporate the temporal dynamics of the population distribution and the spatial heterogeneity of evacuee distribution under various earthquake damage scenarios by using the elliptical attenuation model for seismic intensity to build a dynamic multi-objective optimization model. On the other hand, the previous studies on the shelter location-allocation problem have mainly dealt with the question of how the evacuees can be evacuated as fast as possible and where the candidate shelters can be located as efficiently as possible during an emergency, which can be summarized as the general location-allocation problem. However, a handful of designated shelters exist that have been planned and equipped with relief supplies. Therefore, the main aims for the future planning of disaster emergency shelters should be to supplement the candidate shelters by prioritizing the designated shelters due to the insufficiency of these shelters, which is consistent with the concept of sustainable development. For example, Ma et al. (2019) presented a hierarchical supplement location-allocation optimization model for disaster supply warehouses to achieve resource integration with the demonstration of the Beijing-Tianjin-Hebei region in China. To the best of our knowledge, no study has focused on the application of dynamic multi-objective supplemental location-allocation techniques for the planning of earthquake emergency shelters, which both makes up the research gaps of the uncertainties and designated shelters.

The estimation of spatially and temporally dynamic evacuees, the locations of designated shelters and candidate shelters, and the selection of evacuation routes for the evacuees have a great influence on the emergency transfer and resettlement in disaster relief practice. The number of evacuees is the most important input parameter for location-allocation models of earthquake emergency shelters, and the research on this parameter needs to be further refined. Therefore, in this study, we improved the spatially heterogeneous distribution model of evacuees under different earthquake scenarios by using the elliptical attenuation model for seismic intensity. In addition, the dynamic spatial-temporal variations of population exposure were fully considered to develop a conceptual model of the dynamic spatial-temporal distribution of evacuees based on subway smart card data and building type data. Finally, a scenario-based multi-objective supplemental location-allocation optimization model was developed by considering the designated shelters and the uncertainties of spatial-temporal dynamics with a real-world large-scale case study in the central area of Beijing. Notably, two contributions regarding uncertainties and the designated shelters were considered in this study.

The remainder of this article is organized as follows: Section 2 introduces the research methods of this study. Section 3 presents an application of the central area of Beijing as a case study and analyzes the results. In Section 4, the advantages and disadvantages of the research methods are discussed and the reasons are analyzed. And then, the future research direction is pointed out. Finally, the conclusions are given in section 5.

2. Methods

In this article, we consider the impacts of various objective and subjective factors on the spatial distribution of evacuees at the micro level. The objective factors include the destruction of community buildings at various intensity ranges caused by different earthquake scenarios at specific earthquake magnitudes and epicentre locations. The subjective factors include the willingness of residents to take refuge. The dynamic diurnal changes in the population at the macro level are also considered, and a conceptual model of the dynamic spatial-temporal distribution of evacuees is constructed (Figure 1). Thus, our methods contain two parts: (1) developing the temporally dynamic population distribution model, and (2) modelling the spatially heterogeneous distribution of evacuees.

Figure 1. A conceptual model of the dynamic spatial-temporal distribution of evacuees. Source: the authors.

At the macro level, the total population in a region directly affects the dynamic diurnal variations of the population. In a day, the population shows an inverted U-shaped fluctuation. Similarly, evacuees will also show the characteristics of temporally dynamic distributions.

At the micro level, by setting different earthquake strike scenarios, the exceeding probabilities of different types of buildings reaching different damage states are obtained by using the associated fragility curves. The proportions of different building types are combined to calculate the rate of building damage in each community due to objective factors under different earthquake damage scenarios. Simultaneously, we should fully consider the subjective factors, such as the willingness of the evacuees to remain in their houses under different damage states and the willingness of the evacuees to choose the designated shelters, to simulate the spatially heterogeneous distribution of evacuees.

2.1. Temporally dynamic population distribution model

In research of the dynamic diurnal variations of population in large cities, scholars have found that the population shows an inverted U-shaped fluctuation in a single day. In the same way, similar fluctuations would occur in evacuees. The method proposed by Ma et al. (2017) was used to calculate the dynamic temporal distribution of the population, as shown in Equation (1)(1) $$P_{\mathit{jt}}=P_{j0}+\sum _{n=1}^t\frac{P_{j0}}{P}\times (\frac{\mathit{PE}{X}_n-\mathit{PE}{N}_n}{\delta })\begin{array}{cc}& \end{array}t=1,2,\dots ,20$$(1) . The population parameter in the deterministic location-allocation problem remained constant by deciding on a single period of time. However, the real-world location-allocation problem for earthquake emergency shelters is a decision that has long-term effects; so, the parameters of the system, such as the number of evacuees, construction costs, management costs, and environmental factors, may vary over time, and facility additions can occur at different times. Consequently, a critical decision in the dynamic location-allocation model is not only where but also when to build a facility (Boonmee et al. 2017). (1) $$P_{\mathit{jt}}=P_{j0}+\sum _{n=1}^t\frac{P_{j0}}{P}\times (\frac{\mathit{PE}{X}_n-\mathit{PE}{N}_n}{\delta })\begin{array}{cc}& \end{array}t=1,2,\dots ,20$$(1) where P_jt is the population of community j during period t. For further insight, see Ma et al. (2017) for detailed information on the other variables.

2.2. Spatially heterogeneous distribution model of evacuees

The occurrence of an earthquake is random and uncertain. The epicentral intensities and isoseismal lines vary even for earthquakes with the same magnitude. The damages to community buildings may also differ due to their different structure types and locations. Consequently, the associated number of evacuees is characterized by spatial heterogeneity.

The isoseismic statistics in China indicate that the elliptical attenuation model is more reasonable but more complex than the point source attenuation model for describing the distribution of seismic intensity, especially for large-scale earthquakes (Shi et al. 2011). According to the previous research on the elliptical attenuation relationship, there are many empirical attenuation models at home and abroad. Based on the comparison of multiple attenuation models, an elliptical attenuation model in North China proposed by Cui (2010) is used in this study. Based on the strikes of the faults and the orientations of the maximum principal stress of the earthquake scenarios in the region, the isoseismal lines can be drawn to determine the intensity zones in different communities. In this study, we improved the spatially heterogeneous distribution model of evacuees proposed by Xu et al. (2018) by using the elliptical attenuation model.

The spatiotemporally dynamic evacuation needs in relation to the location-allocation problem for earthquake emergency shelters in this study had a spatiotemporal dynamic change ($N_{\mathit{jt}}^s$) characteristic, which was calculated using Equation (2)(2) $$N_{\mathit{jt}}^s=P_{\mathit{jt}}\times R_s$$(2) . The evacuation rates for each district at the different intensities of the two earthquake scenarios were calculated using Equation (3)(3) $$R_s=\gamma ({\phi }_1\sum _{k=\mathrm{A}}^{\mathrm{D}}{\alpha }^1{}_{\mathit{ks}}{\beta }_{\mathit{kj}}+{\phi }_2\sum _{k=\mathrm{A}}^{\mathrm{D}}{\alpha }^2{}_{\mathit{ks}}{\beta }_{\mathit{kj}}+{\phi }_3\sum _{k=\mathrm{A}}^{\mathrm{D}}{\alpha }^3{}_{\mathit{ks}}{\beta }_{\mathit{kj}})\begin{array}{c}\\ \end{array}k=a, b, c, d\begin{array}{c}\\ \end{array}s=A,B$$(3) based on the building types, fragility curve (Fang et al. 2011) and the empirical equation (Ding et al. 2014) that related seismic intensity and peak ground acceleration. (2) $$N_{\mathit{jt}}^s=P_{\mathit{jt}}\times R_s$$(2) (3) $$R_s=\gamma ({\phi }_1\sum _{k=\mathrm{A}}^{\mathrm{D}}{\alpha }^1{}_{\mathit{ks}}{\beta }_{\mathit{kj}}+{\phi }_2\sum _{k=\mathrm{A}}^{\mathrm{D}}{\alpha }^2{}_{\mathit{ks}}{\beta }_{\mathit{kj}}+{\phi }_3\sum _{k=\mathrm{A}}^{\mathrm{D}}{\alpha }^3{}_{\mathit{ks}}{\beta }_{\mathit{kj}})\begin{array}{c}\\ \end{array}k=a, b, c, d\begin{array}{c}\\ \end{array}s=A,B$$(3) where $N_{\mathit{jt}}^s$ represents the number of evacuees in community j during period t in scenario s, R_s represents the total evacuation rate due to multiple subjective and objective factors in scenario s, different $\phi$ values represent the proportions of people who are willing to evacuate at different house damage states, and $\gamma$ represents the proportion of people who are willing to evacuate to designated shelters (shelter evacuees). We do not consider the non-shelter evacuees who evacuate to hotels or the places of their friends and relatives. For further insight, see Xu et al. (2018) for detailed information on the other variables.

2.3. Supplemental location-allocation optimization model (SLA model)

The probabilities of earthquake occurrence are unknown, so the uncertain parameters are identified using discrete scenarios in our scenario-based location-allocation optimization problem for earthquake emergency shelters. Our scenario-based multi-objective SLA model prioritizes the designated shelters to make full use of them. The model will supplement the surrounding candidate shelters, which can be planned as new earthquake emergency shelters for the future when a designated shelter capacity is insufficient (can no longer serve any other community). (4) $$\min f_1=\sum _{i=1}^N{A}_i{Y}_i^s\begin{array}{cc}& \end{array}\forall j\in J\forall s\in S$$(4) (5) $$\min f_2=\sum d_{\mathit{ji}}{B}_{{}_{\mathit{ji}}}^s\begin{array}{cc}& \end{array}\forall i\in I\forall j\in J\forall s\in S$$(5) subject to (6) $$\sum _{j=1}^M{N}_{\mathit{jt}}^{s}L{B}_{{}_{\mathit{ji}}}^s-A_i{Y}_{{}_i}^s\le 0\begin{array}{cc}& \end{array}\forall i\in I\forall s\in S$$(6) (7) $$d_{\mathit{ji}}{B}_{{}_{\mathit{ji}}}^s-D_j\le 0\begin{array}{cc}& \end{array}\forall i\in I\forall j\in J\forall s\in S$$(7) (8) $$\sum _{i=1}^N{B}_{{}_{\mathit{ji}}}^s{Y}_{{}_i}^s=1\begin{array}{cc}& \end{array}\forall j\in J\forall s\in S$$(8) (9) $$Y_{{}_i}^s=1,\begin{array}{cc}& \forall i\in \overline{\mathit{DS}}\end{array}\forall s\in S$$(9) (10) $$\sum _{j=1}^M(A_i{Y}_{{}_i}^s-N_{\mathit{jt}}^{s}L{B}_{{}_{\mathit{ji}}}^s)<\min N_{j^{\prime }t}^s\begin{array}{cc}& \end{array}\forall i\in \overline{\mathit{DS}}\forall s\in S$$(10) (11) $$B_{{}_{\mathit{ji}}}^s=\{\begin{array}{l}1,\mathit{community}j\mathit{selects}\mathit{shelter}i\mathit{as}\mathit{emergency}\mathit{shelter}\mathit{in}\mathit{scenario}s\\ 0,\mathit{others}\end{array}$$(11) (12) $$Y_{{}_i}^s=\{\begin{array}{l}1,\mathit{shelter}i\mathit{is}\mathit{selected}\mathit{as}\mathit{emergency}\mathit{shelter}\mathit{in}\mathit{scenario}s\\ 0,\mathit{others}\end{array}$$(12) where

I is the set of candidate shelters: I = {1, 2, … i, … N}

J is the set of communities: J = {1,2, … j, … M}

J’ is the remaining communities that cannot be accommodated by the designated shelter

$\overline{\mathit{DS}}$ is the set of designated shelters: $\overline{\mathit{DS}}=\left\{i=1,2,\cdot \cdot \cdot ,i\in I|d\mathit{esignated} \mathit{shelters}\right\}$

S is the set of scenarios: S = {A, B, … s, … Z}

A_i is the area of candidate shelter i

d_ji is the length of the shortest path between community j and shelter i

D_j is the maximal evacuation distance for evacuees in community j

L is the smallest refuge area per capita (1 m²/person)

Our scenario-based multi-objective SLA model first allocates the evacuees into the designated shelters, and then, the model denoted in Equations (4)(4) $$\min f_1=\sum _{i=1}^N{A}_i{Y}_i^s\begin{array}{cc}& \end{array}\forall j\in J\forall s\in S$$(4) and (5)(5) $$\min f_2=\sum d_{\mathit{ji}}{B}_{{}_{\mathit{ji}}}^s\begin{array}{cc}& \end{array}\forall i\in I\forall j\in J\forall s\in S$$(5) selects new shelters and allocates the rest of the evacuees to make location-allocation decisions such that the total shelter area and evacuation distance are minimized, which are both subject to the capacity constraints expressed in Equation (6)(6) $$\sum _{j=1}^M{N}_{\mathit{jt}}^{s}L{B}_{{}_{\mathit{ji}}}^s-A_i{Y}_{{}_i}^s\le 0\begin{array}{cc}& \end{array}\forall i\in I\forall s\in S$$(6) and the distance constraints expressed in Equation (7)(7) $$d_{\mathit{ji}}{B}_{{}_{\mathit{ji}}}^s-D_j\le 0\begin{array}{cc}& \end{array}\forall i\in I\forall j\in J\forall s\in S$$(7) . The capacity constraints require that the areas of the shelters satisfy the daytime and nighttime demands of the evacuees. The distance constraint requires that the distance from community j to a potential shelter i should be less than the distance of the longest evacuation route returned by the algorithm. Equation (8)(8) $$\sum _{i=1}^N{B}_{{}_{\mathit{ji}}}^s{Y}_{{}_i}^s=1\begin{array}{cc}& \end{array}\forall j\in J\forall s\in S$$(8) ensures that each community will select only one shelter. Equation (9)(9) $$Y_{{}_i}^s=1,\begin{array}{cc}& \forall i\in \overline{\mathit{DS}}\end{array}\forall s\in S$$(9) ensures that the designated shelters are selected. Equation (10)(10) $$\sum _{j=1}^M(A_i{Y}_{{}_i}^s-N_{\mathit{jt}}^{s}L{B}_{{}_{\mathit{ji}}}^s)<\min N_{j^{\prime }t}^s\begin{array}{cc}& \end{array}\forall i\in \overline{\mathit{DS}}\forall s\in S$$(10) ensures that the designated shelters can be fully utilized within their service area. Equations (11)(11) $$B_{{}_{\mathit{ji}}}^s=\{\begin{array}{l}1,\mathit{community}j\mathit{selects}\mathit{shelter}i\mathit{as}\mathit{emergency}\mathit{shelter}\mathit{in}\mathit{scenario}s\\ 0,\mathit{others}\end{array}$$(11) and (12)(12) $$Y_{{}_i}^s=\{\begin{array}{l}1,\mathit{shelter}i\mathit{is}\mathit{selected}\mathit{as}\mathit{emergency}\mathit{shelter}\mathit{in}\mathit{scenario}s\\ 0,\mathit{others}\end{array}$$(12) define the decision variables in the model.

2.4. General multi-objective location-allocation optimization model (GLA model)

Most researchers have developed multi-objective models that are GLA models. The GLA model considers all of the designated shelters, parks, green spaces and playgrounds as candidate shelters to accommodate evacuees to have the same priority. Thus, the designated shelters are not the prior shelters to ensure that Equations (9–10) are satisfied, which would be unrealistic and wasteful. In this study, we also developed a GLA model using the same objectives (4-5) and constraints (6-8, 11-12) to compare the solutions obtained by the SLA and GLA models.

The PSO algorithm has been successfully applied to facility location-allocation problems (Bozorgi-Amiri et al. 2012). To effectively solve complex and discrete optimization problems, Zhao et al. (2015) modified the PSO algorithm for the location-allocation problem by applying a von Neumann structure in the former generations, adding global structure in the later generations, and combining with the SA algorithm to help obtain the global optimum. The modified PSO algorithm was used to solve the GLA and SLA models in this study.

3. A case study

3.1. City context and earthquake scenario hypothesis

Beijing is the capital of China and is located in the central part of the Yinshan-Yanshan seismic belt in the major earthquake area of North China, and a number of fault zones are distributed in the surrounding areas. According to historical earthquake situations, it is likely that a moderate or strong earthquake will occur in Beijing. The upper magnitude limit in the latent focal region in Beijing is 8.0, and most of the central area of Beijing is within the upper magnitude limit of 7.0 (Figure 2a) (Zhou 2016). Therefore, four scenarios based on the historical earthquake data with increasing numbers of evacuees were considered: an M7.0 earthquake in South Beijing during the nighttime (Bn) and daytime (Bd) and an M8.0 earthquake in Sanhe-Pinggu during the nighttime (An) and daytime (Ad). In addition, the number of evacuees in community j during the daytime and nighttime in scenario s are set to the calculated maximum number of daytime evacuees and the number of evacuees from the total resident population, respectively. In this study, it is assumed that the intensity of the isoseismic trend of scenario A is consistent with the strike of the fault zone because the M8.0 earthquake in Sanhe-Pinggu is the latest surface rupture event in the fault, and that of scenario B is consistent with the direction of the maximum principal stress because there is no information about typical faults and their activities.

Figure 2. The latent focal region and upper magnitude limits in Beijing. Source: Figure (a) is from the GB 18306-2015 Seismic ground motion parameters zonation map of China, modified by the authors. Figure (b): the authors.

This research considers a study area consisting of the six central districts of Beijing: Dongcheng, Xicheng, Chaoyang, Haidian, Shijingshan and Fengtai. This densely populated region covers an area of 1369 km², and the 6th National Census Data Bulletin recorded a resident population of 11.72 million in 2010. This region accounts for 8% of the area of Beijing, but approximately 61% of the total resident population is concentrated in this region, and population density is approximately seven times the average density of Beijing. The region covers 2039 communities and contains 39 designated earthquake shelters with the capacity to accommodate 1.51 million evacuees (China Earthquake Administration 2015) and has 17 subway lines with 173 subway stations (excluding the airport area) (Figure 2b).

3.2. Data

The data used in this study can be divided into three categories: temporally dynamic population distribution data (including the community population data and the subway smart card data), spatially heterogeneous distribution data of evacuees (including earthquake data, seismic fragility curves, building type data and evacuation willingness), and geospatial data. The proportions of the population who lived in different types of buildings as reported in the 6th National Census in 2010 for the various districts are used to represent the proportions of building type. See Table 1 for details.

Table 1. Data used in this study.

Database	Source	Year	Data Description
Earthquake data	Historical information; books	1057, 1679	Locations and magnitudes of earthquakes and fault trends
Seismic fragility curves (Fang et al. 2011)	Key Laboratory of Environmental Change and Natural Disaster of Ministry of Education, Beijing Normal University (ECND, BNU)	2011	Peak ground acceleration and exceeding probability
Building type data	The 6th National Census Data Bulletin	2010	Proportion of four different types of buildings
Evacuation willingness	A special investigation of resident evacuation conducted by the Committee for Technical Investigation on Countermeasures for Earthquakes in Tokyo	2008	Proportion of residents willing to remain in damaged houses; the proportion of evacuees willing to evacuate to the designated shelters (Zhao et al. 2014)
Community population data	Beijing Municipal Bureau of Civil Affairs	2012	Total population
Subway smart card data	Beihang Interest Group on SmartCity	2013	Bi-hourly population flow at each subway station (Ma et al. 2017)
Geospatial data	ECND, BNU	2014	Location of epicentres, communities, shelters, roads and subway stations

High spatiotemporal resolution evacuees are critical for emergency preparedness and response, including rapid risk assessment, modified evacuation planning and relief delivery so as to improve the regional security. However, population data are usually organized by coarse administrative units or census blocks such as streets, which are too large for identifying evacuation units. Spatially finer data at the community level are more valuable for better planning in disaster risk management. In addition, the dynamics of evacuees on a detailed time scale could serve as a valuable information source (Miyazaki et al. 2015) for disaster risk assessments. Our dynamic hourly population data are then used to estimate the spatially heterogeneous distribution of evacuees under earthquake scenarios at different times to create better preparedness and mitigation plans based on subway smart card data.

There are 39 designated earthquake emergency shelters in the study area. In which, Relic Parks of Yuan Dynasty City and Relic Parks of Imperial City are two banded green spaces and divided into several shelters according to the natural street boundaries. These shelters cannot accommodate the total numbers of evacuees in scenarios B and A, which exceed 1 and 2 million individuals, respectively. Therefore, in accordance with the earthquake emergency shelter (outdoor) programme guidelines (OPG) for the central area of Beijing, 317 open spaces, including parks, green spaces and playgrounds with areas greater than 1000 m² in the Dongcheng and Xicheng districts and 2000 m² in the other districts, were also considered as candidate shelters (Figure 2b). Both types of shelters were treated as temporary shelters, and the shelter capacities that were used in the models were 60% of the total areas excluding surface water areas, as recommended by the OPG.

The best evacuation route for each community to reach its assigned shelter (namely, d_ji) was calculated using the Dijkstra algorithm with the objective of minimizing the total route length. The longest evacuation route distances (D_j) returned by the algorithm for each district are as follows: 6600 m in Fengtai and 7082 m in Haidian, with 6000 m used for all other districts in the GLA model. The D_j values returned by the algorithm for the GLA model are the same with the SLA model when supplementing new candidate shelters. In addition, the D_j values are set to 6000 m in the SLA model for the designated shelters to ensure arrival after 1 hour of walking. Earthquake-induced road damage was ignored because it was assumed that people would evacuate on foot according to the temporary evacuation plan for the city.

3.3. Results

3.3.1. Evacuation rate at the district level

The exceeding probabilities of complete damage, heavy damage and medium damage to buildings can be calculated based on the fragility curve (Fang et al. 2011) and the empirical equation (Ding et al. 2014) that relates seismic intensity and peak ground acceleration. Simultaneously, as shown in Table 2, the forecast results of the average earthquake damage rates for the different intensities of the two earthquake scenarios were derived from the product of exceeding probabilities and the proportions of the four building types. Finally, the evacuation rates for each district at the different intensities of the two earthquake scenarios were calculated using Equation (3)(3) $$R_s=\gamma ({\phi }_1\sum _{k=\mathrm{A}}^{\mathrm{D}}{\alpha }^1{}_{\mathit{ks}}{\beta }_{\mathit{kj}}+{\phi }_2\sum _{k=\mathrm{A}}^{\mathrm{D}}{\alpha }^2{}_{\mathit{ks}}{\beta }_{\mathit{kj}}+{\phi }_3\sum _{k=\mathrm{A}}^{\mathrm{D}}{\alpha }^3{}_{\mathit{ks}}{\beta }_{\mathit{kj}})\begin{array}{c}\\ \end{array}k=a, b, c, d\begin{array}{c}\\ \end{array}s=A,B$$(3) . Table 3 shows the evacuation rates calculated for the different scenarios. In addition, Figure 3 shows the spatial distribution of the evacuation rates.

Figure 3. Spatial distribution of evacuation rate for each district at the different intensities of the two earthquake scenarios. Source: the authors.

Table 2. Forecast results of average earthquake damage rates (%) for the six districts in central Beijing, corresponding to the product of α and β.

Damage		Dongcheng	Xicheng	Chaoyang	Haidian	Fengtai	Shijingshan
VII	basically intact	67.82	71.66	79.30	79.94	77.91	80.75
	slight damage	17.30	15.64	13.02	13.05	13.73	12.74
	moderate damage	9.87	7.77	5.13	4.87	5.72	4.47
	extensive damage	3.55	2.67	1.56	1.44	1.82	1.27
	collapse	1.46	1.13	0.70	0.65	0.82	0.60
VIII	basically intact	45.05	48.30	55.68	55.87	54.10	56.53
	slight damage	23.49	24.06	24.15	24.65	24.39	24.80
	moderate damage	16.15	14.53	12.13	12.10	12.77	11.80
	extensive damage	8.91	7.20	4.94	4.77	5.47	4.45
	collapse	6.40	4.77	2.81	2.56	3.27	2.26
IX	basically intact	22.14	23.73	28.08	27.99	27.05	28.30
	slight damage	22.74	24.13	27.12	27.28	26.53	27.57
	moderate damage	23.12	24.05	24.35	24.96	24.59	25.19
	extensive damage	15.15	13.94	11.93	12.01	12.49	11.79
	collapse	16.85	13.02	8.23	7.71	9.33	6.99

Table 3. Evacuation rates (R_s) calculated using Equation (3)(3) $$R_s=\gamma ({\phi }_1\sum _{k=\mathrm{A}}^{\mathrm{D}}{\alpha }^1{}_{\mathit{ks}}{\beta }_{\mathit{kj}}+{\phi }_2\sum _{k=\mathrm{A}}^{\mathrm{D}}{\alpha }^2{}_{\mathit{ks}}{\beta }_{\mathit{kj}}+{\phi }_3\sum _{k=\mathrm{A}}^{\mathrm{D}}{\alpha }^3{}_{\mathit{ks}}{\beta }_{\mathit{kj}})\begin{array}{c}\\ \end{array}k=a, b, c, d\begin{array}{c}\\ \end{array}s=A,B$$(3) .

District	Rs (%)
	Scenarios Bn and Bd	Scenarios An and Ad
North Dongcheng	6.39	27.92
South Dongcheng	15.00	27.92
North Xicheng	4.93	25.00
South Xicheng	12.34	25.00
West Fengtai	9.71	9.71
East Fengtai	9.71	21.89
North Chaoyang	3.10	20.74
South Chaoyang	8.87	20.74
Haidian	2.91	8.59
Shijingshan	2.64	8.09
Average	5.95	17.35

Table 2 shows that the earthquake damage rates under the different intensities in the different districts have the following consistency rules: With the increase in earthquake strike strength (intensity), the proportions of basically intact and slightly damaged buildings in different districts decreases, while the proportions of moderately, extensively and collapsed buildings increases. Under the conditions of low earthquake strike strength (such as VII intensity), the buildings in each district were mostly intact, followed by slightly damaged, and the proportions of reaching moderately and extensively damaged were both below 10%, while the proportion of collapsed buildings was approximately zero. In contrast, in the case of high earthquake strike strength (such as IX intensity), the differences in the proportions of buildings with different states of damage were small for each district.

As seen from Table 3 and Figure 3, the evacuation rates were lower in scenario B than scenario A, and the average evacuation rate of the communities was 5.95%. South Dongcheng and Shijingshan had the highest and lowest rates of 15.00 and 2.64%, respectively. The evacuation rates were higher in scenario A than scenario B with an average evacuation rate of the communities of 17.35%, and the highest rate was 27.92% in Dongcheng and the lowest was 8.09% in Shijingshan. Compared with the other four districts, the evacuation rates in Dongcheng and Xicheng remained high even in the same intensity zone. The reasons for this result lie in the fact that the proportions of the population within building type A (mainly reinforced concrete structures) and type B (mainly brick-concrete structures) were relatively low, while the proportion in type C (mainly brick-wood structures) buildings was relatively high. Shijingshan had the lowest evacuation rate because it was located in low-intensity zones in both scenarios. On the other hand, type A and B buildings had relatively high proportions with relatively lower exceeding probabilities, while type C had lower proportions with higher exceeding probabilities. Type D buildings (mainly raw-soil and stone structures) were extremely rare in the six districts (Table 4).

Table 4. Proportions of building types in the six districts in central Beijing.

District	Building types (%)
	Type A	Type B	Type C	Type D
Dongcheng	41.79	24.43	33.68	0.10
Xicheng	39.81	34.64	24.30	0.11
Fengtai	47.08	37.48	14.74	0.71
Chaoyang	51.24	35.32	12.72	0.43
Haidian	47.91	40.26	11.56	0.22
Shijingshan	47.31	42.58	9.58	0.36

In ascending order of the sum of the proportions of type A and type B.

3.3.2. Temporal and spatial structures of the dynamic distribution of evacuees

The estimates of the number of evacuees in most location-allocation problems to date have focused on the total population recorded in census data, and consequently, the conventional shelter planning must allocate all of the population of the study area to shelters, leading to the need to build new and redundant shelters. In shelter planning for earthquake emergencies, however, it is usually unnecessary and unrealistic to build many earthquake emergency shelters because many people are either unwilling to evacuate or their buildings are basically intact or only slightly damaged. Consequently, research on the spatially heterogeneous distribution of evacuees under uncertainty is thus imperative for understanding actual evacuation activities. It is also critical to consider the dynamic diurnal variations of the urban population, especially for municipalities due to the frequent movements of the population. Thus, we will be able to begin monitoring and analyzing where and when evacuees allocate to each shelters and which evacuation route they choose based on the spatial-temporal dynamics of evacuees.

The study area can be divided into two intensities (VII and VIII in scenario B, VIII and IX in scenario A) based on the seismic intensity elliptical attenuation model. Figure 4 shows the population density maps of the study area for the two scenarios during nighttime and daytime. On the one hand, the distribution of evacuee density in the six central districts also shows a significant downtown-suburban-outer suburban pattern that is consistent with the associated population densities from the centre to the outer areas, especially during the daytime. On the other hand, there is a marked discrepancy in the evacuee density between the different intensity zones (VII and VIII, VIII and IX), which could be attributed to the spatial heterogeneity of both earthquake intensity and building types. Once the occurrence time, location and magnitude of an earthquake disaster are determined, the spatial-temporal dynamics of evacuees can be rapidly assessed based on our temporally dynamic population distribution model and spatially heterogeneous distribution model of evacuees.

Figure 4. Spatial distribution of the spatial-temporal dynamics of evacuees in the central area of Beijing. The inverse distance weighting (IDW) method was adopted to clarify spatial distribution. Source: the authors.

3.3.3. Comparison of the GLA model and SLA model

Consisting with the current shelter planning of Beijing, all of the modelling results were based on the integrity of the administrative divisions. Consequently, the evacuees cannot evacuate to other districts outside their own. The associated non-inferior solutions of the multi-objective GLA model and SLA model for the six districts in central Beijing are shown in Figure 5. It is not possible to improve each objective value in the alternative non-inferior solution set; consequently, the scatter plots for each scheme shown in Figure 5 maintain a descending order by distance and an ascending order by area. Notably, the solutions significantly differ between scenarios A and B and slightly differ between the daytime and nighttime due to the different distributions of the evacuees. A range of alternative location-allocation schemes exist that demonstrate the inherent trade-off between the conflicting objectives, which could be provided to city planners through the sets of non-inferior solutions based on the multi-objective technique. Decision makers can then apply their judgements to determine the desired balance between the conflicting objectives (Current et al. 1990) while simultaneously considering the differences between the different scenarios. Decision makers with sufficient budgets could select schemes with shorter total evacuation distances (at the bottom right of the plots in Figure 5) and larger shelter areas in the GLA model required in scenario A. Conversely, decision makers with limited budgets could select schemes with smaller total shelter areas (at the top left of the plots) in the SLA model, which would be a more suitable option for realistic emergency shelter planning by preferentially considering the designated shelters.

Figure 5. Non-inferior solutions for all schemes obtained using the GLA model (left) and SLA model (right) under various scenarios for the six districts in central Beijing. Each feasible scheme is represented by one point in the figure and corresponds to a particular spatial allocation result. Source: the authors.

The minimal-area and minimal-distance schemes corresponds to the conflicting objectives between the managers and the evacuees; namely, the minimal total shelter area and the minimal total evacuation distance, are usually used to demonstrate the quality of the solutions for various scenarios. Figures 6 and 7 compare the spatial allocation results, and Tables 5 and 6 present a further comparison of the solutions corresponding to the two schemes obtained using the two models.

Figure 6. Spatial allocation results corresponding to the minimal-area solutions obtained using the GLA model (left) and SLA model (right) for the various scenarios. The red lines are allocation lines for the designated shelters that can accommodate part of the communities for the different scenarios in the SLA model. Source: the authors.

Figure 7. Spatial allocation results corresponding to the minimal-distance solutions obtained using the GLA model (the left) and SLA model (the right) for the various scenarios. The red lines are the same as Figure 6 for different scenarios. Source: the authors.

Table 5. Comparison of solutions of minimal-area scheme for the two models.

District	Supplementary location-allocation model/General location-allocation model
	Bn			Bd			An			Ad
	AEDP	ASA	NSS	AEDP	ASA	NSS	AEDP	ASA	NSS	AEDP	ASA	NSS
Shijingshan	2535/1599	1090/1807+	6/10+	2337/2031	1090/579+	6/8+	2954/1659	1104/2494+	6/9+	2845/1428	1104/3723+	6/10+
Fengtai	1687/2216	4453/5642+	35/22+	1713/2424	4379/6165+	36/21+	1641/1494	7897/8463+	36/40+	1623/1617	5372/7771+	38/33+
Haidian	3005/1737	1013/1844+	20/47+	2990/1699	1021/1838+	21/50+	2758/1598	1490/11406+	27/69+	2579/1674	2876/11804+	31/68+
Chaoyang	2562/2934	8298/8776+	36/21+	2746/2919	7920/8776+	33/21+	2740/2161	8399/12214+	34/41+	2715/1775	9358/12902+	35/59+
Xicheng	1876/1418	1229/3052+	14/19+	2005/2235	1281/3221+	15/4+	2392/2504	2996/3221+	12/4+	2282/2613	3685/3221+	12/4+
Dongcheng	1506/1444	818/789+	6/13+	1485/1097	818/3988+	6/21+	1356/2687	1872/6665+	22/2+	1386/2780	4270/6665+	23/2+

‘⁺’ denotes the results of the GLA model without adding the area and number of the designated shelters. AEDP denotes the average evacuation distance (m) per person, ASA donates the average shelter area (m²) for each community, and NSS denotes the number of selected shelters. Bold part represents the solution about area. The grey entries represent the solution of the GLA model is superior to the SLA model without adding the area and number of the designated shelters.

Table 6. Comparison of solutions of minimal-distance scheme for the two models.

District	Supplementary location-allocation model/General location-allocation model
	Bn			Bd			An			Ad
	AEDP	ASA	NSS	AEDP	ASA	NSS	AEDP	ASA	NSS	AEDP	ASA	NSS
Shijingshan	2095/1047	7492/26331+	7/17+	1992/1038	7492/26331+	7/17+	2262/1042	20787/26331+	11/17+	2246/1033	20610/26331+	10/17+
Fengtai	1431/1318	14903/11237+	50/50+	1474/1296	14708/15062+	49/52+	1398/1454	14031/13540+	49/48+	1377/1381	14361/13987+	48/48+
Haidian	2794/1362	7948/16277+	32/79+	2763/1315	7948/16277+	32/79+	2533/1574	8856/12897+	39/70+	2477/1670	8876/12084+	41/70+
Chaoyang	2377/1506	16856/27777+	54/99+	2329/1447	17513/28508+	56/100+	2537/1912	11483/14820+	41/50+	2526/1755	11951/13145+	40/61+
Xicheng	1838/1143	4102/5864+	17/26+	1966/1197	4102/5734+	17/26+	1527/1296	5079/5508+	18/18+	1623/1571	4822/5161+	22/16+
Dongcheng	1506/783	818/10152+	6/29+	1485/765	818/10138+	6/28+	1282/1171	8278/10316+	22/31+	1357/1272	4417/10383+	24/31+

Bold part represents the solution about distance. The grey entries represent the solution of the SLA model is superior to the GLA model.

As shown in Table 5, Figure 5 and Figure 6, for the minimal-area scheme, the average shelter area (ASA) and the number of selected shelters (NSS) obtained using the SLA model are superior to those obtained using the GLA model, except for scenario Bn in Shijingshan and scenario Ad in Xicheng, as indicated by the dark grey entries. However, the actual ASA in the SLA model is also superior (scenario Bn in Shijingshan) or equal (scenario Ad in Xicheng) to the GLA model when the designated shelter area is added to the selected shelter area. Therefore, the solutions of the SLA model can significantly minimize the shelter area compared with the solutions of the GLA model by prioritizing the designated shelters. In addition, managers can prioritize the designated shelters to make full use of them.

As shown in Table 6, Figure 5 and Figure 7, for the minimal-distance scheme, the average evacuation distance per person (AEDP) achieved using the SLA model is inferior to that achieved using the GLA model. The SLA model generated longer distances when the evacuees first allocated and fully utilized the designated shelters, which inevitably led to the disappearance of the distance advantage in the SLA model compared with the GLA model. Nevertheless, the solutions of the SLA model are slightly inferior to those obtained using the GLA model (scenarios Bn and Bd in Fengtai), even slightly superior to the GLA model (scenarios An and Ad of Fengtai) as indicated by the dark grey entries; because, there are only five designated shelters in the relatively large district of Fengtai. In addition, the five designated shelters with small areas play a small role in the evacuation activities, which leads to the similar solutions in the two models. To sum up, solutions of the GLA model can slightly minimize the evacuation distance compared with the solutions of the SLA model without prior consideration of the designated shelters. In addition, the ASA and NSS significantly increased to slightly decrease the AEDP in the minimal-distance scheme. Therefore, managers can select solutions for the minimal-area scheme with acceptable AEDP and minimal budgets.

4. Discussion

Complex decision problems are frequently encountered in urban planning when focusing on the location-allocation problem for earthquake emergency shelters, and these problems typically involve the consideration of a wide range of incommensurable and conflicting objectives (Coutinho-Rodrigues et al. 2012). Most multi-objective location-allocation models to date do not consider the currently designated shelters as a priority. During the planning for earthquake emergency shelters, however, it is usually necessary to prioritize the designated shelters to avoid wasting resources. The key features that distinguish this research from previous works consist of the supplemental location-allocation optimization model that considers the real-world spatial-temporal dynamic distribution of evacuees and its implementation in the large-scale location-allocation problem of earthquake emergency shelters. As shown in Figures 6 and 7 of the spatial allocation results obtained using the SLA model, the first part evacuees evacuated to the designated shelters, and the rest of the evacuees that could not be accommodated evacuated to other candidate shelters based on the supplemental location-allocation optimization method. However, the post-selected candidate shelters may be better choices for the first part evacuees who were accommodated by the designated shelters. Consequently, the final spatial allocation results of the SLA model remained slightly irrational, as shown in Figures 6 and 7, and the minimal-distance solutions were inferior to those obtained using the GLA model, as shown in Table 6. These results were directly related to the priorities of the designated shelters.

In future studies, in addition to considering the dynamic diurnal variations of the evacuees due to building damage, the shortage of essential services (e.g., water, power) (Li et al. 2017) could further influence the dynamic temporal changes of the evacuees during a longer period after disasters (1 week, 1 month, etc.). Consequently, earthquake shelters could be classified as immediate, short-term, and long-term, which satisfy the different demands during the associated evacuation phases (Chen et al. 2013). Finally, a dynamic hierarchical supplement location-allocation optimization model could be developed to provide a basis for the phased planning of different levels of shelters in Beijing. Nevertheless, the location-allocation problem for earthquake emergency shelters is still a multidisciplinary and complicated systematic project.

5. Conclusion

This article presents a dynamic multi-objective supplemental location-allocation optimization model for earthquake emergency shelters in the case of evacuation via walking that addresses the concerns related to the designated shelters, and this model is believed to be more consistent with the knowledge of shelter planners and managers. This article also incorporates the uncertainties of the temporal dynamics of the population distribution and the spatial heterogeneity of evacuee distribution under various earthquake damage scenarios by using the elliptical attenuation model for seismic intensity into the location-allocation problem. Solutions were derived by applying a modified PSO algorithm to both the SLA model and the conventional GLA model, and the comparison of the results revealed the advantages of the SLA model over the GLA model. A case study in the central area of Beijing was presented to demonstrate the proposed methodology. The solutions of the minimal-area scheme in the SLA model distinctly minimized the shelter area compared with the GLA model by prioritizing the designated shelters. By contrast, the solutions of the minimal-distance scheme in the GLA model, without prior consideration of the designated shelters, slightly minimized the evacuation distance compared with the SLA model. Meanwhile, the ASA and NSS significantly increased to slightly decrease the AEDP in the minimal-distance scheme. Therefore, managers could select the solutions obtained for the minimal-area scheme using the SLA model with acceptable AEDP and minimal budgets by prioritizing the designated shelters. This study aims to make up for the lack of research on supplemental location-allocation planning throughout the world, substantially improve the efficiency and professionalism of disaster reduction and emergency response in China, and provide effective planning strategies and enlightenment for urban planning with respect to public safety in Beijing.

Acknowledgments

The subway fare smart card data were provided by Beihang Interest Group on SmartCity, and we are grateful for their permission to use it in this paper. We also thank the Beijing Municipal Bureau of Civil Affairs for providing access to community population data.

Disclosure statement

No potential conflict of interest was reported by the authors.

Additional information

Funding

This study was funded by the Ministry of Science and Technology, China [2018YFC1508802, 2016YFA0602404]; National Natural Science Foundation of China [41571488, 41201547]; Ministry of Education and State Administration of Foreign Experts Affairs, China [B08008].