A method and system for evaluating seismic risk

Abstract

The application discloses a seismic risk evaluation method and system, relates to the technical field of seismic numerical simulation, and receives a three-dimensional geometric model of a gridded fault and sets initial conditions; a simulation operation equation is constructed to solve a shear stress loading rate, a normal stress loading rate and a state variable, and a common differential equation system is constructed; based on the three-dimensional geometric model and the initial conditions, the common differential equation system is used to realize accelerated simulation operation of the shear stress loading rate, the normal stress loading rate and the state variable, wherein the method of the accelerated simulation operation comprises: using a method combining an H-matrix and MPI to accelerate simulation operation, comprising obtaining a layered matrix based on a Green function under a half-space condition in combination with an adaptive cross approximation algorithm, comprising a normal stiffness H-matrix and a shear stiffness H-matrix, and using MPI parallel operation matrix-vector multiplication. The application introduces a hierarchical matrix and fuses an MPI parallel computing framework, so that the simulation period is greatly shortened.

Description

TECHNICAL FIELD

[0001] The present application relates to the field of seismic simulation, in particular to a seismic risk assessment method and system. BACKGROUND

[0002] Earthquake is one of the most significant natural phenomena in the process of energy release in the earth interior, and its generation mechanism is closely related to fault activity. The sliding process of the fault includes both sudden earthquake rupture and slow aseismic slip, and the two are closely intertwined and interact in time and space. Understanding the interaction mechanism between earthquake slip and aseismic slip has important scientific and practical significance for revealing the earthquake cycle evolution process, assessing seismic risk, and formulating effective disaster reduction strategies.

[0003] In recent years, with the continuous development of earthquake observation methods, such as the popularization of high-precision ground deformation observation technologies such as GPS and InSAR, more and more studies have found that aseismic slip exists widely in many major fault systems and plays a key role in mainshock-aftershock evolution, stress transfer and triggering mechanism. This has prompted the academic community to have a strong interest in the collaborative evolution process of earthquake-aseismic slip (SEAS, Sequences of Earthquakes and Aseismic Slip).

[0004] Seismic hazard assessment aims to quantify the probability of earthquakes and the possible ground shaking they may cause in a specific region within a certain period of time. Traditional seismic hazard assessment methods mostly rely on historical earthquake catalogs, statistical laws (such as Gutenberg-Richter law) and empirical ground motion prediction models. Such methods have certain accuracy in regions with more data, but their effectiveness is limited in regions with sparse seismic activity or potential hazards that cannot be revealed through statistics.

[0005] In contrast, seismic cycle numerical simulation based on physical mechanisms can simulate the long-term sliding evolution process of real fault systems under geological constraints such as fault geometry and friction properties, including the triggering of main earthquakes, aftershock attenuation, aseismic slip and other key processes. In particular, three-dimensional quasi-dynamic simulation can fundamentally reveal the following physical mechanisms closely related to seismic hazard assessment. Numerical simulation can not only be used for scientific research, but also can provide physical constraints and evolution trend prediction for seismic hazard assessment, making up for the shortcomings of traditional statistical methods and improving the scientificity, spatial and temporal resolution and forward-looking of the assessment.

[0006] To gain a deeper understanding of this complex process, researchers have developed various numerical simulation methods. Among them, quasi-dynamic models based on the Boundary Element Method (BEM) have attracted widespread attention due to their high computational efficiency and suitability for long-scale simulations. However, current mainstream simulation codes such as RSQSim and GARNET are mostly two-dimensional simplified models, and they suffer from problems such as closedness, difficulty in expansion, and lack of support for GPUs or MPI, which limit their application on high-performance computing platforms and hinder researchers' exploration under complex three-dimensional tomographic geometry conditions.

[0007] Seismic hazard assessment is one of the core tasks of earthquake engineering and disaster prevention and mitigation research. Traditional methods rely on parameters such as earthquake catalog statistics, earthquake geological surveys, and ground motion attenuation relationships to construct probabilistic models. However, these methods are difficult to fully reflect the physical processes of deep faults, and especially difficult to consider the increased hazard caused by the linkage mechanism of slow slip and local strong earthquakes. Summary of the Invention

