A Structure Economic Loss Optimization Method with the Uncertainty of Ground Motion Amplitude for Chinese Masonry Building

Abstract

In the catastrophe insurance industry, it is impractical for a catastrophe model to simulate millions of sites’ environments in a short time. Hence, the attenuation relation is often adopted to simulate the ground motion on account of calculation speed, and both ground motion expectations and uncertainties must be calculated. Due to the vulnerability curves of our model being based on simulations with a large number of deterministic ground motions, it is necessary but not efficient for loss assessment to analyze all possible ground motion amplitudes and their corresponding loss rates. This paper develops a simplified method to rapidly simulate loss expectations and uncertainties. In this research, Chinese masonry buildings are the focus. The result shows that the modified method gives accurate loss results quickly.

1. Introduction

Since the beginning of the 21st century, the frequency of natural disasters has been on the rise, and with the development of urbanization and industrialization, economic losses and casualties caused by natural disasters have become more and more serious. According to the report issued by the United Nations International Disaster Reduction Agency [1], from 1998 to 2017, climate and geological disasters caused 1.3 million deaths, more than 4.4 billion injuries, and even homelessness, resulting in direct economic losses of up to USD 298 billion, which is 251% higher than those caused by natural disasters in the past 20 years (1978–1997). The seismic intensity of the Wenchuan Earthquake in 2008 reached XI degree, and the building collapse caused serious casualties and infrastructure damage that could not provide normal service functions, which led to the shutdown of enterprises and huge economic losses. From the perspective of industry, the impact of the Great East Japan Earthquake in 2011 was more significant. After the earthquake, the production capacity of the affected areas decreased by 30~50% directly, the recovery time of production capacity was slow, and even 4 years after the earthquake, the industrial production capacity of some affected areas still had not recovered to the pre-earthquake level [2]. The risk management of natural disasters has been paid more and more attention to, and how to scientifically estimate the losses caused by natural disasters has become one of the urgent issues to be solved.

Loss assessment in catastrophe modelling has some differences from that in theoretical research. The loss needs to be analyzed and calculated by using the earthquake risk module and earthquake vulnerability module in combination with the object value [3]. As insurance and reinsurance claims are based on events, for a single object, the output of the earthquake risk module is the estimated result of the ground motion parameters of each simulated event. The earthquake vulnerability is the damage rate of corresponding ground motion parameters, and these curves composed of the economic loss ratio and ground motion parameters are called vulnerability curves in catastrophe modelling [4]. These are different from the structural vulnerability curves in earthquake engineering but are related to structural vulnerability. In seismic engineering, the structural vulnerability curve is usually the statistical analysis of the correlation between the seismic damage investigation data of the structure under a large number of ground motions or the performance data of the structure after elastic–plastic analysis and the ground motion parameters [5]. In catastrophe modelling, the structural vulnerability in earthquake engineering is often referred to as fragility. Vulnerability in the insurance industry represents the weighted sum of the damage ratio and the probability of all different damage levels of a single structure under different earthquake incentives [3].

