Cross-scale broadband seismic wave field simulation method and device

Abstract

The invention discloses a cross-scale broadband seismic wave field simulation method, and the method comprises the following steps: firstly, carrying out the spatial discretization of a research region based on a spectral element method, generating a three-dimensional cross-scale wave velocity structure grid model, constructing GLL integral points in the grid units of the grid model, and integrating a global seismic wave equation based on the integral points; then, converting a seismic wave equation to be solved into a second-order linear ordinary differential equation on a time domain; performing time domain discretization on the second-order linear ordinary differential equation set by adopting a fourth-order symplectic time integral algorithm; local time step form construction is carried out on the fourth-order symplectic time integral algorithm; and finally, in a discrete grid system formed by the grid model, performing multi-step iterative calculation on the physical quantity of the seismic wave field to obtain a three-component seismic oscillation numerical simulation result of the research area. The method improves the calculation efficiency and the precision of a calculation result.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of earthquake simulation, in particular to a cross-scale wide-band seismic wave field simulation method and device. BACKGROUND

[0002] At present, seismic motion simulation helps to study the propagation law of seismic motion, is the basis for studying the internal structure of the earth, and plays an important role in resource exploration and earthquake disaster assessment. In recent years, it has attracted widespread attention in the fields of earthquake engineering, seismology and geophysics. At present, a variety of mature seismic motion numerical simulation methods have been developed, mainly including finite difference method (FDM), finite element method (FEM) and spectral element method (SEM) and the like. Among them, SEM solves the weak form of seismic wave motion equation based on the variational principle, which combines the high precision of spectral method and the grid flexibility of FEM, and the mass matrix obtained by SEM is a diagonal matrix, and the inversion operation of the diagonal matrix has small calculation amount, which can obviously improve the efficiency of seismic wave field simulation.

[0003] However, when dealing with cross-scale models established by high wave velocity structure of the crust and low wave velocity structure of the near-surface, SEM still has a prominent problem: when the model is discretized by using a fixed size of structured hexahedral grid, it is necessary to use the smallest time step to meet the stability of SEM simulation calculation, which greatly reduces the solving efficiency of the method. Moreover, in order to meet the demand of seismic motion input covering the sensitive frequency band of engineering structure, the frequency to be simulated by seismic motion is constantly improved, resulting in smaller and smaller grid size and larger and larger grid quantity, and the time sampling interval is shortened, which inevitably leads to exponential increase of the calculation amount of long-time seismic motion. Therefore, under the condition of simulation time, the increase of simulation frequency will inevitably lead to the increase of total iteration number. At present, the limitation of traditional explicit time integration method in spectral element program gradually appears, and it is urgent to develop a time integration algorithm considering accuracy and efficiency. SUMMARY

[0004] The present application provides a cross-scale wide-band seismic wave field simulation method and device to solve the technical problems in the prior art.

[0005] The technical scheme adopted by the present application to solve the technical problems in the prior art is as follows:

[0006] A cross-scale wide-band seismic wave field simulation method, the method comprising the following steps:

[0007] Step 1, spatially discretizing the research area based on the spectral element method to generate a three-dimensional cross-scale wave velocity structure grid model, constructing GLL integration points in the grid cells of the grid model, and integrating the global seismic wave equation based on the integration points;

[0008] Step 2, the seismic wave equation to be solved is converted into a second-order linear constant differential equation in the time domain;

[0009] Step 3, a fourth-order symplectic time integration algorithm is used to discretize the second-order linear constant differential equation in the time domain;

[0010] Step 4, the fourth-order symplectic time integration algorithm is constructed in a local time step form according to different grid cell sizes and their corresponding maximum stable time steps determined by the wave velocity distribution of the grid model;

[0011] Step 5, in the discrete grid system formed by the grid model, multi-step iterative calculation is performed on the seismic wave field physical quantities to obtain the numerical simulation results of three-component ground motion in the study area.

[0012] Further, in step 1, according to the SEM precision constraint on the grid cell size, combined with the medium wave velocity, simulation frequency and GLL integration point density of the grid cell, the study area is spatially discretized.

[0013] Further, in step 2, the seismic wave equations of all grid cells are assembled to obtain the following global second-order linear constant differential equation in the time domain:

[0014] ;

[0015] In the formula:

[0016] M is the global mass matrix of the grid cell;

[0017] C is the global damping matrix of the grid cell;

[0018] K is the global stiffness matrix of the grid cell;

[0019] F is the global load vector of the grid cell;

[0020] A is the global acceleration vector of the grid cell;

[0021] V is the global velocity vector of the grid cell;

[0022] U is the global displacement vector of the grid cell.

