A two-dimensional three-component surface wave simulation method and system based on finite difference method

Abstract

The application belongs to the technical field of seismic wave forward simulation, and discloses a two-dimensional three-component surface wave simulation method and system based on a finite difference method. The method is based on the discrete processing of continuous medium by the rotated staggered grid finite difference method, and realizes the free surface condition based on the discrete processing strategy; the M-PML is used to process the other absorbing boundary except the free surface in the solving area; the two-dimensional three-component surface wave simulation is carried out through the source loading, and the x, y, z direction seismic wave field containing Rayleigh wave and Love wave is obtained. The application realizes the high-precision two-dimensional three-component surface wave field forward simulation, and realizes the simulation of body wave and surface wave field at the same time, which makes a contribution to the promotion of two-dimensional full wave field forward simulation research, and is beneficial to the promotion of the research on the law of seismic wave propagation.

Description

TECHNICAL FIELD

[0001] The application belongs to the technical field of seismic wave forward modeling, and particularly relates to a two-dimensional three-component surface wave simulation method and system based on a finite difference method. BACKGROUND

[0002] Surface wave forward modeling is an important field in surface wave exploration, is the theoretical basis of surface wave exploration, and is an extremely important component in the field of surface wave exploration. In the early stage of surface wave forward modeling, the dispersion equation of surface wave was used, and related methods include the Thomson-Haskell method, the Schwab-Knopoff method, the delta matrix method, the generalized reflection-transmission coefficient algorithm, and the fast vector method, which are commonly used in the study of surface wave dispersion characteristics in layered media. Based on these methods, researchers have established and improved the surface wave dispersion theory in layered media. However, surface wave simulation based on the dispersion equation cannot directly reflect the propagation process of surface waves in the medium and cannot reflect the relationship between other waves and surface waves in the wave field. Surface wave simulation based on the wave equation can make up for these deficiencies, and thus has become an important field of surface wave forward modeling research.

[0003] Traditional surface wave exploration methods use single-component acquisition, only collecting vertical z-component seismic information. Due to the limited information carried, noise interference, and multi-mode surface wave interference, the inversion results based on single-component have large errors compared to the true situation. Later, researchers found that Rayleigh waves also have seismic responses in the radial x-component record, and their dispersion characteristics are different from those of vertical seismic records, which can supplement the analysis results of the z-component. Therefore, x-component acquisition was started, and x-component and z-component seismic records were used for analysis and research of near-surface velocity structure. In addition, researchers found that Love waves contain different dispersion characteristics from Rayleigh waves, and began to use joint inversion methods of Rayleigh waves and Love waves to overcome the problem of multiple solutions, thereby obtaining higher precision near-surface velocity structure. The seismic response of Love waves can be observed in the tangential y-component, so three-component surface wave exploration technology has been applied to the fields of oil exploration and natural earthquakes.

[0004] Currently, the simulation of Rayleigh waves in the x-component and z-component and Love waves in the y-component in surface wave field numerical simulation is carried out separately, and P-SV wave equation and SH wave equation simulation are carried out respectively to realize Rayleigh wave and Love wave simulation, which cannot directly realize three-component surface wave simulation. With the development of surface wave exploration methods, researchers have begun to develop joint inversion and full waveform inversion methods, and single Rayleigh wave and Love wave simulation methods can no longer meet the demand. SUMMARY

[0005] To overcome the problems in the related art, the embodiment of the present application provides a two-dimensional three-component surface wave simulation method and system based on a finite difference method.The purpose of the present application is to provide a two-dimensional three-component surface wave simulation method based on a finite difference method to realize synchronous simulation of Rayleigh waves and Love waves in a two-dimensional state, and to provide technical support for multi-component surface wave exploration technology, inversion technology and full-wave field simulation technology.

[0006] The technical solution is as follows: a two-dimensional three-component surface wave simulation method based on a finite difference method, comprising:

[0007] S1, performing discrete processing on a continuous medium based on a rotated staggered grid finite difference method, and realizing a free surface condition based on a discrete processing strategy;

[0008] S2, using M-PML to process an absorbing boundary of a solution area except a free surface;

[0009] S3, performing two-dimensional three-component surface wave simulation through source loading to obtain x, y and z direction seismic wave fields containing Rayleigh waves and Love waves.

[0010] Before step S1, the following needs to be performed: based on two-dimensional three-component first-order velocity-stress equations, P-SV wave equations and SH wave equations are solved to obtain equivalence of the two-dimensional three-component first-order velocity-stress equations, the P-SV wave equations and the SH wave equations.

[0011] Further, the equivalence of the two-dimensional three-component first-order velocity-stress equations, the P-SV wave equations and the SH wave equations comprises:

[0012] According to three basic equations of elastic dynamics, a partial derivative of a displacement component with respect to time is obtained to obtain a velocity component, a body force term is omitted, and a three-dimensional first-order velocity-stress elastic wave equation is obtained; when a wave field no longer changes in the y direction, the three-dimensional three-component first-order velocity-stress elastic wave equation degenerates into a two-dimensional three-component first-order velocity-stress elastic wave equation:

[0013]

[0014] wherein, υ i is a velocity component, wherein i=x, y, z; τ ij is a pressure component, wherein i, j=x, y, z; p is a density; C ij is an elastic parameter; wherein i, j=1, 2...6; when a medium is isotropic, an elastic coefficient matrix composed of elastic parameters is as follows:

[0015]

[0016] wherein, l and m are Lame constants;

[0017] Substitute the elastic parameters of the isotropic medium into the two-dimensional three-component first-order velocity-stress equation to obtain:

[0018]

[0019] The P-SV wave equation is:

[0020]

[0021] The SH wave equation is:

[0022]

[0023] Based on the P-SV wave equation, the wave field simulation of Rayleigh waves in the x and z directions is realized, and based on the SH wave equation, the wave field simulation of Love waves in the y direction is realized.

