Numerical simulation method of channel wave in viscoelastic vti medium based on wrk-ssgfd

Abstract

The application discloses a slot wave numerical simulation method for viscoelastic VTI medium based on WRK-SSGFD, first-order velocity-stress equations of the viscoelastic VTI medium are established, the equations are used for high-order staggered grid finite difference discretization in a space domain to obtain semi-discrete difference equations, the semi-discrete equations are converted into the form of ordinary differential equation groups, optimal weight coefficients are solved by minimizing a dispersion error objective function, finally, geological model parameters, boundary conditions and observation system parameters are set, a fourth-order weighted Runge-Kutta method is used for time advancing, wave field values at each time are calculated in sequence until all time steps are completed, and corresponding simulated seismic records are obtained. According to the specific dispersion condition, the optimal weight coefficients are determined, and then the optimal weight is integrated into the fourth-order weighted Runge-Kutta method for simulation calculation, the method effectively suppresses numerical dispersion and improves simulation precision, and provides data support for subsequent research on slot wave propagation rules in the viscoelastic VTI medium.

Description

Technical Field

[0001] This invention relates to a method for simulating channel waves, specifically a numerical simulation method for channel waves in viscoelastic VTI media based on WRK-SSGFD (Runge-Kutta-Staggered finite difference). Background Technology

[0002] Mine channel wave exploration, as a geological exploration technology in coal mines, can effectively provide important geological information for mine layout and disaster prevention. Channel waves have characteristics such as high energy, slow attenuation, and significant dispersion when propagating in coal seams. Through efficient numerical simulation of channel waves, the generation, propagation, and response mechanism of anomalies of channel waves can be effectively analyzed, which is of great significance for guiding field exploration.

[0003] Coal seams, as sedimentary rocks, exhibit significant bedding structures and are a typical vertically isotropic (VTI) medium. Simultaneously, coal seams are rock masses composed of inorganic and organic matter, and compared to the surrounding rocks, they are softer and possess strong viscoelasticity, resulting in a strong absorption and attenuation effect on seismic waves. Therefore, coal seams are media exhibiting both viscoelastic and vertically isotropic characteristics. However, current channel wave research is largely based on the characteristics of coal seams as elastic isotropic media, which does not align with reality. Consequently, there is a significant discrepancy between the theoretical wavefield characteristics and actual wavefield features in current channel wave research.

[0004] Numerical simulation methods for channel waves encompass various approaches, including the Finite Difference Method (FDM), Finite Element Method (FEM), Galerkin Method, and Spectral Element Method. The Finite Difference Method (FDM) is widely used in channel wave simulation due to its simplicity, anisotropy, and high computational efficiency. However, the classical FDTD method is susceptible to numerical dispersion during channel wave propagation simulations, exhibiting significant waveform distortion and accuracy degradation, especially when dealing with high-frequency bands and complex media interfaces. The root cause lies in the fact that when approximating the wave equation with a discretized difference equation, the phase velocity of wave propagation becomes a function of the discrete spatial step size and time step size. Therefore, when the spatial sampling within each wavelength is too small (i.e., the spatial grid is too coarse) or the velocity contrast between subsurface media layers is large, a "tail" appears after the main phase, commonly known as numerical dispersion. The existence of this pseudo-wave makes it difficult to distinguish between real physical waves and pseudo-waves caused by the numerical method; therefore, eliminating numerical dispersion is of great practical significance.

[0005] To mitigate these shortcomings, researchers have successively proposed improvement strategies such as higher-order spatial difference schemes, weighted residual methods, and higher-order time integration schemes. Among them, the higher-order Runge-Kutta time integration method can effectively reduce time discretization errors and maintain the stability and high accuracy of numerical solutions. Some scholars have also approximated higher-order spatial derivatives of the wave equation by using a combination of grid point displacements, particle velocities, and gradient values ​​in the spatial direction, and used third- or fourth-order Runge-Kutta methods for time-progression calculations, thus obtaining the "Runge-Kutta" method for numerical simulation of seismic wave propagation. In addition, numerical tests have further verified the advantages of the Runge-Kutta method in terms of computational efficiency and numerical dispersion under non-uniform time steps, accurately simulating slow P-waves in porous elastic media, and demonstrating its advantages of good stability, high numerical accuracy, and small dispersion errors through theoretical analysis and numerical experiments. However, the existing Runge-Kutta method has not been applied to the field of slot waves. And if it is directly applied to the numerical calculation of slot waves, although it has the advantages of good stability, high numerical accuracy and small dispersion error, it will greatly increase the amount of computation and data storage.

