Earthquake risk evaluation method and system

Abstract

The invention discloses an earthquake risk evaluation method and system, and relates to the technical field of earthquake numerical simulation, and the method comprises the steps: receiving a three-dimensional geometric model of a grid fault, and setting an initial condition; constructing a simulation operation equation to solve a shear stress loading rate, a normal stress loading rate and a state variable, and constructing an ordinary differential equation system; based on the three-dimensional geometric model and the initial conditions, accelerated simulation operation of the shear stress loading rate, the normal stress loading rate and the state variables is achieved through the ordinary differential equation system, and the accelerated simulation operation method comprises the steps that accelerated simulation operation is achieved through the method of combining an H-matrix and MPI, comprising the steps that a layered matrix is obtained based on a Green function under the half-space condition in combination with an adaptive cross-approximation algorithm, the layered matrix comprises a normal stiffness H-matrix and a shear stiffness H-matrix, and MPI parallel operation matrix-vector multiplication is adopted. According to the invention, the hierarchical matrix is introduced and the MPI parallel computing framework is fused, so that the simulation period is greatly shortened.

Description

Technical Field

[0001] The present invention relates to the technical field of earthquake simulation, and in particular to an earthquake risk assessment method and system. Background Art

[0002] Earthquakes are one of the most significant natural phenomena resulting from the release of energy from Earth's interior, and their generation mechanism is closely related to fault activity. Fault slip processes include both sudden seismic ruptures and slow, aseismic slip, which are closely intertwined and influence each other in time and space. Understanding the interaction between seismic and aseismic slip has important scientific and practical implications for revealing the evolution of earthquake cycles, assessing earthquake hazard, and developing effective disaster reduction strategies.

[0003] In recent years, with the continuous development of seismic observation methods, such as the popularization of high-precision surface deformation observation technologies such as GPS and InSAR, an increasing number of studies have found that aseismic slip is widespread across multiple major fault systems and plays a key role in mainshock-aftershock evolution, stress transfer, and triggering mechanisms. This has sparked a strong interest in the co-evolution of earthquakes and aseismic slip (SEAS), a process that occurs when a fault is pulled apart.

[0004] Seismic hazard assessment (SHA) aims to quantify the probability of an earthquake occurring in a specific region within a certain period of time, and the resulting surface shaking. Traditional seismic hazard assessment methods rely primarily on historical earthquake catalogs, statistical laws (such as the Gutenberg-Richter law), and empirical earthquake prediction models. While these methods are somewhat accurate in areas with abundant data, their effectiveness is limited in areas with sparse seismic activity or where potential hazards are difficult to detect statistically.

[0005] In contrast, numerical simulations of earthquake cycles based on physical mechanisms can simulate the long-term slip evolution of real fault systems, including key processes such as mainshock triggering, aftershock attenuation, and aseismic slip, within geological constraints such as fault geometry and frictional properties. In particular, three-dimensional quasi-dynamic simulations can fundamentally reveal the following physical mechanisms closely related to hazard assessment. Numerical simulations are not only useful for scientific research but also provide physical constraints and evolutionary trend predictions for earthquake hazard assessment, addressing the shortcomings of traditional statistical methods and improving the scientific nature, spatial and temporal resolution, and foresight of assessments.

[0006] To gain a deeper understanding of this complex process, researchers have developed a variety of 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-term simulations. However, current mainstream simulation codes, such as RSQSim and GARNET, are mostly simplified two-dimensional models and suffer from limitations such as closed-endedness, difficulty in scalability, and lack of GPU or MPI support. This limits their applicability on high-performance computing platforms and hinders research exploring complex three-dimensional fault geometries.

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

[0008] This application addresses at least one shortcoming of the existing technology and proposes a seismic hazard assessment method based on quasi-dynamic earthquake cycle numerical simulation. This physical mechanism-based earthquake cycle numerical simulation can simulate the long-term slip evolution of a real fault system under geological constraints such as fault geometry and friction properties. The details are as follows: The present invention proposes a method for earthquake risk assessment, comprising the following steps: Receive a gridded 3D fault geometry model and set initial conditions; Construct simulation operation equations to solve shear stress loading rate, normal stress loading rate and state variables, and construct ordinary differential equation systems based on the simulation operation equations; Based on the three-dimensional geometric model and initial conditions, the ordinary differential equation system is used to realize accelerated simulation calculation of shear stress loading rate, normal stress loading rate and state variables, and the earthquake hazard is evaluated according to the simulation calculation results; Among them, the method for accelerating simulation operations includes: using a method combining H-matrix and MPI to accelerate simulation operations, including obtaining a layered matrix based on the Green function under half-space conditions 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.

