Seismic hazard prediction is essential for earthquake preparation in urban areas. Prediction methods based on refined the finite element method (FEM) models and time-history analysis require higher calculation, which is hard to apply universally. In this paper, by applying two simplified models and time-history analysis, the authors were able to develop a refined seismic hazard prediction system for urban buildings using MATLAB and OpenSees. The system is modular in design with a unified data interface to transfer data, where only the macro parameters of buildings are required to complete the computation work. Taking the campus of Tongji University as a test area, the system was applied to predict seismic hazard and compare it with the results of the vulnerability method by applying the seismic damage index to verify its reliability. The results show that the system can meet the needs of accuracy, efficiency, and visualization of seismic hazard prediction in urban areas.

1. Introduction

With the rapid development of the economy and expansion of urbanization in China, along with increased population density and urban building, a sudden unforeseen earthquake would undoubtedly magnify the resulting hazard [1]. In order to improve overall seismic capacity and minimize seismic hazards in urban areas, it is of considerable significance to analyze and study every weakness in earthquake resistance structures and formulate corresponding disaster prevention plans.

Seismic hazard prediction for buildings is an essential work for urban earthquake mitigation [2]. It refers to the quantity distribution of different damage states in a specific area (sample space) under the given seismic intensity; this distribution, called the seismic hazard matrix, is obtained through the calculation and statistics of typical samples [3]. The purpose is to examine the global seismic hazard level of area and the difference of seismic performance among various buildings. Due to the destructive nature of an earthquake, it is impossible to conduct physical experiments during the event; thus, numerical simulations have become an essential means of seismic hazard prediction. The most widely used methods are based on a damage vulnerability matrix [4] and the capacity spectrum [5, 6].

The basic principle of the damage vulnerability matrix is to select various factors that affect the seismic performance of buildings based on historical seismic data and then obtain the damage vulnerability index of the target area according to the actual situation of each factor [7]. This method has been widely applied due to its simplicity and usability [810]. However, the method ignores the duration and the spectral characteristics of ground motions and cannot reflect the damage state of specific buildings [11]. The method based on the capacity spectrum, also known as pushover analysis, uses the earthquake response spectrum as an input to investigate the vulnerability of specific buildings [12]. The excess probability of different damage states is calculated by the intersection of the capacity curve and the demand spectrum (performance point), and then the seismic vulnerability curve expressed by the spectral displacement is constructed to determine the seismic vulnerability. Based on the mature theory and detailed structural parameters, many countries and regions, including America, Canada, China, Israel, and others, have developed seismic damage prediction systems based on the capacity spectrum [1317]. However, this method also has limitations. It is challenging to consider the influence of high-order mode shape on the seismic response and the time-domain characteristics of ground motions, such as the velocity pulse.

In order to avoid these disadvantages, with the in-depth study of performance-based seismic theory [18], the following two improvements have emerged: (1) using a refined model (such as the finite element model, FEM, and discrete element model, DEM) to simulate the structure and (2) applying time-history analysis. Compared with the previous methods, the property of materials and the characteristics of ground motion can be considered more reasonable. Due to the higher accuracy in the theory, there was an extensive application on seismic hazard prediction for a single building. Castaldo et al. [19] performed incremental dynamic analysis (IDA) to evaluate the seismic fragility of a nonlinear hardening and softening structure equipped with friction pendulum isolators. Hajipour et al. [20] constructed FEM models for seven steel buildings with RC shear wall in OpenSees considering the nonlinear geometric effects of materials and employed time-history analysis to attain the seismic performance. Recently, with the rapid improvement of the computational ability of computers, scholars have begun to conduct time-history analysis with refined models in regional seismic hazard prediction gradually. Mahin et al. [21] applied the computational modeling and simulation center (Simcenter) to simulate the entire seismic process and evaluate the seismic loss for over 1.8 million buildings in the San Francisco Bay Area. The Institute of Earthquake Engineering (ERI) [22] developed an integrated earthquake simulation (IES) system and predicted the seismic hazard of Tokyo based on time-history analysis with the fiber beam models and DEM models. Abdurrahman et al. [23] developed a new version in MATLAB, including site response analysis and structural analysis, and applied the system in Istanbul. On the other hand, due to a large number of buildings in urban areas, such systems require extensive computation, although with the emergence of high-performance computing methods such as parallel computation [24], to meet the timeliness of prediction, it needs to be conducted by supercomputers. Moreover, it requires detailed information of buildings in component level to establish the refined models, which takes lots of manpower and time, and it is hard to collect all the data needed for thousands of buildings in the urban area.

