A method for probabilistic seismic hazard assessment considering site effect

Abstract

The application discloses a kind of probability seismic risk evaluation method considering site effect, belong to the field of seismic risk analysis, including: collecting historical earthquake data, measured overburden thickness and soil type, based on GIST model estimation regional overburden thickness, and conversion is V s30 Value to characterize site effect, under the condition of different correlation distance c Smooth estimation of seismic activity parameters a, b;Combination regional tectonic type selects NGA ground motion attenuation model, calculates propagation depth parameter Z 1.0 , Z 2.5 ; Input probability seismic risk calculation engine, simulate the PGA spatial distribution of 50 years exceeding probability 10%; Correlation analysis is carried out between the simulation results under the condition of different correlation distance c and national zoning data, and the optimal correlation distance c is determined by combining the fitting degree of historical earthquake point, and the corresponding PGA distribution is the final evaluation result. Through the application, the probability seismic risk evaluation result more in line with the actual site conditions and historical seismic activity characteristics of research area can be obtained.

Description

Technical Field

[0001] This invention belongs to the field of seismic hazard analysis technology, specifically relating to a probabilistic seismic hazard assessment method that considers site effects. Background Technology

[0002] Seismic hazard assessment is a crucial foundational task for seismic design and site safety evaluation of engineering projects. Its purpose is to quantitatively assess the potential seismic ground motion levels that a project site may experience within its design service life, providing a basis for seismic design of engineering structures, site suitability comparison, and seismic risk prevention and control. Current seismic hazard assessment methods mainly include two categories: deterministic seismic hazard analysis and probabilistic seismic hazard analysis. Deterministic seismic hazard analysis is primarily used for major projects with high seismic fortification requirements and clearly defined controlling earthquake scenarios, such as hydropower and nuclear power plants. Probabilistic seismic hazard analysis, on the other hand, comprehensively considers the uncertainty of earthquake occurrence, magnitude distribution characteristics, spatial distribution of seismic sources, and seismic ground motion attenuation patterns, balancing engineering safety and economic rationality, and has become an important technical means for seismic design in conventional industrial and civil engineering projects.

[0003] The probabilistic seismic hazard analysis method has been gradually developed and improved since its introduction into probability theory in the 1950s. In 1968, Cornell established the theoretical framework for using probabilistic seismic hazard analysis for the assessment of engineering site hazard. Since then, this method has been widely used in seismic zoning, determination of ground motion parameters for engineering sites, and formulation of seismic fortification standards. my country's current national standard GB18306-2015, "Zoning Map of Seismic Ground Motion Parameters in China," is also based on the probabilistic seismic hazard analysis method, using parameters such as peak ground acceleration and characteristic period of response spectrum to characterize regional seismic hazard, providing important support for seismic fortification of engineering projects.

[0004] However, for spoil heaps, construction camps, temporary facilities, and other individual projects along linear engineering projects in western mountainous areas, the site topography, overburden thickness, and soil type often vary significantly over short distances, resulting in a pronounced site effect that directly impacts seismic intensity, spectral characteristics, and duration. Limited by factors such as insufficient measured shear wave velocity, drilling data, and overburden thickness data, traditional probabilistic seismic hazard analysis typically uses bedrock sites or simplified site conditions as a calculation premise. This makes it difficult to fully characterize local site condition differences and their impact on seismic response, easily leading to underestimation or overestimation of site ground motion levels. Consequently, it affects the suitability assessment of engineering sites and the configuration of seismic protection measures, failing to meet the needs of refined comparison and safety evaluation for individual projects.

[0005] Existing related technologies still have certain limitations. For example, patent CN114693066B proposes a seismic hazard analysis method, device, equipment, and storage medium, which realizes seismic hazard analysis through probabilistic seismic hazard decomposition and equivalent earthquake calculation. However, this method does not incorporate site effect parameters into the calculation process, and the evaluation results are difficult to reflect the influence of actual regional site conditions on seismic ground motion response. The paper "Probabilistic Seismic Hazard Analysis Considering Site Effects—Taking Yuxi as an Example" uses the USGS V-seismic hazard analysis method based on topographic slope inversion. S30 The data characterizes the site effect, but does not further consider the influence of overburden thickness and soil type on the site effect; at the same time, the study fixes the relevant distance c=60km, and the parameter values ​​lack regional adaptability verification, which can easily lead to overestimation or underestimation of seismic activity within a specific time window, making it difficult to meet the needs of high-precision and regionalized seismic hazard assessment for single-point projects in complex mountainous areas. Summary of the Invention

[0006] This invention aims to address the problems existing in probabilistic seismic hazard assessment, such as insufficient characterization of site effects, subjective selection of seismic activity parameters, and difficulty in meeting the needs of fine comparison for single-point projects with regional-scale assessment results. It proposes a probabilistic seismic hazard assessment method that considers site effects.

[0007] This invention estimates and corrects the regional overburden thickness based on measured overburden thickness and soil type at specific locations, and converts the overburden thickness into the site shear wave velocity parameter V. S30 Furthermore, the propagation depth parameters required for the seismic motion attenuation model are calculated. By conducting seismic activity smoothing calculations under different correlation distance conditions and comparing and analyzing the PGA simulation results, a regional adaptation optimization method for correlation distance parameters is constructed, thereby obtaining probabilistic seismic hazard assessment results that better reflect the actual site conditions and historical seismic activity characteristics of the study area. To achieve the above objectives, the technical solution of this invention is as follows: A probabilistic seismic hazard assessment method considering site effects includes the following steps: Step 1: Obtain historical earthquake data for the study area, including earthquake time, epicenter coordinates, focal depth, magnitude, and magnitude scale; obtain measured overburden thickness and soil type data for points in the study area through field surveys, geophysical exploration, drilling, or existing data collection. Step 2: Perform magnitude scaling, earthquake catalog removal, and minimum complete magnitude calculation on the historical earthquake data to obtain an independent earthquake catalog that meets the requirements of probabilistic earthquake hazard analysis; based on the measured overburden thickness data, use the GIST model constrained by topographic factors to estimate the overburden thickness of the entire study area to obtain the spatial distribution results of the overburden thickness. Step 3: Determine the quadratic function model parameters between the overburden thickness and shear wave velocity based on the soil type of the study area, and convert the global overburden thickness obtained in Step 2 into shear wave velocity V. S30 To obtain V for characterizing shallow site effects S30 Spatial distribution results; Step 4: Divide the computational grid by combining the scale of the study area, the completeness of the earthquake catalog and the spatial distribution characteristics of historical earthquakes. Use the Gaussian smoothing model as the spatial smoothing model for seismic activity, and calculate the seismic activity distribution of each grid under multiple candidate correlation distances c. Step 5: Based on the preprocessed historical earthquake catalog, statistically analyze the relationship between magnitude and frequency of earthquakes in the study area, establish a magnitude-frequency relationship, and combine it with the grid seismic activity results obtained in Step 4 to determine the seismic activity parameters corresponding to each grid. Step Six: Based on the geological tectonic background, focal depth distribution, and spatiotemporal characteristics of the study area, select a suitable seismic motion attenuation model; based on the V obtained in Step Three... S30 Based on the spatial distribution results, the propagation depth parameters required for the selected seismic motion attenuation model are calculated to characterize the influence of site conditions on the propagation and attenuation process of seismic waves. Step 7: V S30 The spatial distribution, seismic activity parameters, ground motion attenuation model, and propagation depth parameters are input into the probabilistic seismic hazard calculation engine to simulate and calculate the spatial distribution of peak ground acceleration (PGA) in the study area under set exceedance probability and set age conditions under different candidate correlation distance c values; the probabilistic seismic hazard calculation engine is preferably the OpenQuake engine.