[0009] Preferably, the setting of initial conditions includes setting initial stress state, friction parameter configuration and earthquake period calibration, The friction coefficient includes parameters and parameters , the parameters and parameters is an interval adjustable parameter; Combined with the regional characteristic earthquake cycle, the friction parameters a and b are optimized through iteration.

[0010] Preferably, the method of constructing simulation operation equations includes: constructing shear stress and normal stress control equations in combination with the radiation damping assumption to express the elastic stress transfer caused by fault slip; using the rate-state friction law as the boundary condition, combining the control equations and the boundary conditions to construct a group of differential equations; and merging the group of differential equations into a system of ordinary differential equations for solving.

[0011] Preferably, based on the calculation results of the shear stress loading rate, the normal stress loading rate and the state variable, the sliding rate and the sliding amount are calculated, and the time and magnitude of the earthquake are calculated based on the sliding rate and the sliding amount; At the same time, based on the slip rate, slip amount, earthquake occurrence time and magnitude, the occurrence time, magnitude distribution, rupture propagation path, surface deformation pattern and spatial and temporal distribution of high-intensity areas of multiple earthquake events are simulated.

[0012] The method for constructing the gridded three-dimensional geometric model of the fault comprises: Based on the high-precision earthquake positioning results and the trajectory of the surface active fault, the spatial geometry of the fault is determined. Use CAD or GIS software to create a three-dimensional geometric model of the fault; The fault plane is divided into triangular meshes using meshing software.

[0013] Preferably, an adaptive cross-approximation algorithm in an H-matrix framework is used to select some key rows and columns and perform iterative approximation based on the main element selected from the residual to effectively construct a low-rank representation of the matrix and obtain a hierarchical matrix.

[0014] The present invention also provides an earthquake risk assessment system, comprising: An input unit, for receiving a gridded three-dimensional geometric model of the fault and setting initial conditions; The CPU simulation unit constructs simulation operation equations to solve the shear stress loading rate, normal stress loading rate and state variables, and constructs an ordinary differential equation system based on the simulation operation equations; Based on the three-dimensional geometric model and initial conditions, the ordinary differential equation system is used to realize accelerated simulation calculation of shear stress loading rate, normal stress loading rate and state variables; The method for accelerating simulation operations includes: accelerating simulation operations by combining H-matrix and MPI, including obtaining a layered matrix based on the Green's function under half-space conditions 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; An evaluation unit is used to evaluate earthquake hazards based on simulation results.

[0015] This paper proposes a method for earthquake hazard assessment that innovatively introduces a hierarchical matrix (or H-matrix, H-matrix) and applies H-matrix technology to regional-scale earthquake cycle simulation. This effectively compresses the storage size of the inter-fault interaction matrix from the original O(N²) to nearly O(N log N). While ensuring accuracy, it significantly reduces memory consumption, computational complexity, and time, paving the way for high-resolution modeling.

[0016] By combining efficient boundary element calculations with H-matrix acceleration methods, the authors calculated stress changes, slip rate evolution, and rupture pattern evolution at every location on the fault, thereby deriving potential earthquake hazard areas and their probability levels, achieving a transition from "simulating slip behavior" to "assessing damage risk." This method offers the advantages of high precision, high efficiency, and high interpretability, providing important technical support for earthquake prediction and risk decision-making.

[0017] At the same time, by integrating the MPI parallel computing framework, the earthquake simulation tasks can be efficiently run in a multi-core / cluster environment. By dividing the computing tasks according to spatial distribution or time step and using MPI technology for distributed parallel computing, the present invention significantly improves the efficiency of large-area, long-term, and high-spatial-resolution earthquake sequence simulation, realizes the effective utilization of computing resources of more than 100 cores, and greatly shortens the simulation cycle. BRIEF DESCRIPTION OF THE DRAWINGS

