A seismic hazard analysis method based on regional fine geological structure

Abstract

This invention discloses a seismic hazard analysis method based on detailed regional geological structures, including: constructing a detailed three-dimensional geoscientific model; optimizing boundary conditions based on GNSS velocity field data; finite element numerical simulation; and seismic hazard analysis. By constructing a detailed three-dimensional geoscientific model and utilizing the finite element numerical simulation method, combined with historical earthquake data and geological structural characteristics, the uneven strain energy distribution caused by plate movement over a certain period is calculated. The simulated magnitude risk values ​​are weighted by depth and magnitude to reflect the actual seismic damage effects. Finally, the weighted strain energy data of earthquakes at different depths and magnitudes within the study area are summed to obtain the overall seismic hazard distribution of the surface area. This invention overcomes the limitations of existing earthquake observation data, allowing for a broader perspective and time scale to examine seismic activity, and enabling a deeper understanding of the periodicity, frequency, and potential disaster risks of earthquakes.

Description

Technical Field

[0001] This invention belongs to the field of three-dimensional geoscience simulation research. Based on the numerical simulation results of regional fine geological structures, it provides a more comprehensive method for regional seismic hazard analysis. Background Technology

[0002] For a long time, scholars have tried many methods to quantify earthquake risk, such as peak ground acceleration, peak ground velocity, or response spectrum acceleration. Earthquake hazard analysis is mainly divided into deterministic analysis and probabilistic analysis. However, these methods are all based on statistical analysis of recorded earthquake events. In reality, geomagnetic and seismic observatories were only established globally in the late 19th and early 20th centuries, and modern seismic instruments were installed for earthquake observation. Although records from this period provide valuable data, some strong earthquakes may have cycles of thousands of years, and modern seismic observation data alone are insufficient to fully capture these long-term seismic activity patterns. Summary of the Invention

[0003] The purpose of this invention is to provide a seismic hazard analysis method based on regional fine geological structure, so as to break through the limitations of existing seismic observation data, examine seismic activity from a broader perspective and time scale, and provide a scientific basis for formulating effective disaster prevention and mitigation strategies.

[0004] To achieve the above objectives, this invention provides a seismic hazard analysis method based on regional fine geological structures, comprising the following steps:

[0005] S1. Based on the elevation data of the lower boundary of the strata and historical earthquake information of the study area, a three-dimensional geoscientific model containing fault and stratigraphic geometry is constructed, and lithological parameters are assigned to each rock layer in the model.

[0006] S2. Optimize the boundary conditions of the model based on GNSS velocity field data: Collect the measured annual average velocity at GPS stations in and around the study area, use the difference method to obtain the boundary velocity of the model, and invert and verify the obtained boundary velocity;

[0007] S3. The strain energy accumulation of each element in the model at a certain time scale is calculated using the finite element numerical simulation method, and the strain energy density E at each point in the study area is calculated based on the strain energy accumulation of each element. density ;

[0008] S4. Seismic Hazard Analysis: Based on the strain energy at each point within the study area, the magnitude risk value for each point is calculated. The obtained magnitude risk values ​​are then weighted by the focal depth and magnitude to determine the seismic risk weight coefficient for each combination of focal depth and magnitude. Based on the obtained seismic risk weight coefficient, the weighted strain energy of earthquakes at different depths and magnitudes within the study area is summed. Finally, based on the classification of the strain energy release capacity level of the study area, the overall seismic hazard distribution of the study area is obtained.

[0009] When conducting seismic hazard analysis for a target area, step S1 specifically includes the following steps:

[0010] S1.1 Stratigraphic Division: Using the lower boundary elevation data of strata in the global crustal model database, the sedimentary layer, upper crust, middle crust, lower crust and upper mantle of the study area are divided into undulations through grid surface fitting technology;

[0011] S1.2 Fault Structure: Based on historical earthquake information, active faults in the study area are selected, and geometric models of the faults are performed to ensure the rationality of the fault strike and dip angle.

