Probabilistic Seismic Hazard Analysis of a Back Propagation Neural Network Predicting the Peak Ground Acceleration

Abstract

Probabilistic seismic hazard analysis (PSHA) has been recognized as a reasonable method for quantifying seismic threats. Traditionally, this method ignores the effect of the focal depth, in which the ground motion prediction equations (GMPEs) are applied to estimate the probability distribution associated with the possible motion levels induced by the site earthquakes, but it is limited by the unclear geological conditions, which makes it difficult to provide a uniform equation, and these equations cannot express the non-linear relationship under geological conditions. Hence, this paper proposed a method to consider the seismic focal depth for the PSHA with the example of California and used a back propagation neural network (BPNN) to predict the peak ground acceleration (PGA) instead of the GMPEs. Firstly, the measured PGA and unknown PGA seismic data applicable to this method were collected separately. Secondly, the unknown PGA data were supplemented by applying the BPNN based on the measured PGA data. Lastly, based on the full-probability equation, PSHA considering the focal depth was completed and compared with the current California seismic zoning results. The results showed that using the BPNN in the PSHA can ensure computational accuracy and universality, making it more suitable for regions with unclear geological structures and providing the possibility of adding other parameters to be considered for the influence of the PSHA.

1. Introduction

Earthquakes are one of the most influential disasters in the entire life cycle of engineering structures [1]. To avoid damage to engineering structures caused by earthquakes, seismic hazard analysis methods can be used to predict the probability of earthquakes of various strengths in the future and prevent them in their structural designs [2]. Seismic hazard analysis includes the deterministic and probabilistic methods [3]. The deterministic seismic hazard analysis (DSHA) requires the prior determination of the seismic source and magnitude, but neglects the uncertainties in the data and calculations, which could lead to an overestimation of the risk for the seismic source with longer return periods or an underestimation of the risk for the seismic source with shorter return periods. Since Cornell proposed the probabilistic seismic hazard analysis method in 1968, this method has been used in the context of engineering design and in the verification of the performance requirements of structures induced by earthquakes [4]. The probabilistic seismic hazard analysis (PSHA) has been generally recognized as the rational method to quantify seismic threats [5]. The PSHA includes four aspects (as shown in Figure 1): (1) the seismic source characteristics, (2) the relationship between the magnitude and the frequency of the earthquake, (3) the attenuation relationship of the ground motion parameters, and (4) the hazard curve [6]. However, there are three uncertainties in the analysis process: the determination of the potential seismic source area, the estimation of the seismicity parameters, and the relationship between the parameters [7]. Surrounding the uncertain factors, how to improve the PSHA is the focus of scholars’ continuous research.

Researchers have refined the uncertainties in the determination of potential seismic source regions caused by an incomplete understanding of subsurface structures and earthquakes by adjusting the seismic zoning and seismic source models [8,9]. The most important method to determine the potential seismic source area is to gain an in-depth understanding of the geological structural parameters, such as the soil parameters [10] and regional tectonic distribution characteristics [11], along with the deep earth tectonic [12] and seismic zone distributions [13]. In addition, the parameters of the seismic source mechanism can be effectively solved by defining the different seismic source models. The fact that these parameters vary widely in space and time, combined with the effect that different geological formations are interacting systems, makes the accurate analysis of potential seismic source zones very difficult [14,15]. Continuous refinement of seismic zoning can precisely determine the parameters of the area and apply them to calculate the PSHA for each differential zone; however, detailed seismic zoning is not conducive to integration into region-wide, full-probability calculations [16]. Certain parameters, such as the average annual occurrence of earthquakes, earthquake magnitude values, and upper limits of earthquake magnitude, can lead to uncertainties in the valuation of the seismicity parameters. The seismic catalogue obtained from using seismic monitoring techniques can supplement this part of the parameters with a large number of actual measurements, but the number of such monitoring results is still limited by the capacity of the measuring instruments and the number of arrangements. To consider a more complete seismic catalogue, previous studies [17,18] have combined the Monte-Carlo simulations of historical seismicity catalogues based on the measured data but may affect the accuracy of the PSHA. For the uncertainty underlying the determination of the potential seismic source regions and seismic activity parameter valuations, multiple combinations of scenarios may be used to calculate the seismic hazard analysis results of different scenarios while considering the range in the variation of the results of all calculated scenarios to provide the mean and variance of the seismic hazard analysis results.