[0018] In order to more clearly illustrate the embodiments of the present application or the technical solutions in the prior art, the following briefly introduces the drawings required for use in the embodiments or the description of the prior art. Obviously, the drawings described below are only some embodiments of the present application. For ordinary technicians in this field, other drawings can be obtained based on these drawings without paying any creative labor.

[0019] Figure 1 It is a framework diagram of earthquake hazard assessment system; Figure 2 It is a schematic diagram of the structure of the H-matrix and the MPI allocation results during the operation process; Figure 3 This is a diagram comparing the efficiency and memory of the CPU and GPU; Figure 4This is a schematic diagram of the application of PyQuake3D in the seismic hazard assessment of the Anatolian fault in Türkiye. DETAILED DESCRIPTION

[0020] The technical solutions in the embodiments of the present application will be clearly and completely described below in conjunction with the accompanying drawings. Obviously, the described embodiments are only a part of the embodiments of the present application, rather than all of the embodiments. Based on the embodiments in the present application, all other embodiments obtained by those skilled in the art without making creative work are within the scope of protection of this application. The terms "including" and "having" and any variations thereof are intended to cover non-exclusive inclusions, so that a process, method, system, product or device comprising a series of units is not necessarily limited to those units, but may include other units that are not clearly listed or inherent to these processes, methods, products or devices.

[0021] This invention discloses a seismic hazard assessment system, specifically a 3D earthquake cycle simulation system. This system, referred to herein as PyQuake3D, is developed in Python and utilizes the OpenGL graphics library for 3D graphics rendering. After installation, users can invoke PyQuake3D's functions through Python scripts to visualize terrain, faults, and seismic waves. PyQuake3D supports 3D simulation of multiple fault geometries in both full- and half-space conditions. It can handle a variety of complex fault structures, such as curved surfaces, steps, and folds, enhancing its applicability in real-world tectonic scenarios.

[0022] The present invention discloses a method for evaluating earthquake risk, which is applied to the three-dimensional earthquake cycle simulation system, comprising the following steps: Figure 1 : A method for earthquake risk assessment comprises the following steps: Receive a gridded 3D fault geometry model and set initial conditions; Construct simulation operation equations to solve shear stress loading rate, normal stress loading rate and state variables, and construct ordinary differential equation systems based on the simulation operation equations; Based on the three-dimensional geometric model and initial conditions, the ordinary differential equation system is used to realize accelerated simulation calculation of shear stress loading rate, normal stress loading rate and state variables, and the earthquake hazard is evaluated according to the simulation calculation results; The accelerated simulation method includes: using an H-matrix and MPI combined to accelerate the simulation, including using a Green's function under half-space conditions combined with an adaptive cross-approximation algorithm to obtain a layered matrix including a normal stiffness H-matrix and a shear stiffness H-matrix, and using MPI parallel matrix-vector multiplication. The present invention uses an H-matrix acceleration algorithm and MPI parallel technology to simulate earthquake cycles on a large-scale fault system over millennia, significantly outperforming traditional finite differential / finite element methods and breaking through the computing power bottleneck of long-term, multi-cycle simulations.

[0023] Among them, the method of constructing simulation operation equations includes: combining the radiation damping hypothesis to construct shear stress and normal stress control equations to express the elastic stress transfer caused by fault slip; using the rate-state friction law as the boundary condition, combining the control equations and boundary conditions to construct a group of differential equations; merging the group of differential equations into a system of ordinary differential equations for solving.

[0024] The method of constructing a three-dimensional geometric model of the fault and meshing it includes: Determine the spatial geometry of the fault based on high-precision earthquake location results and the trajectory of the active fault on the surface; Use CAD or GIS software to create a three-dimensional geometric model of the fault; The fault plane is divided into triangular meshes using meshing software to obtain the file format .msh that meets the input requirements of the three-dimensional earthquake cycle simulation system.

[0025] Compared with the traditional method of using regular grids or simplified geometric modeling, the present invention combines CAD modeling with high-order grid generation tools such as Gmsh, which can support the refined modeling of arbitrary morphological, branching, and tortuous faults, making the simulated fault geometry closer to the actual geological conditions, thereby improving the accuracy and credibility of the simulation.

