A solid earthquake wave field simulation method based on an elastic multi-relaxation lattice Boltzmann model

Abstract

The application relates to the cross field of geophysical exploration and computational acoustics, in particular to a solid seismic wave field simulation method based on an elastic multi-relaxation lattice Boltzmann model. The method solves the problems of poor numerical stability, strong coupling of physical parameters and inability to finely control numerical dispersion in the prior art. The specific steps are as follows: an elastic MLBM model for solid seismic wave simulation is constructed, including establishing a discrete evolution equation, configuring key parameters, initializing variables; a source function and an elastic force form are defined, and a total external force term is calculated; a balanced state distribution function is designed, and a moment vector thereof is calculated; a collision operation is performed in a moment space; a migration operation is performed in a velocity space, and a boundary condition is applied; macroscopic quantities are calculated, a displacement field is updated, and a strain and stress tensor are calculated; seismic data are recorded, and visualization and result verification are performed. The application constructs a new paradigm for solid seismic wave field simulation, and provides a powerful numerical tool for geophysical exploration, seismic engineering and basic wave theory research.

Description

Technical Field

[0001] This invention relates to the interdisciplinary field of geophysical exploration and computational acoustics, specifically to a method for simulating solid seismic wavefields based on the elastic multi-relaxation lattice Boltzmann model. Background Technology

[0002] Numerical simulation of seismic wavefields is a core computational tool for understanding subsurface structures, aiding resource exploration, and assessing seismic risk. Its essence lies in numerically solving the elastic dynamics wave equations to simulate the propagation, reflection, transmission, and conversion processes of seismic waves (including P-waves and S-waves) in complex geological media. Traditional simulation methods are primarily based on macroscopic continuum frameworks, such as the finite difference method (FDM), finite element method (FEM), and spectral element method (SEM). While these methods are mature and widely used, they still face challenges in handling strongly contrasting interfaces, complex geometric boundaries, and achieving high parallel scalability.

[0003] To address the aforementioned challenges, the Lattice Boltzmann Method (LBM), derived from mesoscopic kinetic theory, offers an alternative numerical simulation approach. Instead of directly solving macroscopic partial differential equations, this method recovers macroscopic dynamic behavior from the bottom up by simulating the two fundamental steps of "collision" and "migration" of virtual particles on a regular lattice. This method possesses inherent advantages such as completely local algorithms, easily handled boundary conditions, and extremely high parallel efficiency, and has been widely applied in computational fluid dynamics. Extending the LBM method to solid mechanics requires constructing an LBM model suitable for elastic media. Its core lies in redesigning the particle equilibrium distribution function so that its macroscopic moments correspond not only to mass and momentum but also to the macroscopic stress tensor; simultaneously, designing collision operators adapted to elastic constitutive relations to recover the Navier-Cauchy equations of motion, i.e., the macroscopic elastic wave equations, during the system's evolution.

[0004] Early Luke-Range Seismic Models (MLBMs) ​​mostly employed a single relaxation time (BGK) model, where the relaxation process of all physical modes is controlled by a single relaxation time parameter. This leads to several drawbacks in solid seismic wavefield simulations: poor numerical stability, rigid coupling of physical parameters, and inability to finely control numerical dispersion. To address the fundamental shortcomings of the BGK model, the multiple relaxation time (MRT) framework was introduced into LBMs. In MLBMs, the collision process occurs in moment space, and moments with different physical meanings (such as conserved moments, stress moments, and higher-order ghost moments) have their own independent relaxation parameters. This method has been proven in fluid simulations to significantly improve numerical stability, optimize dispersion relations, and allow for more flexible boundary treatments, but its application in solid seismic wave simulations remains a research gap.

[0005] In summary, existing technologies still have the following limitations: On the one hand, traditional methods based on macroscopic continuum theory have bottlenecks in parallel efficiency and complex boundary handling; on the other hand, the LBM method, which has natural parallel advantages, cannot meet the requirements of high-precision seismic wavefield simulation due to the instability and flexibility of its traditional BGK model. The multiple relaxation time (MRT) framework has been proven to significantly improve numerical stability, optimize dispersion relations, and support more flexible boundary handling in fluid simulation, but its systematic construction and application in elastic solid wavefield simulation is still a research gap.