[0012] S1.3, Assignment of lithological parameters: Considering the short-term elastic and long-term rheological properties of plate tectonic stress evolution, the constitutive relation of the viscoelastic model is used to simulate the stress evolution at a certain time scale, and the lithological parameters of each stratum are calculated.

[0013] S1.4 Mesh Generation: The 3D geoscientific model is meshed. The size of the mesh unit depends on the geometry of the model. The unit quality must be above 0.8 to be considered an effective calculation unit.

[0014] When conducting seismic hazard analysis for a target area, the lithological parameters mentioned in step S1.3 include density, Young's modulus, Poisson's ratio, and viscosity coefficient; the constitutive equation of the viscoelastic model is:

[0015] K(t)=K (1)

[0016]

[0017] Where K(t) and G(t) are the bulk modulus and shear modulus of the viscoelastic model, respectively, η is the viscosity coefficient, E is Young's modulus, ν is Poisson's ratio, K is the elastic bulk modulus, and G is the elastic shear modulus.

[0018] When conducting seismic hazard analysis for a target area, the Young's modulus E and Poisson's ratio ν are calculated based on the P-wave and S-wave velocity and density parameters of various strata provided by the global crustal model database.

[0019]

[0020] Among them, V s For the P-wave velocity of each stratum, V p Let ρ be the S-wave velocity of each stratum and ρ be the density parameter of each stratum.

[0021] When conducting seismic hazard analysis for a target area, step S2 specifically includes the following steps:

[0022] S2.1, Velocity Field Difference: Collect annual average velocity data from GPS stations within and near the study area, and use MATLAB software to perform difference processing to obtain the boundary velocity of the model;

[0023] S2.2 Reasonableness test: Using the fitted boundary velocity as a known condition, calculate the annual average velocity at GPS observation stations within the study area and compare it with the true value collected in step S2.1 to test the reasonableness of the constraint conditions.

[0024] When conducting seismic hazard analysis for a target area, step S3 specifically includes the following steps:

[0025] S3.1 Preloading process: Load the model to achieve convergence and bring it to equilibrium.

[0026] S3.2 Structural stress simulation: The displacement state of the preloading result is used as the initial state to determine the time scale of structural loading; the annual average velocity is multiplied by the time scale and applied to the nodes around the model as nodal displacement; at the same time, a gravity field is applied, the upper surface of the model is treated as a free boundary, and the bottom boundary is constrained by a normal fixed and horizontal free method.

[0027] S3.3 Post-solution processing: Calculate the volume size V of each element and the cumulative strain energy Es of each element at the end of the solution, and calculate the strain energy density E at each point in the study area. density :

[0028] E density =Es / V (7).

[0029] When conducting seismic hazard analysis for a target area, step S4 specifically includes the following steps:

[0030] S4.1 Calculation of Strain Energy: The formula for calculating the cumulative strain energy Es(H) at depth H within the study area is as follows:

[0031] Es(H)=E density (H)·V0 (8)

[0032] Among them, E density (H) represents the strain energy density at depth H within the study area, and V0 is the unit volume;

[0033] S4.2 Calculate the magnitude: Assume that the accumulated strain energy is completely released, and convert the accumulated strain energy into the corresponding magnitude risk according to the Gutenberg-Richter relation in the following formula (9);

[0034] lg(E)=4.8+1.5·Mw (9)

[0035] The magnitude Mw(H) at depth H is:

[0036]

[0037] S4.3 Calculate the hazard index: Use the weighted coefficient calculation formula (11) to determine the magnitude risk weight coefficient W(H, Mw(H)) under the combination of each source depth H and magnitude Mw;

[0038]

[0039] S4.4 Calculate the weighted cumulative strain energy: Perform a weighted integral on the strain energy at each depth H to obtain the weighted cumulative strain energy per unit area from the upper surface to the lower surface;

[0040]

[0041] Where Hmin is the surface elevation and Hmax is the depth at the bottom of the lower crust;

[0042] S4.5 Seismic Hazard Analysis: The strain energy release capacity of the study area is divided into several levels, thereby classifying the overall seismic hazard level of the study area.

