A seismic risk analysis method considering historical seismic energy release

Abstract

This invention belongs to the field of potential earthquake risk prediction and analysis technology, and discloses a method for earthquake hazard analysis that takes into account the energy release of historical earthquakes. The method includes the following steps: constructing a three-dimensional grid of strain energy in the study area; calculating the energy release and range of historical earthquakes; calculating the cumulative strain energy of each grid; calculating the remaining strain energy and identifying potential seismogenic areas; and performing earthquake hazard analysis. This invention breaks through the existing analysis methods based on numerical simulation results, considering the energy release of historical earthquakes and the accumulation of inter-seismic energy from a new perspective and method, and proposes a new approach to earthquake hazard analysis based on numerical simulation results.

Description

Technical Field

[0001] This invention belongs to the field of potential earthquake risk prediction and analysis technology, specifically relating to an earthquake hazard analysis method that takes into account the energy release of historical earthquakes. Background Technology

[0002] Seismic hazard analysis is a crucial method for assessing the likelihood and intensity of future earthquakes in a given region. Common seismic hazard analysis methods primarily include deterministic and probabilistic approaches, which express seismic hazard through ground motion parameters. However, both methods rely on historical earthquake ground motion monitoring data, which is lacking for older earthquakes, affecting the accuracy of seismic hazard quantification. With the rapid development of science and technology, computer computing power has greatly improved, and finite element numerical simulation has been widely applied in seismic hazard analysis research. However, most numerical simulation methods only simulate the stress and strain of regional faults to determine potential future seismogenic areas, without adequately considering the energy release from historical earthquakes in the region.

[0003] Therefore, a new seismic hazard analysis method that takes into account the energy release of historical earthquakes is needed. Summary of the Invention

[0004] The purpose of this invention is to provide a seismic hazard analysis method that takes into account the energy release of historical earthquakes, in order to solve the problem that most current numerical simulation methods proposed in the background art only simulate the stress and strain of regional faults to determine potential future seismogenic areas, without well considering the energy release of historical earthquakes in the region.

[0005] To achieve the above objectives, the present invention provides a seismic hazard analysis method that takes into account the energy release of historical earthquakes, comprising the following steps:

[0006] S1. Based on the results of numerical simulation of regional fine geological structure, obtain the strain energy density data of each node of fault and stratum. Use interpolation method to interpolate the fault and stratum data respectively, construct a three-dimensional grid of strain energy density in the study area, and then multiply by the grid volume to obtain the three-dimensional grid of strain energy. Calculate the annual average cumulative amount of strain energy in each grid according to the three-dimensional grid of strain energy.

[0007] S2. Obtain and process the historical earthquake catalog of the study area, and calculate the energy released by each earthquake and the extent of energy release.

[0008] S3. Determine the number of years of accumulated strain energy in the grid and calculate the accumulated strain energy of each grid.

[0009] S4. Calculate the three-dimensional grid distribution of residual strain energy in the study area after energy release and energy accumulation, and identify potential seismogenic areas, their source locations, and magnitudes.

[0010] S5. Calculate the peak ground acceleration (PGA) distribution on the Earth's surface after each earthquake based on the location of the earthquake source and the magnitude of the earthquake in the potential seismogenic area. Then, superimpose these values ​​to form the potential PGA distribution on the surface of the entire study area. Classify the regional seismic hazard level according to the PGA.

[0011] In one specific implementation, in step S1, the strain energy density data of faults and stratigraphic nodes are interpolated using an interpolation method to construct a three-dimensional grid of strain energy density for the study area, which is then multiplied by the grid volume to obtain the three-dimensional grid of strain energy. Specifically:

[0012] Cubic spline interpolation, Kriging interpolation, inverse distance weighted interpolation, natural neighborhood interpolation, and Gaussian interpolation were selected for interpolation. The interpolation results were then used for inverse interpolation and compared with the original data. The interpolation method with the smallest error was selected to interpolate the strain energy density data of the strata and fault nodes, construct a three-dimensional strain energy density grid, and then multiply it by the grid volume to obtain the three-dimensional strain energy grid.

[0013] In one specific implementation, step S1 involves calculating the annual average cumulative strain energy for each grid, specifically as follows:

[0014] The strain energy of each grid is calculated by multiplying the obtained three-dimensional strain energy density grid by the grid volume. Then, the overall strain energy distribution is divided by the loading years in the original three-dimensional fine geological structure numerical simulation to obtain the annual average cumulative strain energy of each grid.