[0023] Further, in step 3, the calculation process in each time step is divided into four intermediate steps. In each intermediate step, the displacement, velocity and acceleration updates in the seismic wave field physical quantities follow the mathematical form of the symplectic algorithm, and the acceleration, velocity and displacement are updated alternately by calling the matrix solution of the intermediate step.

[0024] Further, the entire time domain in each time step is uniformly divided into equally spaced intermediate steps; the initial time is calculated, and the solution is recursively solved in units of set time steps; in each time step, the updates of displacement, velocity, and acceleration are as follows:

[0025] ;

[0026] ;

[0027] ;

[0028] ;

[0029] ;

[0030] ;

[0031] ;

[0032] ;

[0033] ;

[0034] In the formula:

[0035] ξ, ψ, λ are calculation coefficients, which are constants;

[0036] Δt represents the time step, which is the discrete time interval for calculation or state update, and the unit is second;

[0037] U(t) represents the displacement at time t;

[0038] V(t) represents the velocity at time t;

[0039] U(t+Δt) represents the displacement at time t+Δt;

[0040] V(t+Δt) represents the velocity at time t+Δt;

[0041] n represents the intermediate step number in the time from t to t+Δt;

[0042] U n represents the displacement at the n-th intermediate step in the time from t to t+Δt;

[0043] V n represents the velocity at the n-th intermediate step in the time from t to t+Δt.

[0044] Further, in step 4, for the multi-size grid elements in the cross-scale wave velocity structure grid model, the maximum stable time step required for each grid element is evaluated using the following formula:

[0045] ;

[0046] wherein:

[0047] Δt cell is the maximum stable time step of the unit; represents the time step requirement to ensure the accuracy of the calculation results of the unit;

[0048] C CFL is the conditional stability number of the fourth-order symplectic time integration algorithm, and is a constant representing the calculation convergence condition of the algorithm;

[0049] x is the edge length of a hexahedral grid unit;

[0050] V P is the P-wave velocity of the grid unit;

[0051] γ is the interval coefficient of the GLL integration points in the grid unit, and is a constant representing the density of the GLL integration points.

[0052] Further, according to different grid sizes determined by the model wave velocity distribution and their corresponding maximum stable time steps Δt cell , a grid grading system based on the size of the time step is established; the grid of the entire calculation domain is divided into different levels, the higher the level, the smaller the corresponding time step, the level of each grid is marked, and the calculation sub-region covered by each level of grid is defined; after the entire calculation domain is divided into sub-regions according to different time steps, the grid nodes on the boundary between different level sub-regions are determined as the basis for data exchange between different level regions.

[0053] Further, in step 5, based on the principle of gradually advancing from the high-level step grid region to the lower-level step grid region, in each global time step, the seismic wave field physical quantity of the nodes in the higher-level step grid region is updated first, and then the initial data is transmitted through the boundary nodes to update the seismic wave field physical quantity of the nodes in the lower-level step grid region, until all the calculation processes of the global calculation time step are completed.

[0054] Further, it is assumed that the step level of the grid unit region of the entire calculation domain is divided into N levels, k is the step level serial number of the grid unit region; the step of the kth level grid unit region is Δt k , the step of the k+1th level grid unit region is Δt k+1 , Δt k =S k Δt k+1 ; S kS k is the ratio of the step length of the kth grid cell region to the step length of the (k+1)th grid cell region, which is a positive integer greater than 1; let the global time step = S 1 x S 2 x…x S k N-1 Δt N =Δt1;Δt N is the step length of the Nth grid cell region;

[0055] According to the following steps A to step B, the seismic wave field physical quantity of all nodes in a global time step is updated from the Nth grid cell region to the first grid cell region:

[0056] Step A, obtain the seismic wave field physical quantity of all nodes in the (k+1)th grid cell region at time t; use Δt k+1 time step to update the seismic wave field physical quantity of all nodes in the (k+1)th grid cell region, and go through S k Δt k+1 time steps; after the update of the Δt k time step, the seismic wave field physical quantity of all nodes in the (k+1)th grid cell region is updated to the state value at t+Δt k time, and save the state value of the seismic wave field physical quantity of all nodes in the (k+1)th grid cell region at t+Δt k time; at the same time, the corresponding region grid node stiffness matrix and boundary matrix are transmitted and updated from the (k+1)th grid cell region to the kth grid cell region through the boundary nodes;

[0057] Step B, obtain the seismic wave field physical quantity of all nodes in the kth grid cell region at time t; use Δt k time step to update the seismic wave field physical quantity of all nodes in the kth grid cell region, and make the seismic wave field physical quantity of all nodes in the kth grid cell region updated to t+Δt k time, save the state value of the seismic wave field physical quantity of all nodes in the kth grid cell region at t+Δt k time.