[0006] Therefore, the research direction of this invention is to provide a new method that can achieve both high accuracy and flexibility in adapting to complex mine structures while maintaining low computational and data storage requirements, thereby effectively improving the accuracy and reliability of the numerical simulation of the channel wave. Summary of the Invention

[0007] To address the problems existing in the prior art, this invention provides a numerical simulation method for groove waves in viscoelastic VTI media based on WRK-SSGFD. The optimal weighting coefficients are determined according to specific dispersion conditions, and then the optimal weights are integrated into the fourth-order weighted Runge-Kutta method for simulation calculation. While maintaining low computational and data storage requirements, this method can effectively suppress numerical dispersion, improve simulation accuracy, and provide data support for subsequent research on the propagation law of groove waves in viscoelastic VTI media.

[0008] To achieve the above objectives, the technical solution adopted by this invention is: a numerical simulation method for groove waves in viscoelastic VTI media based on WRK-SSGFD, comprising the following steps:

[0009] Step 1: Establish the first-order velocity-stress equation for the viscoelastic VTI medium: The first-order velocity-stress equation for the three-dimensional viscoelastic VTI medium is derived from the Kelvin-Voigt constitutive relation, the differential equation of motion, the constitutive equation, and the geometric equation.

[0010] Step 2: Spatial Domain High-Order Staggered Mesh Finite Difference Discretization: Using a staggered mesh and combining it with the equations from Step 1, the velocity components are spatially discretized. The principal stress components are located at the mesh nodes, the shear stress components are located at the mesh center, and the velocity components are located at the mesh boundaries. The spatial derivatives are approximated by a 2M-order optimized finite difference operator to obtain a semi-discrete difference equation.

[0011] Step 3: Time-domain weighted Runge-Kutta integration: Transform the semi-discrete difference equations into a system of ordinary differential equations; use the fourth-order weighted Runge-Kutta method for time progression;

[0012] Step 4: Model Parameter Settings: Establish the geological model and set the corresponding parameters;

[0013] Step 5: Set boundary conditions for the geological model established in Step 4, which will be used for subsequent trough wave numerical simulation.

[0014] Step Six: Setting Observation System Parameters: The source simulation uses a dominant frequency of [missing value]. Delay time The Reich wavelet;

[0015] Step 7: Numerical simulation of trough waves: Based on the parameters set in steps 4 to 6, the fourth-order weighted Runge-Kutta method from step 3 is used to advance the time step, calculating the wave field values ​​at each time step sequentially until all time steps are completed, and obtaining the corresponding simulated seismic records.

[0016] Furthermore, the first-order velocity-stress equation for the three-dimensional viscoelastic VTI medium in step one is as follows:

[0017]

[0018] The first-order velocity-stress equation for a two-dimensional viscoelastic VTI medium considering the source term is:

[0019]

[0020] in, , Let x and z represent the vibration velocities of a point in the x and z directions, respectively. , These represent the principal stresses in the x and z directions at that point, respectively. The shear stress at that point, Represents the elastic coefficient. is the viscosity coefficient.

[0021] Furthermore, the semi-discrete difference equation in step two is specifically as follows:

[0022]

[0023] Furthermore, step three specifically includes:

[0024] The semi-discrete difference equation in step two is expressed as follows:

[0025]

[0026] in, F represents the discrete spatial derivative term. Indicates the source term;

[0027] The fourth-order weighted Runge-Kutta method is used for time advancement: the time step is defined as... ,initialization Solution of time Then calculate the slope vectors for the four stages. Calculated using the following formula:

[0028]

[0029] The node coefficients use standard values:

[0030]

[0031]

[0032] The weighted update solution is:

[0033]