[0015] In one specific implementation, step S2 specifically includes the following steps:

[0016] S2.1 Historical earthquake catalog processing;

[0017] S2.2 Calculation of energy release from historical earthquakes;

[0018] S2.3, Fault identification;

[0019] S2.4 Calculation of fracture length, fracture width and fracture volume;

[0020] S2.5, Distribution of energy release from historical earthquakes.

[0021] In one specific implementation, S2.1, historical earthquake catalog processing, specifically involves: homogenizing the historical earthquake catalog by uniformly converting different magnitude types into surface wave magnitude Ms;

[0022] S2.2 Calculation of energy release from historical earthquakes: Specifically, the energy released by historical earthquakes is calculated based on the Gutenberg-Richter relation.

[0023] S2.3 Fault identification, specifically: based on the three-dimensional coordinates of fault nodes, identify the extent, strike, and dip of faults; use principal component analysis to calculate the strike and dip of each fault; and save the calculation results.

[0024] S2.4 Calculation of rupture length, rupture width and rupture volume: Specifically, the rupture length, rupture width, average slip and rupture volume after historical earthquakes are calculated according to the empirical formula of Wells & Coppersmith (1994).

[0025] In one specific implementation, S2.5, the distribution of historical earthquake energy release, specifically includes:

[0026] First, determine if the hypocenter is on a fault. If it is, construct a cube centered on the hypocenter, along the fault strike with the rupture length as the length, along the dip with the rupture width as the width, and along the vertical direction with the average slip as the height. This ensures that all grid energy within the cube is released, while guaranteeing that the released energy does not exceed the total energy released by the earthquake. If the released energy exceeds the total energy released by the earthquake, the cube is scaled down proportionally until the total energy of the internal grid is less than the energy released by the earthquake. If the released energy does not exceed the total energy released by the earthquake, the remaining energy is distributed to the grid outside the cube using a coordinate system centered on the hypocenter and with the cube's length, width, and height as the x, y, and z axes, respectively. The standard deviations of the 3D Gaussian distribution in the x, y, and z directions are set to RLD / 4, RLD / 6, and RLD / 8, respectively.

[0027] If the epicenter is not on the fault, a coordinate system is established with the epicenter as the center, the mutually perpendicular horizontal directions as the x and y axes, and the vertical direction as the z axis. The standard deviations of the x, y, and z directions are taken as RLD / 4, RLD / 4, and RLD / 6, respectively, where RLD represents the rupture length.

[0028] In one specific implementation, step S3 specifically includes the following steps:

[0029] S3.1 Calculation of the number of years of grid strain energy accumulation: Determine the initial year in which the grid releases energy due to the earthquake. The number of years of grid strain energy accumulation can be obtained by subtracting the initial year from the target year of the earthquake hazard assessment.

[0030] S3.2 Calculation of mesh strain energy accumulation: E ac =E an ×N ac ;

[0031] In the formula, E ac E represents the cumulative strain energy of the mesh. an The annual average cumulative strain energy of the grid, in J and N.ac Years for accumulating grid strain energy.

[0032] In one specific implementation, step S4 specifically includes the following steps:

[0033] S4.1 Calculation of regional residual strain energy: Subtract the strain energy released by historical earthquakes from the three-dimensional strain energy grid obtained by interpolation, and add the inter-seismic accumulated strain energy to obtain the distribution of residual strain energy. The grids that are not affected by earthquakes will no longer accumulate strain energy.

[0034] S4.2 Identification of potential seismogenic zones: Set the strain energy threshold for grid cells, use the seed region growth method to identify grid cells with values ​​above the threshold and form potential seismogenic zones, and calculate the total strain energy within the zone, while ensuring that the total energy does not exceed the energy that can be released by the maximum magnitude.

[0035] S4.3 Determination of the source location in the potential seismogenic zone: First, determine whether there is a fault in the potential seismogenic zone; if there is a fault, take the center point of the grid with the highest strain energy on the fault as the source point; if there is no fault, directly take the center point of the grid with the highest strain energy in the region as the source point.

[0036] S4.4 Calculation of magnitude in potential seismogenic zones: Calculate the possible magnitude of earthquakes in potential seismogenic zones based on the Gutenberg-Richter relation.

[0037] In one specific implementation, step S5 specifically includes the following steps:

[0038] S5.1 Calculation of Peak Ground Acceleration: Based on the study area, select an appropriate seismic attenuation equation to calculate the surface peak ground acceleration distribution after an earthquake occurs in each potential seismogenic zone, and superimpose the peak ground acceleration values ​​of each potential earthquake to form the total peak ground acceleration distribution.

[0039] S5.2 Seismic Hazard Analysis: The peak ground acceleration distribution of the study area is divided into several levels, thereby classifying the overall seismic hazard level of the study area.

[0040] Compared with the prior art, the present invention has the following beneficial effects:

[0041] This invention enables a more scientific and comprehensive assessment of regional seismic hazard. By interpolating the grid node data in the original numerical simulation results, a new three-dimensional distribution of strain energy density is obtained, leading to the distribution of strain energy. Simultaneously, the annual average cumulative strain energy is calculated, considering the scope and manner of energy release from historical earthquakes and post-earthquake energy accumulation, to obtain the distribution of residual strain energy in the region. This further identifies potential seismogenic zones, their source locations, and magnitudes, calculates the PGA distribution, and derives the overall seismic hazard distribution of the regional surface. This invention breaks through the limitations of existing numerical simulation analysis methods, considering the energy release from historical earthquakes and post-earthquake energy accumulation from a new perspective and approach, proposing a new approach to seismic hazard analysis based on numerical simulation results.

[0042] In addition to the objectives, features, and advantages described above, the present invention has other objectives, features, and advantages. The present invention will now be described in further detail. Attached Figure Description

[0043] The accompanying drawings, which form part of this application, are used to provide a further understanding of the invention. The illustrative embodiments of the invention and their descriptions are used to explain the invention and do not constitute an undue limitation of the invention. In the drawings:

[0044] Figure 1 This is a flowchart illustrating one embodiment of the present invention;

[0045] Figure 2 This is a schematic flowchart of a historical earthquake energy release process in one embodiment of the present invention. Detailed Implementation

[0046] The embodiments of the present invention will be described in detail below. The specific embodiments described herein are merely illustrative of the present invention and are not intended to limit the present invention.

[0047] This invention provides a seismic hazard analysis method that takes into account the energy release of historical earthquakes. Based on regional numerical simulation results, a three-dimensional strain energy grid is obtained through interpolation. The total accumulated strain energy is calculated according to the annual average accumulated strain energy of each grid. A historical earthquake catalog is acquired and processed. The energy released by each earthquake is distributed across the grid according to a specific energy release pattern, ultimately yielding the residual strain energy distribution. Seismic hazard analysis is then performed by identifying potential seismogenic zones and calculating their post-seismic surface PGA (Progressive Gain Area). The method includes the following steps:

[0048] S1. Based on the results of numerical simulation of regional fine geological structures, strain energy density data of faults and strata nodes are obtained. Appropriate interpolation methods are used to interpolate the strain energy density data of faults and strata nodes respectively, constructing a three-dimensional grid of strain energy density for the study area. This grid is then multiplied by the grid volume to obtain a three-dimensional grid of strain energy. The annual average cumulative strain energy E of each grid is calculated based on the three-dimensional grid of strain energy.an ;

[0049] S2. Obtain and process the historical earthquake catalog of the study area, and calculate the energy released by each earthquake and the extent of energy release.

[0050] S3. Determine the number of years of accumulated strain energy in the grid and calculate the accumulated strain energy of each grid.

[0051] S4. Calculate the residual strain energy distribution in the study area after energy release and energy accumulation, and identify potential earthquake-generating areas, their source locations, and magnitudes.

[0052] S5. Calculate the distribution of peak ground acceleration (PGA) on the Earth's surface after each earthquake based on the potential earthquake location and magnitude, and superimpose them to form the potential PGA distribution of the entire study area. Then, classify the regional seismic hazard level according to the PGA.

[0053] Example 1

[0054] The present invention provides a seismic hazard analysis method that takes into account the energy release of historical earthquakes, comprising the following steps:

[0055] S1. Constructing a three-dimensional mesh for strain energy in the study area. This aims to construct the three-dimensional distribution of strain energy in faults and strata using numerical simulation results, and includes the following specific steps:

[0056] S1.1 Strain Energy Density Interpolation. Interpolation methods conforming to the geological 3D model were employed, including cubic spline interpolation, Kriging interpolation, inverse distance weighted interpolation, natural neighborhood interpolation, and Gaussian interpolation. The interpolation results were then used for back-interpolation and compared with the original data. The interpolation method with the smallest error was selected to interpolate the strain energy density data of the strata and fault nodes. A 3D grid of strain energy density was constructed in MATLAB, with a grid size of 1km × 1km × 1km.

