A variable-grid spectral element method for seismic wavefield simulation in complex loess plateau region

Abstract

The present application belongs to the technical field of seismic exploration and seismic wave numerical simulation, and relates to a variable grid spectral element method for simulating seismic wave field in complex loess plateau area, comprising the following steps: constructing a complex geological model comprising undulating near-surface and low-speed velocity-decreasing layer; performing variable grid profiling on the model according to the velocity distribution characteristics of underground medium to establish a global conformal quadrilateral spectral element unit model; adopting different longitudinal grid steps for different velocity layers under the condition of keeping consistent horizontal grid step; taking the minimum velocity simulation accuracy as the criterion for time step, performing numerical solution of seismic wave field based on the spectral element method to obtain wave field snapshot and single-shot snapshot. The variable grid spectral element method for simulating seismic wave field in complex loess plateau area adopts irregular grid profiling, avoids the problems of local oversampling of traditional conventional fine grid and insufficient precision of traditional coarse grid, and realizes the balance between simulation precision and calculation efficiency.

Description

Technical Field

[0001] This invention belongs to the field of seismic exploration and seismic wave numerical simulation technology, specifically relating to a variable grid spectral element method for simulating seismic wave fields in complex loess plateau regions. Background Technology

[0002] Complex loess plateau regions are widely distributed in the Loess Plateau and its adjacent areas in my country, and are typical complex surface areas for seismic exploration, engineering seismology, and near-surface structure detection. Complex loess plateau regions are typically characterized by dramatic surface undulations and the development of low-velocity near-surface layers, which can significantly alter the propagation path of seismic waves, causing wavefront distortion, energy scattering, and other problems. This, in turn, affects the accuracy of forward modeling and migration imaging results from seismic data in complex loess plateau regions. Therefore, accurately simulating the propagation process of seismic waves under conditions of undulating near-surface surfaces and low-velocity near-surface layers has become a core problem that needs to be solved.

[0003] Currently, commonly used methods in seismic wavefield simulation studies under complex near-surface conditions include the finite difference method, the finite element method, and the conventional spectral element method. However, the finite difference method, typically based on regular rectangular grids, faces two major challenges for models with significant undulations and dramatic topographic changes, such as those in complex loess plateau regions: 1. The grid discretization of undulating terrain interfaces can only approximate the geometry of the actual terrain as closely as possible, resulting in some errors; 2. Stepped boundaries can generate artificial spurious diffraction waves, affecting the accuracy of numerical simulations. While the finite element method can better adapt to complex geometric boundaries and irregular regions, offering greater flexibility in simulating wavefields of complex geological bodies, traditional finite element methods often suffer from large numerical dispersion, high computational accuracy dependent on grid quality, numerous degrees of freedom, and relatively low computational efficiency. The conventional spectral element method combines the finite element method's ability to handle complex boundaries with the high accuracy and low numerical dispersion of spectral methods, and has received widespread attention in the field of seismic wavefield simulation under complex terrain conditions. However, conventional spectral element methods still have shortcomings in complex loess plateau regions: on the one hand, in order to meet the accuracy requirements of short-wavelength propagation in the near-surface low-velocity deceleration layer, it is often necessary to use finer grid discretization across the entire region, which forces the use of overly dense grids in the deep high-velocity region, resulting in a large amount of unnecessary computational overhead; on the other hand, if local coarse-fine grid splicing or non-coordinated grid processing is directly adopted, wave field propagation errors, artificial scattering, or even numerical instability problems are easily introduced at the junction of coarse and fine grids, which is not conducive to the realization of high-precision wave field simulation under complex models.

[0004] In 2023, Liu Cai et al. published their research on forward modeling of seismic wavefields using a rotating spatiotemporal dual-variable grid finite-difference method in the *Chinese Journal of Geophysics*, pointing out that changing the spatial and temporal grid steps under rotating coordinates can increase computational efficiency. Unlike the finite-difference method, variable-grid strategies based on the spectral element method are relatively unexplored. Therefore, to address the need for seismic wavefield simulation under the combined effects of undulating near-surface morphology and low-velocity deceleration layers in complex loess plateau regions, a wavefield simulation method is needed that can accurately characterize complex surface morphology and near-surface low-velocity structures while balancing numerical accuracy and computational efficiency. This would ensure the accuracy of seismic wavefield simulation under complex near-surface undulations and low-velocity deceleration layers while reducing overall computational costs and improving simulation efficiency. Summary of the Invention

[0005] The purpose of this invention is to provide a variable grid spectral element method for simulating seismic wavefields in complex loess plateau regions, in order to solve the problems of low accuracy and poor computational efficiency in simulating seismic wavefield propagation under the combined effects of undulating near-surface and low-velocity deceleration layers in complex loess plateau regions, as well as the difficulty of conventional uniform grid methods in simultaneously addressing the complex surface characterization and numerical stability.