There are many uncertainties in earthquake risk management, which need to be considered as much as possible when estimating the loss model. Abrahamson [6] and FEMA-P58 [7] indicated that uncertainty may exist in the characteristics of earthquake hazards as well as in the properties of geometric shape and material. Cutfield and Ma [8] studied three fragility curves under different cases, which included two curves without uncertainty and one curve with demand uncertainty, and recommended considering the demand uncertainty when dealing with the post-earthquake data. In addition, De Risi et al. [9] indicated that, when processing fragility curves, the uncertainty of input data should be adequately considered; otherwise, the uncertainty of the output result may be amplified. International scholars have carried out some research on the uncertainty of earthquake ground motion and structural response. In terms of the influence of uncertainties of different factors, Kwon and Elnashai [10], Soleimani [11], and Pan et al. [12] have different thoughts. Kwon and Elnashai [10] indicated that the influence of the randomness of materials on the fragility curve is much less than that of the randomness of ground motion properties for reinforced concrete buildings. Soleimani [11] stated that, compared to the uncertainty of ground motion and materials, the geometric characteristics were the most sensitive for the high box beam concrete bridge. In contrast, Pan et al. [12] found that, in the elastic model, the influence of the uncertainty of the structural geometry on the structural earthquake response is much less than the uncertainty of the material properties, such as damping ratio and yield strength for the transmission tower. The reason for the aforementioned result may be caused by their different structure types. Hence, this study pays more attention to the uncertainty of ground motion parameters. Yamin et al. [13] carried out non-linear time-history analysis based on a 3D structure model and used Monte Carlo simulation to simulate the time and cost of building restoration, ground motion properties, structural response, and building destruction, and finally obtained the fragility of reinforced concrete frame structures. Yazgan [14] provided an empirical method to build a fragility model of the reinforced concrete structure that considers the uncertainty caused by the lack of peak ground motion data. Ansari et al. [15] studied the impact of base types on the fragility of high-rise reinforced concrete buildings, and the seismic uncertainty was considered by adopting a certain number of near-field and far-field earthquakes. Choudhury and Kaushik [16] studied the fragility curve of the steel–concrete frame structure with and without brickwork infill wall, and the open steel–concrete frame structure on the ground floor by analyzing three uncertainties, including the performance difference of non-linear materials, the structural geometric characteristics, and the load situation. Kim et al. [17] analyzed the impact of construction quality defects on the seismic fragility of reinforced concrete frame structures and stated that the seismic fragility of the selected reinforced concrete frame is sensitive to the strength of concrete and the volumetric ratio of transverse reinforcement. Saloustros et al. [18], using the finite element method and the Monte Carlo method, investigated the influence of uncertainty on materials of masonry structure components and present the fragility curve of Santa Maria del Mar church. However, European historical masonry buildings, such as churches and cathedrals, are different from Chinese masonry buildings, which are built with clay bricks instead of stone. In China, many masonry buildings in big cities, such as Beijing and Shanghai, were built from the 1970s to the 1990s to solve the residential problem during urbanization. However, the seismic capacity of a multistory masonry building is poorer than the reinforcement concrete (RC) structure, which means the loss of a masonry building will be more than for RC under a damaging earthquake. Hence, as a more vulnerable structure, the Chinese masonry buildings are selected as calculation examples. Jiang et al. [19] studied four unreinforced masonry structures with different floors, established the finite element model by OpenSees, and discussed the influence of uncertainty of ground motion and structural parameters on fragility by the IDA method and the one-variable quadratic matrix method. It was pointed out that, with the increase of failure degree, the greater the influence of uncertainty on fragility, and, compared with the seismic uncertainty, the structural parameter uncertainty has a more obvious influence on a structure with a lower number of stories, and the structural damping ratio is the most sensitive among all the influence factors.

Although the aforementioned research considered the uncertainties of ground motion and structural response, the method of amplitude adjustment is generally selected to process different known seismic excitations. The amplitude, spectrum characteristics, and duration of each set of seismic excitations are all bound together and fixed. The amplitude adjustment of ground motion means that only the uncertainty of spectrum characteristics and duration is considered in the fragility and vulnerability curve. Potential future earthquakes will not exactly match the fixed combination of amplitude by ground motion prediction equations (GMPEs), spectrum characteristics, and duration of ground motions for historical earthquakes. Most seismic risk assessments are qualification models presented by earthquake intensity parameters and damage probability, and each step of the risk model has uncertainties that should be considered [20]. In addition, Hwang et al. [21] indicated that all the uncertainties need to be considered in the ideal structure model with possible seismic ground motion when seismic vulnerability assessment was processed. Hence, if GMPEs are used to describe the ground motion intensity and estimate the loss, the uncertainties of the calculated amplitude should be considered. For the insurance industry, simplifying the model and saving calculation time is very important, and it is necessary to use a parametric empirical statistical model, such as GMPEs, to get a loss result as soon as possible. All possible ground motion amplitude should be theoretically simulated in this method. However, it is inefficient to take every sample for each site to calculate the damage ratio. Therefore, the uncertainty of ground motion and vulnerability is coupled to generate a modified vulnerability curve, and the results are expressed by the combined mean and variance. The negative effect of this method is that the combined distribution is different from the sampling distribution, but the mean and variance are consistent, and the calculation speed is greatly improved.

The current study aims to overcome the aforementioned issues, such as the lack of consideration of uncertainties and computational efficiency. In this research, PGA (Peak Ground-motion Acceleration) amplitude distribution is applied to create a modified vulnerability curve to simplify loss ratio expectation and its uncertainty to make rapid calculations. The method proposed in this study can not only calculate massive objects’ economic loss and their vulnerability curves in a short time, but also deepen the understanding of how the interaction of uncertainties affects the damage results of objects. It provides a novel idea and a risk-based with the consideration of uncertainties solutions for practical application in the catastrophe insurance industry.

2. Methodology

The attenuation relation model is a parameterized empirical statistical model, which is built by a large amount of data of different magnitudes and distances in a certain area. It has a general rule of seismic response and statistical significance. Specifically, in the insurance industry, the business contract usually calls for loss reports as soon as possible during the renewal season. If using numerical simulation to calculate the ground motion, it usually consumes lots of time. In contrast, attenuation is an efficient and time-saving way to acquire ground motion information. Therefore, this study uses the attenuation relation model to calculate the ground motion parameters.