[0008] Step 8: Compare the degree of fit between the seismic activity distribution and the spatial distribution of historical earthquake points under different candidate correlation distance c values, and conduct correlation analysis between the PGA simulation results under different c values ​​and the national seismic ground motion parameter zoning data; combine the degree of fit of historical earthquake points, Pearson correlation coefficient and Spearman correlation coefficient to determine the optimal correlation distance c in the study area, and use the PGA spatial distribution corresponding to the optimal correlation distance as the final probabilistic seismic hazard assessment result.

[0009] Furthermore, in step two, the unification of earthquake magnitude scaling refers to converting earthquake events expressed in terms of surface wave magnitude and local magnitude into moment magnitude to improve the comparability between earthquake events from different data sources; the removal of clusters from the earthquake catalog refers to using the spatiotemporal window method to identify and remove non-independent events such as foreshocks and aftershocks, while retaining independent mainshock events; the calculation of the minimum complete magnitude refers to using the maximum curvature method to determine the lowest magnitude that can be reliably recorded in the study area within a given time range, that is, simultaneously plotting the magnitude-cumulative frequency curve and the magnitude-non-cumulative frequency curve, and taking the magnitude corresponding to the maximum curvature of the cumulative frequency curve as the minimum complete magnitude.

[0010] Furthermore, the estimation of the total cover thickness of the study area using the GIST model constrained by topographic factors includes: Using a GIS platform, slope units were divided through hydrological process analysis of forward and reverse DEMs, and slope shape was identified and classified based on the longitudinal profile morphology of the slope units. The standard curvature of each grid cell is calculated based on the regional DEM and normalized to obtain the slope curvature contribution Cc. Based on the soil type and the range of its physical and mechanical parameters in the study area, the internal friction angle and critical slope values ​​were determined. The slope grid was then reclassified to obtain the slope contribution A. c ; The measured overburden thickness points were grouped according to slope type, and 80% of the points were randomly selected within each slope type to calculate the relative positional contribution P of the slope. c ; The slope curvature contribution Cc, slope contribution Ac, and slope relative position contribution Pc are multiplied to form a comprehensive control index for overburden thickness, and the correction coefficient K is calculated back based on the measured overburden thickness. c The optimal K corresponding to different slope types c As a correction coefficient for this slope type, the final estimated value of the overburden thickness is obtained.

[0011] Furthermore, in step three, a quadratic function model is used to convert the global overburden thickness into shear wave velocity V. S30 The quadratic function model is as follows:

[0012] In the formula: V s denoted as soil shear wave velocity (m / s); H is the overburden thickness (m); e, f, and g are fitting parameters determined based on the soil type.

[0013] Furthermore, in step four, the Gaussian smoothing model assumes that seismic activity follows a Gaussian decay law in space, and uses the contribution of seismic events within a neighborhood of the grid center with a radius of 3c to the seismic activity of the target grid as the basis for calculating the seismic activity of that grid. The calculation formula is as follows:

[0014] In the formula: The number of earthquakes in the j-th grid that are greater than the minimum intact magnitude; Let i represent the seismic activity of the i-th grid. is the center distance between the i-th and j-th grids, in km; c is the correlation distance.

[0015] Furthermore, in step five, the magnitude-frequency relationship is represented using the Gutenberg-Richter relation:

[0016] In the formula: N is the number of earthquake events with a magnitude greater than or equal to M; M is the magnitude; a is a parameter characterizing the level of regional seismic activity; b is a parameter characterizing the proportion of earthquakes of different magnitudes.

[0017] For each grid cell, the corresponding seismic activity parameters are determined based on the seismic activity results and magnitude-frequency relationship obtained from the Gaussian smoothing model.

[0018] Furthermore, the seismic activity parameter b value of the individual grid is estimated using the maximum likelihood method, and the calculation formula is as follows:

[0019] In the formula: The average magnitude of earthquake events within a neighborhood centered at the grid center with a radius of 3c. The minimum intact magnitude; denoted as the magnitude band width; e represents the grid seismic activity calculated by the Gaussian smoothing model. The seismic activity parameter 'a' value for a single grid is calculated based on the magnitude-frequency relationship.

[0020] Furthermore, in step six, the seismic attenuation model includes ASK14, CY14, and CB14 models, wherein the ASK14 and CY14 models adopt Z... 1.0 The depth at which the shear wave velocity reaches 1.0 km / s is characterized by Z-wave velocity. 2.5 The depth at which the shear wave velocity reaches 2.5 km / s is characterized. Based on V S30 Calculate Z 1.0 and Z 2.5 This allows the site effect to be transmitted not only through V S30The data is incorporated into the seismic hazard calculation and simultaneously fed into the ground motion attenuation model via the propagation depth parameter to improve the completeness of the site effect characterization. 1.0 and Z 2.5 The calculation formula is as follows:

[0021] .