[0058] The application also provides a device for simulating a cross-scale wide-band seismic wave field, which comprises a memory and a processor, wherein the memory is used for storing a computer program; and the processor is used for executing the computer program and realizing the steps of the cross-scale wide-band seismic wave field simulation method as described above when executing the computer program.

[0059] The application has the following advantages and positive effects:

[0060] Improve the calculation efficiency: adopt the LTS form to solve the low efficiency problem caused by using the minimum time step in the global when the traditional SEM is used to process the cross-scale model.

[0061] Results high precision: using SA4 algorithm, the characteristics of structure-preserving symplectic integration makes the calculation results can still maintain high precision in wideband long time simulation of ground motion. BRIEF DESCRIPTION OF DRAWINGS

[0062] Figure 1 It is a cross-scale wideband seismic wave field simulation method workflow schematic diagram of the application.

[0063] Figure 2 It is a schematic diagram of the calculation area and its relative position with the seismogenic fault in the embodiment of the application.

[0064] Figure 3 It is a schematic diagram of three-dimensional cross-scale calculation model grid division in the embodiment of the application.

[0065] Figure 4 It is a schematic diagram of three-component snapshots of seismic wave field at different times in the embodiment of the application.

[0066] In the figure:

[0067] Figure 2 And Figure 4 The abscissa is longitude, and the ordinate is latitude.

[0068] Figure 4 In which, t is different earthquake time.

[0069] TOPO is the abbreviation of Topography, and the Chinese interpretation is topographic elevation data.

[0070] Vel is the abbreviation of Velocity, and the Chinese interpretation is ground motion velocity. DETAILED DESCRIPTION

[0071] The application will be described in detail below with reference to the accompanying drawings and in conjunction with the embodiments, and it should be understood that the preferred embodiments described herein are only used to illustrate and explain the application, and are not used to limit the application.

[0072] The Chinese interpretation of the following English words, abbreviations and phrases is as follows:

[0073] SEM is the abbreviation of Spectral-element method, and the Chinese interpretation is spectral element method.

[0074] GLL is the abbreviation of Gauss-Lobatto-Legendre, and the Chinese interpretation is Gauss-Lobatto-Legendre integral.

[0075] LTS is the abbreviation of Locall Time-Step, and the Chinese interpretation is local time step.

[0076] SA4 is the abbreviation of Fourth-order Symplectic algorithm, which means the fourth-order symplectic time integration algorithm in Chinese.

[0077] Please refer to Figures 1 to 4 A cross-scale wide-band seismic wave field simulation method, the method comprising the following steps:

[0078] Step 1, based on the spectral element method, the spatial domain of the research area is discretized, a three-dimensional cross-scale wave velocity structure grid model is generated, GLL integration points are constructed in the grid cells of the grid model, and the global seismic wave equation is integrated based on these integration point sets.

[0079] Step 2, the seismic wave equation to be solved is converted into a second-order linear constant differential equation in the time domain.

[0080] Step 3, a fourth-order symplectic time integration algorithm is used to discretize the second-order linear constant differential equation in the time domain.

[0081] Step 4, according to the different grid cell sizes determined by the wave velocity distribution of the grid model and the corresponding maximum stable time step, the fourth-order symplectic time integration algorithm is constructed in the form of local time step.

[0082] Step 5, in the discrete grid system composed of the grid model, the seismic wave field physical quantity is calculated by multi-step iteration, and the numerical simulation result of three-component ground motion in the research area is obtained.

[0083] The interpolation basis function selected by the spectral element method is the Lagrange basis function, and the finite element method uses the linear interpolation basis function. Inside the grid cell, the spectral element method uses the gauss lobatto legendre (abbreviated as gll) integration scheme, and the mass matrix and stiffness matrix are converted into a combination of weight coefficients. Due to the orthogonality of the polynomial, the final mass matrix is a diagonal matrix.

[0084] The spectral element method uses Guass-Lobatto-Legendre (GLL) integration to obtain a diagonal mass matrix, which is widely used in seismology due to its high accuracy and efficiency, and can naturally meet the free boundary condition and explicit finite element method. By using the mass concentration method to obtain a diagonal mass matrix, it has the same high efficiency as the spectral element method.

[0085] The symplectic integration method ensures the long-term stability and accuracy of the numerical solution by preserving the symplectic geometric structure of the system (i.e. the area-invariant property in the phase space), which is particularly suitable for conservative systems.