Due to the seismic source, wave propagation, and site effect, the ground motion produced by the same earthquake in different sites or caused by different earthquakes in the same site will be accompanied by the uncertainty of ground motion. In this study, the samples of PGA distribution are adopted instead of the mean value to analyze the influence of loss assessment results. Abrahamson [22,23] put forward that the influence of the uncertainty of the ground motion acceleration spectrum conforms to a lognormal distribution. In addition, Rossetto and Elnashai [24], Lagomarsino and Giovinazzi [25], Bradley and Dhakal [26], Jayaram and Baker [27], and Lallemant et al. [28] stated that a lognormal cumulative distribution function (CDF) is a frequently used method to model seismic structure fragility curves, and it has a significant representation in earthquake risk analysis. Furthermore, according to the model from the Seismic ground motion parameters zonation map of China (GB18306-2015) [29], as shown in Formula (1), the ground motion parameters are also in line with a lognormal distribution:

$$lgY=A+BM+Clg\left(R+D{e}^{EM}\right)$$

(1)

where $Y$ is the ground motion parameter, such as peak ground motion acceleration (PGA) or peak ground motion velocity (PGV), $M$ is the magnitude, $R$ is the epicentral distance, and $A$, $B$, $C$, $D$, and $E$ are regression coefficients.

2.1. Data

The premise of obtaining the vulnerability curve is to have the fragility curve of an object first. The steps of obtaining the fragility curve by the IDA (Incremental Dynamic Analysis) method are usually to build the structural model and calculate the dynamic equation, then use the corresponding calculation and analysis methods and input ground motion parameters; next, calculate the structural response, and finally, according to the relationship between the seismic response and the structural limit under different strength levels, determine the probability of different failure states of the structure under different ground motion parameters [30]. Gattulli et al. [31] indicated that masonry structure is extremely variable in mechanical characteristics due to its features, such as anisotropy, inhomogeneity, and nonlinearity. Hence, the uncertainty of seismic ground motion intensity may cause a greater variation in the loss result for masonry buildings.

The dataset is obtained from the fragility curves of 3–6 floor masonry buildings with different occupancies [32]. In this paper, based on the above dataset, the fragility curve is established by the IDA method with the interlayer elongation coefficient as the index of failure [32]. According to the five elements of floors, design intensity, design code, located area, and building occupancy, a total of 84 fragility curves are formed, and the detailed classification is shown in

Table 1. The classification of the classic Chinese masonry building.

Building Occupancy	Residential						Commercial
Design Code	1989-Code			2001-Code			1989-Code		2001-Code
Floors	1–3	4	4–6	1–3	4	4–6	1–2	3–5	1–2	3–5
Design Intensity	6, 7, 8	9	6, 7, 8	6, 7, 8	9	6, 7, 8	6, 7, 8, 9	6, 7, 8	6, 7, 8, 9	6, 7, 8
Outside Wall Thickness (mm)	240, 370, 490			240, 370, 490			240, 370, 490		240, 370, 490

. Based on the 84 fragility curves, the vulnerability curves are established, then the uncertainty of seismic ground motion intensity is considered, and the calculation method of loss with uncertainty is established to obtain the modified vulnerability curve. The structures and their curves, listed in

Table 2. The information on structure and curve examples.

Structure and Curve	Building Occupancy	Building Age	Floors	Design Intensity	Outside Wall Thickness
0	Residence	89-code	1–3	6	240
8	Residence	01-code	1–3	7	240
9	Residence	01-code	1–3	8	240
10	Residence	01-code	4–6	6	240
16	Residence	89-code	1–3	8	370
19	Residence	89-code	4–6	8	370
24	Residence	01-code	4–6	6	370
26	Residence	01-code	4–6	8	370
27	Residence	89-code	1–3	6	240
51	Commerce	01-code	1–3	8	240
52	Commerce	01-code	5	6	240
54	Commerce	01-code	5	8	240
66	Commerce	01-code	5	6	370
83	Commerce	01-code	4	9	490

, will be shown as examples in this study. In addition, the Chinese seismic design code has a strict limitation on the height and floors of masonry buildings. The building located zone is applied to distinguish the outside wall thickness, and 240 mm, 370 mm, and 490 mm stand for zone III, II, and I, respectively. Figure 1 shows the typical masonry structure in China, where the wall thickness is usually 240 mm and 370 mm, respectively, in Sichuan and Beijing. Figure 1a is a two-floor masonry structure built following 01-code with 7-degree seismic design, which can represent structure 8 in Table 2. The building in Figure 1b represents structure 26 in this study, which is built following the 01-code.

The calculation process is shown in Figure 2. First, each PGA in the original vulnerability curve data is sampled to obtain PGA samples. Next, the PGA samples are transformed into damage ratio samples by the interpolation method. Next, the damage ratio samples are sampled again to get the second damage ratio sample. Finally, the means of second damage ratio samples are calculated to obtain the final damage ratio corresponding to each original PGA.

2.2. Sampling of Ground Motion Parameters