[0008] This application addresses at least one drawback of existing technologies by proposing a seismic hazard assessment method based on quasi-dynamic seismic cycle numerical simulation. This method, based on physical mechanisms, simulates the long-term slip evolution of real fault systems under geological constraints such as fault geometry and friction properties. Details are as follows:

[0009] This invention proposes a method for seismic hazard assessment, comprising the following steps:

[0010] Receive the three-dimensional geometric model of the meshed fault and set the initial conditions;

[0011] Simulation equations are constructed to solve for shear stress loading rate, normal stress loading rate, and state variables. Based on the simulation equations, an ordinary differential equation system is constructed.

[0012] Based on the aforementioned three-dimensional geometric model and initial conditions, the aforementioned ordinary differential equation system is used to accelerate the simulation calculation of shear stress loading rate, normal stress loading rate, and state variables, and the seismic hazard is assessed based on the simulation results.

[0013] The method for accelerating simulation calculations includes: using a combination of H-matrices and MPI to accelerate simulation calculations, including obtaining a layered matrix based on the Green's function under half-space conditions combined with an adaptive cross-approximation algorithm, including the normal stiffness H-matrix and the shear stiffness H-matrix, and using MPI parallel matrix-vector multiplication.

[0014] Preferably, setting the initial conditions includes setting the initial stress state, configuring friction parameters, and calibrating the seismic period.

[0015] The friction coefficient comprises a parameter and a parameter The parameter and the parameter are interval adjustable parameters;

[0016] The friction parameter a and the parameter b are iteratively optimized in combination with a regional characteristic earthquake period.

[0017] Preferably, the method for constructing the simulation operation equation comprises: constructing a shear stress and normal stress control equation expression of elastic stress transfer caused by fault slip in combination with a radiation damping assumption; constructing a differential equation group in combination with the control equation and a boundary condition with the rate-state friction law as the boundary condition; and combining the differential equation group into a common differential equation system for solving.

[0018] Preferably, the sliding velocity and the sliding amount are calculated based on the operation results of the shear stress loading rate, the normal stress loading rate and the state variable, and the earthquake occurrence time and the magnitude are calculated based on the sliding velocity and the sliding amount.

[0019] Meanwhile, the occurrence time, the magnitude distribution, the rupture propagation path, the surface deformation mode of multiple earthquake events and the space-time distribution of high intensity areas are simulated based on the sliding velocity, the sliding amount, the earthquake occurrence time and the magnitude.

[0020] The method for constructing the three-dimensional geometric model of the gridded fault comprises:

[0021] The spatial geometric shape of the fault is determined according to high-precision earthquake positioning results and the trajectory of the surface active fault,

[0022] The three-dimensional geometric model of the fault is established by using CAD or GIS software;

[0023] The fault surface is triangularly meshed by using meshing software.

[0024] Preferably, the adaptive cross-approximation algorithm in the H-matrix framework is used to select part of key rows and columns, and the iteration approximation is performed based on the selected pivotal elements in the residual, so that the low-rank representation of the matrix is effectively constructed, and the layered matrix is obtained.

[0025] The present application further provides a seismic risk assessment system, comprising:

[0026] An input unit is configured to receive the three-dimensional geometric model of the gridded fault and set initial conditions;

[0027] A CPU simulation unit is configured to construct simulation operation equations to solve the shear stress loading rate, the normal stress loading rate and the state variable, and construct a common differential equation system based on the simulation operation equations.

[0028] Based on the three-dimensional geometric model and initial conditions, the acceleration simulation operation of the shear stress loading rate, the normal stress loading rate and the state variable is realized by using the system of ordinary differential equations;

[0029] The method of the acceleration simulation operation comprises: using the method of combining H-matrix and MPI to accelerate the simulation operation, comprising obtaining layered matrix based on the Green function under the condition of half space and combining the adaptive cross approximation algorithm, comprising normal stiffness H-matrix and shear stiffness H-matrix, and using MPI parallel operation matrix-vector multiplication;

[0030] An evaluation unit is configured to evaluate the seismic risk according to the simulation operation result.

