A method for simulating seismic wave based on SBP-SAT block staggered grid
By introducing the SBP-SAT method into the block staggered grid seismic wave simulation, the problems of wasted computational resources and numerical instability in traditional methods are solved, energy conservation and stability coupling are achieved, and the accuracy and efficiency of acoustic wave forward modeling in complex media are improved.
Patent Information
- Authority / Receiving Office
- CN · China
- Patent Type
- Applications(China)
- Current Assignee / Owner
- CHINA UNIV OF MINING & TECH
- Filing Date
- 2026-03-20
- Publication Date
- 2026-06-23
AI Technical Summary
Traditional staggered grid seismic wave simulation methods suffer from computational resource waste and numerical instability when dealing with highly inhomogeneous media. In particular, energy divergence and high-frequency oscillations at the interface make it difficult to achieve energy conservation and stable coupling. Furthermore, perfect matching layers under high-order difference schemes are prone to numerical drift.
A block staggered grid seismic wave simulation method based on SBP-SAT is adopted. By applying SAT boundary conditions to the physical variable equations and synchronously injecting the SAT boundary correction terms into the auxiliary differential equations of PML, SBP energy-compatible interpolation and projection operators are designed, and interface SAT penalty terms are constructed to achieve energy conservation and stability coupling and eliminate numerical drift in long-term simulations.
It significantly improves the accuracy and efficiency of acoustic forward modeling in complex media, ensures the numerical stability of long-term seismic monitoring and multiple wave imaging, reduces the waste of computing resources, and improves computing efficiency.
Smart Images
Figure CN122260434A_ABST
Abstract
Description
Technical Field
[0001] This invention relates to the field of earthquake numerical simulation technology, specifically a block-staggered grid seismic wave simulation method based on summation by parts (SBP) and simultaneous approximation term (SAT). Background Technology
[0002] As seismic exploration advances towards deeper, more complex structures, and wider bandwidths, imaging techniques such as full-waveform inversion and reverse-time migration place extremely high demands on the accuracy and efficiency of forward modeling engines. The staggered-grid finite-difference method, with its advantages of low-frequency dispersion and computational efficiency, remains the mainstream tool for solving acoustic and elastic wave equations. However, when dealing with highly inhomogeneous media (such as those containing low-velocity weathering layers, salt domes, or fault fracture zones), traditional uniform Cartesian grids face severe efficiency bottlenecks. To satisfy the stability conditions (CFL conditions) of local low-velocity bodies and the sampling theorem, the global grid step size is often limited to the minimum wave velocity, leading to severe oversampling in high-velocity regions and resulting in a significant waste of computational resources.
[0003] To address the problem of wasted computational resources, existing techniques often employ non-uniform meshing or block mesh refinement strategies. By deploying meshes of varying densities in different velocity regions, the total number of degrees of freedom is reduced, thereby improving computational efficiency. However, the core challenge of non-uniform meshing lies in the numerical handling of the coarse-to-fine mesh interface. Traditional interpolation methods (such as Lagrange interpolation or overlapping mesh methods) often disrupt the global antisymmetry of the difference operator when processing the interface, leading to energy divergence during long-term simulations and easily causing high-frequency oscillations or even numerical explosions. Especially in staggered mesh architectures, where velocity and stress are defined on different mesh nodes, the exchange of interface fluxes becomes more complex, making it difficult to guarantee energy conservation and obtain strictly stable coupling conditions using traditional interpolation methods.
[0004] Furthermore, the treatment of truncated boundaries in the computational domain is also a key factor limiting the stability of long-term simulations. Currently, complex frequency-shifted perfectly matched layers (CPMLs) and their auxiliary differential equations (ADEs) are commonly used to absorb boundary reflected waves. However, when high-order difference schemes are combined with perfectly matched layers, the lack of closure conditions for the auxiliary differential equations at the boundaries leads to incompatibility between auxiliary and physical variables at the discrete level. This results in non-physical linear or exponential error growth, which is prone to the "PML late instability" phenomenon. This instability is unacceptable in practical applications involving long-term monitoring or multiple-wave imaging.
[0005] Therefore, under the block staggered grid architecture, how to provide a new numerical simulation method for seismic waves that can achieve energy conservation and stable coupling at the interface between coarse and fine grids, and theoretically eliminate the long-term numerical drift problem of perfectly matched layers under high-order schemes, and ultimately significantly improve the accuracy and efficiency of acoustic wave forward modeling in complex media, is a technical problem that urgently needs to be solved in this field. Summary of the Invention
[0006] To address the problems existing in the prior art, this invention provides a block staggered grid seismic wave simulation method based on SBP-SAT. This method, based on the summation-partial integration SBP operator and the simultaneous approximation term SAT technique, not only achieves energy conservation and stable coupling at the interface between non-uniform coarse and fine grids, but also theoretically eliminates the numerical drift problem in long-term simulations by synchronously injecting SAT boundary correction terms into the auxiliary differential equations of the perfectly matched layer (PML), effectively improving the accuracy and efficiency of acoustic wave forward modeling in complex media.
[0007] To achieve the above objectives, the technical solution adopted by this invention is: a block-staggered grid seismic wave simulation method based on SBP-SAT, comprising the following steps: Step 1: Establish the first-order velocity-pressure acoustic wave equation with auxiliary differential equations; introduce velocity-related auxiliary variables. , and stress-related auxiliary variables , A two-dimensional first-order hyperbolic acoustic system with ADE-CPML (auxiliary differential equation complex frequency shift perfect matching layer) and the evolution equation of auxiliary variables were established.
[0008] Step 2: Spatial domain staggered mesh SBP discretization and application of outer boundary SAT conditions; defining the computational domain. , to sound pressure Defined at integer nodes, horizontal velocity and vertical velocity Two-dimensional SBP spatial discrete operators are defined at half-integer nodes to construct a single-block mesh. The free surface boundary conditions are weakly applied using the SAT penalty term, and the SAT boundary correction term is synchronously injected into the physical variable equation and the auxiliary differential equation to form a single-block mesh semi-discrete system.
[0009] Step 3: Construct SAT coupling terms for non-uniform mesh interfaces; divide the computational domain containing the non-uniform medium into coarse and fine mesh blocks with different mesh densities, correct the coarse mesh boundaries using the fine mesh solution variables, and simultaneously correct the fine mesh boundaries using the coarse mesh solution variables, introducing an interface interpolation operator. and interface projection operator Based on the physical requirement of continuous sound pressure and normal velocity, an interface SAT penalty term is constructed to form a two-dimensional semi-discrete system with non-uniform mesh interface conditions.
[0010] Step 4: Simulation parameters and observation system settings; establish a two-dimensional geological model and set the corresponding parameters, including computational domain range, coarse and fine grid parameters, coarse and fine grid PML absorption layer parameters, source parameters, medium parameters, time integration parameters, and observation system parameters.
[0011] Step 5: Time Domain Integration and Wavefield Simulation; Based on the parameters set in Step 4, the semi-discrete system constructed in Steps 1 to 3 is subjected to time integration using the staggered frog-jump method until all time steps are completed, and the simulated seismic record is output and obtained.
[0012] Furthermore, in step one, the two-dimensional first-order hyperbolic acoustic system containing ADE-CPML is as follows: In the formula, , Auxiliary variables related to speed; , An auxiliary variable related to stress; For sound pressure, , They are horizontal velocity and vertical velocity, respectively. For the density of the medium, Bulk modulus , , , All are absorption layer decay functions, which are non-zero only within the PML layer.
[0013] Furthermore, in step two, the semi-discrete system of a single grid includes semi-discrete equations for physical variables and semi-discrete equations for auxiliary variables.
[0014] The semi-discrete equations for the physical variables are: The semi-discrete equation for the auxiliary variable is: In the formula, , , These are the medium coefficient matrices for sound pressure, horizontal velocity, and vertical velocity, respectively, corresponding to... , , Its size is the same as that of the discrete grid nodes of the corresponding system; , , All are two-dimensional norm matrices of the SBP operator; , , , The matrix represents the derivative operators; , , Represents the solution variables in a discrete space; , , , , For boundary SAT penalty terms; This is the penalty coefficient for the SAT term in the auxiliary differential equation, used to control the proportion of the SAT term in the auxiliary differential equation.
[0015] Furthermore, in step two, the boundary SAT penalty term satisfies the condition of no boundary reflected wave. Its expression is: in, , Let be the one-dimensional norm matrix of the SBP operator. , , , These represent the boundary projection operators for the top, bottom, left, and right boundaries, respectively. , , , These represent boundary selection operators for the top, bottom, left, and right boundaries, respectively. , These are unit diagonal matrices of the same dimension as the solution vector P in the vertical and horizontal directions, respectively. This is the penalty coefficient for the SAT term in the solution space, used to control the proportion of the SAT term in the solution space; This represents the matrix tensor product operation.
[0016] Furthermore, in step three, the two-dimensional semi-discrete system with non-uniform mesh interface conditions includes coarse-mesh semi-discrete equations and fine-mesh semi-discrete equations.
[0017] The coarse-grid semi-discrete equation is: The fine-grid semi-discrete equation is as follows: In the formula, subscript C represents the parameter in the coarse mesh, and subscript F represents the parameter in the fine mesh. The definitions of each parameter are the same as in step two.
[0018] Furthermore, in step three, the interface SAT penalty term satisfies the physical requirement of continuity between sound pressure and normal velocity: Where vector express , This indicates element-wise multiplication. Represents the outward unit vector on the interface. This is the interface between coarse and fine grids.
[0019] The expression for the SAT penalty item in the interface is: in, , , , These represent the application applied to the coarse mesh. Variable equations, coarse grid Variable equations, fine mesh Variable equations, fine mesh SAT penalty term on variable equations , , , This represents the corresponding penalty coefficient.
[0020] Furthermore, in step three, the interface interpolation operator and interface projection operator Satisfying SBP energy compatibility requirements: In the formula, The one-dimensional norm matrix of the SBP operator in the coarse mesh. Let be the one-dimensional norm matrix of the SBP operator in a fine mesh.
[0021] Furthermore, in step four, the sampling interval in the time integration parameter... Total time steps The CFL stability condition is satisfied, and its calculation formula is as follows: In the formula, For CFL stability condition parameters, , For finer grid step size, , For coarse grid step size, Maximum wave velocity of the medium; total time steps From the total simulation duration With sampling interval The ratio is determined.
[0022] Furthermore, in step five, the implementation strategy for time integration using the staggered frog-jump method is as follows: update the velocity vector at half a time step, update the sound pressure variable at an integer time step, and calculate the wave field value at each time step in sequence.
[0023] Compared with the prior art, the present invention has the following advantages: 1. This invention innovatively injects SAT boundary correction terms into the auxiliary differential equation of the Perfectly Matched Layer (PML) while applying SAT boundary conditions to the traditional physical variable equations. This eliminates the discretization incompatibility between auxiliary variables and physical variables at the truncation boundary of PML under higher-order difference schemes from the theoretical mathematical framework, and completely solves the problem of "PML late instability" in wavefield simulation. It can meet the application scenarios with extremely high requirements for numerical stability, such as long-term earthquake monitoring and multiple wave imaging.
[0024] 2. For the block staggered mesh architecture, this invention designs interpolation and projection operators that satisfy the SBP energy compatibility condition, and weakly applies the interface flux continuity condition through the SAT penalty term. This coupling method allows the meshes on both sides of the interface to be discontinuous while fully preserving the global antisymmetry of the SBP difference operator, strictly ensuring that the total energy of the system will not grow non-physically at the interface, effectively eliminating the false reflection at the interface between coarse and fine meshes, and realizing high-precision and stable coupling of non-uniform meshes (coarse and fine mesh interfaces).
[0025] 3. The method of this invention allows for the deployment of grids with different densities in different velocity regions of strongly non-uniform media. Fine grids are used in low-speed / complex structure regions to meet stability and sampling requirements, while coarse grids are used in high-speed regions to avoid oversampling. Under the premise of satisfying the local CFL condition and sampling theorem, the total number of system degrees of freedom is effectively reduced, avoiding a huge waste of computational resources. At the same time, it takes into account the dual advantages of low-frequency dispersion of staggered grids and high efficiency of block grids, significantly improving the accuracy and computational efficiency of acoustic wave forward modeling in complex media. Attached Figure Description
[0026] Figure 1 Figure (a) is a schematic diagram of the block staggered mesh discrete method used in this invention, and Figure (b) is a schematic diagram of the distribution of nodes in the block staggered mesh and a schematic diagram of the division of the block mesh and the perfect matching layer (PML) region.
[0027] Figure 2 This is a schematic diagram of the construction of the one-dimensional SBP discrete operator in this invention.
[0028] Figure 3 This is a schematic diagram of the interface transformation operator with a mesh ratio of 2:1 according to the present invention, wherein Figure (a) is the interface projection operator and Figure (b) is the interface interpolation operator.
[0029] Figure 4 The results of this invention are simulations of a coarse-grid seismic source in a homogeneous medium model. Figure (a) is a snapshot of the simulated wavefield, Figure (b) is the total energy curve of the system, Figure (c) is the seismic record, and Figure (d) is the receiver record of a single-channel geophone.
[0030] Figure 5 The results of this invention are simulations of a fine-grid seismic source in a homogeneous medium model. Figure (a) is a snapshot of the simulated wavefield, Figure (b) is the total energy curve of the system, Figure (c) is the seismic record, and Figure (d) is the receiver record of a separate channel geophone. Detailed Implementation
[0031] The present invention will be further described below.
[0032] like Figures 1 to 3 As shown, the specific implementation steps of this embodiment are as follows: Step 1: Establish the first-order velocity-pressure acoustic wave equation including the auxiliary differential equation: Introducing speed-related auxiliary variables , and stress-related auxiliary variables , A two-dimensional first-order hyperbolic acoustic system with an auxiliary differential equation complex frequency shift perfectly matched layer (ADE-CPML) is constructed: The evolution equation for the auxiliary variables is as follows: In the formula, , Auxiliary variables related to speed; , An auxiliary variable related to stress; For sound pressure, , They are horizontal velocity and vertical velocity, respectively. For the density of the medium, Bulk modulus , , , All are absorption layer decay functions, which are non-zero only within the PML layer and zero within the computational domain.
[0033] Step 2: Spatial Domain Staggered Mesh SBP Discretization and Application of Outer Boundary SAT Conditions: Define computational domain , to sound pressure Defined at integer nodes, horizontal velocity and vertical velocity Two-dimensional SBP spatial discrete operators are defined at the half-integer nodes, and constructed on a single mesh. The two-dimensional SBP operator is constructed as follows: In the formula , These are the norm matrix and derivative calculation matrix of the one-dimensional SBP discrete operator, respectively. It is a unit diagonal matrix.
[0034] To achieve reflection-free absorption at the outer boundary of the computational domain and ensure long-term stability, a weak free surface boundary condition is applied using a SAT penalty term. They are then synchronously injected into the physical variable equations and auxiliary differential equations to form a semi-discrete system with a single grid.
[0035] The semi-discrete equations for the physical variables are as follows: The semi-discrete equation for the auxiliary variable is: The expression for the boundary SAT penalty term is: In the above formula, , , The medium coefficient matrix represents, respectively , , Its size is the same as that of the discrete grid nodes of the corresponding system; , Let be the one-dimensional norm matrix of the SBP operator. , , Here is the two-dimensional norm matrix of the SBP operator; , , , The matrix represents the derivative operators; , , Represents the solution variables in a discrete space; , , , , For boundary SAT penalty terms; , The penalty coefficients for the SAT term control the weight of the SAT term in the solution space and the auxiliary differential equation, respectively. , , , These represent the boundary projection operators for the top, bottom, left, and right boundaries, respectively. , , , Select boundary operators to represent the top, bottom, left, and right boundaries respectively; , These are unit diagonal matrices of the same dimension as the solution vector P in the vertical and horizontal directions, respectively. This represents the matrix tensor product operation.
[0036] The optimal value for the penalty parameter is: , , , ; , , , , , .
[0037] Step 3: Construct SAT coupling terms for non-uniform mesh interfaces: The computational domain containing a non-uniform medium is divided into grid blocks of different grid densities, including coarse and fine grid blocks. The fine grid solution variables are used to correct the coarse grid boundaries, and vice versa, introducing an interface interpolation operator. and interface projection operator Based on the physical requirement of continuous sound pressure and normal velocity, an interface SAT penalty term is constructed to form a two-dimensional semi-discrete system with non-uniform mesh interface conditions.
[0038] The coarse-grid semi-discrete equations are as follows: The semi-discrete equations for fine mesh are: The interface SAT penalty term satisfies the physical requirement of continuity between sound pressure and normal velocity: Where vector express , This indicates element-wise multiplication. Represents the outward unit vector on the interface. This is the interface between coarse and fine grids.
[0039] The expression for the SAT penalty item in the interface is: In the above formula, the subscript C represents the parameter in the coarse mesh, and the subscript F represents the parameter in the fine mesh. The definitions of each parameter are the same as in step two. , , , These represent the application applied to the coarse mesh. Variable equations, coarse grid Variable equations, fine mesh Variable equations, fine mesh SAT penalty term on variable equations , , , This represents the corresponding penalty coefficient.
[0040] To ensure the energy stability of the interface coupling, the interpolation operator between coarse and fine meshes... and projection operator SBP energy compatibility conditions must be met: In the formula, The one-dimensional norm matrix of the SBP operator in the coarse mesh. Let be the one-dimensional norm matrix of the SBP operator in a fine mesh.
[0041] The optimal value for the corresponding penalty parameter is: .
[0042] Step 4: Simulation Parameters and Observation System Settings: Establish a two-dimensional geological model and define the computational domain. Set the grid thickness interface to Fine mesh region in the x and y directions The number of variable subgrids is The corresponding grid size is , The coarse grid region in the x and y directions The number of variable subgrids is The corresponding grid size is , The thickness of the absorption layer in the coarse grid region is Reflection coefficient The thickness of the absorption layer in the fine mesh region is Reflection coefficient A Lake wavelet source was selected, and the dominant frequency of the source wavelet was set. Set the medium parameters for the coarse and fine grid regions, and the medium wave velocity in the coarse grid region. ,density Medium wave velocity in fine mesh region ,density Set the total simulation duration. and CFL stability condition parameters Sampling interval Total time steps It can be calculated; set the observation system parameters, including the location of the seismic source. , Detector coordinates , .
[0043] Among them, sampling interval The CFL stability condition is satisfied, and the calculation formula is: In the formula, Maximum wave velocity of the medium; total time steps From the total simulation duration With sampling interval The ratio is determined.
[0044] Step 5: Time-domain integration and wave field simulation: Based on the two-dimensional geological model and corresponding parameters set in step four, the semi-discrete system constructed in steps one to three is subjected to time integration using the staggered frog-jump method until all time steps are completed, and the corresponding simulated seismic records are output and obtained.
[0045] The implementation strategy of the staggered frog-jump time integration algorithm is as follows: update the velocity vector at half a time step, update the pressure variable at an integer time step, and calculate the wave field value at each time step in sequence.
[0046] Effect verification: To verify the effectiveness of the method of the present invention, the following numerical experiments were conducted: Consider a size of In a two-dimensional uniform computational domain, the velocity of sound waves is set to a constant value. , To test the coupling performance of non-uniform meshes, the computational domain was stretched horizontally. The area is divided into two grid blocks: the upper region is discretized with a coarse grid and a grid spacing of [value missing]. , The lower region is discretized using a fine grid with a grid spacing of [value missing]. , The corresponding grid ratio is The thickness of the coarse mesh absorption layer is The thickness of the fine mesh absorption layer is Reflectance coefficient The seismic source uses the dominant frequency. The Reichschild wavelet. The total simulation duration is set to... CFL stability condition parameters The sampling interval is calculated. Total time steps .
[0047] To comprehensively evaluate the interface performance, two source placement scenarios are considered: ① The source is located in a coarse grid region. ① Simulate the process of wave field penetration from coarse to fine grid; ② The seismic source is located in the fine grid region. The simulation simulates the process of the wave field transitioning from a fine grid to a coarse grid. The detector is located at a burial depth of 10m and the channel spacing is 1m.
[0048] The simulation results of coarse-grid and fine-grid seismic sources are as follows: Figure 4 , Figure 5 As shown. By Figure 4 a and Figure 5 The wave field snapshots at 40s, 80s, and 120s in a show that no artificial reflection waves are generated at the interface between the coarse and fine grids or at the surrounding boundaries, regardless of whether the wave enters the fine grid from the coarse grid or the fine grid from the coarse grid. Figure 4 b and Figure 5 b indicates that the system maintains a stable state of energy dissipation during the simulation, with no non-physical energy growth; Figure 4 c and Figure 5 c represents the final simulated earthquake record. Figure 4 d and Figure 5 d represents the extracted records from traces 49 to 54, showing that there are no other noise signals besides the direct wave. The above simulation results fully demonstrate the effectiveness of the method proposed in this invention, which can maintain the numerical stability of seismic wavelength-time simulation under block staggered grid discretization.
[0049] The above description is only a preferred embodiment of the present invention. It should be noted that for those skilled in the art, several improvements and modifications can be made without departing from the principle of the present invention, and these improvements and modifications should also be considered within the scope of protection of the present invention.
Claims
1. A block-staggered grid seismic wave simulation method based on SBP-SAT, characterized in that, Includes the following steps: Step 1: Introduce auxiliary variables related to velocity and pressure to establish the two-dimensional first-order hyperbolic acoustic system of ADE-CPML and the evolution equation of the auxiliary variables; Step 2: Define the computational domain, defining the sound pressure at integer nodes and the horizontal and vertical velocities at half-integer nodes respectively, and construct a two-dimensional SBP spatial discrete operator on a single mesh; apply the SAT penalty term to weakly apply the free surface boundary conditions, and simultaneously inject the SAT boundary correction term into the physical variable equation and the auxiliary differential equation to form a semi-discrete system on a single mesh; Step 3: Divide the computational domain containing the non-uniform medium into coarse and fine grid blocks with different grid densities, introduce interface interpolation operators and interface projection operators, construct interface SAT penalty terms based on the physical requirement of continuous sound pressure and normal velocity, and form a two-dimensional semi-discrete system with non-uniform grid interface conditions. Step 4: Establish a two-dimensional geological model and set the corresponding parameters; Step 5: Based on the parameters set in Step 4, use the staggered frog-jump method to perform time integration on the semi-discrete system constructed in Steps 1 to 3 until all time steps are completed, and output and obtain the simulated seismic record.
2. The block-staggered grid seismic wave simulation method based on SBP-SAT according to claim 1, characterized in that, In step one, the two-dimensional first-order hyperbolic acoustic system containing ADE-CPML is as follows: In the formula, , Auxiliary variables related to speed; , An auxiliary variable related to stress; For sound pressure, , They are horizontal velocity and vertical velocity, respectively. For the density of the medium, Bulk modulus , , , Both are absorption layer attenuation functions.
3. The block-staggered grid seismic wave simulation method based on SBP-SAT according to claim 2, characterized in that, In step two, the semi-discrete system of a single grid includes semi-discrete equations for physical variables and semi-discrete equations for auxiliary variables. The semi-discrete equations for the physical variables are: The semi-discrete equation for the auxiliary variable is: In the formula, , , These are the medium coefficient matrices for sound pressure, horizontal velocity, and vertical velocity, respectively. , , All are two-dimensional norm matrices of the SBP operator; , , , The matrix represents the derivative operators; , , Represents the solution variables in a discrete space; , , , , For boundary SAT penalty terms; This is the penalty coefficient for the SAT term in the auxiliary differential equation, used to control the proportion of the SAT term in the auxiliary differential equation.
4. The block-staggered grid seismic wave simulation method based on SBP-SAT according to claim 3, characterized in that, In step two, the boundary SAT penalty term satisfies the condition of no boundary reflected wave. Its expression is: in, , Let be the one-dimensional norm matrix of the SBP operator. , , , These represent the boundary projection operators for the top, bottom, left, and right boundaries, respectively. , , , These represent boundary selection operators for the top, bottom, left, and right boundaries, respectively. , These are unit diagonal matrices of the same dimension as the solution vector P in the vertical and horizontal directions, respectively. This is the penalty coefficient for the SAT term in the solution space, used to control the proportion of the SAT term in the solution space; This represents the matrix tensor product operation.
5. The block-staggered grid seismic wave simulation method based on SBP-SAT according to claim 4, characterized in that, In step three, the two-dimensional semi-discrete system with non-uniform mesh interface conditions includes coarse mesh semi-discrete equations and fine mesh semi-discrete equations. The coarse-grid semi-discrete equations are: The semi-discrete equations for fine mesh are: In the formula, subscript C represents the parameter in the coarse mesh, and subscript F represents the parameter in the fine mesh. The definitions of each parameter are the same as in step two.
6. The block-staggered grid seismic wave simulation method based on SBP-SAT according to claim 5, characterized in that, In step three, the expression for the interface SAT penalty term is: in, , , , These represent the application applied to the coarse mesh. Variable equations, coarse grid Variable equations, fine mesh Variable equations, fine mesh SAT penalty term on variable equations , , , This represents the corresponding penalty coefficient.
7. The block-staggered grid seismic wave simulation method based on SBP-SAT according to claim 6, characterized in that, In step three, the interface interpolation operator and interface projection operator Satisfying SBP energy compatibility requirements: In the formula, The one-dimensional norm matrix of the SBP operator in the coarse mesh. Let be the one-dimensional norm matrix of the SBP operator in a fine mesh.
8. The block-staggered grid seismic wave simulation method based on SBP-SAT according to claim 1, characterized in that, In step four, the parameters include the computational domain range, coarse and fine grid parameters, coarse and fine grid PML absorbing layer parameters, source parameters, medium parameters, time integration parameters, and observation system parameters.
9. The block-staggered grid seismic wave simulation method based on SBP-SAT according to claim 8, characterized in that, In step four, the sampling interval in the time integration parameter Total time steps The CFL stability condition is satisfied, and its calculation formula is as follows: In the formula, For CFL stability condition parameters, , For finer grid step size, , For coarse grid step size, Maximum wave velocity of the medium; total time steps From the total simulation duration With sampling interval The ratio is determined.
10. The block-staggered grid seismic wave simulation method based on SBP-SAT according to claim 1, characterized in that, In step five, the implementation strategy of time integration using the staggered frog-jump method is as follows: update the velocity vector at half a time step, update the sound pressure variable at an integer time step, and calculate the wave field value at each time step in sequence.