[0024] In step S1, the continuous medium is discretely processed based on the rotated staggered grid finite difference method, including:

[0025] The diagonal direction difference is designed for the rotated staggered grid The xoz coordinate system is used to calculate the difference in the diagonal direction under the xoz coordinate system, and then the linear combination is used to calculate the difference in the xoz coordinate system under the xoz coordinate system.

[0026] Further, based on the discretization strategy of the rotated staggered grid finite difference method, the free surface is placed at the position of the velocity sampling point, and only the density on the free surface is processed for the stress component defined on the same grid point, without adjusting the Lame constant; the adjustment of the density parameter on the free surface is ρ=0.5ρ0, ρ0 is the density of the medium below the free surface; for the wave field component above the free surface, it is directly assigned to zero; after the wave field is initialized, the wave field is updated from the free surface to meet the implicit requirement of the stress condition on the free surface.

[0027] In step S2, the M-PML is used to process the absorbing boundary of the solving region except the free surface, including:

[0028] The absorbing boundary of the solving region except the free surface includes the identified 1, 2 and 3 different absorbing boundaries, and the different regions identified by the identified 1, 2 and 3 different absorbing boundaries are the inlaid layers.

[0029] The wave field propagating into the inlaid layer is divided into x and z two parts according to the propagation direction, and the attenuation coefficient is introduced to attenuate the part perpendicular to the boundary, and the v x component is:

[0030] ​​

[0031] wherein is v x The component propagating along the x direction and the z direction, d x , d z is the attenuation factor in the x direction and the z direction in the inlaid layer; the calculation formula of the attenuation factor is:

[0032]

[0033] wherein, the x direction is the area identified by 1, the z direction is the area identified by 2, p (z / x) is the attenuation factor proportional coefficient in the left and right inlaid layers in the x direction, p (x / z) is the attenuation factor proportional coefficient in the upper and lower inlaid layers in the z direction; the proportional coefficient is taken as p (z / x) = p (x / z) = 0.5; R is a theoretical reflection coefficient, taken as R=0.000001, v pmax is the maximum velocity of the longitudinal wave in the inlaid layer, L is the thickness of the absorbing layer, x and z are the vertical distances of the demarcation between the solving area and the absorbing boundary; for the area identified by 3, the attenuation coefficient is taken as zero, and the attenuation factors calculated in two directions are used for calculation.

[0034] In step S3, the two-dimensional three-component surface wave simulation is performed by loading the seismic source, and x, y, z direction seismic wave fields containing Rayleigh waves and Love waves are obtained, including:

[0035] The v z and v y components at the position of the seismic source are loaded with the simulation seismic source, then the two-dimensional three-component first-order velocity-stress equation is solved by using the rotated staggered grid finite difference method, and after the solving is completed, the three-component surface wave front snapshot and the synthesized seismic record are obtained;

[0036] The simulation seismic source selects different frequencies and different types of wavelets as the simulation seismic source as needed, uses the Ricker wavelet as the simulation wavelet, and simulates the seismic source excitation.

[0037] Another purpose of the present application is to provide a two-dimensional three-component surface wave simulation system based on the finite difference method, which implements the two-dimensional three-component surface wave simulation method based on the finite difference method, and the system comprises:

[0038] The free surface condition implementation module is used for discretely processing the continuous medium based on the rotated staggered grid finite difference method, and realizing the free surface condition based on the discrete processing strategy;

[0039] The absorbing boundary processing module is used for processing the M-PML for the other absorbing boundaries of the solving area except the free surface;

[0040] A three-component seismic wave field obtaining module is configured to obtain x, y, z direction seismic wave fields containing Rayleigh waves and Love waves through two-dimensional three-component surface wave simulation by a seismic source loading.

[0041] Further, the two-dimensional three-component surface wave simulation system based on the finite difference method is loaded on a computer readable storage medium, the computer readable storage medium stores a computer program, and the computer program can realize the functions of the two-dimensional three-component surface wave simulation system based on the finite difference method when executed by a processor.

[0042] Further, the two-dimensional three-component surface wave simulation system based on the finite difference method is applied in the field of oil exploration and natural seismic exploration

[0043] In combination with all the above technical solutions, the present application has the following beneficial effects: the present application realizes high-precision surface wave simulation based on the rotating staggered grid finite difference method, compared with the traditional surface wave simulation method based on the standard staggered grid finite difference method, the processing on the free surface is simpler and easier to realize, and the simulation precision is higher and the difference format stability is better; the currently commonly used surface wave field simulation method is based on the standard staggered grid finite difference method to simulate P-SV wave equation and SH wave equation respectively, so as to realize Rayleigh wave and Love wave simulation respectively. In the traditional surface wave exploration method, the commonly used is single-component surface wave exploration, so the simulation scheme of single-component surface wave has been sufficient to meet the demand. However, with the development of multi-component surface wave exploration, researchers have begun to develop joint inversion and full waveform inversion methods, and the single Rayleigh wave and Love wave simulation method is relatively cumbersome in practical application, and is not convenient for related research and further development. The two-dimensional three-component surface wave simulation method based on the rotating staggered grid difference method provided by the present application can realize surface wave simulation of three directions containing Rayleigh wave and Love wave in a two-dimensional state, which provides a new means for the research of multi-component surface wave exploration method, the research of surface wave multi-component joint inversion method and the research of surface wave propagation law.