[0031] The present application proposes a seismic risk evaluation method, which innovatively introduces hierarchical matrix (or H-matrix), applies H-matrix technology to regional scale seismic period simulation, effectively compresses the storage size of the interaction matrix between faults, reduces it from the original O(N²) to close to O(N log N), greatly reduces the memory consumption, calculation complexity and calculation time under the premise of ensuring accuracy, and provides the possibility for high-resolution modeling.

[0032] Meanwhile, combined with efficient boundary element calculation and H-matrix acceleration method, the stress change, sliding rate evolution and rupture mode evolution of each position on the fault are calculated, and then the potential seismic risk area and its probability level are derived, realizing the conversion from "simulation of sliding behavior" to "evaluation of damage risk". This method has the advantages of high precision, high efficiency and high interpretability, and provides important technical support for earthquake prediction and risk decision-making.

[0033] Meanwhile, the MPI parallel computing framework is fused to realize efficient operation of the seismic simulation task in a multi-core / cluster environment. By dividing the calculation task according to the spatial distribution or time step and using MPI technology for distributed parallel operation, the present application significantly improves the efficiency of large-area, long-time and high-spatial-resolution seismic sequence simulation, realizes effective utilization of more than 100-core computing resources, and greatly shortens the simulation period. BRIEF DESCRIPTION OF DRAWINGS

[0034] In order to more clearly illustrate the technical solutions in the embodiments of the present application or the prior art, the drawings needed to be used in the embodiments or prior art description will be briefly introduced. Obviously, the drawings in the following description are only some embodiments of the present application, and other drawings can be obtained by those skilled in the art without creative labor.

[0035] Figure 1 is a framework diagram of a seismic risk evaluation system;

[0036] Figure 2 is a schematic diagram of the structure of H-matrix in the operation process and the MPI distribution result;

[0037] Figure 3 is a schematic diagram of the efficiency and memory comparison of central processing unit and graphics processing unit;

[0038] Figure 4 is a schematic diagram of the application of PyQuake3D in the seismic risk assessment of Anatolian fault in Turkey. DETAILED DESCRIPTION

[0039] The technical solutions in the embodiments of the present application will be described clearly and completely below with reference to the drawings. Obviously, the described embodiments are only some of the embodiments of the present application, not all the embodiments. Based on the embodiments in the present application, all other embodiments obtained by those skilled in the art without creative work are within the scope of protection of the present application. The terms "comprise" and "have" and any variations thereof are intended to cover non-exclusive inclusion, so that the processes, methods, systems, products or equipment containing a series of units do not have to be limited to those units, but can include other units not clearly listed or inherent to these processes, methods, products or equipment.

[0040] The present application discloses a seismic risk assessment system, which is a three-dimensional seismic cycle simulation system. The present disclosure is called PyQuake3D, which is developed based on Python language and uses OpenGL graphics library for three-dimensional graphics rendering. After installation, users call the functions of PyQuake3D through Python scripts to perform visual analysis of terrain, faults and seismic waves. PyQuake3D supports three-dimensional simulation of multi-fault geometry in full space and half space conditions. It can handle various complex fault structures such as curved surfaces, steps and folds, enhancing its applicability in real tectonic scenarios.

[0041] The present application discloses a seismic risk assessment method, which is applied to the three-dimensional seismic cycle simulation system and includes the following steps, referring to Figure 1 :

[0042] A seismic risk assessment method includes the following steps:

[0043] Receive a three-dimensional geometric model of a gridded fault and set initial conditions;

[0044] Construct simulation operation equations to solve shear stress loading rate, normal stress loading rate and state variables, and construct ordinary differential equation system based on simulation operation equations;

[0045] Based on the aforementioned three-dimensional geometric model and initial conditions, the aforementioned ordinary differential equation system is used to accelerate the simulation calculation of shear stress loading rate, normal stress loading rate, and state variables, and the seismic hazard is assessed based on the simulation results.