[0022] On the other hand, the present invention also proposes a probabilistic seismic hazard assessment device that considers site effects, comprising: The data acquisition module is used to acquire historical earthquake datasets for the study area, as well as measured overburden thickness datasets and soil type datasets at specific locations. The historical earthquake datasets include the time of occurrence, epicenter coordinates, focal depth, magnitude, and magnitude scale. The data preprocessing module is used to perform magnitude scaling, earthquake catalog removal, and minimum complete magnitude calculation on the historical earthquake dataset to obtain an independent earthquake catalog that meets the requirements of probabilistic earthquake hazard analysis. The overburden thickness estimation module is used to estimate the overburden thickness of the entire study area based on measured overburden thickness data and to perform numerical calculations using the GIST model constrained by topographic factors, thereby obtaining spatial distribution data of overburden thickness. The site effect parameter conversion module is used to determine the quadratic function model parameters between the overburden thickness and shear wave velocity according to the soil type in the study area, convert the overburden thickness data of the whole area into shear wave velocity VS30 data, and obtain VS30 spatial distribution data to characterize shallow site effects. The seismic activity calculation module is used to divide the calculation grid by combining the scale of the study area, the completeness of the earthquake catalog and the spatial distribution characteristics of historical earthquakes. It adopts a Gaussian smoothing model as the spatial smoothing model for seismic activity and calculates the seismic activity distribution data of each grid under multiple candidate correlation distances c. The seismic activity parameter determination module is used to establish a magnitude-frequency mathematical model based on the preprocessed independent earthquake catalog through the Gutenberg-Richter relationship, statistically analyze the relationship between magnitude and frequency of occurrence in the study area, and determine the seismic activity parameters corresponding to each grid by combining the grid seismic activity distribution data. The propagation depth parameter calculation module is used to select a ground motion attenuation model based on the geological structure background of the study area, and to calculate the propagation depth parameters required for the selected ground motion attenuation model based on the VS30 spatial distribution data. The hazard simulation calculation module is used to input the VS30 spatial distribution data, seismic activity parameters, ground motion attenuation model and propagation depth parameters into the probabilistic seismic hazard calculation engine to simulate and calculate the spatial distribution data of ground peak acceleration (PGA) in the study area under the set exceedance probability and set age conditions under different candidate correlation distance c values. The optimal parameter determination and evaluation result output module is used to calculate the degree of fit between the distribution of seismic activity and the spatial distribution of historical earthquake points under different candidate correlation distance c values, as well as the correlation between the PGA simulation calculation results under different c values ​​and the nationally published seismic ground motion parameter zoning data. The module comprehensively determines the optimal correlation distance c for the study area and outputs the PGA spatial distribution data corresponding to the optimal correlation distance as the final probabilistic seismic hazard assessment result.

[0023] The present invention also proposes a computer-readable storage medium storing a computer program, which, when executed by a processor, implements the steps in a probabilistic seismic hazard assessment method considering site effects as described above.

[0024] In summary, compared with the prior art, the present invention has the following advantages: 1. Based on the traditional probabilistic seismic hazard assessment framework, this invention considers the combined control of overburden thickness and soil type on site effects, and transforms this into shear wave velocity V. S30 Including seismic hazard calculation overcomes the problem that traditional methods, which only use bedrock site conditions or single site parameters, are unable to reflect the differences in actual site seismic response, thus improving the site adaptability and reliability of the evaluation results.

[0025] 2. This invention constructs a GIST model based on geomorphological features, uses point-based overburden thickness inversion to estimate the spatial distribution of regional overburden thickness, and combines soil type to establish a relationship between overburden thickness and shear wave velocity V. S30 The conversion relationship between these parameters significantly enhances the ability to refine the characterization of site effect parameters. This method can effectively support the selection of site seismic safety suitability for projects such as spoil heaps and construction camps along linear engineering lines, and is particularly suitable for complex mountainous areas with widespread overburden and relatively scarce measured data.

[0026] 3. In the process of seismic activity calculation, this invention sets multiple candidate correlation distances *c* to conduct trial calculations of seismic activity. By comparing the consistency between the distribution of seismic activity and the spatial distribution of historical earthquake points, as well as the correlation between the PGA simulation results and the national seismic ground motion parameter zoning data, the optimal correlation distance *c* and its corresponding probabilistic hazard assessment results are comprehensively determined. This process reduces human intervention in the selection of correlation distance parameters and improves the scientific rigor, objectivity, and regional applicability of parameter selection.

[0027] 4. In processing historical earthquake data, this invention improves the consistency and comparability of earthquake data from different sources by unifying the magnitude scale, removes the influence of non-independent events such as foreshocks and aftershocks by removing clusters from the earthquake catalog, and determines the sample range of earthquake events by calculating the minimum complete magnitude, thereby improving the accuracy of seismic activity analysis and probabilistic seismic hazard assessment results. Attached Figure Description

[0028] Figure 1 The overall flowchart of the probabilistic seismic hazard assessment method considering site effects in this invention is shown below. Figure 2 A spatial distribution map of historical earthquakes in the study area; Figure 3 The graph shows the calculation results of the minimum complete magnitude for the study area; Figure 4 This is a flowchart of the regional overburden thickness inversion estimation based on the GIST model; Figure 5 This is a map showing the general topography and geomorphology of the study area and the elevation distribution characteristics of typical profiles. Figure 6 Distribution map of slope curvature Cc, slope gradient Ac, and slope relative position contribution Pc in the study area; Figure 7 A map showing the distribution of overburden thickness in the study area; Figure 8 For study area V S30 Spatial distribution map; Figure 9 This is a map showing the distribution of seismic activity under different correlation distances c. Figure 10 This is a PGA distribution map of seismic hazard without considering site effects. Figure 11 A PGA distribution map for seismic hazard considering site effects; Figure 12 Box plots showing the seismic activity distribution of all grids and historical earthquake points in the study area under different correlation distances c; Figure 13 Heatmap of Pearson correlation coefficients between PGA results and nationally published data under different correlation distances c; Figure 14 A heatmap showing the Spearman correlation coefficients between PGA results and nationally published data under different correlation distances c. Detailed Implementation

[0029] To make the technical objectives, technical solutions, and beneficial effects of the present invention clearer, the present invention will be further described below in conjunction with preferred embodiments and accompanying drawings. Those skilled in the art should understand that the following detailed description is illustrative rather than restrictive and should not be construed as limiting the scope of protection of the present invention.

[0030] Example 1 See appendix Figure 1 This embodiment proposes a probabilistic seismic hazard assessment method that considers site effects. The implementation process is as follows: Step 1: Collect historical earthquake data for the study area and obtain measured data on overburden thickness and soil type at the study area locations through field surveys.

[0031] The study area selected in this embodiment is located in the transition zone between a basin and a plateau, with a dense distribution of rivers and an overall terrain that slopes from west to east and from north to south. The overburden layer is widespread, faults are well-developed, and there are numerous single-point projects such as spoil heaps and construction camps along the fault lines. The study area extends from 29.6°N–30.5°N and 101.15°E–102.52°E. The topographic and geomorphological overview is as follows: Figure 5 As shown.

[0032] A total of 6,577 historical earthquake events were collected and compiled from databases including the National Earthquake Science Data Center, China Earthquake Networks Center, USGS, and the United States Seismological Research Association. Each earthquake event record includes the time of occurrence, epicentral coordinates, focal depth, magnitude, and magnitude scale, as shown in Table 1. The number of earthquake events under different magnitude scales is further analyzed, as shown in Table 2, where the magnitude scale includes surface wave magnitude. Moment Magnitude Local earthquake magnitude .