[0034] in The weighting coefficients are solved by minimizing the dispersion error objective function, which is a nonlinear optimization problem in mathematical form.

[0035] Furthermore, the objective function for dispersion error in step three is specifically as follows:

[0036]

[0037] The nonlinear optimization problem is then expressed as:

[0038]

[0039] in, Let be the weight vector to be optimized. For the maximum effective wavenumber, This is the theoretical phase velocity (medium-dependent). Numerical phase velocity;

[0040] For the fourth-order Rungekuta form, the analytical form of the numerical phase velocity is:

[0041]

[0042] in For the phase transfer function:

[0043] .

[0044] Furthermore, the optimal weights of the weight coefficients in step three satisfy a symmetric relationship:

[0045]

[0046] And it can be approximated as:

[0047]

[0048] .

[0049] Furthermore, the weighting coefficients in step three should satisfy the following constraints:

[0050]

[0051] Among them, formula (a) is a normalization constraint that makes the Runge-Kutta scheme have a first-order accuracy consistency condition, ensuring that the time discretization process satisfies the law of mass conservation and avoiding non-physical energy growth or decay; formula (b) is a non-negativity constraint that satisfies the maximum principle, prevents numerical oscillations, and ensures that the solution is physically bounded; formula (c) is an extension of the Runge-Kutta scheme with classical CFL conditions, namely a stability constraint.

[0052] Furthermore, the parameters set in step four include: defining the mesh size as... The number of grids is The sampling step size is Sampling time is The spatial domain employs a 2M-order finite difference, and the size of the two-dimensional model is set to [size missing]. The longitudinal wave velocity is The transverse wave velocity is The density is Longitudinal wave quality factor is The shear wave quality factor is The anisotropy parameters are .

[0053] Furthermore, step five specifically involves: adopting a CPML absorbing boundary, setting the absorbing layer thickness to pmln, and the reflection coefficient to R; simultaneously, considering the addition of the absorbing layer pmln, the size of the two-dimensional model should be expanded to... .

[0054] Furthermore, the expression for the Reich wavelet in step six is ​​defined as follows:

[0055]

[0056] The source coordinates are set as follows: The number of earthquake sources is The distance between the guns is Detector coordinates are set to The number of detectors is Lane spacing is .

[0057] Compared with existing technologies, this invention first establishes the first-order velocity-stress equation for viscoelastic VTI media, and then uses the equation for high-order staggered grid finite difference discretization in the spatial domain. The principal stress components are located at grid nodes, the shear stress components at grid centers, and the velocity components at grid boundaries. A 2M-order optimized differential operator is used to approximate the spatial derivative, resulting in a semi-discrete difference equation. The semi-discrete equation is then converted into a system of ordinary differential equations. Weighting coefficients are introduced, and the optimal weighting coefficients are solved by minimizing the dispersion error objective function. Next, geological model parameters, boundary conditions, and observation system parameters are set. The obtained optimal weights are then incorporated into the fourth-order weighted Runge-Kutta method for time-stepping, calculating the wavefield values ​​at each time step sequentially until all time steps are completed, obtaining the corresponding simulated seismic records. This invention determines the optimal weights based on specific dispersion conditions for fourth-order weighted Runge-Kutta method calculations. Experiments demonstrate that, while maintaining low computational and data storage requirements, it effectively suppresses numerical dispersion and improves simulation accuracy, providing data support for subsequent research on channel wave propagation in viscoelastic VTI media. Attached Figure Description

[0058] Figure 1 This is a schematic diagram of the interlaced grid used in this invention;

[0059] Figure 2 This is a schematic diagram of the parameters set in the model of this invention;

[0060] Where (a) is the P-wave velocity, (b) is the S-wave velocity, (c) is the density, (d) is the P-wave quality factor, (e) is the S-wave quality factor, and (f) and (g) are the anisotropy parameters, respectively. , ;

[0061] Figure 3 This is a simulated wavefield snapshot used in this invention;

[0062] Where (a) is a snapshot of the x-component wave field at t=0.05s, (b) is a snapshot of the z-component wave field at t=0.05s, (c) is a snapshot of the x-component wave field at t=0.10s, (d) is a snapshot of the z-component wave field at t=0.10s, (e) is a snapshot of the x-component wave field at t=0.15s, and (f) is a snapshot of the z-component wave field at t=0.15s;