[0043] When conducting seismic hazard analysis for a target area, step S4, depth and magnitude risk weighting, needs to consider that the destructive effect of an earthquake is significantly influenced by both its focal depth and magnitude. The simulated magnitude risk values ​​are then weighted by depth and magnitude. Based on a specific weighting coefficient calculation formula, the risk weight coefficients for each combination of focal depth and magnitude are determined.

[0044] When conducting seismic hazard analysis for a target area, step S5, seismic hazard analysis, requires summing the weighted strain energy data of earthquakes at different depths and magnitudes within the study area based on the obtained seismic risk weight coefficients. By classifying the strain energy release capacity of the study area, the overall seismic hazard distribution of the surface area is obtained.

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

[0046] This invention presents a numerical simulation-based method for regional seismic hazard analysis, aiming to comprehensively assess the seismic hazard of a specific region. By constructing a refined three-dimensional geoscientific model and utilizing the finite element method (FEM) numerical simulation, combined with historical earthquake data and geological structural characteristics, the method calculates the uneven strain energy distribution caused by plate movement over a certain period. The simulated magnitude risk values ​​are weighted by depth and magnitude to reflect the actual seismic damage effects. Finally, the weighted strain energy data of earthquakes at different depths and magnitudes within the study area are summed to obtain the overall seismic hazard distribution of the surface region. This invention overcomes the limitations of existing earthquake observation data, examining seismic activity from a broader perspective and time scale, and enabling a deeper understanding of the periodicity, frequency, and potential disaster risks of earthquakes.

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

[0048] The accompanying drawings are provided to further illustrate embodiments of the present invention and form part of the specification. They are used together with the following detailed description to explain the embodiments of the present invention, but do not constitute a limitation thereof. In the drawings:

[0049] Figure 1 This is a flowchart illustrating an embodiment of the present invention;

[0050] Figure 2 This is a flowchart illustrating the finite element numerical simulation process in one embodiment of the present invention. Detailed Implementation

[0051] The present invention will now be described in detail with reference to the embodiments shown in the accompanying drawings. However, it should be noted that these embodiments are not intended to limit the present invention. Equivalent transformations or substitutions in function, method, or structure made by those skilled in the art based on these embodiments are all within the scope of protection of the present invention.

[0052] In one specific embodiment of the present invention, such as Figure 1 As shown, a seismic hazard analysis method based on strain energy accumulation from regional numerical simulation and weighted calculation is used. A refined three-dimensional finite element geoscientific model is established to highly reproduce the complexity of underground geological structures and the detailed characteristics of interplate interactions. Based on this, annual average velocity field data from GPS is utilized to simulate the long-term trends and dynamic changes of crustal movement, thus providing crucial evidence for regional seismic hazard analysis. The method specifically includes the following steps:

[0053] S1: Constructing a Refined 3D Geoscientific Model. This step aims to construct a refined 3D geoscientific model that includes fault and stratigraphic geometry, and to assign lithological parameters to each rock layer. The model construction in this step includes the following specific steps:

[0054] S1.1: Stratigraphic Division. Using the lower boundary elevation data of the crustal layers in the crust 1.0 global crustal model database, and through grid surface fitting technology, the undulations of the sedimentary layers, upper crust, middle crust, lower crust, and upper mantle in the study area are delineated. Given that earthquakes within the crust cause the most severe damage to the surface, the model depth can be determined based on the crustal thickness of the study area, and the lower boundary of the lower crust can be appropriately extended downwards as the upper mantle layer. The accumulation of strain energy in the upper mantle is not considered within the scope of this invention.

[0055] S1.2: Fault Structure. Based on historical earthquake information, active faults within the study area are selected. Fault data can be obtained from the fault database provided by the Seismic Active Fault Exploration Data Center. For faults lacking dip information, simplified modeling can be performed according to fault type. The dip angle can be set as 30° for thrust faults, 60° for normal and reverse faults, and 90° for strike-slip faults. When modeling fault geometry, the rationality of the fault strike and dip angle must be ensured. For faults with varying dip angles, the curved fault strike can be divided into several segments as horizontal plane baselines. Dip lines on the longitudinal profile are added to control the dip angle of each segment. In areas with significant curvature changes, more segments should be added to ensure that the dip angle does not distort with abrupt changes in strike.