[0033] Table 1. Sources of data for probabilistic seismic hazard assessment

[0034] Table 2. Statistics on the number of earthquakes at different magnitude scales (2009–2022)

[0035] Step 2: Perform magnitude scaling standardization, seismic catalog removal, and minimum complete magnitude calculation on historical earthquake data to obtain an independent earthquake catalog that meets the requirements of probabilistic seismic hazard analysis. The total overburden thickness of the study area is estimated based on field-measured overburden thickness data. Magnitude scaling standardization ensures consistency between earthquake data from different sources in spatiotemporal activity analysis and engineering applications. Since moment magnitude Mw is more closely related to fault rupture parameters and does not have a magnitude saturation problem, this embodiment uses moment magnitude Mw. w As a unified scale, the surface wave magnitude M will be used. s With local magnitude ML The earthquake events represented are converted into moment magnitude M according to equations (1) and (2). w : (1) (2) After unifying the magnitude scale, the moment magnitude, longitude, latitude, and focal depth of each earthquake event were obtained, and a spatial distribution map of historical earthquakes in the study area was drawn based on this data, such as... Figure 2 As shown. By Figure 2 It can be seen that the seismic activity in the study area is significantly controlled by the Xianshuihe fault structure, and the regional earthquakes are mainly of magnitude 0-4, with focal depths mostly concentrated in the range of 0-30 km.

[0036] The earthquake catalog is used to meet the requirement of temporal independence of earthquake events in probabilistic seismic hazard assessment. In this embodiment, the Gardner-Knopoff spatiotemporal window method is used to identify and remove non-independent events such as foreshocks and aftershocks based on the time and spatial windows corresponding to different magnitudes, ultimately obtaining 241 independent mainshock events.

[0037] Minimum intact magnitude M c This represents the lowest reliably recordable magnitude within a given spatiotemporal range. This embodiment employs a simple and efficient maximum curvature method, simultaneously plotting magnitude-cumulative frequency curves and magnitude-non-cumulative frequency curves, and taking the magnitude corresponding to the maximum curvature of the cumulative frequency curve as the minimum intact magnitude M. c Calculation results show that the minimum intact magnitude M in the study area is... c It is 1.7, M c A total of 145 earthquake events were recorded above, and the calculation results are as follows: Figure 3 As shown.

[0038] Secondly, based on the measured overburden thickness data at field points in the study area, the GIST model was used to estimate the overall overburden thickness of the study area. The specific process is as follows: Figure 4 As shown.

[0039] First, using GIS, slope units were delineated through hydrological process analysis of forward and reverse DEMs. This process included DEM depression filling, flow direction extraction, flow calculation, river network extraction and classification, catchment area generation, and merging of positive and negative catchment areas, ultimately obtaining slope units with ridgelines or valleylines as boundaries. The DEM raster resolution used in this embodiment was 12.5m × 12.5m, and the study area contained a total of 12,266,640 raster units.

[0040] Then, slope shape was identified and classified based on the longitudinal profile statistics of the slope units. According to the slope development stage and profile morphology characteristics, the slopes in the study area were divided into three categories: convex slopes, convex-concave slopes, and convex-concave-convex slopes. The standard curvature of each grid was calculated based on the regional DEM and normalized to obtain the slope curvature contribution Cc. Combining the overburden soil type and its physical and mechanical parameter ranges in the study area, the internal friction angle and critical slope value were determined, and the slope grids were reclassified to obtain the slope contribution A. c The field survey points were grouped according to slope type, and 80% of the measured points were randomly selected within each slope type to calculate the relative positional contribution P of the slope. c .

[0041] Finally, the slope curvature contribution Cc, slope contribution Ac, and slope relative position contribution Pc are multiplied to form a comprehensive control index for overburden thickness; the correction coefficient K is then calculated based on the measured overburden thickness. c And assign the corresponding optimal K to different slope types. c Finally, the estimated value of the overburden thickness is obtained.

[0042] The slope curvature contribution C of the study area was obtained according to the above procedure. c Slope contribution A c and relative position contribution P c Spatial distribution, such as Figure 6 As shown, the thickness distribution of the overburden layer across the entire study area is obtained, as shown below. Figure 7 As shown. By Figure 7 The GIST model estimates the total cover thickness in the study area to be 0-63.94 m, generally showing a spatial distribution pattern of thicker cover in valleys and thinner cover at mountain peaks. Specifically, points 1 and 2 are located in valleys with a cover thickness greater than 40 m; points 3 and 4 are located at mountain peaks with a cover thickness less than 20 m; and points 5 and 6 are located at the transition zone between mountain peaks and valleys with a cover thickness of 30-40 m.

[0043] In this embodiment, 502 measured values ​​of overburden thickness were obtained, of which 80% were used for GIST model calibration and estimation, and the remaining 20% ​​were used for model validation. The root mean square error (RMSE) was used to represent the model accuracy. The results show that the RMSE of the overburden thickness estimated by the GIST model is approximately 6.8 m. Furthermore, by combining the engineering survey data from some sites in the study area with the average overburden thickness obtained by traditional exploration methods such as drilling and geophysical exploration, it can be seen that the GIST model estimation results are quite close to the measured results (Table 3), indicating that the model can effectively characterize the spatial distribution characteristics of overburden thickness across the entire study area.

[0044] Table 3 Comparison of overburden thickness obtained by different methods

[0045] Step 3: Convert the overburden thickness in the study area into shear wave velocity V. S30 To construct spatial parameters for characterizing site effects.

[0046] Based on the measured overburden thickness and soil type data at the study sites, the overburden thickness was converted into the corresponding shear wave velocity V. S30 The relationship between the overburden thickness and the shear wave velocity is represented by a quadratic function model: (3) In the formula: V s denoted as soil shear wave velocity in m / s; H is the thickness of the overburden layer in m; e, f, and g are model fitting parameters determined according to the soil type.

[0047] Table 4 shows the model parameters and goodness-of-fit R² corresponding to different soil types. In this embodiment, the overburden soil type in the study area can be generally classified as gravelly soil. Based on the model goodness-of-fit R², e = 206.99, f = 6.5589, and g = -0.026.

[0048] Table 4. Model parameters and goodness of fit for common soil types

[0049] After obtaining the point location or grid scale V S30 After the value is obtained, Kriging interpolation is used to obtain the V value of the study area. S30 Spatial distribution, such as Figure 8 As shown. By Figure 8 It can be seen that V on the north and northeast sides of the study area S30 Higher value, southeast side V S30 The value is relatively low. Based on the topographic transition and active fault distribution in the study area, it is inferred that the thicker overburden on the south side is the cause of the V value in this region. S30 A key reason for the low score.