The uncertainty in the relationship between the seismic parameters originates from the unknown ground motion attenuation relation, which is currently obtained through fitting the empirical formulae to known ground motion-related data, among which include ground motion intensity indicators (e.g., PGA and spectral acceleration [19]), earthquake magnitude, and distance (epicenter distance or hypocentral distance). This formula is known as the GMPE, and currently, GMPEs can be classified into three types [20]: functional, controlling, and predictive. Functional formulas require the determination of functional relationships among the parameters; controlling formulas can estimate the independent effects among the parameters experimentally; and predictive formulas are only used for the calculation of non-deterministic relationships among the data. All physical aspects associated with the earthquake ground motion are not known [21], and even if they are known, they cannot be expressed in the form of a simple equation; the estimated ground motion relationship is the predictive model. After a large number of scholars’ studies on the GMPE, the deviation value of the GMPE in these seismic areas has been reduced year by year [22,23], but it still cannot provide a single definite GMPE [24]. For those areas with low seismicity, the selection of a suitable GMPE for the PSHA is an even more challenging step [25,26].

Undoubtedly, the in-depth analysis of geological formations to solve the problems caused by uncertainty is the fundamental method to improve the PSHA, but facing the problems of many influencing factors and unclear unknown formations is not conducive to the use and design in the engineering field. The effect of the focal depth on the ground motion, which is widely known, has not been systematically considered in the PSHA [27]. This neglect is also present in the GMPE. Therefore, this paper focused on using known data to predict the PGA through the non-linear fitting ability of the BPNN, while the PSHA is based on the full probability formulation considering the effect of the focal depth. The organization of this paper is as follows. Section 2 describes the selection guidelines for the California earthquake catalogue, which includes both known PGA and unknown PGA data, in detail. The BPNN-based PGA prediction method and processes have been described in Section 3, and a preliminary judgment on the accuracy of the fitting results has been made. In Section 4, the PSHA considering the depth of the seismic source has been completed, with the data having been fitted using the BPNN. In Section 5, the resulting hazard curves are compared with the seismic zoning results in ACSE 7–16 to determine the accuracy of the method proposed in this paper. At the end of this paper, the relevant conclusions and recommendations have been summarized in Section 6.

2. Earthquake Catalogs

The BPNN-based PGA predictions need to be supported with as much data as possible, and the data mentioned here include both known PGA data and unknown PGA data. In this context, researchers examining these effects in California [28], as an earthquake-prone region, also started to record and analyze their seismic data earlier in this area; therefore, in this paper, we counted the seismic data in the region, combined with the time of database establishment and the completeness of known data, and determined the time range of data acquisition from the years 1970 to 2021.

2.1. Known PGA Data

The known PGA data were acquired and recorded to the strong motion database through the strong motion station, and the main known strong motion databases include the Next-Generation Attenuation (NGA) database [29], the Center for Engineering Strong Motion Data (CESMD) database [30], and the COSMOS (Consortium of Organizations for Strong-Motion Observation Systems) Virtual Data Center database [31]. Considering the data volume, data time range, and demand parameters, the final data set utilized was the CESMD, jointly established by the United States Geological Survey (USGS) and the California Geological Survey (CGS), and the strong motion data and strong motion station are shown in Figure 2.

The Initial screening of the obtained data was performed. Since the boxed range was controlled by the latitude and longitude, the selected area in Figure 2 is larger than the California area, and the redundant data were removed by combining the actual California range in the data. With the increase in the number of strong seismograph arrangements, more strong motion data were measured for the same strong earthquake. The increase in the amount of data can improve the accuracy of the prediction results, but the propagation of the ground motion is greatly affected by the existence of the earth’s soil, and the more distant the measured locations are, the higher the loss in the confidence rate; therefore, only the data closest to the seismic source were selected to be used in this study. The measured PGAs were provided in the form of two horizontal components and one vertical component, and the largest PGA in the same earthquake was used due to the disparity between the ground motion components, which has not been considered in this paper.