[0063] Figure 4 These are simulated earthquake records obtained after numerical simulation according to the present invention;

[0064] (a) represents the x-component seismic record, and (b) represents the z-component seismic record. Detailed Implementation

[0065] The present invention will be further described below.

[0066] like Figure 1 As shown, the present invention includes the following steps:

[0067] Step 1: Establish the first-order velocity-stress equation for the viscoelastic VTI medium: The first-order velocity-stress equation for the three-dimensional viscoelastic VTI medium is derived from the Kelvin-Voigt constitutive relation, the differential equation of motion, the constitutive equation, and the geometric equation, as follows:

[0068]

[0069] The first-order velocity-stress equation for a two-dimensional viscoelastic VTI medium considering the source term is:

[0070]

[0071] in, , Let x and z represent the vibration velocities of a point in the x and z directions, respectively. , These represent the principal stresses in the x and z directions at that point, respectively. The shear stress at that point, Represents the elastic coefficient. The viscosity coefficient;

[0072] Due to anisotropic parameters and model quality factor The effect on elastic wave propagation can be expressed as the effect on the elastic coefficient. and viscosity coefficient Due to the influence of [the specific factors], the formulas for calculating the elastic coefficient and viscosity coefficient of a two-dimensional viscoelastic VTI medium are as follows:

[0073]

[0074]

[0075] in, The dominant frequency of the wavelet, i.e., the frequency set during simulation. .

[0076] Step 2: Spatial Domain High-Order Staggered Mesh Finite Difference Discretization: Using a staggered mesh and combining it with the equations from Step 1, the velocity components are spatially discretized. The principal stress components are located at the mesh nodes, the shear stress components at the mesh centers, and the velocity components at the mesh boundaries. The specific form is shown in the attached figure. Figure 1 As shown, a 2M-order optimized quantization operator is used to approximate the spatial derivative, resulting in a semi-discrete difference equation, specifically:

[0077]

[0078] The formula for calculating the difference coefficients of the above staggered mesh is:

[0079]

[0080] Step 3: Time-domain weighted Runge-Kutta integration: Transform the semi-discrete difference equations into a system of ordinary differential equations; use the fourth-order weighted Runge-Kutta method for time advancement, specifically:

[0081] The semi-discrete difference equation in step two is expressed as follows:

[0082]

[0083] in, F represents the discrete spatial derivative term. Indicates the source term;

[0084] The fourth-order weighted Runge-Kutta method is used for time advancement: the time step is defined as... ,initialization Solution of time Then calculate the slope vectors for the four stages. Calculated using the following formula:

[0085]

[0086] The node coefficients use standard values:

[0087]

[0088]

[0089] The weighted update solution is:

[0090]

[0091] in The weighting coefficients are solved by minimizing the dispersion error objective function, which is mathematically a nonlinear optimization problem.

[0092] The specific objective function for dispersion error is as follows:

[0093]

[0094] The nonlinear optimization problem is then expressed as:

[0095]

[0096] in, Let be the weight vector to be optimized. For the maximum effective wavenumber, This is the theoretical phase velocity (medium-dependent). Numerical phase velocity;

[0097] For the fourth-order Rungekuta form, the analytical form of the numerical phase velocity is:

[0098]

[0099] in For the phase transfer function:

[0100] .

[0101] Extensive numerical experiments have demonstrated that the optimal weights of the aforementioned weighting coefficients satisfy a symmetric relationship:

[0102]

[0103] And it can be approximated as:

[0104]

[0105] .

[0106] Furthermore, the aforementioned weighting coefficients should satisfy the following constraints:

[0107]

[0108] Among them, formula (a) is a normalization constraint that makes the Runge-Kutta scheme have a first-order accuracy consistency condition, ensuring that the time discretization process satisfies the law of mass conservation and avoiding non-physical energy growth or decay; formula (b) is a non-negativity constraint that satisfies the maximum principle, prevents numerical oscillations, and ensures that the solution is physically bounded; formula (c) is an extension of the Runge-Kutta scheme with classical CFL conditions, namely a stability constraint.