The optimization of ground motion uncertainty originates from quantifying the regular and repeatable parts of random uncertainty and quantifying them in the final risk results. The uncertainty of ground motion is generally reflected by the standard deviation in the attenuation relationship of ground motion. If the standard deviation obtained via ergodic assumption is divided carefully and the identifiable part is quantified, the uncertainty can be reduced, and the calculated result is closer to the real value [33]. In this study, the standard deviation is acquired by the attenuation relation model. However, because the coefficients of the attenuation relation model are diverse, they depend on the geological information of the different areas. Hence, for demonstration, a basic attenuation relation model is adopted in this study. According to the PGA attenuation coefficient and attenuation relationship model (Formula (1)) given in the publicity and implementation textbook of Seismic ground motion parameters zonation map of China [34], the original PGA is sampled. The ground motion parameter PGA conforms to the lognormal distribution, and its standard deviation is a constant that equals 0.236. There are 197 PGA points in the original data, which are arranged from 0.4 m/s² to 20.0 m/s² in steps of 0.1 m/s². Because the Monte Carlo method consumes much computational power, and this algorithm needs to sample twice, the consumption of computational power will increase exponentially. Latin hypercube sampling is adopted to keep the function monotonically increasing, as shown in Figure 3. It is a masonry structure of a residential building with three stories and six degrees of design intensity in site class III with 89-code. Ghotbi and Taciroglu [35] indicated that the Latin hypercube method has a better effect than the Monte Carlo method under the same number of samplings via studying the damage analysis of a normal 4-story non-ductile reinforced concrete structure when selecting seismic ground motion sets.

In addition, the bias of different Latin hypercube sampling times was tested, and the results are shown in

Table 3. The bias of different sampling times.

Sampling Times	Bias
100	0.5%
200	0.07~0.08%
500	0.02~0.03%

. Compared with 100 sampling times, the bias of 500 sampling times is much smaller.

Therefore, the Latin hypercube sampling method is adopted in this algorithm, which can not only reduce the computational power consumption but also satisfy the empirical judgment. In this study, each PGA point is sampled 20 times with truncation in the corresponding lognormal distribution, and the specific truncation method will be introduced in the following. The 3940 PGA samples obtained from the sampling are named ground motion sample 1. It is substituted into the linear relationship between the original PGA and the original damage ratio, and the corresponding 3940 damage ratios are obtained by using the interpolation method, which is named damage ratio sample 1. In the damage ratio sample 1, if the corresponding PGA is larger than 20.0 m/s² (maximum value of the original PGA), the damage ratio takes the value corresponding to 20.0 m/s².

Currently, there are 3940 damage ratio samples, but only 197 corresponding coefficients of variation. Therefore, before sampling the damage ratio, the coefficient of variation should be expanded to 3940. Assuming that the corresponding coefficient of variation follows the linear relationship composed of the initial 197 damage ratios and the coefficient of variation, the corresponding 3940 coefficients of variation can be obtained via the linear interpolation method, namely, the coefficient of variation sample 1. In the process of linear interpolation, the sample value of the damage ratio is less than the original value. In this paper, the original minimum damage ratio is connected with the zero point to interpolate the sample value less than the minimum original minimum damage ratio. It is worth noting that the known condition is the coefficient of variation, so the corresponding standard deviation needs to be calculated via Formula (2), and then linear interpolation must be performed. The results of interpolating $Cv$ first and then calculating $\sigma$ are different from calculating $\sigma$ and then interpolating.

$$\sigma =Cv\times \mu$$

(2)

where $\sigma$ is the standard deviation, $Cv$ is the coefficient of variation, and $\mu$ is the mean value.

According to the attenuation relationship given by the Seismic ground motion parameters zonation map of China (Standardization Administration 2015) [30], the ground motion conforms to the lognormal distribution with a standard deviation equal to 0.236. When the translation between seismic intensity and ground motion parameters is processed, the seismic intensity usually corresponds to a PGA value. In fact, a seismic intensity should correspond to a range of PGA with distribution instead of its mean value. The sample range needs to be truncated because the ground motions cannot tend to infinity or 0. In this study, hence, the correlation between PGA and seismic intensity of site class II in GB 18306-2015 Seismic ground motion parameters zonation map of China [30] was used to determine the truncated range within the corresponding seismic intensity range.

Table 4. The comparison table of PGA and seismic intensity for site class II.

PGA in the Site Class II (m/s2)	0 ≤ αmaxII < 0.392	0.392 ≤ αmaxII < 0.882	0.882 ≤ αmaxII < 1.862	1.862 ≤ αmaxII < 3.724	3.724 ≤ αmaxII < 7.35	7.35 ≤ αmaxII < 14.7	αmaxII ≥ 14.7
Seismic intensity	V and below	VI	VII	VIII	IX	X	XI and above

shows the details of truncation.