[0006] Therefore, it is urgent to construct an elastic multi-relaxation lattice Boltzmann (MRT-LBM) algorithm specifically for solid seismic wavefield simulation. This algorithm will deeply integrate the multi-relaxation time mechanism, which has been proven effective in fluid mechanics, with the elastic LBM framework applicable to solid dynamics, in order to overcome the shortcomings of existing methods in terms of numerical stability, physical parameter decoupling, and fine control of numerical dispersion. Summary of the Invention

[0007] To address the aforementioned technical problems, this invention proposes a solid seismic wavefield simulation method based on the elastic multi-relaxation lattice Boltzmann model, aiming to overcome the shortcomings of existing technologies in terms of numerical stability, physical parameter decoupling, and fine control of numerical dispersion.

[0008] To achieve the above objectives, the present invention provides a method for simulating solid seismic wavefields based on the elastic multi-relaxation lattice Boltzmann model, which specifically includes the following steps:

[0009] (1) Construct an elastic MLBM model for solid seismic wave simulation, including establishing its discrete evolution equation, configuring key parameters, and initializing each field variable;

[0010] (2) Define the source function and the form of elastic force, and calculate the total external force term;

[0011] (3) Design the equilibrium distribution function, calculate its equilibrium moment vector, and perform a collision operation in the moment space;

[0012] (4) Perform the migration operation in the velocity space and apply boundary conditions;

[0013] (5) Calculate macroscopic quantities, update the displacement field, and calculate strain and stress tensor;

[0014] (6) Record earthquake data and visualize and verify the results.

[0015] This invention employs a multi-relaxation-time (MRT) framework as an innovative mechanism. Traditional single-relaxation (BGK) models use only a single relaxation parameter to control all physical modes, resulting in poor numerical stability, tight coupling of physical parameters, and an inability to independently control the dissipation and dispersion behavior of P-waves and S-waves. This invention, by performing a collision operation in moment space and introducing a diagonal relaxation matrix, allows moments with different physical meanings (such as conserved moments, stress moments, and higher-order moments) to have independent relaxation parameters. This mechanism enables the attenuation characteristics, dispersion behavior, and boundary absorption effects of P-waves (P-waves) and S-waves (S-waves) to be finely controlled separately. This significantly improves the physical realism and computational flexibility of the simulation while ensuring numerical stability, making it particularly suitable for simulating seismic wavefields in strongly contrasting media, complex structures, and long-term propagation environments.

[0016] Compared with the prior art, the beneficial effects of the present invention are as follows:

[0017] This invention extends the advanced LBM framework, which is mature in fluid simulation, to solid seismic wavefield simulation. It provides a novel numerical simulation paradigm with different principles, parameter decoupling capabilities, and high parallel efficiency, which addresses the challenges faced by traditional macroscopic methods when dealing with extremely complex structures, multi-scale problems, and multi-physics coupling.

[0018] This invention employs multiple relaxation time frames, allowing the relaxation parameters controlling wavefield dissipation, dispersion, and numerical stability to be independently adjusted. This enables more flexible matching of the physical properties of actual media while ensuring computational stability, and effectively suppresses non-physical oscillations. It is particularly suitable for simulating strongly contrasting media or long-term propagation, exhibiting excellent numerical stability and flexibility. By introducing multiple relaxation time frames, complete decoupling of physical parameters is achieved, allowing independent control of the numerical characteristics of P-waves and S-waves, overcoming the fundamental limitations of traditional single-relaxation models in solid wavefield simulation.

[0019] This invention transforms complex macroscopic boundary conditions such as free surfaces and layered interfaces into local mesoscopic operations on grid point distribution functions. It can spontaneously generate reflected waves, transmitted waves, and wave type transformations that conform to macroscopic theory without explicit connection conditions. Since migration, collision, and boundary processing are all based on pure local grid point operations, there is no need for global matrix solving, and the data access pattern is highly regular. This algorithm has natural inherent parallelism and can be efficiently deployed on high-performance parallel computing platforms such as GPUs or many-core processors. It is suitable for large-scale three-dimensional seismic wavefield simulation. Attached Figure Description

[0020] The invention will now be further described with reference to the accompanying drawings.

[0021] Figure 1 This is a flowchart of the present invention;

[0022] Figure 2This is a schematic diagram of the velocity and density of the two-dimensional homogeneous solid medium in Example 1;