[0109] Step 4: Model Parameter Settings: Establish the geological model and set the corresponding parameters, including: defining the grid size as... The number of grids is The sampling step size is Sampling time is The spatial domain employs a 2M-order finite difference, and the size of the two-dimensional model is set to [size missing]. The longitudinal wave velocity is The transverse wave velocity is The density is Longitudinal wave quality factor is The shear wave quality factor is The anisotropy parameters are .

[0110] Step 5: Boundary Condition Setting: Set boundary conditions for the geological model established in Step 4 for subsequent trough wave numerical simulation. Specifically, use CPML absorbing boundaries, set the absorbing layer thickness to pmln, and the reflection coefficient to R. Also, considering the addition of the pmln absorbing layer, the size of the two-dimensional model should be expanded to [size missing]. .

[0111] Step Six: Setting Observation System Parameters: The source simulation uses a dominant frequency of [missing value]. Delay time The Reich wavelet, specifically expressed as:

[0112]

[0113] The source coordinates are set as follows: The number of earthquake sources is The distance between the guns is Detector coordinates are set to The number of detectors is Lane spacing is .

[0114] Step 7: Numerical simulation of the trough wave: Based on the parameters set in steps 4 to 6, the fourth-order weighted Runge-Kutta method from step 3 is used to advance the time step and calculate the wave field value at each time step in sequence until all time steps are completed, and the corresponding trough wave value is obtained.

[0115] Effect verification:

[0116] To verify the effectiveness of the groove wave simulation in this invention, a numerical simulation experiment was conducted.

[0117] To simplify calculations and demonstrate the channel wave simulation effect, the geological model was set as a rock-coal stratified model, with a survey line located in the exact center of the coal seam. Therefore, all detectors recorded direct waves. The parameters of the geological model are as follows: Figure 2As shown. The two-dimensional model is set with a mesh size of 1×1 and a mesh number of 300×140, where the coal seam centerline is located at z=nz / 2 and the thickness is 20m; the sampling interval is set to dt=0.0001s and the sampling time is set to T=0.2s; the seismic source is a 100Hz Ricker wavelet located at (30, 70); the absorption layer thickness is set to 20, and the reflection coefficient R=1e-6. Numerical simulation is performed using the method of this invention, and the wavefield snapshots during the simulation are shown below. Figure 3 As shown, the simulated earthquake record is as follows Figure 4 As shown. (Through) Figure 4 Earthquake records show that the geological conditions of the geological model can be accurately obtained, thus demonstrating that the present invention can effectively improve the accuracy and reliability of the numerical simulation of channel waves, and provide data support for subsequent research on the propagation law of channel waves in viscoelastic VTI media.

[0118] 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 method for numerical simulation of a channel wave in a viscoelastic VTI medium based on WRK-SSGFD, characterized in that, The method comprises the following steps: Step one, establishing a first-order velocity-stress equation of a viscoelastic VTI medium: a first-order velocity-stress equation of a three-dimensional viscoelastic VTI medium is derived from a Kelvin-Voigt constitutive relation, a motion differential equation, a constitutive equation and a geometric equation; Step two, spatial domain high-order staggered grid finite difference discretization: velocity components are discretized in space by using a staggered grid and combining the equation in step one, wherein main stress components are located at grid nodes, shear stress components are located at grid centers, and velocity components are located at grid boundaries, and 2M-order optimized difference operator approximation is used for spatial derivatives to obtain a semi-discrete difference equation; Step three, time domain weighted Runge-Kutta integration: the semi-discrete difference equation is converted into the form of a system of ordinary differential equations; weight coefficients are optimized by minimizing a dispersion error objective function, and then the optimal weight is integrated into a fourth-order weighted Runge-Kutta method for time marching; Step four, model parameter setting: a geological model is established and corresponding parameters are set; Step five, boundary condition setting: boundary conditions are set for the geological model established in step two, which is used for subsequent channel wave numerical simulation; Step Six: Setting Observation System Parameters: The source simulation uses a dominant frequency of [missing value]. Delay time The Reichswave; Step seven, channel wave numerical simulation: according to the parameters set in steps four to six, the fourth-order weighted Runge-Kutta method in steps three and four is used for time marching, and the wave field value at each time is calculated in sequence until all time steps are completed, and corresponding simulated seismic records are obtained.