[0057] S1.2 Calculation of Annual Average Cumulative Strain Energy. The strain energy density three-dimensional grid obtained in S1.1 is multiplied by the grid volume to calculate the strain energy E of each grid. Then, the overall strain energy distribution is divided by the loading years in the original three-dimensional detailed geological structure numerical simulation to obtain the annual average cumulative strain energy E of each grid. an :

[0058] E an =E / N (1)

[0059] In the formula, E an denoted as the annual average cumulative strain energy of the grid, N is the model loading years, and E is the strain energy of each grid.

[0060] S2. Calculate the energy release and extent of historical earthquakes. Calculate parameters such as the energy released and rupture length of historical earthquakes using empirical formulas to rationally distribute the released seismic energy. This includes the following specific steps:

[0061] S2.1 Obtain and process the historical earthquake catalog. Based on the earthquake catalog downloaded from the USGS official website, filter earthquakes of magnitude 4.0 and above, homogenize the earthquake records, and uniformly convert the magnitudes into surface wave magnitudes Ms;

[0062] S2.2 Calculation of Energy Release from Historical Earthquakes. The energy released by historical earthquakes is calculated according to the Gutenberg-Richter relation in equation (2):

[0063] lg(Es)=1.5×Ms+4.8 (2)

[0064] In the formula, Es is the seismic radiation energy, which can be approximated as the total energy released by the earthquake, and the unit is J. Ms is the surface wave magnitude.

[0065] S2.3 Fault Identification: In MATLAB, the location and extent of faults are identified using the alphaShape function based on fault node data. Principal component analysis is used to determine the best-fit plane of the fault plane, and the strike and dip angle of the fault are calculated. The calculation results are saved, as shown in the following formula:

[0066] (3)

[0067] (4)

[0068] (5)

[0069] In the formula, i represents the i-th fault, j represents the j-th point on a fault, and X j Y j Z j These represent the X, Y, and Z coordinates of the j-th point on a fault, respectively. , and M represents the coordinates of the center point of the i-th fault. j Let P be a matrix containing the coordinates of the j-th point on a fault. i Let C be the coordinate matrix of the center point. i Let λ be the covariance matrix of the i-th fault. i Let v be the eigenvalue of the covariance matrix of the i-th fault. Each fault has 3 eigenvalues. iLet be the eigenvectors corresponding to the eigenvalues ​​of the covariance matrix of the i-th fault, and let be unit vectors. The vector corresponding to the largest eigenvalue represents the fault strike, the eigenvector corresponding to the second largest eigenvalue represents the fault dip direction, and the eigenvector corresponding to the smallest eigenvalue represents the direction perpendicular to the fault plane. The strike and dip angle of each fault can be calculated from these three vectors.

[0070] S2.4 Calculation of rupture length, rupture width, and rupture volume: The rupture length, rupture width, average slip, and rupture volume after historical earthquakes were calculated according to the empirical formulas (2), (3), (4), and (5) of Wells & Coppersmith (1994).

[0071] log(RLD)=0.59×Mw-2.44 (6)

[0072] log(RW) = 0.32 × Mw - 1.01 (7)

[0073] log(AD) = 0.69 × Mw - 4.8 (8)

[0074] V = RLD × RW × AD (9)

[0075] In the formula, RLD, RW, AD and V represent the rupture length, rupture width, average slip and rupture volume, respectively, in kilometers, and Mw represents the moment magnitude.

[0076] S2.5 Historical earthquake energy release and distribution: First, determine whether the hypocenter is on a fault. If the earthquake source is on a fault, a cube is constructed centered on the source point, with the rupture length as the length along the fault strike, the rupture width as the width along the dip direction, and the average slip as the height in the vertical direction. This ensures that all grid energy within the cube is released, while guaranteeing that the released energy within the cube does not exceed the total energy released by the earthquake. If it exceeds the total energy released by the earthquake, the cube is scaled down proportionally until the total energy of the internal grid is less than the energy released by the earthquake. If it does not exceed the total energy released by the earthquake, the remaining energy is distributed to the grid outside the cube according to a three-dimensional Gaussian distribution, with the length, width, and height of the cube as the x, y, and z axes, respectively, using a coordinate system centered on the source point. The standard deviations of the three-dimensional Gaussian distribution in the x, y, and z directions are RLD / 4, RLD / 6, and RLD / 8, respectively. If the source point is not on a fault, a coordinate system is directly constructed centered on the source point, with the mutually perpendicular horizontal directions as the x and y axes and the vertical direction as the z axis, with the standard deviations of the x, y, and z directions being RLD / 4, RLD / 4, and RLD / 6, respectively.