Therefore, based on the above background, the authors of this study developed a refined seismic prediction system suitable for urban buildings. The system uses the Geographic Information System (GIS) to acquire and manage urban building information and visualize prediction results. For different buildings, appropriate simplified models are applied for time-history analysis, which dramatically improves computational efficiency. Based on OpenSees and MATLAB, the automation of hazard prediction and data processing is realized, which reduces manual intervention. Finally, in order to verify the reliability of the system, the campus of Tongji University, Shanghai, was selected as the test area for this study’s new seismic prediction system. It was due to the university’s urban setting and the authors’ affiliation with the Shanghai Institute of Disaster Prevention and Relief located on campus. A comparison between prediction results and results of the vulnerability method showing vulnerability indicates that the system can realize rapid, refined, and visualized seismic hazard prediction in urban areas and can support relevant departments and organizations in disaster prevention planning and emergency rescue decision-making.

Specifically, in Section 2, the frame of the proposed methodology is presented, and the simplified models are introduced briefly. In Section 3, the technical framework and data flow of each module in the prediction system are introduced. Finally, a regional seismic hazard prediction was performed for the test area by applying the proposed system as a case study, and the validation of the system is discussed in Section 4.

2. Methodology

2.1. Methodology Framework

There are three main steps for the seismic hazard prediction based on time-history analysis, as illustrated in Figure 1.

Step 1. Data collection.

It is the essential work of seismic hazard prediction. The most effective way to obtain and manage the building data is to apply the Geographic Information System (GIS). Numerous cities have already implemented their GIS databases, where the macro building parameters (such as height, area, and construction time) could be directly acquired. Otherwise, data could be supplemented by field investigations and CAD drawings.

Step 2. Seismic response analysis.

After obtaining the macro building parameters, the calculation model of buildings in the target area can be established. Two simplified models are utilized in the paper in accordance with building characteristics, which will be introduced in Section 2.2. Subsequently, a nonlinear time-history analysis is conducted to obtain the seismic response of each building under earthquake scenarios in different return periods.

Step 3. Seismic hazard assessment.

Following the seismic response analysis, the limit of each damage state is firstly calculated. Based on the damage limit, the probability of different damage states of each building is determined according to the seismic response of each story. Finally, the prediction results of the target area are input into the GIS database and visualized to illustrate the seismic hazard intuitively for those who lack relevant knowledge.

2.2. Calculation Model

As mentioned in the Introduction, in order to balance calculation accuracy and efficiency, two multiple degrees of freedom (MDOF) models are applied in the system, the parameters of which can be determined by macro building attributes.

For a large number of regular multistory buildings in the urban area, mainly including the RC frame and masonry structure, the deformation characteristic under the lateral load mainly consists of shear deformation. Therefore, a building can be simplified into a nonlinear MDOF shear (NMS) model [25], as illustrated in Figure 2(a). Compared with the single-degree-of-freedom model, the influence of the high-order mode shape on multistory buildings can be considered with the NMS model. The models assume a rigid story and that the mass of each story is concentrated so that each story can be simplified into a mass point, which are connected by a shear spring between them. Notably, the skeleton line of the spring adopts the trilinear skeleton line recommended in the HAZUS report [26], as illustrated in Figure 2(c). Previous studies [27, 28] have demonstrated that the shear spring of the NMS model could reflect the elastoplastic characteristics of the structure, which can simulate the seismic performance of multistory buildings.

For high-rise buildings, where the effects of bending deformation are especially critical, a nonlinear MDOF flexural-shear coupling (NMFS) model was chosen [29], as illustrated in Figure 2(b). The NMFS model discretizes each story into a nonlinear bending spring and a shear spring. The bending spring is utilized to simulate the bending deformation of the shear wall, while the shear spring is used to simulate the shear deformation of the frame. Two springs are connected with a rigid link to ensure horizontal displacement coordination. Similar to the NMS model, the trilinear skeleton line is selected as the skeleton line of springs. In order to verify the accuracy of the simplified models, an NMFS model, an NMS model, and a refined FEM model of a 13-story reinforced concrete frame-shear structure were established and time-history analysis was performed. As shown in Figure 2(d), the result of the NMFS model is in good agreement with that of the FEM model, but the result of the NMS model is far from the above two. It is because the NMS model cannot consider bending deformation. Therefore, the NMFS model is suitable for the seismic hazard prediction of high-rise buildings.