[0023] Figure 3 In a two-dimensional homogeneous solid medium at t=180ms Comparison of component wavefields: (a) High-precision FDM results, (b) Elastic MLBM results, (c) Residuals of the two;

[0024] Figure 4 Comparison of component wave profiles of high-precision FDM and elastic MLBM at t=180ms: (a) depth z=530m, (b) depth z=400m;

[0025] Figure 5 At a depth of z=470m in a two-dimensional homogeneous solid medium Comparison of component common shot gathers: (a) High-precision FDM results, (b) Elastic MLBM results, (c) Residuals of the two;

[0026] Figure 6 For high-precision FDM and elastic MLBM at different detector points Earthquake trace comparison diagram: (a) x=550m, z=470m, (b) x=400m, z=470m;

[0027] Figure 7 This is a schematic diagram of the velocity and density of the two-dimensional layered solid medium in Example 2;

[0028] Figure 8 In a two-dimensional layered solid medium at t=360ms Comparison of component wavefields: (a) High-precision FDM results, (b) Elastic MLBM results, (c) Residuals of the two;

[0029] Figure 9 At a depth of z=1040m in a two-dimensional layered solid medium Comparison of component common shot gathers: (a) High-precision FDM results (b) Elastic MLBM results (c) Residuals of the two;

[0030] Figure 10 Comparison of high-precision FDM and elastic MLBM wavefield results in two-dimensional layered media: (a) at t=360ms, z=1250m Waveform profile, (b) receiver point x=800m, z=1040m Earthquake path. Detailed Implementation

[0031] The present invention will be further described below with reference to the accompanying drawings and specific embodiments.

[0032] like Figure 1As shown, a method for simulating solid seismic wavefields based on the elastic multi-relaxation lattice Boltzmann model includes the following specific implementation steps:

[0033] S1: Model Establishment and Initialization

[0034] An elastic MLBM model for solid seismic wave simulation is constructed, including establishing discrete equations, setting various parameters, and initializing field variables to ensure that the model evolves from physically consistent and undisturbed initial conditions.

[0035] In practical implementation, the first step is to establish the discrete equations for the elastic MLBM:

[0036]

[0037] in, For time ,Location The particle distribution function along the discrete velocity direction i. It is a discrete velocity vector. The transformation matrix is... and These are the particle distribution functions. and equilibrium distribution function In the representation of moment space, It is a diagonal relaxation matrix. It is the form of the external force term in velocity space.

[0038] Based on this, configure three types of key parameters:

[0039] Physical parameters: including reference density Longitudinal wave velocity transverse wave velocity Lamé constant and Equal amounts;

[0040] Numerical parameters: including computational domain size Total number of steps Spatial step size and Time step (CFL conditions must be met) Equal quantities;

[0041] Elastic MLBM parameters: including discrete velocity pairs of the discrete velocity model. and weight Grid speed of sound relaxation matrix Number of lattice cells migrated per unit time Equal amounts.

[0042] Then, the variables are initialized, including Particle distribution function at time t Macro density momentum density velocity field Displacement field strain field and stress field Equal field quantity.

[0043] By establishing discrete evolution equations, configuring physical and numerical parameters, and initializing field variables, a stable and physically consistent numerical model foundation is provided for subsequent wave field simulations, ensuring the numerical stability and physical authenticity of the simulation process.

[0044] S2: Earthquake source excitation and external force calculation

[0045] Seismic waves are generated at a designated hypocenter location, and external force terms are calculated simultaneously. The specific implementation steps are as follows:

[0046] First, determine the source type (point source, explosion source, etc.), location, time function, and detector placement. The source time function uses the Ricker wavelet, the most common seismological component, and its expression is:

[0047]

[0048] in, It is the amplitude factor. It is the dominant frequency of the wavelet, and t is the time variable. It is the wavelet delay time.

[0049] The appropriate elastic force is set based on the above source function as follows:

[0050] ,

[0051] in, and It is Lamé constant, It is a reference density. Yes Directional partial derivative operator It's density.

[0052] Subsequently, the total external force field is calculated based on the elastic force and the focal force at the current time:

[0053] ,

[0054] in, It is the seismic source force. It is elastic force.