[0044] In the traditional seismic wave simulation, it is generally believed that the elastic wave field is the full wave field. However, in the actual situation, the seismic wave field contains not only the body wave in the elastic wave field, but also the surface wave information. Therefore, the seismic wave field containing only the body wave field cannot be called the full wave field, and many researchers have begun to add surface wave simulation to further realize full wave field simulation. However, the current surface wave simulation is limited to two-dimensional x, z component and single y component simulation. The present application realizes high-precision two-dimensional three-component surface wave field forward simulation, and simultaneously realizes body wave and surface wave field simulation, which contributes to the promotion of two-dimensional full wave field forward simulation research and is beneficial to the promotion of seismic wave propagation law research. BRIEF DESCRIPTION OF DRAWINGS

[0045] The accompanying drawings, which are incorporated in and form part of this specification, illustrate embodiments consistent with this disclosure and, together with the description, serve to explain the principles of this disclosure;

[0046] Figure 1 This is a flowchart of a two-dimensional three-component surface wave simulation method based on the finite difference method provided in an embodiment of the present invention;

[0047] Figure 2 This is a schematic diagram of a rotating staggered mesh provided in an embodiment of the present invention;

[0048] Figure 3 This is a boundary condition design diagram provided in an embodiment of the present invention (where the areas marked by 1, 2, and 3 are absorbing boundaries);

[0049] Figure 4 This is a waveform curve of the x-component of the standard staggered grid finite difference method (SSG) provided in Example 1 of the present invention;

[0050] Figure 5 This is a waveform curve of the z-component of the standard staggered mesh finite difference method (SSG) provided in Example 1 of the present invention.

[0051] Figure 6 This is a waveform curve of the z-component of the Rotating Mesh Finite Difference (RSG) method provided in Example 1 of the present invention.

[0052] Figure 7 This is a waveform curve of the x-component of the Rotating Mesh Finite Difference (RSG) method provided in Example 1 of the present invention.

[0053] Figure 8 This is a snapshot of the x-component wavefield at t = 0.15s provided in Example 2 of this embodiment of the invention;

[0054] Figure 9 This is a snapshot of the y-component wavefield at t = 0.15s provided in Example 2 of this embodiment of the invention;

[0055] Figure 10 This is a snapshot of the z-component wavefield at t = 0.15s provided in Example 2 of this embodiment of the invention;

[0056] Figure 11 This is the x-component synthesized seismic record map provided in Example 2 of the present invention;

[0057] Figure 12 This is the y-component synthesized seismic record map provided in Example 2 of the present invention;

[0058] Figure 13 This is the z-component synthesized seismic record map provided in Example 2 of the present invention;

[0059] Figure 14 This is a comparison chart of single-channel surface wave records in the x-component of the present invention (RSG) method and two commonly used standard staggered mesh methods (SSG-SIM, SSG-AEA) when ppw (points per minimum wavelength, the ratio of the minimum wavelength of the surface wave to the spatial step size) = 10, provided in Example 3 of the present invention.

[0060] Figure 15 This is a comparison chart of single-channel surface wave records in the y component of the present invention (RSG) method and two commonly used standard staggered grid methods (SSG-SIM, SSG-AEA) when ppw (points per minimum wavelength, the ratio of the minimum wavelength of the surface wave to the spatial step size) = 10, provided in Example 3 of the present invention.

[0061] Figure 16 This is a comparison chart of single-channel surface wave records in the z-component of the present invention (RSG) method and two commonly used standard staggered grid methods (SSG-SIM, SSG-AEA) when ppw (points per minimum wavelength, the ratio of the minimum wavelength of the surface wave to the spatial step size) = 10, provided in Example 3 of the present invention.

[0062] Figure 17 This is a comparison diagram of single-channel surface wave records of the RSG method, SSG-SIM, and SSG-AEA methods in the x-component when ppw=25, provided in Example 3 of the present invention.

[0063] Figure 18 This is a comparison diagram of single-channel surface wave records of the RSG method, SSG-SIM, and SSG-AEA methods in the y component when ppw=25, provided in Example 3 of the present invention.

[0064] Figure 19 This is a comparison diagram of single-channel surface wave records of the RSG method, SSG-SIM, and SSG-AEA methods in the z component when ppw=25, provided in Example 3 of the present invention.

[0065] Figure 20 This is a comparison diagram of single-channel surface wave records of the RSG method, SSG-SIM, and SSG-AEA methods in the x-component when ppw=50, provided in Example 3 of this invention.

[0066] Figure 21 This is a comparison diagram of single-channel surface wave records of the RSG method, SSG-SIM, and SSG-AEA methods in the y component when ppw=50, provided in Example 3 of the present invention.

[0067] Figure 22 This is a comparison diagram of single-channel surface wave records of the RSG method, SSG-SIM, and SSG-AEA methods in the z component when ppw=50, provided in Example 3 of this invention.

[0068] Figure 23 This is a comparison chart of the L2 norm error of the waveform curves of different methods in the x component provided in Example 3 of the present invention;

[0069] Figure 24 This is a comparison chart of the L2 norm error of the waveform curves of different methods in the y component provided in Example 3 of the present invention;

[0070] Figure 25 This is a comparison chart of the L2 norm errors of the waveform curves of different methods in the z component provided in Example 3 of the present invention;

[0071] Figure 26 This is a synthetic seismic record map of the x component provided in Example 4 of the present invention;

[0072] Figure 27 This is a synthesized seismic record map of the y-component provided in Example 4 of the present invention;

[0073] Figure 28 This is a synthesized seismic record map in the z-component provided in Example 4 of the present invention;

[0074] Figure 29 This is the dispersion energy diagram of the x-component provided in Example 4 of the present invention (where the black dashed line is the theoretical dispersion energy diagram);

[0075] Figure 30 This is the mid-frequency dispersion energy diagram of the y-component provided in Example 4 of the present invention;

[0076] Figure 31 This is the z-component dispersion energy diagram provided in Example 4 of the present invention. Detailed Implementation