The particular calibration methods of the two models can be found in [25, 29].

3. Function Realization

Taking OpenSees 2.5 and MATLAB 2017a as the primary development platforms, the system was divided into three modules: (1) data preprocessing module, (2) structural computing module, and (3) postprocessing module, as shown in Figure 3. The three modules are independent of each other and use a unified data interface to transfer data through files. Adopting such an overall framework has the following advantages: (1) different modules are not subject to the limits of other modules, which facilitates function expansion and upgrade; (2) a unified data interface between different modules ensures consistency of data transfer between modules and facilitates data viewing and processing; and (3) using files for data transfer reduces system memory consumption and improves computational efficiency.

3.1. Data Preprocessing Module

The function of the data preprocessing module is to obtain building attributes from GIS data and calculate parameters of finite element models. The module is mainly composed of three submodules: GIS data-calling, model parameter generation, and data output, as shown in Figure 4.

The GIS data-calling submodule reads the building attributes from GIS data and stores them in local variables, which are called up by MATLAB later. An .xlsx file is used as the exported data format. On the one hand, the .xlsx format is highly compatible and easy to read and has less information loss during the conversion. On the other hand, some fields are often not needed for data preprocessing in GIS attribute tables. The use of Excel also facilitates the modification and management of building attributes. Selected .xlsx files were checked and modified. They generally contain fields such as building code, construction time, and structure type.

The primary function of the model parameter generation submodule is to calculate the resilience parameters required to calibrate the MDOF model. It is composed of several different functions, where the NMS and NMFS functions are used to calibrate the NMS and NMFS models. The NMS function first calculates the shear stiffness of the model and then chooses concrete or masonry function to calculate the resilience parameter according to the structure type. The mainframe of the concrete/masonry function is a loop command that determines the parameters of each story. In the loop body, the base-shear function is first called to calculate the design bearing capacity through the base-shear method. Then, the yield/peak/ultimate function is called to calculate the yield/peak/ultimate capacity and displacement of the skeleton line. Finally, the degradation function is called to determine the hysteresis parameter by the lookup table method [30]. The particular calibration methods of NMS model can be found in [25]. The data flow of the NMS function is shown in Figure 5.

The framework of the NMFS function is more complicated. The elastic parameters of the model should be determined first, which is done by the elastic function. The shear stiffness and bending stiffness can be calculated with formulas (1) and (2). Since these are multivariate nonlinear equations, the Newton iteration method can be used to find their approximate solutions:where is the order period of the structure and is the eigenvalue coefficient related to the order of vibration, which is the function of bending-shear stiffness ratio . The definition of is as follows, where is the story height of the model:

It should be noted that the mode decomposition method is applied to calculate the yield parameters of models. Therefore, first, the “! OpenSees.exe modal-analysis.tcl” command in the NMFS function is run to call OpenSees to perform modal analysis. After obtaining the model shape and period of models, the design bearing capacity can be calculated through the mode-superposition response spectrum method. Finally, the parameters of the skeleton line and degradation are determined. The particular calibration methods of NMFS model can be found in [29]. The data flow of the NMFS function is shown in Figure 6.

The data output submodule summarizes the calculated resilience parameter matrix, generates separate submatrices for each building, and saves them in files for the structural calculation module to read.

3.2. Structural Calculation Module