[0077] S3. Calculate the cumulative strain energy for each grid. The strain energy obtained from the numerical simulation is the result of many years of loading, and there has been energy accumulation over many years. After considering the energy release from historical earthquakes, it is also necessary to consider the re-accumulation of energy and calculate the energy accumulation value for each grid.

[0078] Step S3 includes the following specific steps:

[0079] S3.1 Calculation of the number of years of grid strain energy accumulation. Determine the initial year in which the grid releases energy due to the earthquake. Subtract the initial year from the target year for earthquake hazard assessment to obtain the number of years of grid strain energy accumulation.

[0080] S3.2 Calculation of accumulated strain energy in the mesh. The accumulated strain energy in the mesh can be calculated using formula (10):

[0081] E ac =E an ×N ac (10)

[0082] In the formula, E ac E represents the cumulative strain energy of the mesh. an The annual average cumulative strain energy of the grid, in J and N. ac Years for accumulating grid strain energy.

[0083] S4. Calculate the residual strain energy and identify potential seismic zones. Based on the above steps, the regional residual strain energy distribution can be calculated, and the seed region growth method is used to identify potential seismic zones.

[0084] Step S4 includes the following specific steps:

[0085] S4.1 Calculation of Regional Residual Strain Energy. The residual strain energy distribution is obtained by subtracting the strain energy released by historical earthquakes from the interpolated 3D strain energy grid and then adding the inter-seismic accumulated strain energy. Grids unaffected by earthquakes are not further accumulated for strain energy.

[0086] E r =E t -E s +E ac (11)

[0087] In the formula, E r E represents the residual strain energy. t E is the total strain energy obtained by interpolation of the numerical simulation results. s E represents the energy released by historical earthquakes. ac This represents the cumulative strain energy of the grid, expressed in J.

[0088] S4.2 Identification of Potential Seismogenic Zones. A strain energy threshold is set for each mesh cell. A seed region growth method is used to identify mesh cells exceeding the threshold and form potential seismogenic zones. The total strain energy within each zone is calculated, while ensuring that the total energy does not exceed the energy released by the maximum magnitude earthquake.

[0089] S4.3 Determination of the source location in the potential seismogenic zone. First, determine whether a fault exists within the potential seismogenic zone. If a fault exists, take the center point of the grid with the highest strain energy on the fault as the source point; if no fault exists, directly take the center point of the grid with the highest strain energy in the region as the source point.

[0090] S4.4 Calculation of magnitude of potential seismogenic zone. Calculate the magnitude of possible earthquakes in the potential seismogenic zone according to formula (2).

[0091] S5. Seismic Hazard Analysis. When conducting seismic hazard analysis on the study area, peak ground acceleration (PGA) is considered to quantify seismic hazard. Taking Pakistan as an example, step S5 includes the following steps:

[0092] S5.1 Calculation of Peak Ground Acceleration (PGA). Based on the study area, appropriate seismic attenuation equations (11a) and (11b) are selected to calculate the PGA distribution formed on the Earth's surface after an earthquake in each potential seismogenic zone. The PGA values ​​of each potential earthquake are then superimposed to form the overall PGA distribution:

[0093] a) When Mw ≤ 6.75:

[0094] (11a)

[0095] b) When Mw > 6.75:

[0096] (11b)

[0097] In the formula, PGA is the peak ground acceleration in g, Mw is the moment magnitude, and R... Hyp The distance from the epicenter to the station. R epi H is the epicentral distance, H is the focal depth, and the unit is km; F N F R These are the variables for normal faults and reverse faults, respectively. When the fault is a normal fault, F... N When the fault is a reverse fault, F is 1. R F is 1 when the fault type is other types. N F R All are 0;

[0098] S5.2 Seismic Hazard Analysis. The comprehensive PGA distribution of the study area is divided into several levels, thereby classifying the overall seismic hazard level of the study area.