[0077] To make the above-mentioned objects, features, and advantages of the present invention more apparent and understandable, specific embodiments of the present invention will be described in detail below with reference to the accompanying drawings. Many specific details are set forth in the following description to provide a thorough understanding of the present invention. However, the present invention can be practiced in many other ways different from those described herein, and those skilled in the art can make similar modifications without departing from the spirit of the present invention. Therefore, the present invention is not limited to the specific embodiments disclosed below.

[0078] The innovation of this invention lies in the fact that it combines the P-SV wave equation and the SH wave equation using a two-dimensional three-component first-order velocity-stress equation, and for the first time realizes the simulation of a two-dimensional three-component surface wave including Rayleigh waves and Love waves based on the rotating staggered mesh finite method. This provides convenience for the joint study of multi-component surface waves and is conducive to carrying out full-wavefield simulation.

[0079] Example 1: The two-dimensional three-component surface wave simulation method based on the finite difference method provided in this embodiment of the invention includes: proving the equivalence between the two-dimensional three-component first-order velocity-stress equation and the P-SV wave equation and SH wave equation, and solving the P-SV wave equation and SH wave equation simultaneously based on this; discretizing the continuous medium using the rotated staggered grid finite difference method, and realizing the free surface condition based on this discretization strategy; processing other absorbing boundaries in the solution domain other than the free surface using M-PML; designing a source loading scheme, performing two-dimensional three-component surface wave simulation, and obtaining the seismic wave field in the x, y, and z directions containing Rayleigh waves and Love waves.

[0080] Specifically, the following steps are included:

[0081] Step 1: Based on the two-dimensional three-component first-order velocity-stress equation, simultaneously solve the P-SV wave equation and the SH wave equation to obtain the equivalence between the two-dimensional three-component first-order velocity-stress equation and the P-SV wave equation and the SH wave equation.

[0082] Step 2: Discretize the continuous medium based on the rotating staggered mesh finite difference method, and realize the free surface conditions based on the discretization process;

[0083] Step 3: Use M-PML to process the absorbing boundaries in the solution region other than the free surface;

[0084] Step 4: Perform a two-dimensional three-component surface wave simulation by loading the seismic source to obtain the seismic wave field in the x, y, and z directions, which includes Rayleigh waves and Love waves.

[0085] In this embodiment of the invention, the specific content of step 1 is as follows:

[0086] Based on the two-dimensional three-component first-order velocity-stress wave equation combined with the P-SV wave equation and the SH wave equation, since the Rayleigh wave field simulation in the x and z directions is based on the P-SV wave equation and the Love wave field simulation in the y direction is based on the SH wave equation, the two-dimensional three-component first-order velocity-stress equation can be used to simulate surface wave fields in the x, y, and z directions. The specific steps are as follows.

[0087] By simultaneously solving the three fundamental equations of elastic dynamics, and taking the partial derivative of the displacement component with respect to time, we obtain the velocity component. Neglecting the body force term, we can obtain the three-dimensional first-order velocity-stress elastic wave equation. When the wave field no longer changes in the y-direction, the three-dimensional first-order velocity-stress elastic wave equation degenerates into the two-dimensional first-order velocity-stress elastic wave equation.

[0088]

[0089] Among them, υ i Let τ be the velocity component, where i = x, y, z; ij Let i be the pressure component, and j be the density component, where i, j = x, y, z; ρ be the density component, and C be the density component. ij Let i and j be the elastic parameters, where i and j = 1, 2, ..., 6. The first-order velocity-stress elastic wave equation of the two-dimensional ternary components has been applied to the simulation of two-dimensional ternary seismic waves, providing an excellent approximation of the three-dimensional situation. The elastic coefficient matrix is ​​a commonly used method to describe the elastic properties of a medium, accurately depicting these properties. When the medium is isotropic, the elastic coefficient matrix composed of the elastic parameters is:

[0090]

[0091] Where λ and μ are Lamé constants. Substituting the elastic parameters of the isotropic medium into the two-dimensional three-component first-order velocity-stress equation, we obtain:

[0092]

[0093] It is clear that it contains the P-SV wave equation:

[0094]

[0095] And SH wave equation:

[0096]

[0097] When Rayleigh wave field simulations in the x and z directions are achieved based on the P-SV wave equations, Love wave field simulations in the y direction are achieved based on the SH wave equations. The above has already proven from the principle of elastic waves that the two-dimensional three-component first-order velocity-stress equations can be combined with the P-SV and SH wave equations. Therefore, forward modeling of both the P-SV and SH wave equations can be performed simultaneously based on the two-dimensional three-component first-order velocity-stress equations.

[0098] In this embodiment of the invention, the specific content of step 2 is as follows:

[0099] The continuous medium is discretized using the rotating staggered mesh finite difference method, employing two sets of meshes: a full mesh and a half mesh. Velocity components are defined at points on the full mesh, while stress components are defined at points on the central half mesh. Since density is needed to calculate velocity components and elastic parameters are needed to calculate stress components, density is defined at the same location as the velocity components, and elastic parameters are defined at the same location as the stress components. This avoids interpolation of medium parameters and improves the accuracy of the numerical solution.

[0100] The rotating staggered mesh finite difference method designs two sets of coordinate systems, one of which is along the diagonal direction. The first coordinate system is the xoz coordinate system along the coordinate axes. The spatial difference strategy of the rotated staggered mesh finite difference method is: first calculate the diagonal direction... Difference in coordinate system, and then using The differences in the coordinate system are calculated by linear combination of the differences in the xoz coordinate system along the coordinate axis.

[0101] The discretization strategy based on the rotating staggered grid finite difference method places the free surface at the location of the velocity sampling point. Since the stress components are defined on the same grid point, only the density on the free surface needs to be processed, without adjusting other physical properties. The density parameter on the free surface is adjusted to ρ = 0.5ρ0, where ρ0 is the density of the medium below the free surface. The wave field components above the free surface are directly assigned to zero. After wave field initialization, the wave field is updated starting from the free surface, thus implicitly satisfying the stress condition on the free surface.