[0050] Step 4: Divide the computational grid and use a Gaussian smoothing model to calculate the distribution of seismic activity under different correlation distances c.

[0051] This embodiment employs a Gaussian smoothing model to spatially smooth historical earthquake point sources. This model assumes that seismic activity spatially approximately follows a Gaussian attenuation law. Within the Gaussian kernel smoothing framework, and considering the study area scale and seismic data, a computational grid is first divided according to latitude and longitude of 0.05° × 0.05°, with each grid approximately 5.5km × 5.5km. The number of earthquake events with a minimum complete magnitude or higher within each grid is then counted. Next, a Gaussian kernel function is used to spatially smooth the earthquake events. The contribution of earthquake events within a 3c radius neighborhood centered on the grid center to the target grid characterizes the seismic activity of that grid. The calculation formula is as follows: (4) In the formula: The number of earthquakes in the j-th grid that are greater than the minimum intact magnitude; Let i represent the seismic activity of the i-th grid. is the center distance between the i-th and j-th grids, in km; c is the correlation distance.

[0052] To determine the appropriate correlation distance *c* for the study area, this embodiment conducted trial calculations with *c* values ​​of 2, 3, 4, 5, 10, 20, 30, 40, and 50 km. The results showed that when *c* = 2 km and *c* = 3 km, the smoothed seismic activity field deviated significantly from the spatial distribution of historical earthquakes in the study area, resulting in distorted local smoothing results; therefore, these values ​​were discarded. Finally, correlation distances of *c* = 4, 5, 10, 20, 30, 40, and 50 km were selected as candidate correlation distances, and seismic activity distribution results under different *c* values ​​were obtained, as follows: Figure 9 As shown.

[0053] Step 5: Establish the magnitude-frequency relationship and obtain seismic activity parameters for each grid in the study area.

[0054] Seismic activity typically follows the Gutenberg-Richter magnitude-frequency relationship: (5) In the formula: N is the number of earthquake events with a magnitude greater than or equal to M; M is the magnitude; a is a parameter characterizing the level of regional seismic activity; b is a parameter characterizing the proportion of earthquakes of different magnitudes.

[0055] In probabilistic seismic hazard assessment, N represents the number of seismic events within a neighborhood of 3c centered at the grid center, and M is the minimum intact magnitude of the earthquake in the study area. The b-value is an important parameter characterizing the magnitude-frequency distribution and the level of regional seismic activity. This embodiment uses the maximum likelihood method to calculate the b-value, and the calculation formula is as follows: (6) In the formula: represents the average magnitude of the earthquake events; e represents the seismic activity calculated using the Gaussian smoothing model. The minimum intact magnitude; This refers to the magnitude grading width. Given that the magnitude intervals in actual earthquake catalogs are mostly 0.1, this embodiment... Take 0.2. Mean magnitude. Calculated using the following formula (7): (7) Since this embodiment uses a Gaussian smoothing model to calculate seismic activity at the grid scale, theoretically, the b-value can be calculated grid-by-grid. However, the number of seismic events within a single grid is usually small, making it difficult to obtain stable and reliable b-values. To improve computational stability, this embodiment calculates a unified regional b-value based on a complete seismic catalog of the study area and uses it for subsequent seismic hazard analysis. The calculation results show that the b-value for the study area is 0.638, with a standard deviation of 0.003. Subsequently, the seismic activity of each grid and the Gutenberg-Richte relation are obtained by combining the Gaussian smoothing model, and the corresponding a-value is calculated grid-by-grid.

[0056] Step 6: Select a seismic motion attenuation model and based on V S30 Calculate the propagation depth parameters corresponding to the calculation model.

[0057] Seismic motion attenuation models, also known as seismic motion prediction equations, are a key component of probabilistic seismic hazard assessment. When the magnitude is constant, seismic motion intensity typically decreases with increasing epicentral distance; conversely, when the epicentral distance is constant, seismic motion intensity increases with increasing magnitude. Considering the geological structure and spatiotemporal distribution characteristics of the study area, located on the eastern edge of the Tibetan Plateau, shallow-focus earthquakes are predominant. To account for the cognitive uncertainties surrounding seismic motion attenuation, this embodiment selects three shallow crustal seismic motion attenuation models—CY14, CB14, and ASK14—as candidate models, which can be directly used in the OQ-engine software.

[0058] Site conditions significantly influence seismic response. Bedrock sites typically exhibit smaller vibration amplitudes and shorter durations, while sites with deep overburden may show longer vibration periods and more pronounced seismic amplification effects. The ASK14 and CY14 models use depth Z corresponding to a shear wave velocity Vs reaching 1 km / s. 1.0 The influence of overburden thickness is characterized, while the CB14 model uses the depth Z corresponding to a shear wave velocity Vs reaching 2.5 km / s. 2.5 The influence of the capping layer thickness is characterized. In this embodiment, the shear wave velocity V obtained by fitting the capping layer thickness is utilized. S30 The value is obtained by back-calculating the corresponding Z using equations (8) and (9). 1.0 and Z 2.5 : (8) (9) In the formula: V S30 Z represents the average shear wave velocity over a depth range of up to 30m, expressed in m / s. 1.0 Z represents the depth at which the shear wave velocity reaches 1 km / s, expressed in meters (m). 2.5 The depth at which the shear wave velocity reaches 2.5 km / s, in km.

[0059] Step 7: V S30 The distribution, seismic activity parameters, ground motion attenuation model, and propagation depth parameters are input into the probabilistic seismic hazard calculation engine to simulate and calculate the spatial distribution of the peak ground acceleration (PGA) with a 10% exceedance probability in the study area under different candidate correlation distances c. In this embodiment, the probabilistic seismic hazard calculation engine uses the OpenQuake engine.

[0060] In this embodiment, the upper limit of the focal depth is set to 0 km, and the lower limit is set to 40 km. Given that earthquake focal depths are mainly concentrated in the 0–40 km range, and considering the uncertainty of focal depth, a logic tree model is established with three branches representing focal depths of 5 km, 15 km, and 25 km, with weights of 0.18, 0.70, and 0.12, respectively. The magnitude distribution adopts a truncated Gutenberg-Richter relationship, with a minimum intact magnitude Mc of 1.7 and a maximum magnitude of 6.3, and a magnitude box width of 0.2. Based on previous research results and the actual conditions of the study area, the weights of the ASK14, CY14, and CB14 models are assigned 0.25, 0.25, and 0.50, respectively, to construct a logic tree for the multi-seismic motion attenuation model.

[0061] To account for site effects, each grid center point is assigned a corresponding V. S30 The value is obtained by back-calculating Z using formulas (8) and (9). 1.0 and Z 2.5 This is to reflect the impact of changes in overburden thickness on seismic attenuation.