[0099] This invention presents a numerical simulation method for seismic hazard analysis based on detailed regional geological structures, taking into account historical earthquake energy release. The method interpolates regional strain energy density using nodal data to construct a three-dimensional strain energy density grid. Multiplying this grid by the grid volume yields a three-dimensional strain energy grid. Using a distribution model consistent with historical earthquake energy release, the released energy is allocated to the corresponding grids. Finally, considering the accumulated strain energy in each grid, the remaining strain energy distribution of the regional grid is calculated, and potential seismogenic zones are identified. The post-earthquake ground motion (PGA) of the potential seismogenic zones is calculated using the seismic motion attenuation equation, ultimately achieving an accurate assessment of regional seismic hazard. This method can construct a regional seismic hazard level distribution map, which can then be used for seismic risk analysis.

[0100] The purpose of this invention is to leverage existing numerical simulation results of detailed three-dimensional geological structures in a region. This simulation integrates data on regional stratigraphy and faults to construct a three-dimensional finite element model. Based on the crust 1.0 global crust model, lithological parameters for each layer are determined. The model employs nonlinear frictional contact between faults and strata, and bonded contact between strata. Using surrounding GNSS velocity field data as constraints, the background stress distribution characteristics are simulated and calculated. This invention extracts strain energy density data from each grid node of the numerical simulation results, considers the energy release from historical earthquakes in the region, and predicts the location and magnitude of future earthquakes in the region. This provides new ideas and methods for earthquake prediction and a scientific basis for regional earthquake prevention and disaster reduction.

[0101] The above description, in conjunction with specific preferred embodiments, provides a further detailed explanation of the present invention. It should not be construed that the specific implementation of the present invention is limited to these descriptions. For those skilled in the art, various simple deductions and substitutions can be made without departing from the inventive concept, and all such modifications and substitutions should be considered within the scope of protection of the present invention.

Claims

1. A seismic hazard analysis method that takes into account the energy release of historical earthquakes, characterized in that, Includes the following steps: S1. Based on the results of numerical simulation of regional fine geological structure, obtain the strain energy density data of each node of fault and stratum. Use the interpolation method to interpolate the strain energy density data of each node of fault and stratum respectively, construct a three-dimensional grid of strain energy density in the study area, and then multiply by the grid volume to obtain the three-dimensional grid of strain energy. Calculate the annual average cumulative amount of strain energy of each grid according to the three-dimensional grid of strain energy. S2. Obtain and process the historical earthquake catalog of the study area, and calculate the energy released by each earthquake and the extent of energy release; specifically including the following steps: S2.1 Historical earthquake catalog processing; S2.2 Calculation of energy release from historical earthquakes; S2.3, Fault identification; S2.4 Calculation of fracture length, fracture width and fracture volume; S2.5 Distribution of historical earthquake energy release, specifically as follows: First, determine if the hypocenter is on a fault. If it is, construct a cube centered on the hypocenter, along the fault strike with the rupture length as the length, along the dip with the rupture width as the width, and along the vertical direction with the average slip as the height. This ensures that all grid energy within the cube is released, while guaranteeing that the released energy does not exceed the total energy released by the earthquake. If the released energy exceeds the total energy released by the earthquake, the cube is scaled down proportionally until the total energy of the internal grid is less than the energy released by the earthquake. If the released energy does not exceed the total energy released by the earthquake, the remaining energy is distributed to the grid outside the cube using a coordinate system centered on the hypocenter and with the cube's length, width, and height as the x, y, and z axes, respectively. The standard deviations of the 3D Gaussian distribution in the x, y, and z directions are set to RLD / 4, RLD / 6, and RLD / 8, respectively. If the hypocenter is not on the fault, a coordinate system is established with the hypocenter as the center, the mutually perpendicular horizontal directions as the x and y axes, and the vertical direction as the z axis. The standard deviations of the x, y, and z directions are taken as RLD / 4, RLD / 4, and RLD / 6, respectively, where RLD represents the rupture length. S3. Determine the number of years of accumulated strain energy in the grid and calculate the accumulated strain energy of each grid. S4. Calculate the three-dimensional grid distribution of residual strain energy in the study area after energy release and energy accumulation, and identify potential seismogenic areas, their source locations, and magnitudes. S5. Calculate the peak ground acceleration distribution on the Earth's surface after each earthquake based on the location of the earthquake source and the magnitude of the earthquake in the potential seismogenic area, and superimpose them to form the potential peak ground acceleration distribution on the surface of the entire study area. Classify the regional seismic hazard level according to the peak ground acceleration.