[0056] S1.3: Lithological parameters. Considering the short-term elastic and long-term rheological properties of plate tectonic stress evolution, it is more reasonable to use viscoelastic constitutive relations to simulate the stress distribution formed over a long timescale. The constitutive equation of the viscoelastic model is:

[0057] K(t)=K (1)

[0058]

[0059] Where K(t) and G(t) are the bulk modulus and shear modulus of the viscoelastic model, respectively, K is the elastic bulk modulus, G is the elastic shear modulus, E is Young's modulus, ν is Poisson's ratio, and η is the viscosity coefficient. K and G can be obtained by converting Young's modulus E and Poisson's ratio ν of the linear elastic medium. The P-wave velocity V of various strata is based on the crust 1.0 global crust model database. p S-wave velocity V s With the density parameter ρ, Young's modulus E and Poisson's ratio ν can be calculated.

[0060]

[0061]

[0062] The viscosity coefficient η reflects the resistance of the formation material to flow or deformation, and generally acts on the lower crust and upper mantle. It needs to be flexibly adjusted and optimized according to specific geological conditions and simulation requirements. In addition, the contact with the fault in the model adopts nonlinear frictional contact, and the friction coefficient on the fault plane is taken as the commonly used value of 0.4, with fixed contact between different strata.

[0063] S1.4: Meshing. The model is meshed. The mesh element size depends on the model geometry, and the element quality must be at least 0.8 to be considered an effective computational element. Specifically, while ensuring element quality, the number of elements can be appropriately reduced to improve computational efficiency.

[0064] S2: Optimizing boundary conditions based on GNSS velocity field data. The loading conditions of the model are obtained using the interpolation method based on the measured annual average velocity of GPS stations in and around the study area.

[0065] In a specific implementation process, step S2 includes the following specific steps:

[0066] S2.1: Velocity Field Difference. To best fit the plate tectonics trend of the study area, annual average velocity data from GPS stations within and near the study area were first collected. The griddata function in MATLAB was then used for difference processing to obtain the boundary velocities of the model. Each boundary can be appropriately divided into several parts; for areas with large velocity variations or low regularity, the segmented difference needs to be refined.

[0067] S2.2: Reasonableness Check. Using the boundary velocity obtained from the fitting in step S2.1 as a known condition, the annual average velocity at GPS observation stations within the study area is calculated in reverse and compared with the collected true values ​​to check the reasonableness of the constraints.

[0068] S3: Finite Element Numerical Simulation. This method calculates the strain energy accumulation of each element over a specific time scale using finite element numerical simulation. In practice, when the model is built (before simulation begins), the stress level of each element is zero. At this point, the model is in an unstable state. Preloading is required to bring the model to equilibrium. The equilibrium model is then used as the initial state to calculate the stress changes during the long-term evolution of the model.

[0069] Figure 2 As shown, in a specific implementation process, step S3 includes the following specific steps:

[0070] S3.1: Preloading process. The model is horizontally fixed on all four sides, a gravity field is applied, the bottom surface is fixed in the normal direction, and it is horizontally free. The force convergence criterion adopts the default value of 0.5%, that is, when the ratio of unbalanced force to external force is less than 0.5%, the solution is considered convergent.

[0071] S3.2: Structural Stress Simulation. Using the displacement state from the pre-loading results as the initial state, the time scale of the structural loading is determined. The annual average velocity is multiplied by the time scale and applied as nodal displacements to the nodes around the model. Simultaneously, a gravity field is applied. The upper surface of the model is treated as a free boundary, while the bottom boundary is constrained by a fixed normal direction and horizontal freedom.

[0072] S3.3: Post-processing. Calculate the volume V of each element and the cumulative strain energy Es of each element at the end of time. Divide the cumulative strain energy of each element by its volume to obtain the strain energy density E at each point within the study area. density :