[0006] This invention is achieved through the following technical solution:

[0007] A variable grid spectral element method for simulating seismic wavefields in complex loess plateau regions includes the following steps:

[0008] a. Read in the two-dimensional velocity model and corresponding parameter information of the complex loess plateau region, perform parameter identification, and establish a basic model for two-dimensional seismic wavefield simulation;

[0009] b. Based on velocity differences, according to the two-dimensional seismic wavefield simulation basic model established in step a, set the source, number of time sampling points, time sampling step size, absorbing boundary and detector parameters, which need to meet the CFL stability condition.

[0010] c. The overall area is Based on the magnitude and variation of medium velocity in each region, a variable mesh is generated. While keeping the number of horizontal meshes constant, the vertical mesh step size is adjusted to establish a conformal quadrilateral spectral element model. The mesh step size satisfies the CFL stability condition.

[0011] d. Generate GLL nodes within the spectral element based on the time sampling step size and grid step size;

[0012] e. Obtain the assignment parameter file, and use GLL nodes to assign node values ​​to the longitudinal wave velocity, transverse wave velocity and density in the two-dimensional velocity model to generate the medium parameter file.

[0013] f. Assemble the mass matrix and stiffness matrix based on the spectral element model, apply source loading conditions and boundary conditions, and use the time-domain stepwise integration method to iteratively solve the wavefield to obtain the two-dimensional seismic wavefield propagation results.

[0014] g. Update the displacement field, output the detector seismic record and seismic wave field propagation results obtained in step f, and analyze the simulation results to characterize the seismic wave propagation law under the complex loess plateau undulating surface and non-uniform medium conditions.

[0015] Further, in step a, the parameter information includes the model's lateral length, longitudinal length, undulating surface parameters, medium density parameters, and medium wave velocity parameters. The undulating surface parameters are used to describe the surface elevation variation characteristics of the complex loess plateau region, and the medium wave velocity parameters include longitudinal wave velocity and transverse wave velocity, which are used to characterize the underground velocity stratification structure.

[0016] Further, in step b, the source is selected as the Ricker wavelet to determine the excitation location and mode of the seismic wave. The number of time sampling points and the time sampling step size are used to determine the time domain calculation range. The absorbing boundary is used to attenuate the artificial reflected waves at the boundary of the calculation area. The detector parameters are the seismic response information received during the propagation of the wave field. The source function is:

[0017] ;

[0018] In the formula, M is the seismic moment tensor, a second-order tensor describing the intensity and properties of the earthquake source; Indicates the spatial distribution of earthquake sources. The location of the earthquake source is S(t); S(t) is the source time function, which describes the change of the seismic wave over time.

[0019] The stability conditions for CFL are as follows:

[0020] ;

[0021] In the formula, V is the distance between two adjacent GLL integration points; p represents the P-wave velocity in the medium; C is the CFL condition number, which gradually decreases as the wave velocity ratio increases.

[0022] Furthermore, the time step is determined based on the lowest velocity and the smallest grid size within the computational region.

[0023] Further, in step c, the variable mesh generation is based on the medium velocity distribution. A smaller longitudinal mesh step size is used in low-velocity regions and regions with large velocity variations, while a larger longitudinal mesh step size is used in high-velocity regions and regions with small velocity variations. The established spectral elements are conformal quadrilateral elements, and each element is connected through a common boundary. In the elastic medium, the propagation behavior of seismic waves is described by the linear elastic wave equation:

[0024] ;

[0025] In the formula, ρ is the density of the medium, u is the displacement, σ is the stress, and f is the source. This represents the spatial gradient.

[0026] Using the weighted residual principle, the differential equation is transformed into an integral equation to solve the wave equation; a test function is introduced. In the overall region Integrating within:

[0027] ;

[0028] By performing integration by parts, we obtain the weak form of the wave equation:

[0029] .

[0030] In the formula, It is the unit vector of the outward normal on the free surface boundary conditions.

[0031] Furthermore, in step c, the variable mesh partitioning divides the computational domain into layers based on the velocity distribution of the subsurface medium. The low-velocity layer uses a smaller vertical mesh step size, while the high-velocity layer uses a larger vertical mesh step size. The variable mesh partitioning maintains a consistent horizontal mesh step size throughout the entire computational domain, adjusting only the vertical mesh step size. Each spectral element is directly connected through a common boundary, and no transition region is set between coarse and fine meshes. The ratio of the longer side to the shorter side of each spectral element is no greater than 4. For complex near-surface areas, in step c, the mesh step size is set based on the maximum near-surface layer thickness, and the mesh naturally refines the local density in areas with severe topographic relief.

[0032] Further, in step d, the GLL nodes are the interpolation and integration nodes in the spectral element unit, used to implement higher-order interpolation and numerical integration in the spectral element method. The time step and spatial discrete step size should satisfy the CFL stability condition; the GLL nodes are the n+1 roots of the following equation:

[0033] ;

[0034] In the formula, for The first derivative of a Legendre polynomial of order 1. The order Lagrange interpolation basis functions are expressed as follows:

[0035] .

[0036] In the formula, Lagrange interpolation basis functions This indicates that the basis function corresponds to the i-th node. x is the independent variable.

[0037] Further, in step e, the medium parameter file includes P-wave velocity, S-wave velocity, and density values ​​at each node, used to characterize the medium distribution characteristics at different locations within the computational domain; for numerical integration on the elements, GLL integration is used for calculation.

[0038] ;

[0039] Where, ω α ω β α, β = 0, and n is the weight coefficient of the GLL integral; it involves a weighted summation of the values ​​of the approximation function at the nodes of the GLL integral, where the weight coefficients are constants independent of each element.

[0040] Further, in step f, the mass matrix and stiffness matrix are assembled from the geometric information, nodal information and medium parameters of the spectral element. The boundary conditions include the free surface boundary conditions of the undulating ground surface and the absorbing boundary conditions of the perfectly matched layer. The wave field iterative solution is used to obtain the displacement field and velocity field distribution results at each time step.

[0041] For the dielectric unit Ω e The solution is:

[0042] ;

[0043] ;

[0044] in, , , , The time iteration result is:

[0045] ;

[0046] in, For the quality matrix, Let K be the velocity vector for the next time step, K be the global stiffness matrix, and F be the source term.

[0047] Further, in step g, the geophone seismic record includes information on the displacement response and velocity response at each geophone location as a function of time.

[0048] Compared with the prior art, the beneficial effects of the present invention are:

[0049] 1. This method fully utilizes the technical advantages of the spectral element method in high-precision fitting of complex boundaries and high-precision solution of wave field propagation. It can realistically characterize the undulating surface morphology and near-surface velocity structure characteristics of complex loess plateau areas, thereby realizing the fine simulation of seismic wave propagation process under complex surface conditions.

[0050] 2. This method implements variable grid partitioning based on the velocity distribution characteristics of the underground medium. A smaller vertical grid step size is used in the low-velocity deceleration layer region and a larger vertical grid step size is used in the high-velocity region. This can effectively balance the simulation accuracy of the short-wavelength propagation characteristics of the low-velocity layer with the computational efficiency requirements of the deep region.

[0051] 3. This method constructs a global conformal quadrilateral spectral element model while maintaining a consistent horizontal grid step size, enabling the natural transfer of wave field information between grid regions of different scales. This reduces the artificial scattering, wave field distortion, and numerical instability problems caused by improper transition between coarse and fine grids in traditional local refinement methods.

[0052] 4. This method can significantly outperform the traditional coarse-grid method in simulation results while maintaining high simulation accuracy, and can obtain wave field simulation results that are highly consistent with the traditional fine-grid method.

[0053] 5. This method can significantly reduce computational costs and improve computational efficiency while maintaining high-precision simulation results, thus achieving a better balance between simulation accuracy and computational efficiency;

[0054] 6. This method can provide reliable technical support for forward modeling of seismic wavefields, seismic record analysis, research on wavefield propagation mechanisms, and subsequent seismic data processing and interpretation under the conditions of undulating near-surface and low-velocity deceleration layers in complex loess plateau areas. It has high theoretical significance and engineering application value. Attached Figure Description

[0055] To more clearly illustrate the technical solutions of the embodiments of the present invention, the accompanying drawings used in the embodiments will be briefly introduced below. It should be understood that the following drawings only show some embodiments of the present invention and should not be regarded as a limitation on the scope. For those skilled in the art, other related drawings can be obtained based on these drawings without creative effort.

[0056] Figure 1 This is a flowchart illustrating the steps of the variable grid spectral element method for simulating seismic wavefields in complex loess plateau regions according to the present invention.

[0057] Figure 2This is a detailed schematic diagram illustrating the implementation process of a specific embodiment of the variable grid spectral element method for simulating seismic wavefields in complex loess plateau regions according to the present invention. In this diagram, (a) is the embodiment model, (b) is a schematic diagram of the variable grid, (c) is a schematic diagram of the coarse grid, and (d) is a schematic diagram of the fine grid.

[0058] Figure 3 The following are examples of seismic record analysis diagrams, where (a) is a wavefield snapshot of the traditional coarse-grid spectral method, (b) is a wavefield snapshot of the traditional fine-grid spectral method, and (c) is a wavefield snapshot of the variable-grid spectral method.

[0059] Figure 4 The following are examples of single-track record analysis diagrams: (a) is a comparison diagram of single-track records with traditional coarse mesh, fine mesh and variable mesh; (b) is an error analysis diagram of traditional coarse mesh, fine mesh and variable mesh; and (c) is a magnified view of (a).