In this algorithm, for each ground motion parameter P, the algorithm finds seismic intensity level IP using Table 4 and uses [IP − 1, IP + 1] as the initial range. The left endpoint of the new range is the left endpoint of the minimum value of the three seismic intensities, which is recorded as the temporary left endpoint. The right endpoint of the new range is the left endpoint of the maximum value of the three seismic intensities, which is recorded as the temporary right endpoint. For example, the initial range of 6-level seismic intensity is the range of the entirety of the 5-level, 6-level, and 7-level. Next, according to the original ground motion parameter data, its quantile is obtained in its corresponding seismic intensity range (an intensity range). The final truncated range of the original ground motion parameter data is calculated by using the temporary right endpoint, temporary left endpoint, and quantile. The calculation process is shown in Formulas (3)–(5), which are equivalent to the range of two seismic intensities.

$$Q_i=(PG{A}_i-PG{A}_{i,j})/(PG{A}_{i,j+1}-PG{A}_{i,j})$$

(3)

$$Lowe{r}_i=PG{A}_{i,j-1}+Q_i(PG{A}_{i,j}-PG{A}_{i,j-1})$$

(4)

$$Uppe{r}_i=PG{A}_{i,j+1}+Q_i(PG{A}_{i,j+2}-PG{A}_{i,j+1})$$

(5)

where $Q_i$ is the quantile of $PG{A}_i$ in its corresponding seismic intensity range, $PG{A}_i$ is the ith element in the original PGA,$PG{A}_{i,j}$ is $PG{A}_i$ corresponding to the left end of the seismic intensity range, $PG{A}_{i,j+1}$ is $PG{A}_i$ corresponding to the right end of the seismic intensity range, $Lowe{r}_i$ is the lower limit of the final truncation range, $Uppe{r}_i$ is the upper limit of the final truncation range, $PG{A}_{i,j-1}$ is the temporary left endpoint, and $PG{A}_{i,j+2}$ is the temporary right endpoint.

When $PG{A}_i$ is in the range corresponding to seismic intensity V and below, or XI and above, the truncation range will be different from Formulas (3)–(5). When $PG{A}_i$ is in the range of seismic intensity V and below, due to it is the minimum seismic intensity range, it is unable to obtain a smaller level, so its truncation lower limit is fixed equal to 0, and the truncation upper limit still follows the algorithm of Formula (5). When $PG{A}_i$ is in the range of seismic intensity XI and above because it is the largest seismic intensity range, it is unable to obtain a higher level, so its upper truncation limit is fixed at 20 (the maximum PGA of the calculation example of this algorithm is 20 m/s²), and the lower truncation limit is still calculated according to Formula (4).

2.3. Damage Ratio Sampling

Most of the previous studies used lognormal distribution to describe the distribution characteristics of building damage ratio [22,23,24,25,26,27,28,36], and this would lead to more than one sample, but this is not logical, because the damage rate cannot exceed 100% of the replacement value. Some scholars pointed out that if the cost of cleaning up the ruins is superimposed, it will exceed the replacement cost. However, if this factor is not considered in the statistical process, then it is not logical to consider it at this time. If truncated sampling is used, the distribution characteristics will lead to a change in mean and standard deviation. Therefore, the β distribution is used to describe the distribution of structural loss. The mean value and standard deviation of the distribution are kept unchanged and, as known conditions, the coefficients α and β of β distribution are calculated by substituting Formulas (6) and (7).

$$\alpha =\frac{{\mu }^2{\sigma }^4-\mu {\sigma }^4+{\mu }^2{\left(1-\mu \right)}^2}{\left(1-\mu \right){\sigma }^4}$$

(6)

$$\beta =\frac{\mu {\sigma }^4-{\sigma }^4+\mu {\left(1-\mu \right)}^2}{{\sigma }^2}$$

(7)

where μ is the mean value of the structural damage ratio and σ is the standard deviation of the structural damage ratio.

If the α and β corresponding to 3940 damage ratios are known, the second Latin hypercube sampling can be carried out, and the sampling frequency is the same as 20. After sampling, the sample was expanded to 78,800, equivalent to 400 damage ratio samples derived from each original PGA point, which results in damage ratio sample 2.

The 197 final damage ratios are obtained via calculating the mean of corresponding damage ratio sample 2, which is, namely, damage ratio sample 3 with the uncertainty of ground motion intensity. Finally, according to damage ratio sample 2, the coefficient of variation of damage ratio sample 3 considering the coupling effects of various uncertainties is calculated, namely, the coefficient of variation sample 3.

3. Result and Discussion

Vulnerability Curve Analysis

To more intuitively understand the twice sampling process and the distribution of PGA and damage ratio samples, the vulnerability curves obtained from each group of Latin hypercube samples are displayed in three-dimensional (3D) form. Figure 4 shows the 3D vulnerability diagrams of structure 0 and structure 54. Structure 0 is a 3-story masonry structure of a residential building with 6-degree design intensity and site zone III according to 89-code, and structure 54 is a 6-story masonry structure of an office building with 8-degree design intensity and site zone III according to 01-code. The 89-code and 01-code mean the Code for seismic design of buildings in China, which are published in 1989 (GB50011-1989) and 2001 (GB50011-2001), respectively, and the key difference between the above standards is shown in