[0086] Fourth-order symplectic time integration algorithm is mainly used in numerical calculation scenarios that need to keep system energy conservation, especially in the fields of dynamics system, celestial mechanics, material science, etc. Fourth-order symplectic algorithm effectively solves the problem of energy loss in long-time calculation of traditional numerical integration method by keeping the symplectic structure (i.e. long-term stability of Hamilton structure) of the system. Its precision and stability are better than low-order algorithms, and it is suitable for processing highly nonlinear systems.

[0087] According to the different grid cell sizes determined by the model wave velocity distribution and their corresponding maximum stable time steps, a grid cell hierarchical system based on time step size is established. Different grid cell sizes have different iteration time steps, which can improve the time efficiency of calculation.

[0088] Preferably, in step 1, the spatial discretization of the study area can be performed according to the constraint of SEM precision on grid cell size, combined with medium wave velocity, simulation frequency and grid cell GLL integration point density.

[0089] Preferably, in step 2, the seismic wave equation of all grid cells can be assembled to obtain the following time-domain second-order linear ordinary differential equation globally:

[0090] ;

[0091] In the formula:

[0092] M is the global mass matrix of the grid cell;

[0093] C is the global damping matrix of the grid cell;

[0094] K is the global stiffness matrix of the grid cell;

[0095] F is the global load vector of the grid cell;

[0096] A is the global acceleration vector of the grid cell;

[0097] V is the global velocity vector of the grid cell;

[0098] U is the global displacement vector of the grid cell.

[0099] Preferably, in step 3, the calculation process within each time step can be divided into four intermediate steps. In each intermediate step, the displacement, velocity and acceleration updates in the seismic wave field physical quantity follow the mathematical form of symplectic algorithm, and the matrix is called to solve the intermediate step acceleration, velocity and displacement alternating step update.

[0100] Preferably, the entire time domain within each time step can be evenly divided into equidistant intermediate steps; the calculation is recursively solved from the initial time with a set time step; in each time step, the update of displacement, velocity and acceleration can be shown as follows:

[0101] ;

[0102] ;

[0103] ;

[0104] ;

[0105] ;

[0106] ;

[0107] ;

[0108] ;

[0109] ;

[0110] In the formula:

[0111] ξ, ψ, λ are calculation coefficients, which are constants;

[0112] Δt represents the time step, which is the discrete time interval for calculation or state update, and the unit is second;

[0113] U(t) represents the displacement at time t;

[0114] V(t) represents the velocity at time t;

[0115] U(t+Δt) represents the displacement at time t+Δt;

[0116] V(t+Δt) represents the velocity at time t+Δt;

[0117] n represents the intermediate step number within t to t+Δt time;

[0118] U n represents the displacement at the n-th intermediate step within t to t+Δt time;

[0119] V n represents the velocity at the n-th intermediate step within t to t+Δt time.

[0120] Preferably, in step 4, for multi-size grid elements in the cross-scale wave velocity structure grid model, the following formula can be used to evaluate the maximum stable time step required for each grid element:

[0121] ;

[0122] wherein:

[0123] Δt cell is the maximum stable time step of the element; represents the time step requirement to ensure the accuracy of the calculation results of the element;

[0124] C CFL is the conditional stability number of the fourth-order symplectic time integration algorithm, and is a constant representing the calculation convergence condition of the algorithm;

[0125] x is the edge length of a hexahedral grid element;

[0126] V P is the P-wave velocity of the grid element;

[0127] γ is the interval coefficient of the GLL integration points in the grid element, and is a constant representing the density of the GLL integration points.

[0128] Preferably, a grid hierarchical system based on the size of the time step can be established according to different grid sizes determined by the model wave velocity distribution and their corresponding maximum stable time steps Δt cell , the grid of the entire calculation domain is divided into different levels, the higher the level, the smaller the time step, the level to which each grid belongs is marked, and the calculation sub-area covered by each level of grid is defined; after the entire calculation domain is divided into sub-areas according to different time steps, the grid nodes on the boundary between different level sub-areas are determined as the basis for data exchange between different level areas.

[0129] Preferably, in step 5, based on the principle of gradually advancing from the high-level step grid area to the lower-level step grid area, in each global time step, the seismic wave field physical quantity of the nodes in the calculation of the higher-level step grid area is updated first, and then the initial data is transmitted through the boundary nodes to update the seismic wave field physical quantity of the nodes in the calculation of the lower-level step grid area, until all the calculation processes of the global calculation time step are completed.