[0046] The accelerated simulation method includes: using a combination of H-matrix and MPI to accelerate simulation calculations, including obtaining a layered matrix (including the normal stiffness H-matrix and the shear stiffness H-matrix) based on the Green's function under half-space conditions combined with an adaptive cross-approximation algorithm, and employing MPI parallel matrix-vector multiplication. This invention, using an H-matrix acceleration algorithm and MPI parallel technology, can achieve large-scale fault system seismic period simulation on a millennium scale, significantly outperforming traditional finite differential / finite element methods and breaking through the computational bottleneck of long-term, multi-period simulations.

[0047] The methods for constructing simulation equations include: constructing shear stress and normal stress control equations based on the radiation damping assumption to express the elastic stress transfer caused by fault slip; using the rate-state friction law as boundary conditions, constructing a system of differential equations by combining the control equations and boundary conditions; and combining the differential equations into a system of ordinary differential equations for solving.

[0048] The methods for constructing and meshing a three-dimensional geometric model of a fault include:

[0049] Based on the high-precision seismic location results and the trajectory of active surface faults, the spatial geometry of the faults is determined;

[0050] Use CAD or GIS software to create a three-dimensional geometric model of the fault;

[0051] Mesh generation software was used to divide the fault plane into triangular meshes, resulting in a file format (.msh) that meets the input requirements of the three-dimensional seismic period simulation system.

[0052] Compared to traditional methods that use regular meshes or simplified geometric modeling, this invention combines CAD modeling with high-order mesh generation tools such as Gmsh, which can support refined modeling of faults of arbitrary shape, branches, and tortuous faults. This makes the simulated fault geometry closer to the actual geological conditions, thereby improving the accuracy and reliability of the simulation.

[0053] Setting initial conditions includes setting the initial stress state, configuring friction parameters and calibrating the seismic period, and the friction coefficient includes parameters. and parameters As a preferred option, the parameters and parameters The parameters are adjustable within an interval, and the friction parameters a and b are iteratively optimized by combining regional seismic cycles.

[0054] Wherein, the initial stress condition and the friction coefficient are input according to fault related information, mainly from the observation, inversion result of predecessors, or experience formula. Unlike the traditional manual setting regional stress direction, the present application determines the principal stress direction and stress ratio by using regional focal mechanism solution and GPS data, and automatically projects the regional stress tensor to the fault plane by PyQuake3D, which simplifies the parameter setting process and improves the physical rationality of stress boundary condition setting.

[0055] Based on the differential differential equation, a constant differential equation group with dimension 4N is constructed , and the H-matrix algorithm Runge-Kutta method is used to solve the constant differential equation.

[0056] Wherein, the control equation of shear stress and normal stress is shown as formula (1) and formula (2),

[0057] (1)

[0058] (2)

[0059] Wherein, the third term in formula (1) is a radiation damping term, and the introduction of the term can avoid the divergence of slip rate caused by instability in the quasi-static model;

[0060] Indicates the slip rate of tectonic loading;

[0061] Indicates the stiffness modulus;

[0062] S is the S wave velocity, that is, the propagation speed of seismic transverse wave in medium;

[0063] The slip amount of the jth fault unit;

[0064] Is the stiffness matrix, which is used to quantify the stress caused by the unit slip of the jth fault unit on the ith fault unit, wherein, Is the normal stiffness matrix, Is the shear stiffness matrix;

[0065] Is the slip amount of the ith fault unit, which represents the relative displacement of the unit;

[0066] Is the initial shear stress, which represents the shear stress state of the fault when it is not disturbed;

[0067] Is the initial normal stress, which represents the normal stress state of the fault when it is not disturbed.

[0068] The friction coefficient expression of the rate-and-state friction law is expressed in the form of a regularized aging law, as shown in equations (3) and (4) below:

[0069] (3)

[0070] (4)

[0071] wherein, is the characteristic slip distance, is the reference slip rate, is the reference friction coefficient, is the state variable, is the slip rate, and parameters and are used to describe the direct effect and evolution effect of the shear resistance during the rate jump and thereafter, respectively.

[0072] Through a series of operations of the chain rule, combined with the control equations and boundary conditions, a solvable difference differential equation system is constructed, as shown in equations (5)-(8) below:

[0073] (5)

[0074] (6)

[0075] (7)

[0076] (8)

[0077] Based on equations (5)-(8), a four-dimensional ordinary differential equation system is constructed, as shown in equations (9) and (10) below:

[0078] (9)

[0079] (10)

[0080] wherein the superscripts 1, 2 in equations (5) and (6) represent the shear stress or slip velocity in the strike-slip and dip-slip directions, respectively, and are the external shear stress and normal stress loading rates, is the state variable after the rupture mode conversion, and the conversion formula is as follows:

[0081]

[0082] The 5th order Runge-Kutta method with adaptive step size is used to solve the system of ordinary differential equations, and the MPI parallel operation matrix-vector multiplication is used. Specifically, it includes:

[0083] Based on the Green function and the adaptive cross approximation algorithm, a hierarchical matrix, H-matrix, is constructed. The principle is as follows:

[0084] The three-dimensional geometric model of the fault is discretized into a plurality of elements, each element is controlled by a set of coupled ordinary differential equations, and the ordinary differential equations describe the quasi-dynamic rupture evolution. The interaction between the elements is adjusted by a pre-calculated stress Green function, which defines how the slip on the source element causes stress changes on all other receiving elements. The H-matrix divides a large dense matrix into a plurality of smaller sub-blocks through a hierarchical (tree-like) structure, and each MPI process is responsible for processing a plurality of sub-blocks. These sub-blocks are processed as (1) low-rank sub-blocks: stored in the form of U and V, used to represent local low-rank matrices; or (2) full matrix sub-blocks: directly store all data of the sub-block. In the process of constructing low-rank sub-blocks, each process independently calls algorithms such as adaptive cross approximation to reduce inter-process communication. In the MPI parallel operation matrix-vector multiplication phase, each process locally calculates the product part of the sub-block it is responsible for and the vector, and then uses the MPI_Reduce operation to combine the partial results of all processes into a complete vector y. As shown in Figure 3 , the H-matrix structure and the parallel distribution of sub-blocks (sub-matrices) are shown.

[0085] Further, the present application discloses a comparative analysis of the storage and calculation efficiency of H-matrix (also known as hierarchical matrix or H-matrix) and NumPy, wherein NumPy is a dense matrix (Dense Matrix) calculation library of Python, which is a widely used open source scientific calculation library in Python language, and can efficiently handle large-scale multi-dimensional array and matrix operations, and provide rich mathematical function support. While the H-matrix technology divides a large-scale dense matrix into blocks and constructs a low-rank approximation, the storage and calculation complexity is reduced from to N *log( N ). With the increase of the number of elements, the performance of H-matrix in matrix-vector multiplication operation is obviously better than that of NumPy. As shown in Figure 2 , wherein Figure 2 (a) of FIG. 5 is a comparison diagram of the time required for the central processing parallelism to achieve graph parallelism and MPI acceleration in a matrix-vector volume of 500 steps; Figure 2 (b) is a comparison diagram of the memory cost of the dense matrix and the H-matrix.

[0086] In terms of parallel computing, MPI uses 14 cores for acceleration, which improves the computational efficiency of matrix-vector multiplication by about two times compared with using only 2 cores. However, since NumPy itself has implemented some parallel optimization at the bottom layer, the acceleration effect is not significant under different core numbers, such as Figure 2 (a) of FIG. 1.

[0087] Although the performance of GPU is much better than that of CPU in terms of matrix-vector operation efficiency, the performance gap can reach nearly 1000 times. However, since GPU computing needs to be based on a complete dense matrix, its storage resource consumption is relatively higher, such as Figure 2 (b) of FIG. 1. In addition, the hardware cost of GPU is usually higher than that of CPU. Therefore, the combination of H-matrix and MPI proposed in the present application is more advantageous for large-scale numerical simulation tasks such as seismic numerical simulation, while GPU is more suitable for processing small-scale models.