[0062] For ease of comparison, this embodiment also calculates the probabilistic seismic hazard without considering site effects, i.e., bedrock site conditions, and takes c=50km. The results are shown in [reference needed]. Figure 10 .Depend on Figure 10 It can be seen that, without considering the site effect, the seismic hazard PGA of the study area is distributed in the range of 0.29-0.337g.

[0063] According to relevant standards, this embodiment selects a seismic protection standard with a 50-year exceedance probability of 10% for general engineering projects. Figure 11 The distribution of the 50-year exceedance probability (PGA) for the study area at different correlation distances *c*, considering site effects, is presented. The results show that, after considering site effects, the PGA distribution of the seismic hazard in the study area ranges from 0.32 to 0.36g, an overall increase of approximately 0.03g compared to not considering site effects. Its spatial distribution characteristics generally exhibit a banded pattern from west to east, resembling "low-medium-high-medium," which is largely consistent with the basic peak ground acceleration distribution law of the study area published in GB18306-2015 "Zoning of Seismic Ground Motion Parameters in China," indicating the rationality of the calculation results of this invention.

[0064] Step 8: Analyze the fit between seismic activity results and the spatial distribution of historical earthquake locations under different correlation distances c. Conduct correlation analysis between the PGA simulation results and the data published by the state, and comprehensively determine the optimal correlation distance c and the final probabilistic seismic hazard assessment results for the study area.

[0065] Under different candidate correlation distances *c*, seismic activity values ​​were extracted for all grids in the study area and the grids containing historical earthquake points, and box plots were drawn, as shown below. Figure 12 As shown. From Figure 12 It can be seen that when the correlation distance c=4km, the spatial distribution of historical earthquake points has a high degree of spatial overlap with areas of high seismic activity, and the standard deviation of seismic activity in the grid containing the historical earthquake points is significantly higher than that of all grids in the study area. This indicates that the seismic activity calculation results under this condition can better reflect the spatial distribution characteristics of historical earthquakes. In contrast, when c=10km to 50km, the distribution range of seismic activity values ​​in both types of grids narrows significantly, the maximum and minimum values ​​are close, and outliers appear in some cases, making it difficult to fully reflect the spatial differences in actual seismic activity in the study area.

[0066] Furthermore, Pearson correlation analysis and Spearman correlation analysis were used to evaluate the correlation between the PGA simulation results and the nationally published seismic motion parameter zoning data under different correlation distances c.

[0067] Pearson correlation analysis measures the degree of linear correlation between two variables. Its correlation coefficient ranges from -1 to 1, where 1 indicates a perfect positive correlation, -1 indicates a perfect negative correlation, and 0 indicates no linear correlation. The formula for calculating the Pearson correlation coefficient is as follows: (10) In the formula: r is the Pearson correlation coefficient; X and Y are two data samples; the numerator is the sample covariance, and the denominator is the product of the standard deviations of the two samples.

[0068] This embodiment performs Pearson correlation analysis on the PGA simulation results under different correlation distances c and the national seismic motion parameter zoning results. The results are as follows: Figure 13 As shown in the figure, the results indicate that when c=4km, the Pearson correlation coefficient between the PGA simulation results and the nationally published zoning results reaches a maximum of 0.5882; in contrast, when c=40km, the Pearson correlation coefficient is only 0.301. This result suggests that the PGA simulation results under the c=4km condition are in better agreement with the nationally published zoning results and better reflect the spatial distribution characteristics of PGA in the study area.

[0069] Spearman correlation analysis is a rank-based nonparametric correlation analysis method suitable for evaluating monotonic correlations between variables. It has good applicability to non-normal distributions and nonlinear relationships. The formula for calculating the Spearman correlation coefficient is as follows: (11) In the formula: r s is the Spearman correlation coefficient; n is the number of samples; d is the difference set formed by the rank differences of corresponding elements in the sample set.

[0070] This embodiment performs Spearman correlation analysis on the PGA simulation results under different correlation distances c and the national seismic motion parameter zoning results. The results are as follows: Figure 14 As shown in the figure, the results indicate that the Spearman correlation coefficient between the PGA simulation results and the nationally published zoning results is highest at c=4km, reaching 0.55; in contrast, the Spearman correlation coefficient is only 0.38 when c=30km. This result further demonstrates that the PGA simulation results under the c=4km condition have a higher consistency with the nationally published zoning results and can more accurately reflect the actual spatial distribution characteristics of PGA in the study area.

[0071] Example 2 Based on the same inventive concept, this embodiment proposes a probabilistic seismic hazard assessment device considering site effects, including a data acquisition module, a data preprocessing module, a cover layer thickness estimation module, a site effect parameter conversion module, a seismic activity calculation module, a seismic activity parameter determination module, a propagation depth parameter calculation module, a hazard simulation calculation module, and an optimal parameter determination and evaluation result output module. The functions of each module correspond completely to those in the aforementioned method embodiment (steps one to eight), and will not be elaborated further. Among them: The data acquisition module is used to acquire historical earthquake datasets for the study area, as well as measured overburden thickness datasets and soil type datasets for points in the study area; the historical earthquake datasets include the time of occurrence, epicenter coordinates, focal depth, magnitude, and magnitude scale.

[0072] The data preprocessing module is used to perform magnitude scaling, earthquake catalog removal, and minimum complete magnitude calculation on the historical earthquake dataset to obtain an independent earthquake catalog that meets the requirements of probabilistic earthquake hazard analysis.

[0073] The overburden thickness estimation module is used to estimate the overburden thickness of the entire study area based on measured overburden thickness data and the GIST model constrained by topographic factors, thereby obtaining spatial distribution data of overburden thickness.

[0074] The site effect parameter conversion module is used to determine the quadratic function model parameters between overburden thickness and shear wave velocity based on the soil type in the study area, converting the overburden thickness data across the entire area into shear wave velocity V. S30 Data were used to obtain V for characterizing shallow site effects. S30 Spatial distribution data.

[0075] The seismic activity calculation module is used to divide the calculation grid by combining the scale of the study area, the completeness of the earthquake catalog and the spatial distribution characteristics of historical earthquakes. A Gaussian smoothing model is used as the spatial smoothing model for seismic activity, and the seismic activity distribution data of each grid is calculated under multiple candidate correlation distances c.

[0076] The seismic activity parameter determination module is used to establish a magnitude-frequency mathematical model based on the preprocessed independent earthquake catalog through the Gutenberg-Richter relationship, statistically analyze the relationship between magnitude and frequency of occurrence in the study area, and determine the seismic activity parameters corresponding to each grid by combining the grid seismic activity distribution data.

[0077] The propagation depth parameter calculation module is used to select a seismic motion attenuation model based on the geological structure background of the study area, and based on the V... S30 Spatial distribution data is used to calculate the propagation depth parameters required for the selected seismic attenuation model.