Table 5. The design parameters of 89-code and 01-code.

Standard	89-Code	01-Code
Brick Intensity	MU7.5 (compressive strength ≥ 7.5 MPa)	MU10.0 (compressive strength ≥ 10.0 MPa)
Concrete Strength	C15 (compressive strength ≥ 15.0 MPa)	C20 (compressive strength ≥ 20.0 MPa)
Roof Dead Load (including self-weight) (kN/m2)	6.0	6.0
Roof Live Load (kN/m2)	1.5	2.0
Floor Dead Load (including dead weight) (kN/m2)	3.3	3.3
Floor Live Load (kN/m2)	1.5	2.0

. In this study, the wall thickness of the masonry structure is 240 mm in site zone III, 370 mm in site zone II, and 490 mm in site zone I. The x-axis of the 3D vulnerability diagram is the original PGA in m/s², the y-axis is ground motion sample 1 from the first sampling in m/s², and the z-axis is the damage ratio. Each original PGA on the x-axis has 20 PGA samples on the y-axis, and 3940 points on the X-Y plane have 20 damage ratio samples on the z-axis. In the 3D coordinate system, there are 78,800 (197 × 20 × 20) points, representing 197 original PGA points sampled separately 20 times, and then the corresponding damage ratio (damage ratio sample 1) of 3940 PGA samples (ground motion sample 1) was sampled, respectively, with 20 sampling frequency, and 78,800 damage ratios obtained from the final sampling are distributed in three-dimensional space.

It can be seen in the 3D diagram that the vulnerability curve becomes a vulnerability surface, which can help to understand the sampling process. Its shape becomes wider along the boundary of the range and finally ends at the end of the plane composed of PGA and PGA samples, namely, points (20, 20). The scatter color in Figure 4 gradually changes from blue to yellow, which means that the value of the damage ratio gradually increases from small to large. The surface color changes from light to dark, which means that the value of the damage ratio changes from low to high. As Figure 4a,c shows, the number of yellow scattered points is much denser than the number in Figure 4b,d. Especially when the value of PGA is small, the number of light color points in Figure 4b,d is much more than curve 27, and, due to this, in the two-dimensional (2D) vulnerability curve, as shown in Figure 5, curve 52 is much steeper than curve 27. It means that a higher damage ratio has been reached when the seismic excitation has not reached a certain intensity. In addition, the area of the surface is negatively related to the loss of the structure. The larger the area of the surface is, the more vulnerable the structure is. It is worth noting that, as shown in Figure 4 and Figure 5, the edge of the surface formed by the damage ratio points and damage ratio curve is not completely smooth, which is caused by the change in the range of truncated sampling. Namely, the algorithm of Formulas (3) and (5) lead to this result.

In this paper, some typical comparison results are selected. As shown in Figure 4, the red line is the vulnerability curve without considering the uncertainty, and the green line is the vulnerability curve considering the uncertainty. It should be noted that the original vulnerability curve starts at the point that PGA and the damage ratio are equal to 0.4 m/s² and 3%, respectively. However, this study makes the original minimum damage ratio extend to 0 (i.e., the point of PGA = 0, damage ratio = 0), so the modified curve starts from (0, 0). From vulnerability curve 0 to 83, a total of 84 vulnerability curves describe different masonry structures. It can be seen that the vulnerability curve of masonry structures with the uncertainty of ground motion intensity is lower than the original vulnerability curve in most PGA ranges. Only when the PGA is small is the damage ratio with the uncertainty of ground motion intensity higher or nearly overlaps with the loss ratio without the uncertainty. Ioannou et al. [37] used the building damage database of the 1980 Irpinia earthquake as the data source via the Bayesian framework to explore whether the uncertainty of ground motion intensity can significantly change the shape of the fragility curve, and the results show that the complex model with considering more uncertainty will not change the shape of the fragility curve. It can be seen that the vulnerability curve obtained in this study also has no obvious change in shape, but the damage ratio before and after modification is different, especially when the PGA value range is in the middle, namely, in the range of moderate earthquake.

Curve 0 is the vulnerability curve of a 3-story residential masonry structure with 6-degree design intensity and site zone III according to 89-code, and curve 54 is the vulnerability curve of a 6-story office masonry structure with 8-degree design intensity and site zone III according to 01-code, as shown in Figure 6a,c. Under the same earthquake excitation, the stricter the design code of the building structure, the smaller the structure’s seismic response and the smaller the damage ratio.