[0088] In practical applications, even combined with parallel computing technology, it is still very time-consuming to obtain a high-precision complete dense matrix. To solve this problem, the adaptive cross approximation algorithm in the H-matrix framework is used to select some key rows and columns, and based on the selected pivot elements in the residual, iterative approximation is performed to effectively construct a low-rank representation of the matrix, forming a low-rank layered matrix, thereby avoiding the overall calculation of the entire matrix.

[0089] The present application also discloses a seismic risk assessment system, comprising:

[0090] an input unit configured to receive a three-dimensional geometric model of a gridded fault and set initial conditions;

[0091] a CPU simulation unit configured to construct a simulation operation equation to solve a shear stress loading rate, a normal stress loading rate and a state variable, and construct an ordinary differential equation system based on the simulation operation equation;

[0092] based on the three-dimensional geometric model and the initial conditions, the ordinary differential equation system is used to realize accelerated simulation operation of the shear stress loading rate, the normal stress loading rate and the state variable;

[0093] wherein the method of accelerated simulation operation comprises: using an H-matrix and MPI combined method to accelerate simulation operation, including obtaining a layered matrix based on a Green function under a half-space condition combined with an adaptive cross approximation algorithm, including a normal stiffness H-matrix and a shear stiffness H-matrix, and using MPI parallel operation matrix-vector multiplication;

[0094] an evaluation unit configured to evaluate seismic risk according to simulation operation results.

[0095] Further, the seismic risk assessment system further comprises a GPU simulation unit,

[0096] a GPU simulation unit configured to construct a dense stiffness matrix based on a Green function, including a dense shear stress stiffness matrix and a normal stress stiffness matrix;

[0097] based on the three-dimensional geometric model and initial conditions, the constant differential equation system is used to realize accelerated simulation operation of shear stress loading rate, normal stress loading rate and state variable;

[0098] wherein, the simulation operation equation: the shear stress and normal stress control equation expression caused by elastic stress transfer due to fault slip is constructed by combining the radiation damping assumption; the solvable differential equation group is constructed by combining the control equation and the boundary condition with the rate-state friction law as the boundary condition; the solvable differential equation group is converted into an N-dimensional constant differential equation system;

[0099] the accelerated simulation operation method: the constant differential equation is solved by using the 5th order Runge-Kutta method, and the GPU parallel acceleration is used;

[0100] The CPU simulation unit and the GPU simulation unit are selectively applied to simulation operation.

[0101] Further, the application of PyQuake3D in the earthquake risk assessment of Anatolian fault in Turkey is disclosed as an application example, which proves the advantages of the method and system disclosed in the application.

[0102] (1) Constructing a three-dimensional geometric model of the fault and meshing

[0103] First, the spatial geometry of the fault is determined according to the high-precision earthquake location results and the trajectory of the surface active fault. A three-dimensional geometric model of the fault is established using CAD or GIS software (such as AutoCAD, QGIS). Further, the fault surface is triangularly meshed by using the open-source meshing software Gmsh, so that the model has the ability to adapt to any complex geometry, and finally the.msh file format conforming to the input requirements of PyQuake3D is exported.

[0104] (2) Setting the initial stress state

[0105] The stress parameters are selected according to the focal mechanism solution, GPS observation data, regional stress field and other information, combined with the principal stress direction and stress ratio. PyQuake3D internally integrates stress projection function, which can automatically project the three-dimensional regional principal stress tensor onto the fault surface to generate the initial shear stress and normal stress on the fault surface, providing a reasonable starting point for subsequent slip evolution simulation.

[0106] (3) Friction parameter configuration and earthquake period calibration

[0107] The friction parameters a and b are crucial for controlling the seismic period. They are usually set within a reasonable range based on experimental research and previous experience. For example, a = 0.005-0.015, b = 0.01-0.025. Through trial and error, the parameter combination is adjusted to simulate different seismic periods, which are then compared with the characteristic seismic period recorded in history to optimize the model parameters and make the simulation results as close to the actual observations as possible.

[0108] (4) Conducting millennial-scale seismic period simulation and disaster assessment

[0109] Based on the set three-dimensional geometric model of the fault, the initial stress state and the friction parameters, PyQuake3D is used to simulate the long-time (millennial-scale) seismic period evolution of the target fault.