2. The method of numerical simulation of the slot wave in the viscoelastic VTI medium based on the WRK-SSGFD according to claim 1, characterized in that, The first-order velocity-stress equation of the three-dimensional viscoelastic VTI medium in step one is specifically: The first-order velocity-stress equation of a two-dimensional viscoelastic VTI medium considering a source term is: wherein, , respectively denote the vibration velocity of a point x in the x direction and the z direction, , respectively denote the principal stress of the point x in the x direction and the z direction, is the shear stress of the point x, denotes the elastic coefficient, is the viscous coefficient.

3. The method of claim 1, wherein the WRK-SSGFD-based viscoelastic VTI medium's channel wave numerical modeling method is characterized by, The semi-discrete difference equation in step two is specifically:

4. The method of claim 1, wherein the WRK-SSGFD-based viscoelastic VTI medium's channel wave numerical modeling method is characterized by, Step three is specifically: The semi-discrete difference equation in step two is expressed as wherein F represents a discrete spatial derivative term, represents a source term; The fourth-order weighted Runge-Kutta method is used for time advancement: the time step is defined as... ,initialization Solution of time Then calculate the slope vectors for the four stages. Calculated using the following formula: Wherein the node coefficient adopts a standard value: The weighted updated solution is: wherein are weight coefficients, solved by minimizing a dispersion error objective function, which is mathematically formulated as a nonlinear optimization problem.

5. The method of numerical simulation of the slot wave in the viscoelastic VTI medium based on the WRK-SSGFD according to claim 4, characterized in that, The dispersion error objective function in step three is specifically: Then the nonlinear optimization problem is expressed as: wherein, is the weight vector to be optimized, is the maximum effective wave number, is the theoretical phase velocity, is the numerical phase velocity; For the fourth-order Runge-Kutta form, the analytical form of the numerical phase velocity is: wherein is the phase transfer function: 。 6. The method of claim 4, wherein the WRK-SSGFD-based viscoelastic VTI medium's channel wave numerical modeling method is characterized by, The optimal weight of the weight coefficient in step three satisfies a symmetry relationship: And is approximately expressed as: 。 7. The method of claim 6, wherein the WRK-SSGFD-based viscoelastic VTI medium's channel wave numerical modeling method is characterized by, The weight coefficient in step three should satisfy the following constraint conditions: Wherein, formula (a) is a normalization constraint, which makes the Runge-Kutta format have a first-order accuracy consistency condition, ensures that the time discretization process satisfies the law of conservation of mass, avoids non-physical growth or decay of energy; formula (b) is a non-negativity constraint, which satisfies the maximum principle, prevents numerical oscillation, and ensures that the solution is physically bounded; formula (c) is an extension of the classical CFL condition of the Runge-Kutta format, that is, a stability constraint.

8. The method of slot wave numerical simulation of viscoelastic VTI medium based on WRK-SSGFD according to claim 1, characterized in that, The parameter setting in the fourth step includes: defining the grid size as , the grid number as , the sampling step as , and the sampling time as ; the spatial domain adopts 2M-order difference, the two-dimensional model size is set as , the longitudinal wave velocity is , the transverse wave velocity is , the density is , the longitudinal wave quality factor is , the transverse wave quality factor is , and the anisotropy parameter is .

9. The method of slot wave numerical simulation of a viscoelastic VTI medium based on the WRK-SSGFD according to claim 1, characterized in that, The step five is specifically: using CPML absorbing boundary, the thickness of the absorbing layer is set as pmln, and the reflection coefficient is R; meanwhile, considering the addition of the absorbing layer pmln, the size of the two-dimensional model should be expanded to .

10. The method of slot wave numerical simulation of a viscoelastic VTI medium based on the WRK-SSGFD according to claim 1, characterized in that, The Ricker wavelet expression in step six is defined as: Wherein, the source coordinate is set as , the source number is , the shot interval is ; the detector coordinate is set as , the detector number is , and the trace interval is .

Images