[0130] Preferably, the step length level of the grid element area of the entire calculation domain can be divided into N levels, k is the step length level number of the grid element area; the step length of the kth level grid element area is Δt k , the step length of the k+1th level grid element area is Δt k+1 , Δt k =S k Δt k+1 ; S kS k+1 is the step of the k+1 level grid cell region, which is a positive integer greater than 1; let the global time step = S 1 x S 2 x … x S k+1 N-1 Δt N =Δt1;Δt N S N is the step of the N level grid cell region.

[0131] The seismic wave field physical quantity of all nodes in a global time step can be updated from the N level grid cell region to the 1 level grid cell region according to the following steps A to B:

[0132] Step A, obtaining the seismic wave field physical quantity of all nodes in the k+1 level grid cell region at time t; using Δt k+1 S k+1 time steps to update the seismic wave field physical quantity of all nodes in the k+1 level grid cell region, and the k+1 level grid cell region experiences S k Δt k+1 time steps, the seismic wave field physical quantity of all nodes in the k+1 level grid cell region is updated to the state value at time t+Δt k t+Δt k , and the state value of the seismic wave field physical quantity of all nodes in the k+1 level grid cell region at time t+Δt k is saved; at the same time, the corresponding region grid node stiffness matrix and boundary matrix are transmitted and updated from the k+1 level grid cell region to the k level grid cell region through the boundary nodes.

[0133] Step B, obtaining the seismic wave field physical quantity of all nodes in the k level grid cell region at time t; using Δt k S k time steps to update the seismic wave field physical quantity of all nodes in the k level grid cell region, so that the seismic wave field physical quantity of all nodes in the k level grid cell region is updated to time t+Δt k , and the state value of the seismic wave field physical quantity of all nodes in the k level grid cell region at time t+Δt k is saved.

[0134] The application also provides a device of the cross-scale wide frequency band seismic wave field simulation method, which comprises a memory and a processor, the memory is used for storing a computer program; the processor is used for executing the computer program and realizing the steps of the cross-scale wide frequency band seismic wave field simulation method as described above when executing the computer program.

[0135] The working principle and working process of the application will be further described below with reference to a preferred embodiment of the application:

[0136] As Figure 1 shown, a cross-scale wide frequency band seismic wave field simulation method, the method comprising the following steps:

[0137] Firstly, based on the spectral element method, the research area is spatially discrete, a three-dimensional cross-scale wave speed structure grid model is generated, GLL integration points are constructed in the grid cells of the grid model, and the global seismic wave equation is integrated based on these integration points; then, the seismic wave equation to be solved is converted into a second-order linear constant differential equation in the time domain; then, a fourth-order symplectic time integration algorithm is used to discretize the second-order linear constant differential equation in the time domain; and the fourth-order symplectic time integration algorithm is constructed in the form of local time step; finally, in the discrete grid system formed by the grid model, the seismic wave field physical quantity is calculated by multi-step iteration, and the three-component ground motion numerical simulation result of the research area is obtained.

[0138] The cross-scale wide frequency band seismic wave field simulation method provided by the application integrates LTS, SA4 and various algorithms, as shown in the figure, the generated grid model has been spatially discrete based on SEM solving, GLL integration points are constructed in the grid model to convert the wave equation to be solved into a second-order linear constant differential equation in the time domain; the SA4 algorithm is used to discretize the second-order linear constant differential equation in the time domain; then, the SA4 algorithm is constructed in the form of LTS, and the stress, strain, acceleration, velocity and displacement are calculated by multi-step iteration on the entire discrete grid system, and the three-component ground motion numerical simulation result of the region is obtained. The specific method steps are as follows: Figure 1 (1) Grid generation and equation construction:

[0139] Firstly, according to the requirement of SEM accuracy on grid size, combined with medium wave speed, simulation frequency and unit GLL integration point density, the generated three-dimensional cross-scale wave speed structure grid model of the research area is constructed, and then GLL integration points are constructed in the grid model, and the global equation is integrated based on these integration points, so that the wave equation to be solved is converted into a second-order linear constant differential equation in the time domain about the whole displacement vector U, please refer to the following formula:

[0140]

[0141] ;

[0142] In the formula:

[0143] M is the global mass matrix of the grid cell;

[0144] C is the global damping matrix of the grid cell;

[0145] K is the global stiffness matrix of the grid cell;

[0146] F is the global load vector of the grid cell;

[0147] A is the global acceleration vector of the grid cell;

[0148] ​V is the global velocity vector of the mesh element;

[0149] U is the global displacement vector of the mesh element.

[0150] V is the first-order time derivative of U, and A is the second-order time derivative of U.

[0151] This step simplifies the complex wave problem into a standard solvable form, making it easier to construct a fourth-order explicit time integration algorithm and directly use numerical calculation methods to solve the time-domain evolution results of the overall displacement vector U.