[0026] Among them, setting the initial conditions includes setting the initial stress state, friction parameter configuration and earthquake cycle calibration, and the friction coefficient includes parameters and parameters As a preference, the parameter and parameters The friction parameters a and b are interval adjustable parameters, and are optimized through iteration in combination with the regional characteristic earthquake cycle.

[0027] Initial stress conditions and friction coefficients are input based on fault-related information, primarily derived from previous observations, inversion results, or empirical formulas. Unlike traditional manual setting of regional stress directions, this method utilizes regional focal mechanism solutions and GPS data to determine principal stress directions and stress ratios. PyQuake3D then automatically projects the regional stress tensor onto the fault plane, simplifying the parameter setting process and improving the physical rationality of stress boundary condition settings.

[0028] Constructing a 4N-dimensional system of ordinary differential equations based on differential equations , using the H-matrix algorithm Runge-Kutta method to solve ordinary differential equations.

[0029] Among them, the control equations of the shear stress and normal stress are shown in formula (1) and formula (2), (1) (2) Among them, the third term in formula (1) is the radiation damping term, which can be introduced to avoid the divergence of slip rate caused by instability in the quasi-static model; represents the slip rate of tectonic loading; represents the stiffness modulus; is the S-wave velocity, i.e. the propagation speed of seismic shear waves in the medium; is the slip of the jth fault unit; is the stiffness matrix that quantifies the stress induced on the i-th fault element by a unit slip on the j-th fault element, where is the normal stiffness matrix, is the shear stiffness matrix; is the slip of the i-th fault unit, which represents the relative displacement of the unit; is the initial shear stress, which represents the shear stress state of the fault when it is undisturbed; is the initial normal stress, which represents the normal stress state of the fault when it is undisturbed.

[0030] The friction coefficient expression of the rate-state friction law adopts the regularized aging law expression form, as shown in the following formulas (3) and (4): (3) (4) in, is the characteristic slip distance, is the reference slip rate, is the reference friction coefficient, is the state variable, is the sliding rate, parameter and parameters They are used to describe the direct effect and evolution effect of shear resistance during and after rate mutation, respectively.

[0031] Through a series of operations of the chain rule, the control equations and boundary conditions are combined to construct a solvable differential equation system. The solvable differential equation system is as follows: (5) (6) (7) (8) Based on formula (5)-formula (8), a system of ordinary differential equations with a dimension of 4N is constructed as follows: formula (9) and formula (10): (9) (10) The superscripts 1 and 2 in formula (5) and formula (6) represent the shear stress or sliding velocity in the strike-slip and dip-slip directions, respectively. ,as well as are the external shear stress and normal stress loading rates, is the state variable after the rupture mode is transformed. The transformation formula is as follows:

[0032] A 5th-order Runge-Kutta method with adaptive step size is used to solve the system of ordinary differential equations, and MPI parallel operation matrix-vector multiplication is used. Specifically, it includes: A hierarchical matrix is ​​constructed based on Green's function and adaptive cross approximation algorithm: H-matrix. The principle is as follows: The three-dimensional geometric model of the fault is discretized into several units, each of which is controlled by a set of coupled ordinary differential equations that describe the quasi-dynamic fracture evolution. The interaction between the units is regulated by a pre-calculated stress Green's function, which defines how slip on the source unit causes stress changes on all other receiving units. The H-matrix divides a large dense matrix into multiple smaller sub-blocks through a hierarchical (tree-like) structure, and each MPI process is responsible for processing several sub-blocks. These sub-blocks are processed as (1) low-rank sub-blocks: stored in the form of U and V to represent local low-rank matrices; or (2) full-matrix sub-blocks: all data of the sub-block is directly stored. 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 stage, each process locally calculates the product part of the sub-block it is responsible for and the vector, and then summarizes the partial results of all processes into the complete vector y through the MPI_Reduce operation. Figure 3 As shown, the H-matrix structure and the parallel allocation of sub-blocks (sub-matrices) are demonstrated.