[0110] The simulation results of PyQuake3D include: based on the calculation results of the shear stress loading rate, the normal stress loading rate and the state variable, the slip rate and the slip amount are calculated, and based on the slip rate and the slip amount, the earthquake occurrence time, the magnitude and other rich information are calculated; at the same time, based on the slip rate, the slip amount, the earthquake occurrence time and the magnitude, the occurrence time, the magnitude distribution, the rupture propagation path, the surface deformation mode and the spatio-temporal distribution of high intensity area of multiple earthquake events are simulated. These results can be used to quantify the risk level that may be brought by future earthquakes, and provide scientific basis for seismic hazard zoning, disaster reduction and prevention strategy formulation and infrastructure planning.

[0111] The simulation results are shown in Figure 4 . Figure 4 (a) of which is a seismic rupture process and aftershock distribution map, showing the co-seismic slip spatial distribution of six typical earthquake events, with color from dark blue (0 m) to bright yellow (10 m) representing the slip amount. Each subgraph represents an earthquake event, the red arrow indicates the rupture propagation direction, and the red star indicates the rupture starting point (focus), and the yellow star indicates the main shock location. The fault geometry is complex, including multiple bending and branching structures, showing that the earthquake rupture may propagate along different paths. Figure 4 (b) of which is a seismic period and magnitude information map. The bar chart in the figure represents the interval period of each event, with units of years (yr) or minutes (minutes), and the blue column corresponds to the repeating period of the static slip event, and the red bar represents the main shock magnitude (Mw). The middle text annotation indicates that a certain event is a slow earthquake occurring on a fault branch. The earthquake magnitude is mainly concentrated in the Mw7.4 to Mw8.3 interval, showing that the fault has the ability to produce large earthquakes and there is diversity in rupture recurrence. The period length in the figure varies greatly, from tens of minutes to more than two hundred years, revealing the irregularity between strong earthquakes. Figure 4The (c) is the moment rate change graph in the same earthquake time sequence. Each graph shows the moment rate evolution with time in an earthquake event, which can be used to analyze the speed of energy release of the earthquake. Different curve shapes reflect different rupture duration and energy release mode of different earthquakes, for example: the first event has long duration and high peak, which represents a typical strong earthquake; the fourth earthquake has only tens of seconds of duration, which represents a fast rupture event; the sixth earthquake has multiple energy release peaks, which may represent multi-segment rupture or composite event.

[0112] Based on the above disclosure, the present application can complete the simulation of the millennium-level earthquake cycle of complex structures such as fault segmentation, branch fault, fault jump, etc. within a reasonable time

[0113] Also disclosed is a computer readable storage medium storing a computer program, which makes a computer execute to realize the above-mentioned earthquake risk assessment method.

[0114] For example, the computer program can be divided into one or more modules / units, one or more modules / units are stored in the memory and executed by the processor, and the I / O interface transmission of data is completed by the input interface and the output interface to complete the present application. One or more modules / units can be a series of computer program instruction segments capable of completing a specific function, which are used to describe the execution process of the computer program in the computer device.

[0115] The above is only a specific embodiment of the present application, but the protection scope of the present application is not limited thereto. Any change or replacement within the technical scope disclosed by the present application should be covered within the protection scope of the present application. Therefore, the protection scope of the present application should be subject to the protection scope of the claims.

Claims

1. A method for evaluating seismic hazard, characterized by, The method comprises the following steps: receiving a three-dimensional geometric model of a gridded fault and setting initial conditions; constructing simulation operation equations to solve shear stress loading rate, normal stress loading rate and state variables, and constructing an ordinary differential equation system based on the simulation operation equations; based on the three-dimensional geometric model and the initial conditions, using the ordinary differential equation system to realize accelerated simulation operation of the shear stress loading rate, the normal stress loading rate and the state variables, and evaluating seismic risk according to the simulation operation result; wherein the method for accelerated simulation operation comprises: using a method combining H-matrix and MPI to accelerate simulation operation, including obtaining a layered matrix based on Green function under half-space condition combined with an adaptive cross approximation algorithm, including normal stiffness H-matrix and shear stiffness H-matrix, and using MPI parallel operation matrix-vector multiplication.