The structural calculation module based on OpenSees uses nonlinear time-history analysis to calculate the seismic response of the NMS or NMFS model under different seismic intensities. The module is the core calculation and analysis part of the system, including the following four submodules:(1)The information reading submodule reads the model parameters of each building and ground motion data generated by the data preprocessing module through files according to the number of buildings. After that, parameters are stored in the matrices and then assigned to global variables by the “lassign” command.(2)The modeling submodule retrieves model parameters from memory; generates the nodes, elements, and boundary conditions of the calculation model according to the modeling rules; and defines the mass, stiffness matrix, and material recovery force characteristics (skeleton line and hysteresis model) of each node. Since the calculation model assumes that the stiffness and mass of each story are evenly distributed along with the height, determining the loop body controls each parameter in the modeling submodule.(3)The analysis submodule defines the static analysis conditions of the self-heavy load and the time-history analysis conditions of different ground motion periods and determines the selection of the degree of freedom in the calculation of the constraint function, iterative algorithm, control algorithm, and system equation. The dynamic analysis uses the Newmark implicit integration method. Each seismic wave is analyzed in 5000 steps, and each step is 0.01 s.(4)The recording submodule records the seismic response of structural models, such as node displacement, inter-story displacement angle, and mode shapes of each order, and saves it as a file through the unified data interface.

For the calculation of a single model, the above four modules are connected in a series. When the calculation of a single model is completed, the entire memory storing the model variables and temporary data is released, and the information reading submodule is returned to the next model until all calculations are completed. The framework of the structural calculation module is shown in Figure 7.

In the modeling submodule, both of the simplified models apply the uniaxial hysteretic material to simulate the inter-story resilience characteristics of the structure and the two-node link element to connect the concentrated mass nodes of each story. Uniaxial hysteretic materials are often used to construct constitutive models of uniaxial trilinear hysteresis materials. The inter-story force-displacement relationships and hysteresis parameters of the hysteretic material are defined in Figure 8 [31]. The two-node link element is a widely used joint element, which can directly use the uniaxial stress-strain material relationship model as the force-displacement relationship model of the structure [32]. However, this element has many features that are difficult to grasp. Two points require special attention: (1) for a non-zero length element, the local x-axis points are from node i to node j, and (2) the local degrees of freedom 2 and 3 (shear and flexural) of the two-node link elements are coupled, and the shear distance controls the stiffness matrix. The shear distance is defined as the ratio of the distance between the shear center point and the node i to the unit length. In the default state, the shear distance is defined as the middle of the element, as shown in Figure 9. For the meanings of the variables in Figures 8 and 9, refer to [31, 32].

3.3. Data Postprocessing Module

The primary function of the data postprocessing module is to process the enormous seismic response data generated by the structural calculation module and determine the seismic damage state of the structure according to the calculated damage limit. The module consists of three submodules: response processing, damage limit calculation, and seismic damage result. The framework of the module is shown in Figure 10.

The response processing submodule first reads the seismic response of each building and then calculates the response envelope of each story under all seismic conditions. The response parameter required for damage discrimination is divided into two types: (1) interstory displacement angle and (2) interstory slope . The recorder command in OpenSees can directly record the interstory displacement angle, and the interstory slope needs to be calculated from the displacement parameters of each story, as shown in formulas (4) and (5):where is the displacement, is the height, and is the slope of story (when , ).

The damage limit calculation submodule needs to call different functions according to structural types and calculation models when determining the damage limit of each story. Among them, the NMS RC and NMS M functions are used to calculate the damage limit of the RC frame structure and multistory masonry structure in the NMS model, respectively. The NMFS RC function is used to calculate the damage limit of the RC frame shear/cylinder structure in the NMFS model. Two main methods are used to determine the damage limit through force-based and deformation-based failure criteria [33]. One is to directly read the prewritten file and look up the table according to structure types to obtain the limit of displacement. The second is to call the key points of the building capacity curve calculated by the data preprocessing module from memory. The data flow of the damage limit calculation submodule is shown in Figure 11.

The seismic damage submodule determines the damage state of each story according to the response envelope and damage limit and calculates the regional seismic hazard. This study applies the seismic hazard matrix proposed by Yin to characterize the seismic performance level of an area [34]. The seismic damage matrix describes the probability of a certain level of damage occurring from an earthquake, which can be expressed aswhere indicates the probability of level damage when an earthquake of intensity is experienced by a building in a region and indicates the probability of an earthquake of intensity occurring over a while.

4. Case Study

4.1. Overview of the Test Area

To verify the reliability of the system, the authors of this study selected the campus of Tongji University, Shanghai, which has a complete computer-aided design (CAD) planning map and GIS data for the sections of the test area, to conduct seismic hazard prediction. Missing or insufficient data of building materials were obtained through field surveys and by reviewing building design drawings in the university archives. The collected building data were stored and managed using ArcGIS 10.4. Through GIS data and field surveys, information on a total of 421 buildings was obtained. There is a variety of building types in the test area, most of which are typical multistory masonry houses in urban areas, RC frame office buildings, and RC frame shear/cylinder high-rise buildings, as shown in Figure 12.