[0033] Furthermore, the present invention discloses a comparative analysis of the storage and computational efficiency of H-matrix (also known as hierarchical matrix or H-matrix) and NumPy. NumPy is a dense matrix computing library for Python. It is an open source scientific computing library widely used in the Python language. It can efficiently handle large-scale multidimensional arrays and matrix operations and provides rich mathematical function support. H-matrix technology reduces storage and computational complexity from 100% to 100% by partitioning large-scale dense matrices and constructing low-rank approximations. Reduce to N *log( N ). As the number of units increases, H-matrix performs significantly better than NumPy in matrix-vector multiplication operations. Figure 2 As shown, Figure 2 (a) is a schematic diagram comparing the time required to achieve graph parallelism and MPI-accelerated central processing parallelism for a matrix-load volume of 500 steps; Figure 2 (b) is a schematic diagram comparing the memory costs of dense matrices and H matrices.

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

[0035] Although GPUs outperform CPUs in terms of matrix-vector computing efficiency, with the performance gap reaching nearly a thousand times, GPU computing requires full dense matrices, so its storage resource consumption is relatively high. Figure 2 (b). Furthermore, GPU hardware costs are generally higher than CPUs. Therefore, the combination of H-matrix and MPI proposed in this invention is more advantageous for large-scale numerical simulation tasks, such as earthquake numerical simulation, while GPUs are more suitable for processing small-scale models.

[0036] In practical applications, even with the use of parallel computing techniques, obtaining a complete dense matrix with high precision is often very time-consuming. To address this issue, the present invention utilizes an adaptive cross-approximation algorithm within the H-matrix framework. By selecting a few key rows and columns and performing iterative approximation based on the selected pivot elements from the residual, the algorithm effectively constructs a low-rank representation of the matrix, forming a low-rank hierarchical matrix, thus avoiding the need for a full computation of the entire matrix.

[0037] The present invention also discloses an earthquake risk assessment system, comprising: An input unit, for receiving a gridded three-dimensional geometric model of the fault and setting initial conditions; The CPU simulation unit constructs simulation operation equations to solve the shear stress loading rate, normal stress loading rate and state variables, and constructs an ordinary differential equation system based on the simulation operation equations; Based on the three-dimensional geometric model and initial conditions, the ordinary differential equation system is used to realize accelerated simulation calculation of shear stress loading rate, normal stress loading rate and state variables; The method for accelerating simulation operations includes: accelerating simulation operations by combining H-matrix and MPI, including obtaining a layered matrix based on the Green's function under half-space conditions 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; An evaluation unit is used to evaluate earthquake hazards based on simulation results.

[0038] Furthermore, the earthquake hazard assessment system also includes a GPU simulation unit. GPU simulation unit, used to construct dense stiffness matrices based on Green's function, including dense shear stress stiffness matrix and normal stress stiffness matrix; Based on the three-dimensional geometric model and initial conditions, the ordinary differential equation system is used to realize accelerated simulation calculation of shear stress loading rate, normal stress loading rate and state variables; Among them, the simulation operation equations: combine the radiation damping assumption to construct the shear stress and normal stress control equations to express the elastic stress transfer caused by fault slip; use the rate-state friction law as the boundary condition, combine the control equations and boundary conditions to construct a solvable differential equation system; convert the solvable differential equation system into an N-dimensional ordinary differential equation system; The accelerated simulation calculation method: uses the 5th order Runge-Kutta method to solve ordinary differential equations and accelerates them in parallel on the GPU; The CPU simulation unit and the GPU simulation unit can be selectively applied to simulation operations.

[0039] Furthermore, an example of the application of PyQuake3D in the seismic hazard assessment of the Anatolian fault in Türkiye is disclosed to demonstrate the advantages of the method and system disclosed in the present invention.

[0040] (1) Constructing a three-dimensional geometric model of the fault and meshing it First, the spatial geometry of the fault is determined based on high-precision earthquake location results and the trajectory of active surface faults. A 3D geometric model of the fault is constructed using CAD or GIS software (such as AutoCAD or QGIS). Using the open-source meshing software Gmsh, the fault plane is triangulated, enabling the model to adapt to arbitrarily complex geometries. Finally, the model is exported in the .msh file format, which complies with PyQuake3D input requirements.

[0041] (2) Setting the initial stress state Stress parameters are selected based on information such as focal mechanism solutions, GPS observations, and the regional in-situ stress field, combined with principal stress directions and stress ratios. PyQuake3D integrates a stress projection function that automatically and accurately projects the three-dimensional regional principal stress tensor onto the fault plane, generating initial shear and normal stresses on the fault plane and providing a reasonable starting point for subsequent slip evolution simulations.