By controlling the unique variable, it is better to study the influence of different factors on the results of the vulnerability curve. Curve 24 is the vulnerability curve of the 6-story residential building structure with 6 seismic degree design intensity and site zone II according to 01-code, and curve 26 is the vulnerability curve of the 6-story residential building structure with 8-degree design intensity and site zone II according to 01-code. The only difference between curve 24 and curve 26 is seismic design intensity. It can be seen from Figure 6c,d that, in the case of the same height, site class, and design code, the seismic intensity is negatively correlated with its damage ratio. By observing curves 16 and 19, as Figure 6e,f show, they are, respectively, the 3-story residential masonry structure vulnerability curves with 8-degree design intensity and site zone II according to 89-code, and the 6-story with 8-degree design intensity and site zone II according to 89-code. It is found that for masonry structures, under the same seismic intensity, design code, and site zone, the higher the height, the higher the damage rate will be. Curve 26 is the vulnerability curve of the 6-story residential masonry building with 6-degree design intensity and site zone II according to 01-code. Comparing curve 19 and curve 26, as shown in Figure 6d,f, it can be seen that the newer the design code, the lower the damage ratio under the conditions of the same building height, seismic intensity, and region. Curve 24 and curve 10 take the site zone as a single variable, and other conditions are the same, which are 6-story masonry housings with, 6-degree design intensity according to 01-code. Comparing Figure 6c,g, it can be seen that under the conditions of different site zone and the same building height, seismic intensity, and design code, the higher the number of the site zone, the greater the loss, which means that the thicker the wall, the smaller the loss will be if other conditions remain unchanged. Curve 66 is a 6-story office building with 6-degree design intensity and site zone II according to 01-code. Compared with curve 24, the difference is only functional use. Through comparison, it can be seen that under the same PGA, the damage ratio of the office building is smaller than that of the residential building, as shown in Figure 6c,h. However, it can also be seen that no matter which single variable is changed, the trend of the curve is almost unchanged, and the resulting damage ratio does not change significantly.

4. Discussion

To better show the effect of our combined uncertainty method, the case of the Jiuzhaigou earthquake in 2017 is chosen for verification. On 7 August 2013, the earthquake occurred in Jiuzhaigou county, Sichuan, China, and the magnitude reached 7.0. This study uses the earthquake damage data collected by Zhang et al. [38] in the Jiuzhaigou earthquake. There is a total of eleven investigation areas, eight of which are located in the IX-degree area, two of which are located in the VIII-degree area, and the last one is in the VII-degree area. Detailed information, such as the geographic location and seismic parameters, is listed in

Table 6. Detailed information of the investigation sites.

	Epicenter	Investigation Point
Site number	0	1	2	3	4	5	6	7	8	9	10	11
Intensity	IX	VIII	IX	IX	IX	IX	IX	VIII	VIII	VII	IX	IX
Longitude	103.82	103.78	103.80	103.83	103.85	103.92	103.88	103.92	103.93	103.98	103.92	103.90
Latitude	33.20	33.32	33.30	33.30	33.30	33.28	33.28	33.27	33.27	33.32	33.22	33.20
Distance to epicenter r(km)	0	13	11	12	13	12	14	16	17	24	14	12
PGA (m/s2)		2.9	3.5	4.7	4.6	3.7	3.9	3.3	3.5	2.5	4.6	5.7

. According to the geographic information of the investigation areas, the relevant PGA can be calculated, and the calculation method adopts attenuation relation with Formula 1. In addition, all the investigation sites are located in site class II in China, as the information in

Table 7. The relationship between equivalent shear wave velocity and overburden thickness of site class II in China.

Site Class II
Equivalent Shear Wave Velocity of Soil (m/s2)	Overburden Thickness (m)
$500 \ge v_{se}$ > 250	≥5
$250 \ge v_{se}$ > 150	3–50
$v_{se}\le$ 150	3–15

shows.

According to the investigation information of local constructions, the vulnerability curves most likely to represent the masonry in this area are selected empirically in this study. Theoretically, local building characteristics should be adopted, but it is difficult to get all detailed information on all local building types. The two constructions chosen in the Jiuzhaigou area are built after 2001 adopted 01-code, the seismic intensity is 8-degree, the wall thickness of the buildings is 240 mm, and, because of the hilly county, the constructions usually are low-rise. They represent residential buildings and commercial buildings, respectively, with 01-code, 8-degree design intensity and fewer than three floors. During the sampling, there are two groups of PGA, which are the revised group and the original group. The revised group is sampled 10,000 times and follows truncation log-normal distribution to simulate the precise ground motion of the different buildings’ location. The resulting PGA samples are inserted into the original vulnerability curve and revised vulnerability curve, respectively, to get the two damage ratio samples. Then, based on

Table 8. The relationship between building damage state and damage ratio.

Damage State	Basic Intact	Slight Damage	Moderate Damage	Heavy Damage	Ruin
Damage Ratio	[0,10)	[10,30)	[30,55)	[55,85)	[85,100)

, the damage ratio can be transformed to the damage state, which is the earthquake damage matrix.

The earthquake damage matrix collected by Zhang et al. [38], the revised matrix, and the actual statistical matrix are listed in