[0102] In this embodiment of the invention, the specific content of step 3 is as follows:

[0103] The actual subsurface medium is a semi-infinite space medium, but the computer storage space used in forward modeling is limited, and efficiency considerations prevent the simulation of seismic waves in infinite space. In order to absorb boundary reflections caused by artificially truncated boundaries and improve the accuracy of the simulated wavefield in the solution domain, absorbing boundary techniques are used to process the non-free surface boundaries.

[0104] The boundary is handled using the M-PML technique, which is an improvement on the traditional PML technique. M-PML is a further development of the traditional PML method. It attenuates the split wave field in the mosaic layer in both the x and z directions based on the traditional PML.

[0105] Where d x d z p represents the attenuation factors in the x and z directions of the mosaic layer. (z / x) p is the attenuation factor scaling factor in the x-direction, i.e., in the left and right mosaic layers. (x / z) It represents the attenuation factor scaling factor in the z-direction, i.e., in the upper and lower mosaic layers.

[0106] In this embodiment of the invention, the specific content of step 4 is as follows:

[0107] Simulating the seismic source is a crucial component of seismic wave forward modeling. Currently, the Ricker wavelet is commonly used to simulate the source wavelet, thus achieving source simulation. This invention uses the commonly used Ricker wavelet as the simulation wavelet to simulate source excitation, and its mathematical expression is as follows:

[0108]

[0109] Where (x0, z0) is the location of the earthquake source, f p t0 is the dominant frequency of the seismic source, and t0 is the wavelet delay. Let v be the spatial attenuation function, where α is the spatial attenuation factor, representing the rate of attenuation of the seismic wavelet in space. At the source location, v... z and v y The simulated seismic source is loaded onto the components, and then the two-dimensional three-component first-order velocity-stress equation is solved using the rotating staggered mesh finite difference method. After the solution is completed, surface wavefront snapshots and synthetic seismic records of the three components can be obtained.

[0110] Example 2, as another embodiment of the present invention, shows that step S1 in Example 1 is a theoretical proof and derivation, and this step does not need to be performed every time in the specific implementation of the present invention. Figure 1 As shown, the two-dimensional three-component surface wave simulation method based on the finite difference method provided in this embodiment of the invention includes:

[0111] S1, based on the rotating staggered mesh finite difference method, discretizes the continuous medium and realizes the free surface condition based on the discretization strategy:

[0112] A schematic diagram of a rotating staggered mesh is shown below. Figure 2 As shown, the velocity component is defined at the full grid point, the stress component is defined at the central half grid point, and the density is defined at the same location as the velocity component. The elastic parameter is defined at the same location as the stress component. The rotated staggered mesh design uses two coordinate systems, one diagonally... The second is the xoz coordinate system along the coordinate axes, such as... Figure 2 As shown. The spatial difference strategy of the rotating staggered mesh finite difference method is: first calculate the diagonal direction. Difference in coordinate system, and then using The differences in the coordinate system are calculated by linear combination of the differences in the xoz coordinate system. That is, the differences in the coordinate axis directions are calculated using the differences in the diagonal directions.

[0113] Place the free surface Figure 2As shown, the free surface is located at the same location as the velocity sampling point. Since the stress component sampling points are located at the same grid point, sampling is not performed on the free surface, and no special processing is required for the stress components on the free surface. However, the density and velocity components are defined at the same grid point, and the stress components and elastic modulus are defined at the same grid point. Therefore, only the density on the free surface needs to be processed, and no adjustment to other physical properties is required.

[0114] The density parameter on the free surface is adjusted to ρ = 0.5ρ0, where ρ0 is the density of the medium below the free surface. The wave field components above the free surface are directly assigned zero. After wave field initialization, the wave field is updated starting from the free surface, thus implicitly satisfying the stress condition on the free surface.

[0115] S2, M-PML is used to process the absorbing boundaries in the solution domain other than the free surface:

[0116] To achieve surface wave simulation, this invention employs a combination of free surface boundary conditions and absorbing boundary conditions to handle the boundary, as shown in the structure... Figure 3 As shown in the figure, the regions marked by 1, 2, and 3 are different absorption boundaries.

[0117] The free surface is treated in the manner described in step S1. Figure 3 The absorption boundaries marked in sections 1, 2, and 3 are processed using the M-PML technique. Then, an absorption attenuation layer of finite thickness is embedded outside the region. Figure 3 The regions marked 1, 2, and 3 in the diagram constitute the mosaic layer. The wave field propagating into the mosaic layer is divided into two parts, x and z, according to the propagation direction. An attenuation coefficient is introduced to attenuate the portion perpendicular to the boundary, using v x For example, the components:

[0118]

[0119] in For v x The components propagating along the x and z directions, d x d z Let be the attenuation factors in the x and z directions of the mosaic layer. The formula for calculating the attenuation factor is:

[0120]

[0121] Wherein, the x-direction is... Figure 3 The area marked in Figure 1, in the z-direction, is... Figure 3 The area marked in Figure 2, p (z / x) p (x / z)The attenuation factor scaling factors in the x and z directions are respectively determined based on actual experiments. This invention suggests using p based on experimental results. (z / x) =p (x / z) =0.5. R is the theoretical reflection coefficient, which is taken as R = 0.000001 in this paper, v p max Let L be the maximum longitudinal wave velocity within the mosaic layer, L be the thickness of the absorbing layer, and x and z be the perpendicular distances between the solution domain and the absorbing boundary. Figure 2 For the region marked by 3, the attenuation coefficient is set to zero, and the attenuation factor calculated in both directions is used for calculation.

[0122] S3, a two-dimensional three-component surface wave simulation was performed using source loading to obtain the seismic wavefield in the x, y, and z directions, which includes Rayleigh and Love waves:

[0123] At the same time, at the epicenter location v z and v y Simulated seismic sources are loaded onto the components, and different frequencies and types of wavelets can be selected as simulated seismic sources according to experimental needs. This invention uses a Ricker wavelet with a dominant frequency of 40Hz and a delay of 30ms as the seismic source. The two-dimensional three-component first-order velocity-stress equation is solved using the rotated staggered mesh finite difference method. After the solution is completed, surface wavefront snapshots and synthetic seismic records of the three components can be obtained.

[0124] Example 3: The two-dimensional three-component surface wave simulation system based on the finite difference method provided in this embodiment of the invention includes:

[0125] The free surface condition implementation module is used to discretize the continuous medium based on the rotating staggered mesh finite difference method, and to implement the free surface condition based on the discretization strategy.

[0126] The absorbing boundary processing module is used to process other absorbing boundaries in the solution domain, excluding free surfaces, using M-PML.

[0127] The three-component seismic wavefield acquisition module is used to perform two-dimensional three-component surface wave simulation through source loading, and obtain seismic wavefields in the x, y, and z directions containing Rayleigh waves and Love waves.

[0128] In the above embodiments, the descriptions of each embodiment have different focuses. For parts that are not described in detail or recorded in a certain embodiment, please refer to the relevant descriptions of other embodiments.

[0129] The information interaction and execution process between the above-mentioned devices / units are based on the same concept as the method embodiments of the present invention. For details on their specific functions and technical effects, please refer to the method embodiments section, and they will not be repeated here.

[0130] Those skilled in the art will clearly understand that, for the sake of convenience and brevity, the above-described division of functional units and modules is merely an example. In practical applications, the above functions can be assigned to different functional units and modules as needed, that is, the internal structure of the device can be divided into different functional units or modules to complete all or part of the functions described above. The functional units and modules in the embodiments can be integrated into one processing unit, or each unit can exist physically separately, or two or more units can be integrated into one unit. The integrated unit can be implemented in hardware or as a software functional unit. Furthermore, the specific names of the functional units and modules are only for easy differentiation and are not intended to limit the scope of protection of this invention. The specific working process of the units and modules in the above system can be referred to the corresponding process in the foregoing method embodiments.

[0131] This invention also provides a computer device comprising: at least one processor, a memory, and a computer program stored in the memory and executable on the at least one processor, wherein the processor executes the computer program to implement the steps in any of the above method embodiments.

[0132] This invention also provides a computer-readable storage medium storing a computer program that, when executed by a processor, can implement the steps described in the various method embodiments above.

[0133] This invention also provides an information data processing terminal, which, when executed on an electronic device, provides a user input interface to implement the steps described in the above method embodiments. The information data processing terminal is not limited to mobile phones, computers, or switches.

[0134] This invention also provides a server that, when executed on an electronic device, provides a user input interface to implement the steps described in the above method embodiments.

[0135] This invention provides a computer program product that, when run on an electronic device, enables the electronic device to implement the steps described in the various method embodiments above.

[0136] If the integrated unit is implemented as a software functional unit and sold or used as an independent product, it can be stored in a computer-readable storage medium. Based on this understanding, all or part of the processes in the methods of the above embodiments of this application can be implemented by a computer program instructing related hardware. The computer program can be stored in a computer-readable storage medium, and when executed by a processor, it can implement the steps of the various method embodiments described above. The computer program includes computer program code, which can be in the form of source code, object code, executable files, or certain intermediate forms. The computer-readable medium can include at least: any entity or device capable of carrying the computer program code to a photographing device / terminal device, a recording medium, a computer memory, a read-only memory (ROM), a random access memory (RAM), an electrical carrier signal, a telecommunication signal, and a software distribution medium. Examples include USB flash drives, portable hard drives, magnetic disks, or optical disks.

[0137] To further illustrate the effects of the embodiments of the present invention, the following experiments were conducted.

[0138] Example 1:

[0139] To verify that the difference scheme of the rotated staggered mesh finite difference method used in this invention has better stability than the standard staggered mesh finite difference method, this example uses the P-SV wave equation simulation as an example, and designs a homogeneous medium model with a P-wave velocity of 300 m / s, a S-wave velocity of 1200 m / s, and a density of 2000 kg / m³. 3 The simulation parameters are: time step Δt = 0.2 ms, spatial step Δh = 1 m. The value of v at (200 m, 200 m) is... z A Ricker wavelet with a dominant frequency of 40Hz and a delay of 30ms is applied to the component. Numerical simulations of the P-SV wave equation are performed using a second-order time, 16th-order space rotating staggered grid difference scheme and a standard staggered grid difference scheme, respectively. The v at (199m, 240m) is obtained. x Components and v z The seismic record consists of components. The simulated waveforms obtained using the standard staggered grid method are shown below. Figures 4-5 As shown, the simulated waveform obtained according to the finite difference method of rotated staggered mesh is as follows: Figures 6-7 As shown. Figures 4-5 The waveform curve in the figure is distorted due to the instability of the difference format and cannot reflect the true waveform information, while Figure = Figures 6-7 The waveform curves remain stable. This indicates that the rotated staggered mesh finite difference method has better stability than the standard staggered mesh finite difference method.

[0140] Example 2:

[0141] A simple, uniform half-space model is set up with a longitudinal wave velocity of 260 m / s, a transverse wave velocity of 100 m / s, and a density of 1500 kg / m³. 3 A source wavelet is applied to a free surface. The simulation region is 80m × 40m. The calculation time step is Δt = 0.1ms, and the spatial step is Δh = 0.1m. The result is v. x v y and v z A snapshot of the wave field at t = 0.15s, as shown below. Figures 8-10 As shown in the figure, the wave field snapshot clearly shows various waveforms, including Rayleigh waves (R), shear head waves (Head), longitudinal waves (P), as well as SV waves and SH waves. Their morphologies conform to the principles of exploration geophysics, proving that the present invention can achieve forward modeling of the P-SV wave equation and the SH wave equation simultaneously.