2. The method of seismic hazard assessment according to claim 1, wherein, The setting of the initial conditions comprises setting an initial stress state, friction parameter configuration and seismic period calibration, The friction parameters comprise the parameters and the parameters , the parameters and the parameters are interval adjustable parameters; combining with regional characteristic seismic periods, iteratively optimizing friction parameters a and b.

3. The method of claim 1, wherein the step of determining the seismic hazard comprises the step of: The method for constructing simulation operation equations comprises: constructing shear stress and normal stress control equation expressions caused by elastic stress transfer due to fault slip combined with radiation damping assumption; constructing a differential equation group combined with control equations and boundary conditions taking rate-state friction law as boundary conditions; and combining the differential equation group into an ordinary differential equation system for solving. ​ 4. The method of claim 1, wherein the step of determining the seismic hazard comprises the step of: Based on the operation results of the shear stress loading rate, the normal stress loading rate and the state variables, the sliding rate and the sliding amount are calculated, and the occurrence time and magnitude of an earthquake are calculated based on the sliding rate and the sliding amount; ​ Meanwhile, based on the sliding rate, the sliding amount, the occurrence time and the magnitude of an earthquake, the occurrence time, the magnitude distribution, the rupture propagation path, the surface deformation mode of multiple earthquake events, and the spatial and temporal distribution of high intensity areas are simulated.

5. The method of seismic hazard assessment according to claim 1, wherein, The method for constructing the three-dimensional geometric model of the gridded fault comprises: determining the spatial geometric shape of the fault according to high-precision earthquake positioning results and the trajectory of the surface active fault, establishing a three-dimensional geometric model of the fault using CAD or GIS software; using a meshing software to perform triangular mesh partitioning on the fault surface.

6. The method of seismic hazard assessment according to claim 1, wherein, The accelerated simulation operation comprises: using an adaptive cross approximation algorithm in the H-matrix framework to select part of key rows and columns, and iteratively approximating based on the selected pivotal elements in the residual, to effectively construct a low-rank representation of the matrix and obtain a layered matrix.

7. A seismic hazard assessment system, characterized by, It comprises: an input unit for receiving a three-dimensional geometric model of a gridded fault and setting initial conditions; a CPU simulation unit for constructing simulation operation equations to solve shear stress loading rate, normal stress loading rate and state variables, and constructing an ordinary differential equation system based on the simulation operation equations; based on the three-dimensional geometric model and the initial conditions, using the ordinary differential equation system to realize accelerated simulation operation of the shear stress loading rate, the normal stress loading rate and the state variables; and The method for accelerating the simulation operation comprises: accelerating the simulation operation by using a method combining an H-matrix and MPI, comprising obtaining a layered matrix based on a Green function under a half-space condition and combining an adaptive cross approximation algorithm, comprising a normal stiffness H-matrix and a shear stiffness H-matrix, and performing matrix-vector multiplication by using MPI parallel operation; The evaluation unit is configured to evaluate the seismic risk according to the simulation operation result.

8. A seismic hazard assessment system as in claim 7, wherein, The GPU simulation unit is further configured to construct a dense stiffness matrix based on the Green function, the dense stiffness matrix comprising a dense shear stress stiffness matrix and a normal stress stiffness matrix. The GPU simulation unit is further configured to accelerate the simulation operation of the shear stress loading rate, the normal stress loading rate and the state variable by using the system of ordinary differential equations based on the three-dimensional geometric model and the initial condition. The simulation operation equation is constructed by combining a radiation damping assumption to construct shear stress and normal stress control equations to express elastic stress transfer caused by fault slip. A rate-state friction law is used as a boundary condition, the differential equation set is constructed by combining the control equation and the boundary condition, and the differential equation set is constructed into a solvable system of ordinary differential equations. The method for accelerating the simulation operation is based on an H-matrix to construct a low-rank matrix and uses MPI parallel acceleration; meanwhile, an option of GPU parallel acceleration is provided. The CPU simulation unit and the GPU simulation unit are selectively applied to the simulation operation. ​

Images