The GMPE generally calculates the PGA using the magnitude and distance (epicenter distance or hypocentral distance), but in fact, the epicenter distance is the distance between the ground location of the seismic source and the seismic station, and the hypocentral distance is the distance between the seismic source and the seismic station, so the calculated PGA is affected by the focal depth. To obtain more accurate PGA calculation results, the magnitude, epicentral distance, focal depth, and PGA were obtained in this paper. In total, 107 known PGA data points were obtained through combining the above selection principles. The relationship between the magnitude and the epicenter distance of these data is shown in Figure 3.

2.2. Unknown PGA Data

In addition to the strong motion database, the earthquake database mainly records the magnitude, epicenter distance, and focal depth. There is no PGA (unknown PGA), but the PGA data is a necessary part of PSHA. The unknown PGA database can complete the earthquake catalogues to reduce the problem of overestimating the annual PGA exceedance probability due to the strong motion data. According to the basic requirements for strong motion data selection, this part of the data needs to ensure at least three kinds of information with the magnitude, epicenter distance, and focal depth, which can be provided by the United States Geological Survey (USGS) [32], and the data acquisition process is shown in Figure 4.

The initial screening of the obtained data was performed. Consistent with the screening of the known strong motion data, the data for the non-California regions and the presence of the required parameters (magnitude, epicenter distance, and focal depth) were removed. Since the seismic hazard is calculated for destructive earthquakes, the magnitude of the destructive earthquakes is generally taken to be about 4.5 or higher and considering that the minimum value of earthquake magnitude in the strong motion data is 4.1, some data with an earthquake magnitude of 4.0 or higher have been retained. A total of 2645 data points were obtained from the California area. The frequency distribution histogram (Figure 5) was obtained with a magnitude difference of 0.1. The data displayed in Figure 5 are concentrated in lower magnitudes, which is consistent with the occurrence pattern of earthquakes.

3. BPNN-Based PGA Prediction

The application of neural networks to compute multi-parameter fitting problems has been widely used in structural engineering [33], where a non-linear computational power can effectively reduce the time to analyze the mechanism while ensuring computational accuracy. Lin Jyh-Woei [34] used the BPNN to predict the probability of earthquakes in Taiwan, a method that can be commercialized at a relatively low cost and with minimal resources and equipment using only the earthquake catalogues.

Combined with the statistical data, it is easy to find that the PGA, as an important part of the strong motion data, is affected by the propagation speed and the location of the strong motion station and has not yet achieved the purpose of obtaining effective real motion data for each earthquake; the amount of recorded seismic data is much larger than the strong motion data, and the supplemental PGA, through seismic data, has also become an important part of the PSHA.

Previously, the GMPE was used to supplement the PGA, and this method generally fitted the seismic magnitude separately according to the polar difference of one (e.g., Richter scales four, five, and six). The fitted equations also differed due to the distinct geological considerations, and the goodness-of-fit was low. The GMPE is an empirical prediction relationship, as shown in Equation (1) [35] and Equation (2) [36], that represents the ground motion parameters (including the PGA, spectral acceleration at a given frequency, peak velocity, or peak displacement) as a function of other parameters (such as the earthquake magnitude and epicentral distance):

$$\ln PGA=C_1+C_{2}m+C_3\ln R+C_{4}R$$

(1)

$$\ln PGA=C_1+C_{2}m-C_3\ln \left(R+C_4\exp \left(C_{5}m\right)\right)$$

(2)

where m is the earthquake magnitude; R is the epicentral distance or hypocentral distance; and C₁, C₂, C₃, C₄, and C₅ are constants that are independent of m and R.