[0073] E denstity =Es / V (7)

[0074] S4: Seismic Hazard Analysis. When conducting seismic hazard analysis for the study area, given that the destructive effects of an earthquake are significantly influenced by both its focal depth and magnitude, it is necessary to consider the hazard indices corresponding to strain energies at different depths and magnitudes for weighted processing. Since the energy release range of each earthquake cannot be quantitatively calculated, in practice, the strain energy per unit volume V0 (i.e., a 1×1×1m cube) is used as the object; that is, the strain energy density at any point multiplied by the unit volume V0 is taken as the strain energy at that point.

[0075] In a specific implementation process, step S4 includes the following steps:

[0076] S4.1: Calculate the strain energy. The cumulative strain energy Es(H) at depth H (km) can be expressed as the strain energy density E density The product of the volume and the unit volume V0 is:

[0077] Es(H)=E density (H)·V0 (8)

[0078] Among them, E density (H) represents the strain energy density at depth H within the study area;

[0079] S4.2: Calculate the magnitude. Assume that the accumulated strain energy is fully released and convert the strain energy into the corresponding magnitude risk according to the Gutenberg-Richter relation (Equation (9)).

[0080] lg(E)=4.8+1.5·Mw (9)

[0081] The magnitude Mw(H) at depth H(km) is:

[0082]

[0083] S4.3: Calculate the risk index. Based on the weighted coefficient calculation formula (11), determine the risk weight coefficient W(H, Mw(H)) under the combination of each source depth H and magnitude Mw.

[0084]

[0085] In fact, the weighting formula used in this embodiment is derived from the PGA attenuation prediction model and is only a preferred calculation method.

[0086] S4.4: Calculate the weighted cumulative strain energy. To calculate the weighted cumulative strain energy per unit area from the top surface to the bottom surface, it is necessary to perform a weighted integral over the strain energy at each depth H:

[0087]

[0088] Where Hmin is the surface elevation and Hmax is the depth at the bottom of the lower crust.

[0089] S4.5: Seismic Hazard Analysis. The strain energy release capacity of the study area can be divided into several levels, thereby classifying the seismic hazard level.

[0090] This invention discloses a seismic hazard analysis method based on detailed regional geological structures. By constructing a refined three-dimensional geoscientific model, combining it with a global crustal model database and seismic fault data, and optimizing boundary conditions using GNSS velocity field data, the method then uses finite element numerical simulation to calculate the strain energy accumulation of each element at a specific time scale, ultimately achieving an accurate assessment of seismic hazard. This enables not only the creation of seismic risk zoning maps but also a deeper understanding of the periodicity, frequency, and potential disaster risks of earthquakes.

[0091] The above description is merely a preferred embodiment of the present invention and is not intended to limit the invention. Various modifications and variations can be made to the present invention by those skilled in the art. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the scope of protection of the present invention.

Claims