[0055] The total external force term in the discrete equations of elastic MLBM is expressed as follows:

[0056] ,

[0057] in, These are the weights of the discrete velocity direction i. Discrete velocity vector in directional components, It is the total external force. directional components, It is the square of the speed of sound in the grid.

[0058] By defining the source function and the form of elastic force, the total external force term is calculated and discretized, providing a physical driving source for wave field simulation and ensuring that seismic waves can be correctly excited and propagated in the model.

[0059] S3: Equilibrium State Construction and Moment Space Collision

[0060] To enable the elastic MLBM to recover the wave equation of the elastic solid, the equilibrium distribution function was redesigned. This invention employs an extended form of the equilibrium distribution function:

[0061]

[0062] in, It is the weight of the discrete velocity direction i. It is density, It is the equilibrium state of momentum. Components of direction and It is the component of the discrete velocity vector in different directions. It is the speed of sound of the grid. It is the stress tensor component, It is the Kronecker symbol.

[0063] The momentum term is adopted in the following corrected form:

[0064]

[0065] in, It is momentum in directional components, It is the time step. It is the total external force. The directional component.

[0066] This dynamic term is compatible with the non-zero stress response of solid media caused by shear deformation.

[0067] The constructed equilibrium distribution function consists of three parts: the density term... Momentum term With key stress terms This ensures that its second moment matches the macroscopic stress, thus embedding a linear elastic constitutive relation at the mesoscopic level. The corresponding equilibrium state is updated in moment space as follows:

[0068]

[0069] in, It's density. and These are the components of momentum in different directions. It is a reference density. It is the time step. and It represents the components of the total external force in different directions. It is the square of the speed of sound in the grid. , and The components of the stress tensor in different directions.

[0070] Based on this, the velocity space particle distribution function is transformed using spatial moment transformation. Transform to moment space And carry out a physically meaningful collision process in that space:

[0071]

[0072] in, It is the moment vector after the collision. It is the moment vector before the collision. It is the equilibrium moment vector. It is the time step. It is the identity matrix. The relaxation matrix represents the external force term in moment space. The parameters can be adjusted independently within the range of (0,2).

[0073] The innovation of this step lies in the introduction of a multiple relaxation time (MRT) framework. This addresses the problem in traditional BGK models that use a single relaxation time, leading to strong coupling of all physical processes and difficulty in independently controlling numerical dissipation and dispersion. This invention constructs a diagonal relaxation matrix... Each relaxation parameter can be set independently, corresponding to the relaxation rate of mass, momentum, stress, and their higher-order moments, respectively. For example, by adjusting the relaxation parameters related to stress moment, the numerical attenuation and dispersion characteristics of P-waves (P-waves) and S-waves (S-waves) can be independently controlled; by adjusting the parameters related to higher-order "ghost moments," numerical stability and boundary absorption effects can be further optimized. This parameter decoupling mechanism allows the invention to flexibly adapt to the physical properties of different geological media, exhibiting superior numerical stability and fidelity in complex wavefield simulations. This multi-relaxation mechanism allows the attenuation characteristics, dispersion behavior, and boundary absorption effects of P-waves and S-waves to be controlled separately, ensuring numerical stability while flexibly adapting to the physical properties of complex geological media. After the collision is completed, an inverse transformation is performed... Return to velocity space.

[0074] Therefore, by designing an equilibrium distribution function that includes stress terms and performing a collision operation with multiple relaxation parameters in the moment space, a rigorous recovery of the elastic wave equation can be achieved. At the same time, the attenuation and dispersion characteristics of longitudinal and transverse waves can be independently controlled, thereby improving numerical stability and physical flexibility.

[0075] S4: Velocity Space Migration and Application of Boundary Conditions

[0076] In velocity space, the distribution function updated after the collision is migrated to adjacent grid points along its discrete velocity direction:

[0077]

[0078] in, It is the particle distribution function vector after migration. It is the particle distribution function vector after a collision in moment space.

[0079] Subsequently, corresponding boundary conditions were applied according to the simulation requirements: free surface conditions (stress of 0) were achieved by using non-equilibrium extrapolation; absorbing boundary conditions (very weak truncation boundary reflection) were achieved by loading a perfectly matched layer or a sponge layer; and fixed boundary conditions (displacement of 0) were achieved by using a bounce scheme.