As a traditional PGA prediction method, the GMPE varies greatly under different regions and is influenced by the geological soil and other conditions; however, there is still no uniform model. This paper makes use of the better non-linear computational ability of the BPNN to compensate for the computational errors caused by the lack of understanding of geology in the engineering field and to improve the problem that multiple GPMEs need to be considered for supplementing PGAs across different regions. Due to the complex non-linear correlation between the input and the output quantities, the accuracy of the previous primary function fitting was deemed to be poor, while the neural network, as a tool to describe the correlation between the multiple quantities, has a strong non-linear fitting ability, and the above mapping relationship has been considered by the neural network. The process of predicting the PGA using the BPNN is as follows:

(1)

The known PGA data are randomly split into a training set, a validation set, and a test set in the ratio of 7:1.5:1.5 to avoid errors in the prediction results caused by large differences in the order of magnitude using data normalization.

(2)

The BPNN (Figure 6) is constructed with three elements (a magnitude of Richter scale, focal depth, and the epicentral distance) in the input layer and one element (PGA) in the output layer. To avoid underfitting and overfitting, the number of hidden layers is set to two, each layer has four neurons, and the number of neurons is obtained by considering the number of elements in the input and output layers in combination with several attempts, where the transfer function (transig) and the training function (trainlm) are non-linear functions [37].

(3)

The data set is trained with the determinants of the coefficients that are greater than 0.95 as the index to obtain the fitted network, where the regression fitting results are shown in Figure 7.

(4)

After saving the network, the normalized unknown PGA data are inputted into the network. The predicted data are obtained and inverse normalized to obtain the final predicted PGA.

The regression fitting results of the California data are shown in Figure 7. There is randomness in the calculation process of the BPNN, and a 0.95 goodness-of-fit value of the total data was selected as the calculation index in this study, although, this index can still be stricter. As the total data was divided into different data sets, the fit superiority between each data set differs, and the superiority of the total data set is easier to reach the index when the training data set accounts for the most and its superiority is higher.

4. PSHA

4.1. The PSHA Method

The PSHA is obtained based on the full probability formula [38]; to consider the influence of the focal depth on the ground motion, this paper combined the existing data available for regression analysis to select the best-fit formula. Adding the probability density function about the depth of the earthquake, in this regard, the annual average transcendence probability of a seismic event that meets the PGA greater than a given value is shown in Equations (3)–(8):

$$\begin{array}{l}P\left(PGA>x\right)=\lambda \left(m>m_{\min }\right)\\ {\displaystyle \iiint P\left(PGA>x|m,R,D\right)f\left(m\right)f\left(R\right)f\left(D\right)dmdRdD}\end{array}$$

(3)

$$\lambda \left(m>m_{\min }\right)=\left(A_D/A\right)\lambda$$

(4)

$$f(m)=\frac{b\ln \left(10\right){10}^{-b\left(m-m_{\min }\right)}}{1-{10}^{-b\left(m_{\max }-m_{\min }\right)}}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em} m_{\min }<m<m_{\max }$$

(5)

$${\log }_{10}\left(N_m/T\right)=a-bm$$

(6)

$$f\left(R\right)=\frac{d}{dR}F\left(R\right)=\frac{d}{dR}\left(\frac{\pi R^2}{\pi R_t^2}\right)=\frac{2R}{R_t^2} 0<R<R_t$$

(7)

$$f\left(D\right)=a_1{e}^{\left(a_{2}D\right)}+a_3{e}^{\left(a_{4}D\right)}$$

(8)

where $P\left(PGA>x|m,R,D\right)$ denotes the conditional probability of the PGA exceeding x given m, R, and D; $\lambda \left(m>m_{\min }\right)$ is the average annual occurrence rate of earthquake events with a magnitude greater than the minimum magnitude [39], based on the California earthquake catalogue [40]; $A_D$ is the area of the analyzed region; $A$ is the area of the analyzed region and its 300 km radius; $\lambda$ is the average annual occurrence rate; $f(m)$ is the probability density function of the magnitude; $m_{\min }$ and $m_{\max }$ are the minimum and maximum earthquakes that may occur in the seismic source, respectively [41]; $b$ denotes the Gutenberg–Richter coefficient [42]; $N_m$ is the annual number of seismic events in a given magnitude range; $T$ is the earthquake catalogue period; $f\left(R\right)$ is the probability density function of the epicentral distance; $R$ is the epicentral distance [43]; $R_t$ denotes the limit of the epicentral distance; $f\left(D\right)$ is the probability density function of the focal depth; $D$ is the focal depth; and $a_1,a_2,a_3,a_4$ represent the fitting coefficient.