[0060] Figure 5 The following is a comparison of earthquake records in a complex loess plateau construction area. (a) is the actual earthquake record in a complex loess plateau construction area, and (b) is the earthquake record simulated by the variable grid spectral element method. Detailed Implementation

[0061] The present invention will now be described in further detail with reference to the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are merely illustrative of the invention and not intended to limit it. Furthermore, it should be noted that, for ease of description, the accompanying drawings show only the parts relevant to the present invention, and not all of the structures.

[0062] This invention addresses the challenge of conventional uniform grid wavefield simulation methods in complex loess plateau regions, characterized by dramatic near-surface undulations and the development of low-velocity deceleration layers. These methods struggle to simultaneously achieve detailed representation of the complex surface, accurate description of low-velocity layer propagation characteristics, and computational efficiency. To address this, a variable-grid spectral element method for seismic wavefield simulation in complex loess plateau regions is proposed. This method uses variable grid partitioning based on the velocity distribution characteristics of the subsurface medium. While maintaining a consistent horizontal grid step size, the vertical grid step size is adjusted, employing a smaller vertical grid for low-velocity layers and a larger vertical grid for high-velocity layers. A globally conformal quadrilateral spectral element model is established. This approach ensures the accuracy of seismic wavefield simulation under complex near-surface undulations and low-velocity deceleration layer conditions while reducing overall computational costs and improving simulation efficiency.

[0063] like Figure 1 As shown, the present invention provides a variable grid spectral element method for simulating seismic wavefields in complex loess plateau regions, comprising the following steps:

[0064] a. Read in the two-dimensional velocity model and corresponding parameter information of the complex loess plateau region, identify the parameters, and establish a basic model for two-dimensional seismic wave field simulation. The parameter information includes the model's lateral length, longitudinal length, undulating surface parameters, medium density parameters, and medium wave velocity parameters. The undulating surface parameters are used to describe the surface elevation variation characteristics of the complex loess plateau region, and the medium wave velocity parameters include P-wave velocity and S-wave velocity, which are used to characterize the underground velocity stratification structure.

[0065] b. Based on velocity differences, and using the two-dimensional seismic wavefield simulation model established in step a, set the source, number of time sampling points, time sampling step size, absorbing boundary, and detector parameters. The source uses a Ricker wavelet to determine the seismic wave excitation location and mode. The number of time sampling points and the time sampling step size determine the time-domain calculation range and must satisfy the CFL stability condition. The absorbing boundary attenuates artificial reflections at the boundary of the calculation area. The detector parameters receive the seismic response information during wavefield propagation. The source function is:

[0066] ;

[0067] In the formula, M is the seismic moment tensor, a second-order tensor describing the intensity and properties of the earthquake source. Indicates the spatial distribution of earthquake sources. The location of the seismic source is given by S(t). S(t) is the source time function, describing the change of the seismic wave over time. The stability conditions of the CFL are as follows:

[0068] ;

[0069] In the formula, V is the distance between two adjacent GLL integration points; p represents the P-wave velocity in the medium; C is the CFL condition number, which gradually decreases as the wave velocity ratio increases.

[0070] c. The overall area is A variable mesh was used to determine the medium velocity magnitude and variation differences in each region. While maintaining a constant number of horizontal meshes, the vertical mesh step size was adjusted to establish a conformal quadrilateral spectral element model. The variable mesh was based on the medium velocity distribution, using a smaller vertical mesh step size in low-velocity regions and regions with large velocity variations, and a larger vertical mesh step size in high-velocity regions and regions with small velocity variations. The established spectral elements were conformal quadrilateral elements, connected by a common boundary to avoid artificial boundary scattering caused by non-conformal connections. The mesh step size satisfied the CFL stability condition. In elastic media, the propagation behavior of seismic waves is described by the linear elastic wave equation:

[0071] ;

[0072] In the formula, ρ is the density of the medium, u is the displacement, σ is the stress, and f is the source. This represents the spatial gradient.

[0073] Using the weighted residual principle, the differential equation is transformed into an integral equation, and the wave equation is solved. A test function is introduced. In the overall region Integrating within:

[0074] ;

[0075] By performing integration by parts, we obtain the weak form of the wave equation:

[0076] ;

[0077] In the formula, It is the unit vector of the outward normal on the free surface boundary conditions.

[0078] d. Generate GLL nodes within the spectral element based on the time sampling step size and grid step size. GLL nodes are interpolation and integration nodes within the spectral element, used to implement higher-order interpolation and numerical integration in the spectral element method. The time step size and spatial discrete step size should satisfy the CFL stability condition. The GLL nodes are the n+1 roots of the following equation:

[0079] ;

[0080] In the formula, for The first derivative of a Legendre polynomial of order 1. The order Lagrange interpolation basis functions are expressed as follows:

[0081] ;

[0082] In the formula, Lagrange interpolation basis functions This indicates that the basis function corresponds to the i-th node. x is the independent variable;

[0083] e. Obtain the assigned parameter file. Using GLL nodes, assign node values ​​to the P-wave velocity, S-wave velocity, and density in the two-dimensional velocity model to generate a medium parameter file. This file includes the P-wave velocity, S-wave velocity, and density values ​​at each node, characterizing the medium distribution features at different locations within the computational domain and providing a parameter basis for subsequent spectral element construction and wavefield solution. Numerical integration on the elements is calculated using GLL integration.

[0084] ;

[0085] In the formula, ω α ω β α, β = 0, n is the weight coefficient of the GLL integral, and any function in two-dimensional space can be approximately represented as: , , This involves a weighted summation of the approximation function's values ​​at the GLL integration nodes, where the weighting coefficients are constants independent of each element.

[0086] f. Set boundary conditions, solve the weak form of the elastic wave equation, obtain the cyclic mass matrix, and perform CFL stability condition judgment. Assemble the mass and stiffness matrices based on the spectral element model, apply source loading conditions and boundary conditions, and use a time-domain stepwise integration method to iteratively solve the wavefield, obtaining the two-dimensional seismic wavefield propagation results. The mass and stiffness matrices are assembled from the spectral element geometric information, nodal information, and medium parameters. Boundary conditions include undulating free surface boundary conditions and perfectly matched layer absorbing boundary conditions. The wavefield iterative solution is used to obtain the displacement and velocity field distributions at each time step. For the medium element Ω... e The solution is:

[0087] ;

[0088] ;

[0089] in, , , , The time iteration result is:

[0090] ;

[0091] in, For the quality matrix, Let K be the velocity vector for the next time step, K be the global stiffness matrix, and F be the source term.

[0092] The loop iterates and performs a CFL stability check as time t increases from 0 to its maximum value t. max The steps performed simultaneously in the loop will terminate if the stability condition is not met during the loop, and steps b-e will be repeated until the condition is met.

[0093] g. Update the displacement field u(t+dt), output the geophone seismic record and seismic wave field propagation results obtained in step f, and analyze the simulation results. The geophone seismic record includes information on the displacement response and velocity response at each geophone location as a function of time. The analysis is used to characterize the seismic wave propagation law under the complex loess plateau undulating surface and non-uniform medium conditions.

[0094] In step b, the time step is determined based on the lowest velocity and the smallest grid size within the computational region to meet the accuracy requirements for wave field propagation calculation in low-velocity regions.

[0095] In step c, the variable mesh partitioning is based on the velocity distribution of the underground medium to divide the computational domain into layers. The low-velocity layer uses a smaller vertical mesh step size, while the high-velocity layer uses a larger vertical mesh step size, so as to reduce the overall spatial discretization cost while ensuring the accuracy of the wave field simulation.

[0096] In step c, the variable mesh partitioning maintains a consistent horizontal mesh step size across the entire computational domain, adjusting only the vertical mesh step size.

[0097] In step c, all the established spectral elements are conformal quadrilateral elements, and the spectral elements are directly connected through a common boundary. No transition region is set between coarse and fine meshes, so that the wave field information can be naturally transmitted between the spectral elements and artificial scattering can be avoided.

[0098] In step c, the ratio of the longer side to the shorter side of the spectral element is no greater than 4 to avoid excessive mesh distortion and to ensure that the Jacobian matrix is ​​not 0 during the element mapping process.

[0099] For complex near-surface areas, step c sets the grid step size based on the maximum near-surface layer thickness. The grid is naturally locally refined in areas with severe topographic undulations to improve the accuracy of wave field simulation under complex terrain conditions.

[0100] Example 1:

[0101] a. such as Figure 2 As shown in (a), the two-dimensional velocity model and corresponding parameter information are read in. The horizontal length is 2000m, the vertical length is 1500m, the velocity of the first layer is 800m / s, the velocity of the second layer is 2000m / s, the velocity of the third layer is 3000m / s, the velocity of the fourth layer is 4000m / s, and the density of the first layer is 1300kg / m³. 3 The density of the second layer is 1600 kg / m³. 3 The density of the third layer is 2000 kg / m³. 3 The density of the fourth layer is 2300 kg / m³. 3 .

[0102] b. Set the seismic source. Here, the Ricker wavelet is used with a dominant frequency of 20Hz. The shot is fired from the center of the ground. The seismic source location is (1000, 1200). The number of time sampling points is 50,000, the time sampling step is 0.05ms, the thickness of the coarse grid perfect matching layer is 10, the thickness of the boundary between the variable grid and the fine grid is 30, and there are 2,000 detectors. See Table 1. The grid step size of the coarse grid is 40m × 30m, the grid step size of the fine grid is 4m × 3m, the grid step size of the first layer of the variable grid is 4m × 3m, the grid step size of the second layer is 4m × 5m, the grid step size of the third layer is 4m × 6m, and the grid step size of the fourth layer is 4m × 8m.