[0078] The hazard simulation calculation module is used to calculate the V S30 The spatial distribution data, seismic activity parameters, ground motion attenuation model, and propagation depth parameters are input into the probabilistic seismic hazard calculation engine to simulate and calculate the spatial distribution data of ground peak acceleration (PGA) in the study area under different candidate correlation distance c values, with set exceedance probabilities and set time limits.

[0079] The optimal parameter determination and evaluation result output module is used to calculate the degree of fit between the seismic activity distribution and the spatial distribution of historical earthquake points under different candidate correlation distance c values, as well as the correlation between the PGA simulation calculation results under different c values ​​and the national published seismic motion parameter zoning data. The optimal correlation distance c for the study area is determined in a comprehensive manner, and the PGA spatial distribution data corresponding to the optimal correlation distance is output as the final probabilistic seismic hazard assessment result.

[0080] Example 3 This invention also provides a computer-readable storage medium storing a computer program, which, when executed by a processor, implements the steps of a probabilistic seismic hazard assessment method considering site effects as described in Embodiment 1 above.

[0081] Specifically, the computer-readable storage medium can be any tangible medium capable of carrying and storing computer program code, including but not limited to: USB flash drives, portable hard drives, read-only memory (ROM), random access memory (RAM), magnetic disks, optical disks (such as CD-ROM and DVD-ROM), flash cards, solid-state drives (SSDs), or server storage and other non-volatile storage media.

[0082] When the computer program runs on a computing device such as a personal computer, server, workstation, or dedicated seismic data processing terminal, the processor calls and executes the instructions in the computer program to perform the following operations in sequence: Acquire historical earthquake datasets, measured overburden thickness datasets, and soil type datasets for the study area through data interfaces or input devices; The historical earthquake dataset was subjected to magnitude scaling unification, earthquake catalog removal, and minimum complete magnitude calculation to form an independent earthquake catalog. The spatial distribution of overburden thickness in the entire study area was estimated based on the GIST model. Convert the overall capping layer thickness into shear wave velocity V. S30 Spatial distribution data is used to calculate the propagation depth parameter Z required for the selected seismic attenuation model. 1.0 or Z 2.5 ; The Gaussian smoothing model was used to calculate the grid seismic activity distribution under multiple candidate correlation distances c, and the seismic activity parameters a and b of each grid were determined by combining the Gutenberg-Richter relationship. Input the above parameters into a probabilistic seismic hazard calculation engine (such as the OpenQuake engine) to simulate and obtain the spatial distribution of the peak ground acceleration (PGA) corresponding to different c values ​​under the conditions of set exceedance probability and set age. By combining the spatial fit between the distribution of seismic activity and historical earthquake locations, the Pearson correlation coefficient and Spearman correlation coefficient between PGA simulation results and the nationally published seismic ground motion parameter zoning data, the optimal correlation distance c is determined, and the corresponding PGA spatial distribution is output as the final probabilistic seismic hazard assessment result.

[0083] It should be noted that the specific implementation details of each of the above steps have been described in detail in Embodiment 1, and will not be repeated in this embodiment. Those skilled in the art will understand that all or part of the process of the method of Embodiment 1 can be implemented by a computer program instructing related hardware (such as a processor, memory, and communication interface). The computer program can be stored in the computer-readable storage medium described in this embodiment, and when the program is executed, it can execute the process including the steps of Embodiment 1 or any combination thereof.

[0084] The computer-readable storage medium described in this embodiment can be sold or used as a standalone software product, or it can be built into a seismic hazard analysis system, an engineering seismic fortification auxiliary decision-making platform, or a regional seismic safety evaluation device, and provided to end users in the form of a pre-installed or upgraded installation package.

[0085] Furthermore, the computer-readable storage medium described in this embodiment is not limited to a standalone operating environment, but can also be applied to a distributed computing system or cloud computing platform, where multiple computing nodes collaboratively execute the computer program to jointly complete the probabilistic seismic hazard assessment task for the study area.

Claims

1. A probabilistic seismic hazard assessment method considering site effects, characterized in that, Includes the following steps: Step 1: Obtain historical earthquake data for the study area, as well as measured overburden thickness and soil type data at the study area locations; Step 2: Preprocess the historical earthquake data to obtain an independent earthquake catalog that meets the requirements of probabilistic earthquake hazard analysis; Based on measured overburden thickness data, the overburden thickness of the entire study area was estimated using the GIST model constrained by topographic factors, and the spatial distribution results of overburden thickness were obtained. Step 3: Determine the quadratic function model parameters between the overburden thickness and shear wave velocity based on the soil type of the study area, and convert the total overburden thickness into shear wave velocity V. S30 To obtain V for characterizing shallow site effects S30 Spatial distribution; Step 4: Divide the computational grid by combining the scale of the study area, the completeness of the earthquake catalog and the spatial distribution characteristics of historical earthquakes. Use the Gaussian smoothing model as the spatial smoothing model for seismic activity, and calculate the seismic activity distribution of each grid under multiple candidate correlation distances c. Step 5: Establish a magnitude-frequency relationship based on the independent earthquake catalog, and determine the seismic activity parameters corresponding to each grid by combining the grid seismic activity distribution; Step 6: Select a seismic motion attenuation model based on the V S30 Spatial distribution calculation of the propagation depth parameters required for the selected ground motion attenuation model; Step 7: V S30 The seismic hazard calculation engine inputs spatial distribution, seismic activity parameters, ground motion attenuation model and propagation depth parameters to simulate and calculate the spatial distribution of peak ground acceleration (PGA) in the study area under set exceedance probability and set age conditions under different candidate correlation distance c values. Step 8: Based on the degree of fit between the distribution of seismic activity and the spatial distribution of historical earthquake points under different candidate correlation distance c values, as well as the correlation between the PGA simulation results and the nationally published seismic ground motion parameter zoning data, determine the optimal correlation distance for the study area, and use the PGA spatial distribution corresponding to the optimal correlation distance as the final probabilistic seismic hazard assessment result.

2. The probabilistic seismic hazard assessment method considering site effects according to claim 1, characterized in that, The historical earthquake data includes the time of occurrence, epicenter coordinates, focal depth, magnitude, and magnitude scale; The historical earthquake data is preprocessed, including: magnitude scaling, removal of clusters from the earthquake catalog, and calculation of the minimum complete magnitude. The unified magnitude scaling refers to converting earthquake events expressed as surface wave magnitude and local magnitude into moment magnitude; the earthquake catalog cluster removal refers to using the spatiotemporal window method to identify and remove non-independent events such as foreshocks and aftershocks, while retaining independent mainshock events; the minimum complete magnitude calculation refers to using the maximum curvature method to determine the lowest magnitude that can be reliably recorded in the study area within a given time range, that is, simultaneously plotting the magnitude-cumulative frequency curve and the magnitude-non-cumulative frequency curve, and taking the magnitude corresponding to the maximum curvature of the cumulative frequency curve as the minimum complete magnitude.