[0080] The propagation of the wave field between grid points is achieved through migration operations. Combined with boundary handling methods such as non-equilibrium extrapolation, perfectly matched layers or sponge layers, and bounce schemes, the wave field is simulated under actual physical conditions such as free surfaces, absorbing boundaries, and fixed boundaries to ensure the correct behavior of the wave field under complex boundaries.

[0081] S5: Macroscopic quantity calculation, displacement field update and stress calculation

[0082] Macroscopic quantities with clear physical meaning are extracted from the particle distribution function, and strain-stress relationships are explicitly constructed based on displacement, thereby achieving a complete description of the elastic wave field.

[0083] First, the density is calculated sequentially based on the particle distribution function after migration. Momentum density and speed Macroscopic quantities:

[0084]

[0085]

[0086]

[0087] in, It is the particle distribution function. It is a discrete velocity vector in A vector of direction.

[0088] Based on this, the displacement field is updated using time-domain integration. :

[0089]

[0090] in, It is a displacement field. It is a velocity field. This is the time step. The displacement field must be explicitly stored for strain calculation.

[0091] Subsequently, the displacement gradient is calculated using the central difference method, and the strain tensor is obtained. :

[0092]

[0093] in, and Yes and Directional partial derivative operator, and Displacement components.

[0094] by Taking the component as an example, its discrete form is:

[0095]

[0096] in, It is a location displacement in the x direction, It is the spatial step size.

[0097] Finally, based on the generalized Hooke's law for isotropic linear elastic media, the Cauchy stress tensor is calculated from the strain tensor:

[0098]

[0099] in, and It is Lamé's constant. and These are the components of the strain tensor in different directions. This is the Kronecker notation. To further adapt to the dimensionless LBM frame, the stress is converted to a dimensionless form:

[0100]

[0101] in, It is the component of the Cauchy stress tensor. It is a reference density. and It is Lamé constant, and These are the components of the strain tensor in different directions. It is the Kronecker symbol.

[0102] By extracting macroscopic quantities such as density, momentum, and velocity from the distribution function after migration and boundary processing, updating the displacement field, and calculating strain and stress tensors, a physical quantity mapping from mesoscopic particle evolution to macroscopic continuum mechanics is achieved, providing a direct data source for seismic records.

[0103] S6: Seismic Data Recording, Visualization, and Result Verification

[0104] After the simulation, velocity and stress fields were recorded at preset receiver points to form multi-component seismic gathers. Systematic verification was performed by comparing analytical solutions or high-precision numerical solutions: quantitative analysis of P / S wave travel time errors, amplitude attenuation, and phase dispersion was conducted to verify the viscoelastic response of the medium; qualitative examination of wavefront morphology, polarization characteristics, and converted waves generated by complex interfaces was performed in the wavefield snapshots. Finally, the verified data were synthesized into seismic profiles, common-shot gathers, and wavefield snapshots, providing a reliable basis for subsequent geophysical interpretation.

[0105] S7: Specific Implementation

[0106] (1) Example of two-dimensional homogeneous solid medium