The test area has experienced several developments and transformations throughout its history. The construction times of different buildings are diverse, and the distribution of structure types varies from generation to generation. The system developed in this study considers the influence of when the construction took place. The time of construction determines the design specifications adopted by the structure, which affects the value of the bearing capacity and the damage limit of the model. The distribution of construction time for different structural types in the test area is shown in Figure 13. This figure shows many of RC structures built since the 1980s, while the proportion of masonry structures gradually decreased.

4.2. Selection of Ground Motion

Collecting ground motion data is a prerequisite for regional seismic hazard prediction. Due to the lack of detailed geological data on the test area, the ground motion parameter plot of Shanghai was used as the basis of ground motion data, as shown in Figure 14 [35]. The response spectrum applied for the test area is shown in Figure 15. For the meanings of variables, refer to [36]. According to the contour map of the surface peak acceleration provided by the plot and the regulations for seismic design of Shanghai buildings (J10284-2013) [36], the highest surface peak acceleration in the test area with a 50-year probability of surpassing 63%, 10%, and 3% was scaled to peak ground-level intensities of 35 gal, 100 gal, and 200 gal; the site feature cycle takes 0.9 s, 0.9 s, and 1.1 s; and the maximum earthquake influence coefficient takes values of 0.09, 0.23, and 0.45.

The predominant period of the input ground motion should align as closely as possible with the characteristic period of the test area, and the response spectrum should also align well with the design response spectrum within a specified period. According to the design response spectrum and considering the site conditions in Shanghai, 28 ground motions (14 for frequency and moderate earthquakes and 14 for rare earthquakes) were selected from the strong earthquake database of the Pacific Earthquake Engineering Research Center (PEER) [37] and the KiK-net Digital Strong Earthquake Recording System of the National Research Institute for Earth Science and Disaster Resilience (NIED) [38]. A comparison of the average response spectrum for the 28 ground motions and the recorded response spectra for this project is shown in Figure 16. The average response spectrum of the selected ground motions is in good agreement within the period 0 to 6 s. Since the natural period of most buildings in the test area was within 6 s, the 28 ground motions were selected to conduct the elastoplastic time-history analysis to obtain the seismic response and damage state of each building.

4.3. Prediction Results

By using the system developed in this study, 421 buildings in the test area were investigated based on seismic hazard prediction under frequent earthquakes (intensity VI 0.05 g), moderate earthquakes (intensity VII 0.1 g), and rare earthquakes (intensity VIII, 0.2 g). A total of 17,682 calculations with a laptop (i7 7700 HQ, 32 g memory) took about 58 minutes. The seismic hazard prediction results of the test area under different earthquake levels were visualized on the ArcScene module in ArcGIS 10.4, as shown in Figure 17.

The damage states of different structural types are shown in Figure 18. From the prediction results, the seismic capacity of various structural types in the test area was determined to be uneven due to construction time, design level, and structural problems. The seismic capacity of the RC frame and RC frame-shear/cylinder structure with a seismic design is the best. Most buildings were intact under moderate earthquake prediction evaluation and moderately damaged under rare earthquake prediction evaluation. Because the construction time of multistory masonry buildings represents the oldest buildings and their seismic capacity is slightly worse than that of RC structures, about half of the masonry buildings showed moderate damage under rare earthquakes. However, in general, the damage assessment of buildings in the test area met the target levels of the three-level seismic design.

4.4. Discussion

In order to verify the reliability of the prediction results in this system, the results were compared with results based on the vulnerability matrix [7]. Comparison data were selected from the disaster prevention planning report on Weifang Street in the Pudong New Area of Shanghai. Ouyang et al. [39] conducted a sample survey of the target area and applied a semi-experiential and semi-theoretical method to obtain the seismic damage matrices of different structural types in Weifang Street under different construction years. Due to the sample of buildings, this section only compares the results of a reinforced concrete frame and multistory masonry structures, as shown in Figures 19 and 20. In figures, THA means the results based on time-history analysis, and VM means the results based on a vulnerability matrix.