[0042] (3) Friction parameter configuration and earthquake period calibration Friction parameters a and b are crucial for controlling the earthquake period. They are typically set within a reasonable range based on experimental research and previous experience. For example, a = 0.005–0.015 and b = 0.01–0.025 are used. Through trial and error, parameter combinations are adjusted to simulate different earthquake periods. These are then compared with characteristic earthquake periods recorded in historical records to optimize model parameters and ensure that the simulation results are as close to actual observations as possible.

[0043] (4) Conduct millennial-scale earthquake cycle simulation and disaster assessment Based on the set three-dimensional geometric model, initial stress state and friction parameters of the fault, PyQuake3D is used to simulate the long-term (millennial scale) earthquake cycle evolution of the target fault.

[0044] The PyQuake3D simulation results include: calculations of slip rate and slip amount based on the shear stress loading rate, normal stress loading rate, and state variables; and calculations of earthquake occurrence time, magnitude, and other rich information based on the slip rate and slip amount. Furthermore, the simulations use slip rate, slip amount, earthquake occurrence time, and magnitude to simulate the occurrence time, magnitude distribution, rupture propagation paths, surface deformation patterns, and the spatial and temporal distribution of high-intensity zones of multiple earthquakes. These results can be used to quantify the potential risk level of future earthquakes, providing a scientific basis for earthquake hazard zoning, disaster mitigation and prevention strategy development, and infrastructure planning.

[0045] The simulation results are as follows Figure 4 shown. Figure 4 (a) shows the earthquake rupture process and aftershock distribution, showing the spatial distribution of coseismic slip for six typical earthquake events. Colors range from dark blue (0 m) to bright yellow (10 m) to indicate slip magnitude. Each sub-graph represents an earthquake event. Red arrows indicate the direction of rupture propagation, red asterisks indicate the rupture origin (source), and yellow asterisks indicate the location of the main shock. The fault geometry is complex, including multiple bends and branching structures, indicating that earthquake ruptures can propagate along different paths. Figure 4 (b) shows a graph of earthquake period and magnitude. The bars in the graph represent the interval between each event, expressed in years (yr) or minutes (minutes). The blue bars correspond to the recurrence period of static slip events, and the red bars represent the mainshock magnitude (Mw). The text in the middle indicates that one of the events was a slow earthquake occurring on a fault branch. The magnitudes of the earthquakes are mostly concentrated in the Mw7.4 to Mw8.3 range, demonstrating the fault's capacity to produce large earthquakes and the diversity of rupture recurrence. The period lengths in the graph vary widely, ranging from tens of minutes to over 200 years, revealing the irregularity between strong earthquakes. Figure 4 (c) shows the moment rate evolution during the co-seismic time course. Each plot shows the temporal evolution of the moment rate during a seismic event, which can be used to analyze the rate of energy release. The different curve shapes reflect the rupture duration and energy release pattern of different earthquakes. For example, the first event lasts for a long time and has a high peak, representing a typical strong earthquake; the fourth earthquake lasts only tens of seconds, indicating a fast rupture event; and the sixth earthquake has multiple energy release peaks, possibly indicating a multi-stage rupture or complex event.

[0046] Based on the above disclosure, the present invention disclosed by the present invention can complete the millennial earthquake cycle simulation of complex structures such as multiple fault segments, branch faults, fault jumps, etc. within a reasonable time. Also disclosed is a computer-readable storage medium storing a computer program, which enables a computer to implement the above-mentioned earthquake hazard assessment method when executed.

[0047] Exemplarily, a computer program may be divided into one or more modules / units, one or more modules / units being stored in a memory and executed by a processor, and the input interface and output interface completing the I / O interface transmission of data to complete the present invention. One or more modules / units may be a series of computer program instruction segments capable of completing specific functions, and the instruction segments are used to describe the execution process of the computer program in a computer device.

[0048] The above description is merely a specific embodiment of the present invention, but the scope of protection of the present invention is not limited thereto. Any changes or substitutions within the technical scope disclosed by the present invention shall be covered by the scope of protection of the present invention. Therefore, the scope of protection of the present invention shall be subject to the scope of protection of the claims.

Claims