[0152] (2) Time integration method:

[0153] The fourth-order explicit SA4 algorithm is used to discretize the above second-order linear ordinary differential equation set in the time domain. The specific process is as follows: the entire time domain is uniformly divided into equal time steps (the length of the interval is the time step size Δt). The calculation starts from the initial time (t=0) and recursively solves in units of selected time step size (Δt). In each time step, the update of field variables (displacement, velocity, acceleration) strictly follows the mathematical form of the symplectic algorithm, and the specific numerical calculation process is as follows:

[0154] ;

[0155] ;

[0156] ;

[0157] ;

[0158] ;

[0159] ;

[0160] ;

[0161] ;

[0162] ;

[0163] In the formula:

[0164] ξ, ψ, λ are calculation coefficients, which are constants;

[0165] Δt represents the time step size, which is the discrete time interval for calculation or state update, and the unit is second;

[0166] U(t) represents the displacement at time t;

[0167] V(t) represents the velocity at time t;

[0168] U(t+Δt) represents the displacement at t+Δt;

[0169] V(t+Δt) represents the velocity at t+Δt;

[0170] n represents the intermediate step number within t to t+Δt;

[0171] U n represents the displacement of the nth intermediate step within t to t+Δt;

[0172] V n represents the velocity of the nth intermediate step within t to t+Δt.

[0173] U n and V n are the displacement and velocity of the intermediate step, n=1,2,3,4; the subscript number n represents the nth intermediate step within t to t+Δt. The SA4 method is divided into 4 intermediate steps in one time step, and the core of each intermediate step is to call the matrix to solve the acceleration of the intermediate step (i.e. ), and an alternating step update strategy of velocity and displacement is adopted, which avoids the storage dependence of a large amount of intermediate data in single-step calculation, thereby reducing the memory and parallel reading requirements of the solving program, and providing an efficient and accurate time integration method for cross-scale wide-band seismic wave field simulation.

[0174] (3) Constructing a local time step method:

[0175] Based on the above SA4 algorithm, an LTS form is further constructed to obtain a more efficient cross-scale model wide-band seismic wave field simulation method. The specific construction process is as follows: first, for the multi-size grid elements generated by the multi-scale wave velocity model, the maximum stable time step required for each grid element is evaluated using the following formula.

[0176] ;

[0177] In the formula:

[0178] Δt cell is the maximum stable time step of the element; represents the time step requirement to ensure the accuracy of the calculation results of the element;

[0179] C CFL is the conditional stability number of the fourth-order symplectic time integration algorithm, which is a constant representing the calculation convergence condition of the algorithm;

[0180] x is the edge length of the hexahedral grid element;

[0181] V P is the P-wave velocity of the grid element;

[0182] γ is the interval coefficient of GLL integration points in the grid element, which is a constant representing the density of GLL integration points.

[0183] Through this calculation, it can be prepared for subsequent time step grading according to grid size.

[0184] Then, according to the different grid sizes determined by the model wave velocity distribution and their corresponding maximum stable time steps Δt cell , a grid grading system based on time step size is established. The grid of the entire calculation domain is divided into different levels (the higher the level, the smaller the time step), and each grid belonging to the level is marked, and the calculation sub-area covered by each level grid is defined. After the sub-area division of the entire calculation domain according to different time steps is realized, the boundaries between different level sub-areas are identified, and the specific grid nodes on these area boundaries are clearly marked as the basis for data exchange between different level areas.

[0185] Finally, a solution process within a time step needs to be constructed to realize the recursive solution of the model from the initial time (t=0) to the end of the simulation. The core principle of this process is: in each time step, first complete the calculation of the high-level step (small time step) grid area, and then perform the calculation of the low-level step (large time step) grid area. For simplicity, this description will take a case of a two-level grid system (N=2) as an example to elaborate the solution process. For a multi-level grid system, the calculation process between two-level grids and three-level grids can be constructed according to the correlation logic between two-level grids and one-level grids. The specific process of the two-level grid system is as follows:

[0186] Let the global time step = Δt, S1=2, Δt2=1 / 2*Δt, Δt1=Δt.

[0187] Step S1, obtain the displacement and velocity values of all nodes in the one-level grid area at the current time, i.e. t time, and update the displacement and velocity values of the nodes at t time through the boundary grid nodes to the two-level grid.