1. A seismic hazard analysis method based on regional fine geological structures, characterized in that, Includes the following steps: S1. Based on the elevation data of the lower boundary of the strata and historical earthquake information of the study area, a three-dimensional geoscientific model containing fault and stratigraphic geometry is constructed, and lithological parameters are assigned to each rock layer in the model. S2. Optimize the boundary conditions of the model based on GNSS velocity field data: Collect the measured annual average velocity at GPS stations in and around the study area, use the difference method to obtain the boundary velocity of the model, and invert and verify the obtained boundary velocity; S3. The strain energy accumulation of each element in the model at a certain time scale is calculated using the finite element numerical simulation method, and the strain energy density at each point in the study area is calculated based on the strain energy accumulation of each element. ; Specifically, the following steps are included: S3.1 Preloading process: Load the model to achieve convergence and bring it to equilibrium. S3.2 Structural stress simulation: The displacement state of the preloading result is used as the initial state to determine the time scale of structural loading; the annual average velocity is multiplied by the time scale and applied to the nodes around the model as nodal displacement; at the same time, a gravity field is applied, the upper surface of the model is treated as a free boundary, and the bottom boundary is constrained by a normal fixed and horizontal free method. S3.3 Post-solution processing: Calculate the volume of each element. and the cumulative strain energy of each element at the end of the solution. And calculate the strain energy density at each point within the study area. : (7) S4. Seismic Hazard Analysis: Based on the strain energy at each point within the study area, the corresponding magnitude risk value is calculated. The obtained magnitude risk values ​​are then weighted by focal depth and magnitude to determine the seismic risk weight coefficient for each combination of focal depth and magnitude. Based on the obtained seismic risk weight coefficient, the weighted strain energy of earthquakes at different depths and magnitudes within the study area is summed. Finally, based on the classification of the strain energy release capacity level of the study area, the overall seismic hazard distribution of the study area is determined. This includes the following steps: S4.1 Calculate strain energy: depth within the study area Cumulative strain energy at the point The calculation formula is: (8) in, For the depth within the study area Strain energy density at the point, per unit volume; S4.2 Calculate the magnitude: Assume that the accumulated strain energy is completely released, and convert the accumulated strain energy into the corresponding magnitude risk according to the Gutenberg-Richter relation in the following formula (9); (9) Then depth The corresponding magnitude for: (10) S4.3 Calculation of Hazard Index: The depth of each seismic source is determined using the weighted coefficient calculation formula (11). With magnitude Combined magnitude risk weighting coefficient ; (11) S4.4 Calculate the weighted cumulative strain energy: for each depth The strain energy is weighted and integrated to obtain the cumulative weighted strain energy per unit area from the top surface to the bottom surface; (12) in, It is the surface elevation. It is the depth at the bottom of the lower crust; S4.5 Seismic Hazard Analysis: The strain energy release capacity of the study area is divided into several levels, thereby classifying the overall seismic hazard level of the study area.

2. The method according to claim 1, characterized in that, Step S1 specifically includes the following steps: S1.1 Stratigraphic Division: Using the lower boundary elevation data of strata in the global crustal model database, the sedimentary layer, upper crust, middle crust, lower crust and upper mantle of the study area are divided into undulations through grid surface fitting technology; S1.2 Fault Structure: Based on historical earthquake information, active faults in the study area are selected, and geometric models of the faults are performed to ensure the rationality of the fault strike and dip angle. S1.3, Assignment of lithological parameters: Considering the short-term elastic and long-term rheological properties of plate tectonic stress evolution, the constitutive relation of the viscoelastic model is used to simulate the stress evolution at a certain time scale, and the lithological parameters of each stratum are calculated. S1.4 Mesh Generation: The 3D geoscientific model is meshed. The size of the mesh unit depends on the geometry of the model. The unit quality must be above 0.8 to be considered an effective calculation unit.

3. The method according to claim 2, characterized in that, The lithological parameters mentioned in step S1.3 include density, Young's modulus, Poisson's ratio, and viscosity coefficient; the constitutive equation of the viscoelastic model is: (1) (2) (3) (4) in, and These are the bulk modulus and shear modulus of the viscoelastic model, respectively. The viscosity coefficient, For Young's modulus, Poisson's ratio, For elastic bulk modulus, It is the elastic shear modulus.

4. The method according to claim 3, characterized in that, The Young's modulus is calculated based on the P-wave and S-wave velocity and density parameters of various strata provided by the global crustal model database. Compared to Poisson : (5) (6) in, For the P-wave velocity of each stratum, For the S-wave velocity of each stratum, These are the density parameters for each stratum.

5. The method according to claim 1, characterized in that, Step S2 specifically includes the following steps: S2.1, Velocity Field Difference: Collect annual average velocity data from GPS stations within and near the study area, and use MATLAB software to perform difference processing to obtain the boundary velocity of the model; S2.2 Reasonableness test: Using the fitted boundary velocity as a known condition, calculate the annual average velocity at GPS observation stations within the study area and compare it with the true value collected in step S2.1 to test the reasonableness of the constraint conditions.

Images