1. A method for earthquake risk assessment, characterized in that: The steps include: Receive a gridded 3D fault geometry model and set initial conditions; Construct simulation operation equations to solve shear stress loading rate, normal stress loading rate and state variables, and construct ordinary differential equation systems based on the simulation operation equations; Based on the three-dimensional geometric model and initial conditions, the ordinary differential equation system is used to realize accelerated simulation calculation of shear stress loading rate, normal stress loading rate and state variables, and the earthquake hazard is evaluated according to the simulation calculation results; Among them, the method for accelerating simulation operations includes: using a method combining H-matrix and MPI to accelerate simulation operations, including obtaining a layered matrix based on the Green function under half-space conditions 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.

2. A method for earthquake risk assessment according to claim 1, characterized in that: The setting of initial conditions includes setting initial stress state, friction parameter configuration and earthquake period calibration, The friction coefficient includes parameters and parameters , the parameters and parameters is an interval adjustable parameter; Combined with the regional characteristic earthquake cycle, the friction parameters a and b are optimized through iteration.

3. The earthquake risk assessment method according to claim 1, wherein: The method for constructing simulation operation equations includes: constructing shear stress and normal stress control equations in combination with the radiation damping hypothesis to express the elastic stress transfer caused by fault slip; using the rate-state friction law as a boundary condition, constructing a differential equation group in combination with the control equations and the boundary conditions; and merging the differential equation group into an ordinary differential equation system for solving.

4. A method for earthquake risk assessment according to claim 1, characterized in that: Calculating the slip rate and slip amount based on the calculation results of the shear stress loading rate, the normal stress loading rate, and the state variables, and calculating the earthquake occurrence time and magnitude based on the slip rate and slip amount; At the same time, based on the slip rate, slip amount, earthquake occurrence time and magnitude, the occurrence time, magnitude distribution, rupture propagation path, surface deformation pattern and spatial and temporal distribution of high-intensity areas of multiple earthquake events are simulated.

5. The earthquake risk assessment method according to claim 1, wherein: The method for constructing the gridded three-dimensional geometric model of the fault comprises: Based on the high-precision earthquake positioning results and the trajectory of the surface active fault, the spatial geometry of the fault is determined. Use CAD or GIS software to create a three-dimensional geometric model of the fault; The fault plane is divided into triangular meshes using meshing software.

6. A method for earthquake risk assessment according to claim 1, characterized in that: The accelerated simulation operation includes: The adaptive cross approximation algorithm in the H-matrix framework is used to select some key rows and columns, and iterative approximation is performed based on the main elements selected from the residual to effectively construct a low-rank representation of the matrix and obtain a hierarchical matrix.

7. An earthquake risk assessment system, characterized in that: include: An input unit, for receiving a gridded three-dimensional geometric model of the fault and setting initial conditions; The CPU simulation unit constructs simulation operation equations to solve the shear stress loading rate, normal stress loading rate and state variables, and constructs an ordinary differential equation system based on the simulation operation equations; Based on the three-dimensional geometric model and initial conditions, the ordinary differential equation system is used to realize accelerated simulation calculation of shear stress loading rate, normal stress loading rate and state variables; The method for accelerating simulation operations includes: accelerating simulation operations by combining H-matrix and MPI, including obtaining a layered matrix based on the Green's function under half-space conditions 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; An evaluation unit is used to evaluate earthquake hazards based on simulation results.

8. The earthquake risk assessment system according to claim 7, wherein: Also includes GPU simulation unit, GPU simulation unit, used to construct dense stiffness matrices based on Green's function, including dense shear stress stiffness matrix and normal stress stiffness matrix; Based on the three-dimensional geometric model and initial conditions, the ordinary differential equation system is used to realize accelerated simulation calculation of shear stress loading rate, normal stress loading rate and state variables; Among them, the simulation operation equations: Combined with the radiation damping assumption, the shear stress and normal stress control equations are constructed to express the elastic stress transfer caused by fault slip; Taking the rate-state friction law as the boundary condition, combining the control equation and the boundary condition to construct a differential equation system, and constructing the differential equation system into a solvable ordinary differential equation system; The accelerated simulation operation method: constructs a low-rank matrix based on the H-matrix and uses MPI parallel acceleration; and also provides the option of GPU parallel acceleration; The CPU simulation unit and the GPU simulation unit can be selectively applied to simulation operations.

Images