[0188] Step S2, the one-level grid area needs to be calculated first, and the displacement and velocity values of all nodes in the two-level grid area at t time are obtained, and the displacement and velocity values of all nodes in the two-level grid area are updated using Δt / 2 time step according to formula (2). After two Δt / 2 time steps, the velocity and displacement values of all nodes in the two-level grid area are updated to t+Δt, and the acceleration, velocity, and displacement results at this time are saved. At the same time, the stiffness matrix K and the boundary matrix C of the relevant regional grid nodes are transmitted from the two-level grid to the one-level grid through the boundary nodes and updated. Complete the calculation of the two-level grid at t+Δt.

[0189] Step S3, time domain iteration calculation is carried out on the primary grid area, displacement and velocity values of all nodes in the primary grid area at t time are obtained, all node velocity and displacement values in the primary grid area are updated using formula (2) using Δt time step, and the velocity and displacement values of all nodes in the primary grid area are updated to t+Δt time. The acceleration, velocity and displacement results of the primary grid area at this time are saved and combined with the results in the secondary area to output the t+Δt time data of the entire model, and the solving process of the entire model from t time to t+Δt time, i.e. a time step, is completed.

[0190] A mountain once had a 6.2 magnitude earthquake. Through the comparison of the damage degree of houses, the number of casualties and the distribution of peak ground acceleration of different earthquakes of the same magnitude in history, it is concluded that this earthquake shows the characteristics of "same earthquake disaster" compared with the same level earthquake in the past. In order to analyze the main reason for the serious earthquake damage caused by this earthquake, the method of the present application is used to simulate the wideband seismic response of the surrounding area when the earthquake occurs.

[0191] The method of the present application is used to simulate the seismic wave field, the simulation calculation area of this time is shown in Figure 2 According to the target area, the calculation model along the east-west direction is about 111 km long, the north-south direction is about 108 km wide, and the longitudinal direction is about 35 km. The three-dimensional cross-scale velocity structure model of the mountain is discretized in multiple scales, and the model grid size is expanded from 50 m at the surface to 200 m at the bottom of the model, as shown in Figure 3 . Figure 4 The snapshots of the seismic wave field at different times of the earthquake are shown in the figure. Near the ridge and canyon and other complex topography, the waveforms show obvious deformation and distortion, which indicates that near the complex topography, the seismic wave can cause multiple reflections, which may cause significant peak amplification of the seismic ground motion in the Jishishan area.

[0192] The above-described embodiments are only used to illustrate the technical ideas and characteristics of the present application, and the purpose is to enable those skilled in the art to understand the content of the present application and implement it, and cannot be limited to the patent scope of the present application by the above-described embodiments, i.e. any equivalent changes or modifications made according to the spirit disclosed by the present application still fall within the patent scope of the present application.

Claims

1. A method of cross-scale broadband seismic wavefield simulation, characterized in that, The method comprises the following steps: Step 1, spatially discretize the research area based on the spectral element method, generate a three-dimensional cross-scale wave velocity structure grid model, construct GLL integration points in the grid cells of the grid model, and integrate the global seismic wave equation based on the integration points; Step 2, convert the seismic wave equation to be solved into a second-order linear constant differential equation in the time domain; Step 3, use a fourth-order symplectic time integration algorithm to discretize the second-order linear constant differential equation in the time domain; Step 4, according to the different grid cell sizes determined by the wave velocity distribution of the grid model and the corresponding maximum stable time step, construct the fourth-order symplectic time integration algorithm in the form of local time step; Step 5, in the discrete grid system formed by the grid model, multi-step iterative calculation is performed on the seismic wave field physical quantity to obtain the numerical simulation result of three-component ground motion in the research area.

2. The cross-scale broadband seismic wavefield simulation method of claim 1, wherein, In step 1, according to the constraint of SEM accuracy on the size of the grid cell, combined with the medium wave velocity, the simulation frequency and the GLL integration point density of the grid cell, the research area is spatially discretized.

3. The cross-scale broadband seismic wavefield simulation method of claim 1, wherein, In step 2, the seismic wave equation of all grid cells is assembled to obtain the following second-order linear constant differential equation in the time domain: ； In the formula: M is the global mass matrix of the grid cell; C is the global damping matrix of the grid cell; K is the global stiffness matrix of the grid cell; F is the global load vector of the grid cell; A is the global acceleration vector of the grid cell; V is the global velocity vector of the grid cell; U is the global displacement vector of the grid cell.

4. The cross-scale broadband seismic wavefield simulation method of claim 3, wherein, In step 3, the calculation process in each time step is divided into four intermediate steps. In each intermediate step, the displacement, velocity and acceleration updates of the seismic wave field physical quantity follow the mathematical form of the symplectic algorithm, and the acceleration, velocity and displacement are updated alternately by calling the matrix solution of the intermediate step.