[0107] To verify the baseline accuracy of the method of this invention, a two-dimensional homogeneous solid medium model was first tested, and the results of the eighth-order finite difference method (FDM), which has a small time sampling step and high spatial discretization accuracy, were used as a reference solution for comparison. The model size was 1000m × 1000m, and a 1m × 1m grid was used for discretization. The medium parameters were set as follows: P-wave velocity 2400m / s, S-wave velocity 1380m / s, and density 1200kg / m³. The seismic source was located at the center of the model (e.g., Figure 2As shown in the figure, the dominant frequency 30Hz Ricker wavelet acts on the vertical velocity component. The total simulation time is 180ms, the time step of the elastic MLBM method of this invention is 0.5ms, and the time step of high-precision FDM is 0.05ms.

[0108] Figure 3 The horizontal velocity components calculated by the high-precision FDM algorithm and the elastic MLBM algorithm of this invention at time t=180ms are shown. Wavefield snapshots and their residuals. It can be seen that the wavefront profiles of the P-wave and S-wave simulated by the method of this invention are clear and smooth, and highly coincide with the reference solution in terms of energy distribution and spatial position. The residual amplitude is much smaller than the signal amplitude, indicating that this method can accurately decouple and propagate longitudinal and transverse waves in a homogeneous medium. Figure 4 The results further demonstrate the comparison of waveform profiles extracted at depths of 530m and 400m at t=180ms. The results show that the waveform profiles obtained by the two methods maintain a high degree of consistency in travel time, amplitude envelope, and waveform phase, intuitively confirming the fidelity of the method of this invention in terms of kinematic and dynamic characteristics. Figure 5 A comparison of the common shot gather at a depth of z=470m is shown. The gather simulated by the method of this invention matches well with the high-precision FDM gather in terms of the first arrival time of P-waves and S-waves, the phase axis morphology, and the amplitude decay trend with offset. Moreover, the residual gather has weak energy, which proves the reliability of the method in simulating long-distance propagation of wave fields. Figure 6 Velocity vibration maps from two representative receivers were selected for point-to-point comparison. Throughout the recording time, the waveform output by the method of this invention almost completely overlapped with the reference solution, further quantitatively verifying its extremely high local accuracy. Therefore, in the benchmark test of a homogeneous medium, the method of this invention passed the comprehensive verification of the high-precision reference solution, and its accuracy, fidelity, and physical correctness in wavefield simulation were fully demonstrated.

[0109] (2) Example of two-dimensional layered solid medium

[0110] To verify the ability of the method of this invention to handle non-homogeneous media and physical interfaces, a horizontal two-layer solid medium model was further designed. The model size is 2000m × 2000m, with a grid step size of 1m × 1m. With a depth of z = 1200m as the boundary, the upper medium has a P-wave velocity of 2200m / s, a S-wave velocity of 1270m / s, and a density of 1000kg / m³; the lower medium has a P-wave velocity of 3000m / s, a S-wave velocity of 1730m / s, and a density of 1200kg / m³. The seismic source is located in the upper medium and is a vertical velocity point source with a dominant frequency of 35Hz. The total simulation time is 360ms.

[0111] Figure 8A comparison of wavefield snapshots at t=360ms is shown. It is evident that in the upper medium, the direct P-waves and S-waves simulated by the two methods are perfectly matched. After passing through the interface, the overall wavefield structures of the generated reflected and transmitted waves are highly similar, although there are slight differences in local amplitude and phase. This stems from the different underlying principles of the interface processing mechanisms of the two algorithms. Importantly, through the results of the method of this invention (… Figure 8 (b) Detailed annotations were performed, identifying all theoretical seismic phases, including direct P-waves, S-waves, reflected waves (RPP, RPS, RSS, RSP), and transmitted waves (TPP, TPS, TSS, TSP), where RPS, RSP, TPS, and TSP are converted waves. This fully demonstrates that the method of this invention can completely and accurately simulate the complex physical processes of wave reflection, transmission, and wave mode conversion at the interface. Figure 9 A comparison of common shot gathers at a depth of z=1040m was performed. The results of the two algorithms showed a high degree of consistency in the travel time of the phase axes of the direct wave from the upper medium and the reflected wave from the interface, further verifying the kinematic accuracy of the method of this invention. Figure 10 Further comparisons between deep wave profiles and single-point seismograms are presented, showing that the arrival times and waveform characteristics of the main seismic phases are in good agreement.

[0112] This example demonstrates that the elastic MLBM algorithm proposed in this invention possesses the ability to accurately simulate seismic wave propagation in non-homogeneous solid media and the complex interactions between waves and geological interfaces. The fundamental difference between this algorithm and high-precision FDM in terms of interface response details lies in their methodological principles: this invention is based on mesoscopic particle evolution, with interface effects implicitly implemented through lattice rules; while high-precision FDM is based on macroscopic equations, requiring explicit application of continuity conditions. This mesoscopic characteristic makes the method of this invention more physically intuitive and has greater parallel potential when simulating complex interfaces.

[0113] The embodiments of the present invention have been described in detail above, but the content described is only a preferred embodiment of the present invention and should not be considered as limiting the scope of the present invention. All equivalent changes and improvements made within the scope of the present invention should still fall within the patent coverage of the present invention.

Claims

1. A method for simulating solid seismic wavefields based on the elastic multi-relaxation lattice Boltzmann model, characterized in that, Includes the following steps: S1: Construct an elastic MLBM model for solid seismic wave simulation, including establishing its discrete evolution equations, configuring key parameters, and initializing field variables; S2: Define the source function and elastic force form, and calculate the total external force terms; S3: Design the equilibrium distribution function, calculate its equilibrium moment vector, and perform a collision operation in the moment space; S4: Perform a migration operation in the velocity space and apply boundary conditions; S5: Calculate macroscopic quantities, update the displacement field, and calculate strain and stress tensor; S6: Record earthquake data and visualize and verify the results.

2. The method for simulating solid seismic wavefields based on the elastic multi-relaxation lattice Boltzmann model according to claim 1, characterized in that, The discrete evolution equation in S1 is as follows: in, For time ,Location The particle distribution function along the discrete velocity direction i. It is a discrete velocity vector. The transformation matrix is... and These are the particle distribution functions. and equilibrium distribution function In the representation of moment space, It is a diagonal relaxation matrix. It is the form of the external force term in velocity space.

3. The method for simulating solid seismic wavefields based on the elastic multi-relaxation lattice Boltzmann model according to claim 1, characterized in that, The key parameters in S1 include: physical parameters, numerical parameters, and elastic MLBM parameters.

4. The method for simulating solid seismic wavefields based on the elastic multi-relaxation lattice Boltzmann model according to claim 1, characterized in that, The formula for calculating the total external force term in S2 is as follows: ， in, It is the seismic source force. It is elastic force; The expression for the elastic force is as follows: ， in, and It is Lamé constant, It is a reference density. Yes Directional partial derivative operator It is density; The expression for the focal force is as follows: in, It is the amplitude factor. It is the dominant frequency of the wavelet, and t is the time variable. It is the wavelet delay time; The total external force term is applied in the elastic MLBM discrete equation in the following form: ， in, These are the weights of the discrete velocity direction i. Discrete velocity vector in directional components, It is the total external force. directional components, It is the square of the speed of sound in the grid.

5. The method for simulating solid seismic wavefields based on the elastic multi-relaxation lattice Boltzmann model according to claim 1, characterized in that, The equilibrium distribution function in S3 is as follows: in, It is the weight of the discrete velocity direction i. It is density, It is the equilibrium state of momentum. Components of direction and It is the component of the discrete velocity vector in different directions. It is the speed of sound of the grid. It is the stress tensor component, It is the Kronecker symbol.

6. The method for simulating solid seismic wavefields based on the elastic multi-relaxation lattice Boltzmann model according to claim 1, characterized in that, The formula for calculating the equilibrium moment vector in S3 is as follows: in, It's density. and These are the components of momentum in different directions. It is a reference density. It is the time step. and It represents the components of the total external force in different directions. It is the square of the speed of sound in the grid. , and The components of the stress tensor in different directions.

7. The method for simulating solid seismic wavefields based on the elastic multi-relaxation lattice Boltzmann model according to claim 1, characterized in that, The collision operation in S3 is implemented by the following equation: in, It is the moment vector after the collision. It is the moment vector before the collision. It is the equilibrium moment vector. It is the time step. It is the identity matrix. Let be the representation of the external force term in moment space. Represents the relaxation matrix; the relaxation matrix The value of is in the range of 0 to 2, and is used to control the attenuation characteristics, dispersion behavior and boundary absorption effect of longitudinal and transverse waves respectively; after the collision is completed, the updated moment vector is converted back to the velocity space through inverse moment transformation.

8. The method for simulating solid seismic wavefields based on the elastic multi-relaxation lattice Boltzmann model according to claim 1, characterized in that, The migration operation in S4 is implemented by the following equation: in, It is the particle distribution function vector after migration. It is the particle distribution function vector after a collision in moment space.

9. The method for simulating solid seismic wavefields based on the elastic multi-relaxation lattice Boltzmann model according to claim 1, characterized in that, The boundary conditions in S4 include: applying a non-equilibrium extrapolation method to the free surface grid points, and applying a perfect matching layer or a sponge absorption layer to the truncation boundary of the model.

10. A method for simulating solid seismic wavefields based on an elastic multi-relaxation lattice Boltzmann model according to claim 1, characterized in that, The calculation of macroscopic quantities in S5 includes density. ,momentum and speed They are given by the following formulas: in, It is the particle distribution function. It is a discrete velocity vector in The directional component.

Images