[0142] To observe the morphology of seismic records in different components, the receiver line was positioned on a free surface with a minimum offset of 10m, 70 receiver channels, and a channel spacing of 1m. Three-component seismic records were received with a recording interval of 0.1ms and a recording length of 1s, yielding v. x v y and v z Composite seismic records of components, such as Figures 11-13 As shown. Where v x Components and v z The component is a Rayleigh shot gather record, showing low-velocity, high-energy Rayleigh waves and body waves. The body waves are much weaker than Rayleigh waves and their waveforms are not obvious in the seismic record. y The component contains only SH waves and no P-waves. Its morphological characteristics conform to the principles of exploration geophysics, further proving the effectiveness of the invention.

[0143] Example 3:

[0144] To demonstrate that the method of this invention has higher accuracy than traditional surface wave simulation methods, a uniform half-space model was designed with a longitudinal wave velocity of 260 m / s, a transverse wave velocity of 100 m / s, and a density of 1500 kg / m³. 3 Since the accuracy of surface wave simulation is closely related to the accuracy of mesh generation, mesh generation parameters with different accuracies were designed. The parameter ppw (points perminimum wavelength, the ratio of the minimum wavelength of the surface wave to the spatial step size) was introduced to measure the accuracy of mesh generation. Three sets of simulation parameters with different ppw values ​​were designed: ppw = 10, Δt = 0.10 ms, Δh = 0.125 m; ppw = 10, Δt = 0.04 ms, Δh = 0.050 m; ppw = 10, Δt = 0.02 ms, Δh = 0.025 m.

[0145] Simulations were performed using three methods: SSG-SIM, SSG-AEA, and the RSG method proposed in this invention. The focal depth was 1m, and the receiver was located on a free surface, resulting in simulated seismic records. Single-channel records with a simulated shot concentration offset of 30m were extracted and compared with analytical solutions obtained based on the Cagniard-De Hoop technique. A comparison was made when ppw = 10. Figures 14-16 As shown, the comparison when ppw=25 is as follows: Figures 17-19 As shown, the comparison when ppw=50 is as follows: Figures 20-22 As shown in the figure, the comparison clearly shows that the method proposed in this invention has higher accuracy than traditional surface wave simulation methods, and maintains high accuracy even with lower meshing accuracy. To quantitatively compare the simulation accuracy of different methods, L2 norm error is introduced to quantify the error of the numerical solution. The L2 norm error between the simulated waveform curve and the theoretical value is calculated, and the comparison results are shown in the figure. Figures 23-25 As shown in the figure Figures 23-25 As can be seen, the surface wave simulation method proposed in this invention has higher accuracy than the traditional surface wave simulation method.

[0146] Example 4:

[0147] Based on the actual situation of surface wave exploration, a common four-layered velocity-increasing model is established. The first layer has a P-wave velocity of 490 m / s, a S-wave velocity of 200 m / s, and a density of 2000 kg / m³. 3 The thickness is 5m; the longitudinal wave velocity of the second layer is 750m / s, the transverse wave velocity is 300m / s, and the density is 2000kg / m³. 3 The thickness is 6m; the longitudinal wave velocity of the third layer is 980m / s, the transverse wave velocity is 400m / s, and the density is 2000kg / m³. 3 The thickness is 6m; the longitudinal wave velocity of the fourth layer is 1230m / s, the transverse wave velocity is 500m / s, and the density is 2000kg / m³. 3 The thickness is infinite. A seismic source is loaded onto a free surface, with a time step Δt = 0.1 ms and a spatial step Δh = 0.25 m. The receiver line is placed on the free surface, with a minimum offset of 10 m, 100 receiver channels, a channel spacing of 1 m, and a recording duration of 1 s. The simulation yields vx. 、 v y and v z Component-based composite gun gathering records, such as Figures 26-28 As shown. Figures 26-28 Chinese v x and v z The component is a Rayleigh record, v y The component is the Love wave record, v x and v zThe component shot shows well-developed and very clear Rayleigh waves with obvious dispersion, exhibiting a broom-like pattern, consistent with the characteristics of high energy and low velocity. It is also accompanied by significant body wave development, consistent with exploration geophysical principles. y The component shot concentration shows a broom-like Love wave, with clear development of the in-phase axis and obvious dispersion phenomenon.

[0148] To test whether the dispersion characteristics of the synthetic seismic records obtained in this invention match the theory, high-resolution LRT was used to extract the dispersion energy maps of each component shot gather records, and these maps were compared with the dispersion curves calculated using the fast vector transfer algorithm and the generalized reflection-transmission coefficient algorithm, as shown in Figure 29-. Figure 31 As shown in the figure, the energy distribution of the dispersion energy map corresponding to each component matches the theoretical dispersion curve well in both the basic and higher orders, proving that the synthetic seismic record simulated by this invention can guarantee the correctness of the surface wave dispersion characteristics.

[0149] The above description is merely a preferred embodiment of the present invention, but the scope of protection of the present invention is not limited thereto. Any modifications, equivalent substitutions, and improvements made by those skilled in the art within the scope of the technology disclosed in the present invention, and within the spirit and principles of the present invention, should be covered within the scope of protection of the present invention.

Claims