[0103] c. For example Figure 2 (b) Figure 2 (c) and Figure 2 As shown in (d) (the actual mesh density is too large, so the image is a schematic diagram for ease of display), a variable mesh is generated based on the magnitude and variation of the medium velocity in each region. The vertical mesh step size is adjusted while keeping the number of horizontal meshes constant, and a conformal quadrilateral spectral element model is established. The variable mesh generation is based on the medium velocity distribution. A smaller vertical mesh step size is used in low-velocity regions and regions with large velocity variations, while a larger vertical mesh step size is used in high-velocity regions and regions with small velocity variations. The established spectral element is a conformal quadrilateral element, and each element is connected through a common boundary to avoid artificial boundary scattering caused by non-conformal connections.

[0104] d. Generate GLL nodes within the spectral element based on the time sampling step size and spatial discrete step size. GLL nodes are interpolation and integration nodes within the spectral element, used to implement higher-order interpolation and numerical integration in the spectral element method. The time step size, spatial discrete step size, medium velocity, and wavelength should satisfy the CFL stability condition.

[0105] e. Using GLL nodes, assign node values ​​to the P-wave velocity, S-wave velocity, and density in the two-dimensional velocity model to generate a medium parameter file. The medium parameter file includes the P-wave velocity, S-wave velocity, and density values ​​at each node, which are used to characterize the medium distribution characteristics at different locations in the computational region and provide a parameter basis for subsequent spectral element construction and wave field solution.

[0106] f. Assemble the mass and stiffness matrices based on the spectral element model, apply source loading conditions and boundary conditions, and use the time-domain stepwise integration method to iteratively solve the wavefield to obtain the two-dimensional seismic wavefield propagation results. The mass and stiffness matrices are assembled from the geometric information, nodal information and medium parameters of the spectral elements. The boundary conditions include free surface boundary conditions of undulating ground and perfectly matched layer absorbing boundary conditions. The wavefield iterative solution is used to obtain the displacement and velocity field distribution results at each time step.

[0107] g. Output the geophone seismic records and seismic wave field propagation results obtained in step f, and analyze the simulation results. The geophone seismic records include information on the displacement response and velocity response at each geophone location as a function of time. The analysis is used to characterize the seismic wave propagation law under the complex loess plateau undulating surface and non-uniform medium conditions.

[0108] Figure 3 and Figure 4 The results of the comparative analysis of 0.5s wavefield snapshots and single-channel records in the example model using the method of the present invention are shown. The results of this embodiment are analyzed below in conjunction with the accompanying drawings.

[0109] As shown in Figure 3(a), the wave field snapshot obtained by the traditional coarse grid spectral element method is relatively insufficient in characterizing local wave field details under complex undulating near-surface and low-velocity deceleration layer conditions. Especially in the near-surface low-velocity layer region, the wavefront morphology and energy distribution are distorted to a certain extent. This indicates that when the grid scale is large, it is difficult to fully characterize the short-wavelength propagation characteristics in low-velocity media, thus adversely affecting the simulation accuracy.

[0110] like Figure 3 As shown in (b), the wave field snapshot obtained by the traditional fine-grid spectral element method can reflect the wave field propagation process under complex undulating surface and low-velocity deceleration layer conditions in a relatively fine manner. The wavefront morphology is more complete and the wave field energy distribution is more reasonable, which can be used as a reference standard for high-precision simulation results.

[0111] like Figure 3 As shown in (c), the wavefield snapshot obtained by the variable grid spectral element method of this invention and Figure 3 The wavefield results of the traditional fine-grid spectral element method shown in (b) are highly consistent, indicating that the method of the present invention can accurately recover the seismic wave propagation process under complex undulating near-surface and low-velocity deceleration layer conditions, and has high simulation accuracy.

[0112] like Figure 4 Figure (a) shows a comparison of single-channel records obtained by the traditional coarse-grid spectral element method, the traditional fine-grid spectral element method, and the variable-grid spectral element method of this invention. As can be seen from the figure, the single-channel records obtained by the variable-grid spectral element method of this invention are largely consistent with the results of the traditional fine-grid spectral element method in terms of arrival time, amplitude, and waveform characteristics, indicating that the method of this invention can obtain seismic records with high precision similar to those obtained by the fine-grid method. In contrast, there are significant differences between the records obtained by the traditional coarse-grid spectral element method and the fine-grid results, reflecting its insufficient accuracy under complex model conditions.