To facilitate the comparative analysis, the seismic damage index was used to indicate the seismic performance of buildings of different ages assessed under the vulnerability matrix used in [39] and the time-history analysis method used in this paper. The seismic damage index is a dimensionless index ranging from 0 to 1. The smaller the index value, the lighter the damage [40]. The average seismic damage index of an area can be calculated aswhere indicates the regional average seismic damage index under certain intensity, indicates the seismic damage index corresponding to level damage under certain seismic intensity, and indicates the probability of level damage level under certain seismic intensity. The division of damage states and corresponding seismic damage indices are shown in Table 1.

A comparison of seismic damage indices of the vulnerability matrix and time-history analysis is shown in Figure 21. The seismic damage index of buildings built in the same era increases with increased seismic intensity; the older the construction time, the lower the earthquake damage index under the same seismic intensity. Due to differences in building properties in the two areas (Weifang Street in the Pudong New Area and the Tongji campus), a certain degree of dispersion was found in the results, but the trend of the seismic damage index proved consistent; moreover, in newer buildings, a smaller difference between the two methods verified the accuracy of the results to some extent.

Notably, seismic hazard prediction based on a vulnerability matrix developed in [39] used seismic intensity as an evaluation index, which is too rough. The system developed in this study is based on a nonlinear time-history analysis, which can fully consider the time-domain and frequency-domain characteristics of ground motion. Besides, although construction age is taken into account to improve the accuracy in the vulnerability-matrix-based approach, the results are based on statistical methods that cannot reflect the damage state of a specific building. Therefore, there is a tendency to overestimate the damage of older buildings, especially under high seismic intensity. The system applies multiple degrees of freedom (MDOF), finite element models to simulate buildings, which can obtain the damage of each story and predict seismic hazard more accurately. Moreover, the seismic hazard results are more refined, which can provide a data foundation for future economic loss analysis and seismic resilience evaluation. In summary, the method used in this paper is better able to take into account the building specificities of the target area than the vulnerability matrix method. For predicting seismic damage in areas where there is a lack of historical seismic data, the proposed method is preferred.

5. Conclusions

In this paper, two simplified models using a time-history analysis method were applied. Based on MATLAB and OpenSees, a refined seismic hazard prediction system for urban buildings was developed. The campus of Tongji University, Shanghai, was chosen as the test area to demonstrate the application effects of the system. Conclusions concerning the results of this study are as follows:(1)The system adopts a unified data interface, which can be seamlessly connected between the modules based on MATLAB and OpenSees and realizes the automatic calculation of seismic hazard prediction.(2)The system applies different simplified models for different structures. Computation work can be done only by the macro parameters of buildings in GIS data. While this ensures calculation accuracy, calculation efficiency is improved, and the manual operation and hardware requirements of the prediction are reduced.(3)The prediction results of this time-history analysis system were compared with the results of the vulnerability method, and the trend of seismic hazard predictions was consistent in the results of both, which verified the reliability of the system. At the same time, the system can obtain the specific seismic response and damage state of each story of a building, which achieves the refined prediction necessary to assess the results of seismic hazards accurately.

With the continuous expansion of the city scale and higher requirements for disaster risk response, there is a broad application prospect and research demand for seismic hazard prediction in urban areas based on the time-history analysis. The following issues should be studied in greater depth in the future:(1)The calculation accuracy of the method applied in the paper is theoretically higher than the traditional methods, but the applicability of the method and the rationality of the parameters are to be tested. Therefore, it is necessary to compare the prediction results with the actual seismic damage in the future.(2)Although the simplified models used in this paper can greatly reduce the computation time required for modeling and time-history analysis, it is still necessary to adopt methods such as parallel computing and machine learning to reduce the computation time to meet the timeliness of seismic hazard prediction for urban-scale in the future.(3)Due to the lack of time and detailed building data, the vulnerability matrix of the test area was not obtained, so the prediction results were compared to the results of another area. In the future, more study cases in the same or similar areas are needed to verify the system’s reliability further.

Data Availability

The data used to support the findings of this study are included in the article.

Conflicts of Interest

The authors declare no conflicts of interest.


This work was supported by the National Key Research and Development Program of China (nos. 2016YFC0800209 and 2017YFC0803300) and the National Natural Science Foundation of China (no. 51178351).