4.2. Cumulative Probability Distribution of the PGA

It is necessary to determine the probability distribution of the PGA in the PSHA and judge the accuracy of the BPNN in supplementing the PGA data through the probability distribution results that were previously obtained through fitting. In this section, all the California data were divided into two categories: Sample A only contains the measured PGA and Sample B contains the measured data and the BPNN prediction data.

Four commonly used asymmetric probability models, namely the Gumbel, Frechet, Weibull, and generalized Pareto distribution (GPD), were selected to model the PGA probability distribution for analysis. The cumulative distribution function (CDF) of the above four models is shown in Equations (9)–(12), and the values of the parameters in the probability distribution function through regression analysis are listed in

Table 1. Regression results of the CDF of the PGA.

Type	CDF	a	b	c	R2
Sample A	Gumbel	122.4	145.5	——	0.9797
	Frechet	−157.6	268.8	2.179	0.9942
	Weibull	10.17	223.5	0.9294	0.9971
	GPD	6.825	215.6	0.09101	0.9969
Sample B	Gumbel	46.7	42.6	——	0.9724
	Frechet	−3.254 × 104	3.259 × 104	765.3	0.9723
	Weibull	55.5	136	2.958	0.9755
	GPD	6.944	73.03	2.676 × 10−7	0.9591

.

$$\mathrm{Gumbel}: F\left(x\right)=\exp \left[-\exp \left(-\frac{x-a}{b}\right)\right]$$

(9)

$$\mathrm{Frechet}: F\left(x\right)=\exp \left[-{\left(-\frac{x-a}{b}\right)}^{-c}\right]$$

(10)

$$\mathrm{Weibull}: F\left(x\right)=1-\exp \left[-{\left(\frac{x-a}{b}\right)}^c\right]$$

(11)

$$\mathrm{GPD}: F\left(x\right)=1-{\left[1+c\left(x-a\right)/b\right]}^{-\frac{1}{c}}$$

(12)

The closer that the determinants of the coefficients (R²) in Table 1 are to one, the better their goodness-of-fit values. Weibull displayed the best goodness-of-fit value in both samples. BPNN prediction had less impact on the regularity of the data, which still conforms to the original probability distribution. The CDFs of the two samples are plotted in Figure 8.

It is not difficult to interpret that inflexion points appeared in the regions with more concentrated PGA values. Combining Figure 3 and Figure 5, the known PGAs were often larger ground motion events, while the more frequent earthquakes were relatively concentrated at lower magnitudes. This indicates that the PGA results obtained via the BPNN prediction will increase the number of medium and low PGAs, and the change of this data volume eventually affects the CDF of the PGAs.

4.3. Logic Tree

In the PSHA, the uncertainties underlying the magnitude, focal depth, epicenter distance, b-value, PGA data prediction, and the CDF of the PGA have significant effects on the accuracy of the results. Section 3 adjusted the PGA data prediction method through the use of the BPNN, while improving the uncertainty of the magnitude and the uncertainty of the focal depth through the regression analysis. The logic tree method can be used to refine the uncertainty caused by other factors.

Branch I (Figure 9) is used to consider the uncertainty of the b-value, whose value can be determined from the seismic data. According to Roshan [44], three values of the b-value, namely b − ε, b, and b + ε, as calculated from the maximum likelihood method, are considered. For the current database, the b-value that was estimated using the maximum likelihood method was calculated to be 0.9125 (±0.0718). A 75% weightage was assigned to the estimated value that was estimated from the database and the bounds on the b-value following the accounting error were each assigned a weightage of 12.5%.