1. A two-dimensional three-component surface wave simulation method based on the finite difference method, characterized in that, The method includes: S1, based on the rotating staggered mesh finite difference method, the continuous medium is discretized, and the free surface condition is realized based on the discretization strategy; S2, M-PML is used to process the absorbing boundaries of the solution region except for the free surface; S3, by performing two-dimensional three-component surface wave simulation through source loading, the seismic wave field in the x, y, and z directions containing Rayleigh waves and Love waves is obtained; Before step S1, the following steps are required: Based on the two-dimensional three-component first-order velocity-stress equation, the P-SV wave equation and the SH wave equation are combined to obtain the equivalence between the two-dimensional three-component first-order velocity-stress equation and the P-SV wave equation and the SH wave equation. Substituting the elastic parameters of the isotropic medium into the two-dimensional three-component first-order velocity-stress equation, we obtain: ； The P-SV wave equation is: ； The SH wave equation is: ； Wavefield simulations of Rayleigh waves in the x and z directions were achieved based on the P-SV wave equations, and wavefield simulations of Love waves in the y direction were achieved based on the SH wave equations; among them, Let Lamé constant be denoted by . In step S1, based on the discretization strategy of the rotated staggered mesh finite difference method, the free surface is placed at the location of the velocity sampling point. For stress components defined at the same mesh point, only the density on the free surface is processed, without adjusting the Lamé constant; the density parameter on the free surface is adjusted to... , The density of the medium below the free surface is given; the wave field components above the free surface are directly assigned to zero; after wave field initialization, the wave field is updated starting from the free surface to satisfy the implicit stress condition requirements on the free surface. In step S2, the absorbing boundaries of the solution region, excluding the free surface, are processed using M-PML, including: The absorbing boundaries of the solution domain, excluding the free surface, include the three different absorbing boundaries marked 1, 2, and 3. The different regions marked by the three different absorbing boundaries are called mosaic layers. The wave field propagating into the mosaic layer is divided into two parts, x and z, according to the propagation direction. An attenuation coefficient is introduced to attenuate the part perpendicular to the boundary. The components are: ； in for along direction and The component of directional propagation, where t is time. In the mosaic layer direction and Directional attenuation factor; Let be the velocity component, where ; The pressure component, where ; For density; the formula for calculating the attenuation factor is: ； In this context, the x-direction represents the region marked by 1, and the z-direction represents the region marked by 2. for The attenuation factor scaling factor in the left and right mosaic layers. for The attenuation factor scaling factor in the upper and lower tiling layers; the scaling factor is taken as... ; Let be the theoretical reflection coefficient, and take . , The maximum velocity of the longitudinal wave within the mosaic layer, Let x be the thickness of the absorption layer, and z be the projections of the minimum distance from a point within the absorption boundary region to the boundary of the solution region in the x and z directions, respectively. For the region marked by 3, the attenuation coefficient is set to zero, and the attenuation factor calculated in both directions is used for calculation.

2. The two-dimensional three-component surface wave simulation method based on the finite difference method according to claim 1, characterized in that, Obtaining the equivalence of the two-dimensional three-component first-order velocity-stress equation with the P-SV wave equation and the SH wave equation includes: By simultaneously solving the three fundamental equations of elastic dynamics, and taking the partial derivative of the displacement component with respect to time, we obtain the velocity component. Neglecting the body force term, we obtain the three-dimensional first-order velocity-stress elastic wave equation. When the wave field no longer changes in the y-direction, the three-dimensional first-order velocity-stress elastic wave equation degenerates into the two-dimensional first-order velocity-stress elastic wave equation. ； in, Let be the velocity component, where ; The pressure component, where ; For density, Let be the elastic parameter, where When the medium is isotropic, the elastic coefficient matrix composed of elastic parameters is: ； in, Let be Lamé's constant.

3. The two-dimensional three-component surface wave simulation method based on the finite difference method according to claim 1, characterized in that, In step S1, the continuous medium is discretized based on the rotated staggered mesh finite difference method, including: For the diagonal direction of the rotational staggered mesh design Coordinate system and coordinate axis directions Coordinate system, calculate the diagonal direction Difference in coordinate system, then utilize Linear combination of differences in a coordinate system is used to calculate the coordinate axis directions. Difference in coordinate system.

4. The two-dimensional three-component surface wave simulation method based on the finite difference method according to claim 1, characterized in that, In step S3, the two-dimensional three-component surface wave simulation performed by source loading to obtain the seismic wavefield in the x, y, and z directions, including Rayleigh waves and Love waves, includes: At the epicenter and The simulated seismic source is loaded onto the component, and then the two-dimensional three-component first-order velocity-stress equation is solved using the rotating staggered mesh finite difference method. After the solution is completed, surface wavefront snapshots and synthetic seismic records of the three components are obtained. The simulated seismic source is selected as a wavelet of different frequencies and types as needed, and the Rick wavelet is used as the simulated wavelet to simulate the excitation of the seismic source.

5. A two-dimensional three-component surface wave simulation system based on the finite difference method, characterized in that, The system implements the two-dimensional three-component surface wave simulation method based on the finite difference method as described in any one of claims 1-4, and the system includes: The free surface condition implementation module is used to discretize the continuous medium based on the rotating staggered mesh finite difference method, and to implement the free surface condition based on the discretization strategy. The absorbing boundary processing module is used to process other absorbing boundaries in the solution domain, excluding free surfaces, using M-PML. The three-component seismic wavefield acquisition module is used to perform two-dimensional three-component surface wave simulation through source loading, and obtain seismic wavefields in the x, y, and z directions containing Rayleigh waves and Love waves.

6. The two-dimensional three-component surface wave simulation system based on the finite difference method according to claim 5, characterized in that, The two-dimensional three-component surface wave simulation system based on the finite difference method is mounted on a computer-readable storage medium, which stores a computer program. When the computer program is executed by a processor, it can realize the function of the two-dimensional three-component surface wave simulation system based on the finite difference method.

7. The two-dimensional three-component surface wave simulation system based on the finite difference method according to claim 5, characterized in that, Application of the two-dimensional three-component surface wave simulation system based on the finite difference method in the fields of petroleum exploration and natural earthquakes.

Images