[0113] like Figure 4As shown in (b) of the table, the relative and absolute error analysis results of the three methods are presented. The accuracy and efficiency analysis of the examples under different mesh partitions is shown in Table 1. Quantitative analysis shows that the normalized root mean square error (NRMSE) of the variable mesh spectral element method of this invention compared to the traditional fine mesh spectral element method is 0.103%, indicating that the method of this invention has extremely high consistency with the high-precision results of the fine mesh method; the NRMSE of the traditional coarse mesh spectral element method compared to the traditional fine mesh spectral element method is 4.720%. These results indicate that although the traditional coarse mesh spectral element method can reduce the computational load, its accuracy is significantly insufficient, while the method of this invention can significantly improve the reliability of the simulation results.

[0114] like Figure 5 The image shows a comparison between actual seismic records from a certain work area and seismic records simulated using the variable grid spectral element method. For example... Figure 5 (a) in the image represents an actual, complex loess plateau work area's earthquake record. Figure 5 (b) in the figure shows the seismic record of the model in this work area when the variable grid spectral element method is applied. The first arrival wave and the seismic record in the figure are highly similar in morphology, which can be considered as good simulation effect. Figure 5 In (a), the signal-to-noise ratio of the actual work area is low. Figure 5 The results obtained from simulation (b) provide a certain reference for analyzing actual data.

[0115] Table 1

[0116]

[0117] Further analysis of the computational efficiency of the three methods reveals that the computation time for the traditional coarse-grid spectral element method is 986.74 s, for the traditional fine-grid spectral element method it is 21202.80 s, and for the variable-grid spectral element method of this invention it is 14650.23 s. Compared to the traditional fine-grid spectral element method, the computation time of the variable-grid spectral element method of this invention is reduced by 6552.57 s, and its computation time is only 69.09% of that of the traditional fine-grid spectral element method. This demonstrates that the method of this invention, while maintaining approximately the same simulation accuracy as the fine-grid spectral element method, can significantly reduce computational costs and improve computational efficiency.

[0118] comprehensive Figure 3 and Figure 4 The results show that the variable-grid spectral element method proposed in this invention, under the conditions of undulating near-surface and low-velocity deceleration layers in complex loess plateau regions, can obtain high-precision wavefield simulation results that are highly consistent with the traditional fine-grid spectral element method, while significantly improving computational efficiency. Compared with the traditional coarse-grid spectral element method, the method of this invention significantly reduces simulation errors. Therefore, the method of this invention can effectively achieve a balance between accuracy and computational efficiency in seismic wavefield simulation under complex surface conditions, and has good application effects and engineering promotion value.

[0119] It should be noted that those skilled in the art will understand that the present invention is not limited to the specific embodiments described herein, and various obvious changes, readjustments, and substitutions can be made by those skilled in the art without departing from the scope of protection of the present invention. Therefore, although the present invention has been described in detail through the above embodiments, the present invention is not limited to the above embodiments, and may include more other equivalent embodiments without departing from the concept of the present invention, and the scope of the present invention is determined by the scope of the appended claims.

Claims

1. A variable grid spectral element method for simulating seismic wavefields in complex loess plateau regions, characterized in that, Includes the following steps: a. Read in the two-dimensional velocity model and corresponding parameter information of the complex loess plateau region, perform parameter identification, and establish a basic model for two-dimensional seismic wavefield simulation; b. Based on velocity differences, according to the two-dimensional seismic wavefield simulation basic model established in step a, set the source, number of time sampling points, time sampling step size, absorbing boundary and detector parameters, which need to meet the CFL stability condition. c. The overall area is Based on the magnitude and variation of medium velocity in each region, a variable mesh is generated. While keeping the number of horizontal meshes constant, the vertical mesh step size is adjusted to establish a conformal quadrilateral spectral element model. The mesh step size satisfies the CFL stability condition. d. Generate GLL nodes within the spectral element based on the time sampling step size and grid step size; e. Obtain the assignment parameter file, and use GLL nodes to assign node values ​​to the longitudinal wave velocity, transverse wave velocity and density in the two-dimensional velocity model to generate the medium parameter file. f. Assemble the mass matrix and stiffness matrix based on the spectral element model, apply source loading conditions and boundary conditions, and use the time-domain stepwise integration method to iteratively solve the wavefield to obtain the two-dimensional seismic wavefield propagation results. g. Update the displacement field, output the detector seismic record and seismic wave field propagation results obtained in step f, and analyze the simulation results to characterize the seismic wave propagation law under the complex loess plateau undulating surface and non-uniform medium conditions.

2. The variable grid spectral element method for simulating seismic wavefields in complex loess plateau regions according to claim 1, characterized in that, Step a, the parameter information includes the model's lateral length, longitudinal length, undulating surface parameters, medium density parameters, and medium wave velocity parameters. The undulating surface parameters are used to describe the surface elevation variation characteristics of the complex loess plateau region, and the medium wave velocity parameters include longitudinal wave velocity and transverse wave velocity, which are used to characterize the underground velocity stratification structure.