Branch Ⅱ is used to consider the uncertainty of the CDF. Ayman Shama [35] took equal weights for each model. Villani Manuela [45] combined a stochastic model weight of 35% and an empirical model weight of 65% to fit the ground motion results. DU W [46] assigned the 10 published GMPEs to three zones: the subduction zone, the strike-slip fault, and the deformation background zone, with different weights. The models that were selected in this study took generality into account, meaning that each model had a good fit, but there was also a certain gap, and using equal weights will pay too much attention to the models with a lower fit. Thus, this paper assigned weights according to the weighted average of the fit coefficients in Table 1. The branch Ⅲ of the logic tree is the final weight of each item.

5. Results

After considering the influence of the uncertainty factors on the PSHA using the logic tree, the annual average occurrence of seismic events regarding the PGA can be determined according to Equation (3). In Section 4.2, the Sample B database had a similar fitting result as Sample A. Meanwhile, in order to determine the advancement in the BPNN prediction results, this paper plotted the sample B database in Figure 10 and compared it with the results in the specification.

Assuming that the number of seismic events in a given time obeys the Poisson distribution [47], the probability of k occurrences of a seismic event can be expressed as:

$$P\left(X=k\right)=\frac{{\lambda }^k{e}^{-\lambda }}{k!}$$

(13)

The probability that no earthquake occurs in a given time can be expressed by taking a value of zero for k. Therefore, the probability that an earthquake event satisfying the PGA is greater than a given limit, x, which occurs at least once a year, t, can be expressed as:

$$P\left(\mathrm{PGA}>x \mathrm{in} t \mathrm{years}\right)=1-e^{-P\left(\mathrm{PGA}>x\right)t}$$

(14)

The inverse of Equation (14) yields the formula for the average annual transcendence rate of a certain intensity seismic event:

$$P\left(\mathrm{PGA}>x\right)=\frac{\ln \left[1-P\left(\mathrm{PGA}>x \mathrm{in} t \mathrm{years}\right)\right]}{t}$$

(15)

Comparing the PGAs for the California zone and ASCE 7-16 [48], the probability of exceeding the ground motion in 50 years was found to be equal to 10%. With t = 50, the average annual occurrence corresponding to a given intensity seismic event was determined to be 2.10 × 10⁻³. Based on the seismic hazard curves determined in this section, the PGA value of the corresponding average annual occurrence seismic event was calculated as 0.345 g. ASCE 7-16 provides two different frequencies of spectral acceleration across the California area, where the short-period PGA was found to be larger compared to the long-period PGA. Hence, the 0.2 s spectral acceleration for the main seismic zone located in California was selected and converted for the Class D sites, resulting in about 0.353 g. The error between the two was about 2.27%. The calculation results obtained in this study underestimate the seismic hazard of the region, probably due to the fact that the data used in this paper covered the whole region of California, which contains some areas with the lower data, and the data used as the comparison terms are the results of the denser seismic areas in the California region. Moreover, the certain gap observed in the initial data caused the final error.

6. Conclusions

In this paper, the BPNN prediction PGA approach has been proposed instead of the GMPE, which takes into account the effect of the focal depth and is based on the full probability formula, while adding the focal depth to the analysis of the PSHA. Following the calculation with data from the California area as an example, the results are compared with the ASCE results, and the main conclusions are drawn as follows:

(1)

The results of the seismic hazard curves show that the probability of exceeding the ground motion equal to 10% in 50 years for the California region is in close agreement with the data in ASCE 7-16, with an error of 2.27%.

(2)

The BPNN method of predicting the PGA is convenient for the PSHA in areas with unknown geological conditions, which can be only performed using the seismic catalogues alone and can ensure the accuracy and generalizability of the analysis results while having the ability to take into account the non-linear relationships that exist between multiple parameters.

(3)

Based on the full probability formula, adding the effect of the focal depth on the annual average transcendental probability formula for the seismic events that satisfy the PGA greater than a given value can improve the PSHA.

(4)

When determining the weights of the CDFs through the logic tree, the use of weighting the goodness-of-fit coefficients enables more emphasis to be placed on the CDFs with better fitting results, but without biasing the results due to overemphasis.

The choice of using the BPNN to predict the PGA in this paper raises the possibility of adding more influencing factors to the PSHA. This method can be used in more areas that have more earthquake data, with a better accuracy and generalizability, making a tentative step for future scholars’ research.