3. The probabilistic seismic hazard assessment method considering site effects according to claim 1, characterized in that, The estimation of the total overburden thickness in the study area using the GIST model constrained by topographic factors includes: Using a GIS platform, slope units were divided through hydrological process analysis of forward and reverse DEMs, and slope shape was identified and classified based on the longitudinal profile morphology of the slope units. The standard curvature of each grid cell is calculated based on the regional DEM and normalized to obtain the slope curvature contribution Cc. Based on the soil type and the range of its physical and mechanical parameters in the study area, the internal friction angle and critical slope values ​​were determined. The slope grid was then reclassified to obtain the slope contribution A. c ; The measured overburden thickness points were grouped according to slope type, and 80% of the points were randomly selected within each slope type to calculate the relative positional contribution P of the slope. c ; The slope curvature contribution Cc, slope contribution Ac, and slope relative position contribution Pc are multiplied to form a comprehensive control index for overburden thickness, and the correction coefficient K is calculated back based on the measured overburden thickness. c The optimal K corresponding to different slope types c As a correction coefficient for this slope type, the final estimated value of the overburden thickness is obtained.

4. The probabilistic seismic hazard assessment method considering site effects according to claim 1, characterized in that, In step three, a quadratic function model is used to convert the global overburden thickness into shear wave velocity V. S30 The quadratic function model is as follows: In the formula: V s denoted as soil shear wave velocity (m / s); H is the thickness of the overburden layer (m); e, f, and g are fitting parameters determined based on the soil type.

5. The probabilistic seismic hazard assessment method considering site effects according to claim 1, characterized in that, In step four, the Gaussian smoothing model uses the contribution of seismic events within a neighborhood of 3c, centered on the grid center, to the seismic activity of the target grid as the basis for calculation. The calculation formula is as follows: In the formula: The number of earthquakes in the j-th grid that are greater than the minimum intact magnitude; Let i represent the seismic activity of the i-th grid. is the center distance between the i-th and j-th grids, in km; c is the correlation distance.

6. The probabilistic seismic hazard assessment method considering site effects according to claim 5, characterized in that, In step five, the magnitude-frequency relationship is represented using the Gutenberg-Richter relation: In the formula: N is the number of earthquake events with a magnitude greater than or equal to M; M is the magnitude; a is a parameter characterizing the level of regional seismic activity; b is a parameter characterizing the proportion of earthquakes of different magnitudes.

7. The probabilistic seismic hazard assessment method considering site effects according to claim 6, characterized in that, The seismic activity parameter b for a single grid is estimated using the maximum likelihood method, and the calculation formula is as follows: In the formula: The average magnitude of earthquake events within a neighborhood centered at the grid center with a radius of 3c. The minimum intact magnitude; denoted as the magnitude band width; e represents the grid seismic activity calculated by the Gaussian smoothing model. The seismic activity parameter 'a' value for a single grid is calculated based on the magnitude-frequency relationship.

8. The probabilistic seismic hazard assessment method considering site effects according to claim 1, characterized in that, The seismic attenuation model is selected from ASK14, CY14, and CB14 models, among which ASK14 and CY14 models adopt Z-axis... 1.0 The depth at which the shear wave velocity reaches 1.0 km / s is characterized by Z-wave velocity. 2.5 The depth at which the shear wave velocity reaches 2.5 km / s; based on V S30 Calculate Z 1.0 and Z 2.5 The calculation formulas are as follows: 。 9. A probabilistic seismic hazard assessment device considering site effects, characterized in that, include: The data acquisition module is used to acquire historical earthquake datasets for the study area, as well as measured overburden thickness datasets and soil type datasets at specific locations. The historical earthquake datasets include the time of occurrence, epicenter coordinates, focal depth, magnitude, and magnitude scale. The data preprocessing module is used to perform magnitude scaling, earthquake catalog removal, and minimum complete magnitude calculation on the historical earthquake dataset to obtain an independent earthquake catalog that meets the requirements of probabilistic earthquake hazard analysis. The overburden thickness estimation module is used to estimate the overburden thickness of the entire study area based on measured overburden thickness data and to perform numerical calculations using the GIST model constrained by topographic factors, thereby obtaining spatial distribution data of overburden thickness. The site effect parameter conversion module is used to determine the quadratic function model parameters between overburden thickness and shear wave velocity based on the soil type in the study area, and to convert the overburden thickness data of the entire area into shear wave velocity V. S30 Data were used to obtain V for characterizing shallow site effects. S30 Spatial distribution data; The seismic activity calculation module is used to divide the calculation grid by combining the scale of the study area, the completeness of the earthquake catalog and the spatial distribution characteristics of historical earthquakes. It adopts a Gaussian smoothing model as the spatial smoothing model for seismic activity and calculates the seismic activity distribution data of each grid under multiple candidate correlation distances c. The seismic activity parameter determination module is used to establish a magnitude-frequency mathematical model based on the preprocessed independent earthquake catalog through the Gutenberg-Richter relationship, statistically analyze the relationship between magnitude and frequency of occurrence in the study area, and determine the seismic activity parameters corresponding to each grid by combining the grid seismic activity distribution data. The propagation depth parameter calculation module is used to select a seismic attenuation model based on the geological structure background of the study area, and based on the V... S30 Spatial distribution data is used to calculate the propagation depth parameters required for the selected seismic attenuation model. The hazard simulation calculation module is used to calculate the V S30 Spatial distribution data, seismic activity parameters, ground motion attenuation models, and propagation depth parameters are input into a probabilistic seismic hazard calculation engine to simulate and calculate the spatial distribution data of peak ground acceleration (PGA) in the study area under different candidate correlation distance c values, with set exceedance probabilities and set time limits. The optimal parameter determination and evaluation result output module is used to calculate the degree of fit between the distribution of seismic activity and the spatial distribution of historical earthquake points under different candidate correlation distance c values, as well as the correlation between the PGA simulation calculation results under different c values ​​and the nationally published seismic ground motion parameter zoning data. The module comprehensively determines the optimal correlation distance c for the study area and outputs the PGA spatial distribution data corresponding to the optimal correlation distance as the final probabilistic seismic hazard assessment result.

10. A computer-readable storage medium, characterized in that, The computer-readable storage medium stores a computer program that, when executed by a processor, implements the steps of a probabilistic seismic hazard assessment method considering site effects as described in any one of claims 1 to 8.

Images