3. The variable grid spectral element method for simulating seismic wavefields in complex loess plateau regions according to claim 1, characterized in that, Step b: The source is selected as the Ricker wavelet to determine the excitation location and mode of the seismic wave. The number of time sampling points and the time sampling step size are used to determine the time domain calculation range. The absorbing boundary is used to attenuate artificial reflected waves at the boundary of the calculation area. The detector parameters are used to receive seismic response information during the propagation of the received wavefield. The source function is: ； In the formula, M is the seismic moment tensor, a second-order tensor describing the intensity and properties of the earthquake source; Indicates the spatial distribution of earthquake sources. The location of the earthquake source is S(t); S(t) is the source time function, which describes the change of the seismic wave over time. The stability conditions for CFL are as follows: ； In the formula, V is the distance between two adjacent GLL integration points; p represents the P-wave velocity in the medium; C is the CFL condition number, which gradually decreases as the wave velocity ratio increases.

4. The variable grid spectral element method for simulating seismic wavefields in complex loess plateau regions according to claim 3, characterized in that: The time step is determined based on the lowest velocity and the smallest grid size within the computational region.

5. The variable grid spectral element method for simulating seismic wavefields in complex loess plateau regions according to claim 1, characterized in that, Step c, the variable mesh generation is based on the medium velocity distribution. A smaller longitudinal mesh step size is used in low-velocity regions and regions with large velocity variations, while a larger longitudinal mesh step size is used in high-velocity regions and regions with small velocity variations. The established spectral elements are conformal quadrilateral elements, and each element is connected through a common boundary. In elastic media, the propagation behavior of seismic waves is described by the linear elastic wave equation: ； In the formula, ρ is the density of the medium, u is the displacement, σ is the stress, and f is the source. For spatial gradient; Using the weighted residual principle, the differential equation is transformed into an integral equation to solve the wave equation; a test function is introduced. In the overall region Integrating within: ； By performing integration by parts, we obtain the weak form of the wave equation: ； In the formula, It is the unit vector of the outward normal on the free surface boundary conditions.

6. The variable grid spectral element method for simulating seismic wavefields in complex loess plateau regions according to claim 5, characterized in that: Step c: The variable mesh partitioning divides the computational domain into layers based on the velocity distribution of the subsurface medium. The low-velocity layer uses a smaller vertical mesh step size, while the high-velocity layer uses a larger vertical mesh step size. The variable mesh partitioning maintains a consistent horizontal mesh step size throughout the entire computational domain, adjusting only the vertical mesh step size. Each spectral element is directly connected through a common boundary, and no transition region is set between coarse and fine meshes. The ratio of the longer side to the shorter side of the spectral element is no greater than 4; for complex near-surface areas, the grid step size is set based on the maximum near-surface layer thickness, and the grid is naturally locally refined in areas with severe topographic relief.

7. The variable grid spectral element method for simulating seismic wavefields in complex loess plateau regions according to claim 1, characterized in that, Step d: GLL nodes are interpolation and integration nodes in the spectral element unit, used to implement higher-order interpolation and numerical integration in the spectral element method. The time step and spatial discrete step size should satisfy the CFL stability condition; GLL nodes are the n+1 roots of the following equation: ； In the formula, for The first derivative of a Legendre polynomial of order 1, for The order Lagrange interpolation basis functions are expressed as follows: ； In the formula, Lagrange interpolation basis functions This indicates that the basis function corresponds to the i-th node. x is the independent variable.

8. The variable grid spectral element method for simulating seismic wavefields in complex loess plateau regions according to claim 1, characterized in that, Step e: The medium parameter file includes P-wave velocity, S-wave velocity, and density values ​​at each node, used to characterize the medium distribution characteristics at different locations within the computational domain; for numerical integration on the elements, GLL integration is used for calculation. ； In the formula, ω α ω β α, β = 0, and n is the weight coefficient of the GLL integral.

9. The variable grid spectral element method for simulating seismic wavefields in complex loess plateau regions according to claim 1, characterized in that, Step f: The mass matrix and stiffness matrix are assembled from the geometric information, nodal information and medium parameters of the spectral element. The boundary conditions include the free surface boundary conditions of the undulating ground surface and the absorbing boundary conditions of the perfectly matched layer. The wave field iterative solution is used to obtain the displacement field and velocity field distribution results at each time step. For the dielectric unit Ω e The solution is: ； ； in, , , , The time iteration result is: ； in, For the quality matrix, Let K be the velocity vector for the next time step, K be the global stiffness matrix, and F be the source term.

10. The variable grid spectral element method for simulating seismic wavefields in complex loess plateau regions according to claim 1, characterized in that, Step g: The geophone seismic record includes information on the displacement response and velocity response at each geophone location over time.

Images