2. The seismic hazard analysis method according to claim 1, characterized in that, In step S1, the strain energy density data of faults and stratigraphic nodes are interpolated using an interpolation method to construct a three-dimensional grid of strain energy density for the study area. This grid is then multiplied by the grid volume to obtain the three-dimensional grid of strain energy. Specifically: Cubic spline interpolation, Kriging interpolation, inverse distance weighted interpolation, natural neighborhood interpolation, and Gaussian interpolation were selected for interpolation. The interpolation results were then used for inverse interpolation and compared with the original data. The interpolation method with the smallest error was selected to interpolate the strain energy density data of the strata and fault nodes, and a three-dimensional strain energy density grid was constructed. Multiplying the grid volume, a three-dimensional strain energy grid was obtained.

3. The seismic hazard analysis method according to claim 1, characterized in that, In step S1, the annual average cumulative strain energy of each grid is calculated, specifically as follows: The strain energy of each grid is calculated by multiplying the obtained three-dimensional strain energy density grid by the grid volume. Then, the overall strain energy distribution is divided by the loading years in the original three-dimensional fine geological structure numerical simulation to obtain the annual average cumulative strain energy of each grid.

4. The seismic hazard analysis method according to claim 1, characterized in that, S2.1 Historical earthquake catalog processing, specifically: homogenize the historical earthquake catalog and uniformly convert different magnitude types into surface wave magnitude Ms; S2.2 Calculation of energy release from historical earthquakes: Specifically, the energy released by historical earthquakes is calculated based on the Gutenberg-Richter relation. S2.3 Fault identification, specifically: based on the three-dimensional coordinates of fault nodes, identify the extent, strike, and dip of faults; use principal component analysis to calculate the strike and dip of each fault; and save the calculation results. S2.4 Calculation of rupture length, rupture width and rupture volume: Specifically, the rupture length, rupture width, average slip and rupture volume after historical earthquakes are calculated according to the empirical formula of Wells & Coppersmith, 1994.

5. The seismic hazard analysis method according to claim 1, characterized in that, Step S3 specifically includes the following steps: S3.1 Calculation of the number of years of grid strain energy accumulation: Determine the initial year in which the grid releases energy due to the earthquake. The number of years of grid strain energy accumulation can be obtained by subtracting the initial year from the target year of the earthquake hazard assessment. S3.2 Calculation of mesh strain energy accumulation: E ac =E an ×N ac ; In the formula, E ac E represents the cumulative strain energy of the mesh. an The annual average cumulative strain energy of the grid, in J and N. ac Years for accumulating grid strain energy.

6. The seismic hazard analysis method according to claim 1, characterized in that, Step S4 specifically includes the following steps: S4.1 Calculation of regional residual strain energy: Subtract the strain energy released by historical earthquakes from the three-dimensional strain energy grid obtained by interpolation, and add the inter-seismic accumulated strain energy to obtain the distribution of residual strain energy. The grids that are not affected by earthquakes will no longer accumulate strain energy. S4.2 Identification of potential seismogenic zones: Set the strain energy threshold for grid cells, use the seed region growth method to identify grid cells with values ​​above the threshold and form potential seismogenic zones, and calculate the total strain energy within the zone, while ensuring that the total energy does not exceed the energy that can be released by the maximum magnitude. S4.3 Determination of the source location in the potential seismogenic zone: First, determine whether there is a fault in the potential seismogenic zone; if there is a fault, take the center point of the grid with the highest strain energy on the fault as the source point; if there is no fault, directly take the center point of the grid with the highest strain energy in the region as the source point. S4.4 Calculation of magnitude in potential seismogenic zones: Calculate the possible magnitude of earthquakes in potential seismogenic zones based on the Gutenberg-Richter relation.

7. The seismic hazard analysis method according to claim 1, characterized in that, Step S5 specifically includes the following steps: S5.1 Calculation of Peak Ground Acceleration: Based on the study area, select an appropriate seismic attenuation equation to calculate the surface peak ground acceleration distribution after an earthquake occurs in each potential seismogenic zone, and superimpose the peak ground acceleration values ​​of each potential earthquake to form the total peak ground acceleration distribution. S5.2 Seismic Hazard Analysis: The peak ground acceleration distribution of the study area is divided into several levels, thereby classifying the overall seismic hazard level of the study area.

Images