5. The cross-scale broadband seismic wavefield simulation method of claim 4, wherein, The entire time domain in each time step is uniformly divided into equidistant intermediate steps; the calculation is recursively solved from the initial time with a set time step; in each time step, the displacement, velocity and acceleration updates are as follows: ； ； ； ； ； ； ； ； ； In the formula: ξ, ψ, λ are calculation coefficients, which are constants; Δt represents the time step, which is the discrete time interval for calculation or state update, and the unit is second; U(t) represents the displacement at time t; V(t) represents the velocity at time t; U(t+Δt) represents the displacement at time t+Δt; V(t+Δt) represents the velocity at time t+Δt; n represents the intermediate step number in the time from t to t+Δt; U n denotes the displacement of the n-th intermediate step in the time interval t to t+Δt; V n Vn(t) represents the velocity of the nth intermediate step in the time interval t to t+Δt.

6. The cross-scale broadband seismic wavefield simulation method of claim 1, wherein, In step 4, for the multi-size grid cells in the cross-scale wave velocity structure grid model, the maximum stable time step required by each grid cell is evaluated using the following formula: ； In the formula: Δt cell is the maximum stable time step for the cell; represents the time step requirement to guarantee the accuracy of the calculation results of the cell; C CFL C is the condition number of the fourth-order symplectic time integration algorithm, and is a constant representing the computational convergence condition of the algorithm. x is the edge length of the hexahedral grid cell; V P P is the P-wave velocity for the grid cell; γ is the interval coefficient of the GLL integration points in the grid cell, which is a constant representing the density of the GLL integration points.

7. The cross-scale broadband seismic wavefield simulation method of claim 6, wherein, Different grid sizes and their corresponding maximum stable time steps Δt determined by the model wave speed distribution cell , a grid hierarchy system based on time step size is established; the grids of the entire calculation domain are divided into different levels, the higher the level, the smaller the corresponding time step, and the level to which each grid belongs is marked to define the calculation sub-area covered by each level grid; after the entire calculation domain is divided into sub-areas according to different time steps, the grid nodes on the boundaries between sub-areas of different levels are determined as the basis for data exchange between sub-areas of different levels.

8. The cross-scale broadband seismic wavefield simulation method of claim 7, wherein, In step 5, based on the principle of gradually advancing from high-level step grid area to lower-level step grid area, in each global time step, first update the seismic wave field physical quantity of the node of the higher-level step grid area, and then transfer the initial data through the boundary node to update the seismic wave field physical quantity of the node of the lower-level step grid area, until all the calculation processes of the global time step are completed.

9. The cross-scale broadband seismic wavefield simulation method of claim 8, wherein: The step level of the grid cell region of the whole calculation domain is divided into N levels, k is the step level serial number of the grid cell region; the step of the kth level grid cell region is Δt k , the step of the k+1th level grid cell region is Δt k+1 , Δt k =S k Δt k+1 ; S k is the ratio of the step of the kth level grid cell region to the step of the k+1th level grid cell region, and is a positive integer greater than 1; Let global time step = S1 x S2 x... x SN N-1 Δt N =Δt1;Δt N is the step size for the Nth level grid cell region; According to the following steps A to B, the seismic wave field physical quantity of all nodes in a global time step is updated from the grid cell area of the Nth level to the grid cell area of the 1st level: Step A, obtain the seismic wave field physical quantity of all nodes in the k+1 level grid cell region at time t; use Δt k+1 time step to update the seismic wave field physical quantity of all nodes in the k+1 level grid cell region, experience S k Δt k+1 time step updates, the seismic wave field physical quantity of all nodes in the k+1 level grid cell region is updated to the state value at t+Δt k time, save the seismic wave field physical quantity state value of all nodes in the k+1 level grid cell region at t+Δt k time; at the same time, the stiffness matrix and the boundary matrix of the corresponding region grid nodes are transmitted and updated from the k+1 level grid cell region to the k level grid cell region through the boundary nodes; Step B, obtaining the seismic wave field physical quantity of all nodes in the kth level grid cell region at time t; using Δt k updating the seismic wave field physical quantity of all nodes in the kth level grid cell region by time step, so that the seismic wave field physical quantity of all nodes in the kth level grid cell region is updated to t+Δt k time, saving the state value of the seismic wave field physical quantity of all nodes in the kth level grid cell region at t+Δt k time.

10. An apparatus of a cross-scale broadband seismic wavefield simulation method, comprising a memory and a processor, characterized in that, The memory is configured to store a computer program; and the processor is configured to execute the computer program and implement the cross-scale wide-band seismic wave field simulation method steps in any one of claims 1 to 9 when executing the computer program.

Images