Table 9. The comparison of the earthquake damage matrixes.

	Basic Intact	Slight Damage	Moderate Damage	Heavy Damage	Ruin
Statistical matrix	0	0.141361	0.602094	0.204188	0.052356
Original residence matrix	0.012231	0.194815	0.350715	0.404741	0.037497
Revised residence matrix	0.006984	0.219444	0.40487	0.362677	0.006025
Original business matrix	0.029034	0.23363	0.324603	0.366655	0.046079
Revised business matrix	0.017947	0.259816	0.381534	0.333405	0.007297

, in which the numbers represent the proportion of buildings in each state of damage to the total.

In Figure 6, 01-8-R-3-Re-trunc stands for revised curve with truncation sampling, and 01-8-R-3-Ori-norm means original curve with random sampling. In Figure 7a, compared with the statistical data, the basic intact and moderate damage of the revised curve is much better than the original curve, and the rest of the damage situations are very similar. In addition, for business construction, the performance of the revised curve for basic intact and moderate damage is superior to the original curve, as shown in Figure 7b. Generally, the revised curve can provide a more conservative result than the original curve, and, without sacrificing the accuracy of the rest of the damaged state, a more reasonable performance is achieved, especially in the highest proportion of damaged state, the moderate damage.

However, in the remote mountainous area, there are many residential constructions built by local villagers, which cannot satisfy the seismic requirement. Therefore, in this research, the curve of the masonry structure with the best earthquake resistance and the masonry structure with the worst earthquake resistance are put together to study, as shown in Figure 7. The 89-6-R-3-Re-trunc (Curve 0) stands for the revised curve of residential masonry buildings with 89-code, 6-degree design intensity, 6-floor, and 240 mm wall-thickness via truncation sampling, which is the worst of all the curves. The 01-9-B-4-Re-trunc (Curve 83) means the revised curve of business masonry buildings with 01-code, 9-degree design intensity, 4-floor, and 490 mm wall-thickness via truncation sampling, which has the best seismic performance of the curves. In Figure 8, except for basic intact and moderate damage, the rest of the damage states are within the best performance curve and the worst performance curve. In terms of basic intact, the percentage of basic intact in statistical data is zero, because areas with notable damage are usually selected for earthquake damage surveys. Hence, the statistical magnitude of the intact building may have bias, which can explain why both calculated curves are above the statistical data in the damage state of basic intact. The actual results range between the best and worst curves, which depends on what is known about the local masonry buildings. Therefore, the curves are selected empirically in this paper, and the damage result obtained is acceptable, which can indirectly provide a comparatively accurate reference, if detailed building information could be gained, for the prediction of earthquake damage in the insurance industry.

5. Conclusions

It is necessary to quantify the independent uncertainty of the ground motion intensity in the process of calculating the loss expectations. In this paper, a novel method has been purposed for rapid calculation and the amplitude distribution selection was discussed. To verify the effect of the method, the damage to masonry buildings in the Jiuzhaigou earthquake was discussed.

The conclusions of this study are divided into the following three points:

• Through this method, PGA and damage ratio are no longer required to be sampled during calculation, and a fast calculation is supported via the revised curves. Although the distribution of damage ratio has been changed after combination, it still meets the needs of the insurance industry.
• Whether uncertainty is considered or not, under the same seismic excitation input, the higher the design intensity level, the lower the damage ratio. Therefore, the strict seismic design code can not only reduce the casualties of personnel in the earthquake but also reduce the seismic economic losses of the structure and ensure the residual value of the structure.
• Between the vulnerability curves with and without considering the uncertainty of ground motion intensity, there was no evident difference in shape, but due to the truncation range set, there are two nodes that make the curve divided into three sections, which can be regarded as a continuous piecewise linear function. In the range of frequent earthquakes or basic seismic intensity, the corrected damage ratio is higher than the original damage ratio. Next, approximately in the range of rare earthquakes, the original damage ratio exceeds the corrected damage ratio first, and the difference between the two curves gradually increases from small to large, then from large to small, and finally, the difference between the two curves is further reduced to almost overlapping.

In loss estimation of catastrophe insurance, amplitude uncertainty is particularly important. For possible earthquakes in the future, this paper provides a novel view and calculation for building vulnerability that fully considers the uncertainty of the three elements of ground motion, the uncertainty of vulnerability, and the coupling uncertainty between them. In addition, during processing the translation between seismic intensity and PGA, the distribution is adopted instead of PGA mean value, and the uncertainty in this process is considered as well. Furthermore, it can efficiently calculate the vulnerability of a large number of structures, which provides a fast and relatively accurately loss result. In the damage ratio sampling, β distribution is used to replace the lognormal distribution, and the modified vulnerability curve of the building structure is estimated. The method proposed in this study solves the problem that putting the sample mean into a complex calculation gives a different result from taking the sample mean after a complex calculation. However, the question of the percentage of moderate damage taking too much has not been solved and still needs to be studied from another perspective.