Two-dimensional SHTE mode seismic wave field time domain numerical simulation method based on high-order difference

Abstract

The invention belongs to the field of geophysical exploration, and relates to a two-dimensional SHTE mode seismic wave field time domain numerical simulation method based on high-order difference, which comprises the following steps of: subdividing a seismic wave calculation region by adopting a staggered grid, and discretizing an eight-order precision difference space of a seismic wave variable based on a time domain finite difference method to obtain a discrete format of a time domain partial derivative term; region division: discretizing the spatial difference of the electromagnetic wave variable to obtain a discrete format of a spatial partial derivative term; substituting the discrete format into a time domain control equation of a two-dimensional SHTE mode wave to respectively obtain a finite difference iteration format of a seismic wave vector and a finite difference iteration format of an electromagnetic wave vector, and loading a seismic source; iterating once by using the finite difference iteration format of the seismic wave vector, and then iterating the finite difference iteration format of the electromagnetic wave vector for multiple times; and repeating until the number of iterations of the seismic waves is reached. According to the invention, two-dimensional SHTE mode seismic wave field time domain numerical simulation based on a high-order difference method can be realized.

Description

Technical Field

[0001] The present invention belongs to the field of geophysical exploration, and specifically is a time-domain numerical simulation method of a two-dimensional SHTE mode seismoelectric wave field based on high-order differences. Background Art

[0002] The seismoelectric effect refers to the phenomenon that in a two-phase porous medium, a double electric layer is formed at the interface between the solid skeleton and the porous fluid, and the propagation of seismic waves in the porous medium causes the movement of charged ions in the porous fluid, thereby inducing the generation of an electromagnetic wave field. The seismoelectric detection method is a geophysical detection method based on the seismoelectric effect. It not only has the advantages of high resolution of traditional seismic exploration methods, but can also distinguish interfaces with electrical differences. Therefore, it can be used to detect complex oil and gas reservoirs that have no obvious difference in mechanical properties from the surrounding formations but obvious differences in electrical properties.

[0003] Research on seismoelectric wave field forward modeling can guide seismoelectric detection methods to achieve more effective seismoelectric detection. At present, the research on seismoelectric wave field forward modeling methods at home and abroad is mainly based on analytical methods and numerical simulation methods. For analytical solutions, Pride and Haartsen (1996) solved the solution of the Pride equation in the whole space and gave the Green's function of solid displacement, solid-flow displacement and electric field generated by point force source and current source. Gao and Hu (2010) derived the full-space Green's function of displacement, electric field and magnetic field generated by seismic double-dipole force source. Haartsen and Pride (1997) developed an algorithm to simulate the seismoelectric wave field generated by point source in layered medium. Compared with analytical methods, numerical simulation methods are more suitable for forward modeling of complex underground models with inhomogeneous medium parameters and irregular interfaces. Haine and Pride (2006) used the time domain finite difference method to simulate the seismoelectric response of inhomogeneous porous media. When simulating the electromagnetic response, the displacement current and conduction current terms in Maxwell's equations were ignored, and the Poisson equation was used to simulate the electromagnetic field. However, the quasi-static approximation method cannot effectively simulate the accompanying electric field generated by shear waves and produces numerical errors. Therefore, it is necessary to perform time-domain numerical simulation based on Maxwell's full equation to ensure the simulation accuracy. Gao et al. (2019) used the frequency domain finite difference method to realize the seismic electrical simulation of the two-dimensional SHTE mode. However, the results calculated in the frequency domain require the time-frequency conversion method to convert the time domain signal to the frequency domain. In order to ensure the simulation accuracy, it is necessary to select a wide frequency band for time-frequency conversion to ensure numerical accuracy, which results in high computational cost and low efficiency. Using high-precision finite differences in the direct time domain to simulate the seismoelectric wave field can efficiently and accurately obtain the seismoelectric response of complex underground porous media.

[0004] The Chinese patent publication number CN116184490A discloses a forward modeling method for the seismoelectric wave field of VTI porous media. The seismoelectric wave field of the VTI porous media horizontal layering model is solved using the global matrix method.

[0005] Chinese Patent Publication No. CN117111174A discloses an interface detection method for seismic-electric logging. When solving the problem, the electric field is regarded as a quasi-static field, thereby calculating the electric field induced by the acoustic wave.

[0006] In summary, the current analytical methods mainly calculate the full space and layered models, and numerical simulation methods are required for models of inhomogeneous media and irregular interfaces. The frequency domain seismoelectric wave field simulation method based on Maxwell's full equation requires a large computational cost to ensure the accuracy of time-frequency conversion calculations, and the use of quasi-static approximation methods in numerical simulations in the time domain will cause certain numerical errors. Therefore, the study of two-dimensional SHTE mode seismoelectric wave field simulation based on the full equation can overcome the simulation error caused by the quasi-static approximation process. At the same time, the high-order difference method can reduce the error in numerical approximation, realize high-precision simulation of the seismic electrical response of irregular anomalies in underground porous media, and guide instrument design and observation. Summary of the invention

[0007] The technical problem to be solved by the present invention is to provide a time-domain numerical simulation method of a two-dimensional SHTE mode seismoelectric wave field based on high-order differences. Aiming at the problem that the quasi-static approximation method may cause numerical simulation errors, in order to realize the simulation of seismic electrical response based on the Maxwell full equation, the two-dimensional SHTE mode seismoelectric wave field propagation equation is obtained as the control equation by ignoring the feedback term of electromagnetic waves to seismic waves. In order to improve the simulation accuracy of seismic electrical response, based on the spatial high-order difference approximation method, the spatial high-order difference coefficients are solved, and the differential form of the two-dimensional SHTE mode seismoelectric wave field propagation equation is given. The boundary conditions of seismic waves and electromagnetic waves are set respectively, and the seismic wave and electromagnetic wave control equations are alternately iterated after the earthquake source is loaded, and finally the high-precision simulation of the two-dimensional SHTE mode seismoelectric wave field is realized.

[0008] The present invention is achieved in this way.

[0009] A time-domain numerical simulation method of two-dimensional SHTE mode seismoelectric wave field based on high-order differences, the method includes:

[0010] S1. By neglecting the feedback term L(ω)E of the electromagnetic wave field to the seismic wave field in the Pride equation of the two-dimensional SHTE mode, the time domain control equation of the two-dimensional SHTE mode wave is derived;

[0011] S2. The seismic wave calculation area is divided by staggered grids, and the eighth-order precision differential spatial discretization of the seismic wave variables based on the time-domain finite difference method is used to obtain the discrete format of the time-domain partial derivatives; the electromagnetic wave calculation area is divided by Yee grids, and the spatial differential discretization of the electromagnetic wave variables based on the spatial high-order precision differential method is used to obtain the discrete format of the spatial partial derivatives;

[0012] S3, substituting the discrete format into the time domain control equation of the two-dimensional SHTE mode wave to obtain the finite difference iteration format of the seismic wave vector and the finite difference iteration format of the electromagnetic wave vector, setting the PML boundary for the seismic wave area and the electromagnetic wave simulation area, and loading the seismic source;

[0013] S4, iterating once using the finite difference iterative format of the seismic wave vector and then iterating multiple times using the finite difference iterative format of the electromagnetic wave vector;

[0014] S5. Repeat S4 until the number of seismic wave iterations ends, and display the calculation results.

[0015] Furthermore, in S1, the time domain control equation of the two-dimensional SHTE mode wave is expressed as:

[0016]

[0017] In formula (1), F represents the loaded source, E y represents the y-component electric field strength, H x ,H z They represent the magnetic field strength of the x component and the z component respectively, L is the seismic-electric coupling coefficient, σ(t) is the conductivity of the underground medium, ε is the dielectric constant of the medium, μ is the magnetic permeability of the medium, and v y represents the solid particle velocity in the y direction, q y represents the relative particle velocity of the solid flow in the y direction, τ xy With τ zy is the component of the stress tensor in the rectangular coordinate system, ρ is the equivalent density of the porous medium, ρ f is the pore fluid density, ρ m is the additional mass density, η is the viscosity coefficient of the pore fluid, κ 0 is the permeability of the porous medium, and G represents the shear modulus.

[0018] Furthermore, in S2, the eighth-order precision differential space discretization form of the seismic wave variable is:

[0019]

[0020] v yis the velocity component of the solid phase particles of the seismic wave, i and j represent the grid numbers in the x and z directions respectively, n represents the number of the selected Taylor expansion node, Δx is the side length of the unit grid in the x direction, is the differential coefficient of the spatial derivative with eighth-order accuracy, and the solution is as follows:

[0021]

[0022] The spatial difference discrete form of electromagnetic wave variables is:

[0023]

[0024] Further: The finite difference iteration format of the seismic wave velocity vector in S3 is:

[0025]

[0026] Load the source into F in equation (5) y Term, k represents the number of seismic wave time steps, seismic wave velocity vector Defined at the half-time node, the vector consisting of stress components Defined at an integer number of time nodes, D x With D z They represent the high-order difference operators in the x-direction and the z-direction respectively. The specific form is shown in formula (2). Δt represents the time step of the seismic wave.

[0027] The iterative form of the seismic wave stress component vector is:

[0028]

[0029] The finite difference iterative format for the electromagnetic wave vector is shown in equation (8):

[0030]

[0031] Δ tEM Represents one time step of an electromagnetic wave.

[0032] Furthermore, in S4, N EM is the number of electromagnetic iteration updates, and the number of updates is calculated as follows:

[0033]

[0034] In formula (9), Δt is a time step of the seismic wave. After one seismic wave iteration through formula (5) and formula (7), formula (8) is used to perform N EM Electromagnetic wave iterations.

[0035] Compared with the prior art, the present invention has the following beneficial effects:

[0036] The present invention can overcome the shear wave associated electric field simulation error caused by the use of quasi-static approximation in time domain seismoelectric wave field simulation in current research methods, and realizes the time domain numerical simulation of two-dimensional SHTE mode seismoelectric wave field based on high-order difference method. BRIEF DESCRIPTION OF THE DRAWINGS

[0037] Figure 1 It is a schematic diagram of a time-domain numerical simulation method of a two-dimensional SHTE mode seismoelectric wave field provided by an embodiment of the present invention;

[0038] Figure 2 It is a comparison between the numerical simulation results of the two-dimensional SHTE mode under the porous medium layered model provided by the embodiment of the present invention and the analytical solution;

[0039] Figure 3 It is a wave field snapshot of the solid phase particle velocity under the porous medium layered model provided by the embodiment of the present invention;

[0040] Figure 4 It is a wavelength snapshot of the electric field component Ey under the porous medium layered model provided by an embodiment of the present invention. DETAILED DESCRIPTION

[0041] In order to make the purpose, technical solution and advantages of the present invention more clearly understood, the present invention is further described in detail below in conjunction with the embodiments. It should be understood that the specific embodiments described herein are only used to explain the present invention and are not used to limit the present invention.

[0042] A time-domain numerical simulation method of two-dimensional SHTE mode seismoelectric wave field based on high-order differences includes

[0043] S1. By neglecting the feedback term L(ω)E of the electromagnetic wave field to the seismic wave field in the Pride equation of the two-dimensional SHTE mode, the time domain control equation of the two-dimensional SHTE mode wave is derived;

[0044] S2. The seismic wave calculation area is divided by staggered grids, and the eighth-order precision differential spatial discretization of the seismic wave variables based on the time-domain finite difference method is used to obtain the discrete format of the time-domain partial derivatives; the electromagnetic wave calculation area is divided by Yee grids, and the spatial differential discretization of the electromagnetic wave variables based on the spatial high-order precision differential method is used to obtain the discrete format of the spatial partial derivatives;

[0045] S3, substituting the discrete format into the time domain control equation of the two-dimensional SHTE mode wave to obtain the finite difference iteration format of the seismic wave vector and the finite difference iteration format of the electromagnetic wave vector, setting the PML boundary for the seismic wave area and the electromagnetic wave simulation area, and loading the seismic source;

[0046] S4, iterating once using the finite difference iterative format of the seismic wave vector and then iterating multiple times using the finite difference iterative format of the electromagnetic wave vector;

[0047] S5. Repeat S4 until the number of seismic wave iterations ends, and display the calculation results.

[0048] In S1, the time domain control equation of the two-dimensional SHTE mode wave is expressed as:

[0049]

[0050] In formula (1), F represents the loaded source, E y represents the y-component electric field strength, H x ,H z They represent the magnetic field strength of the x component and the z component respectively, L is the seismic-electric coupling coefficient, σ(t) is the conductivity of the underground medium, ε is the dielectric constant of the medium, μ is the magnetic permeability of the medium, and v y represents the solid particle velocity in the y direction, q y represents the relative particle velocity of the solid flow in the y direction, τ xy With τ zy is the component of the stress tensor in the rectangular coordinate system, ρ is the equivalent density of the porous medium, ρ f is the pore fluid density, is the additional mass density, η is the viscosity coefficient of the pore fluid, κ 0 is the permeability of the porous medium, and G represents the shear modulus.

[0051] In S2, the eighth-order precision difference space discretization form of the seismic wave variable is:

[0052]

[0053] v y is the velocity component of the solid phase particles of the seismic wave, i and j represent the grid numbers in the x and z directions respectively, n represents the number of the Taylor expansion node, Δx is the side length of the unit grid in the x direction, is the differential coefficient of the spatial derivative with eighth-order accuracy, and the solution is as follows:

[0054]

[0055] The spatial difference discrete form of electromagnetic wave variables is:

[0056]

[0057] The finite difference iteration format of the seismic wave velocity vector in S3 is:

[0058]

[0059] Load the source into F in equation (5) y Term, k represents the number of seismic wave time steps, seismic wave velocity vector Defined at the half-time node, the vector consisting of stress components Defined at an integer number of time nodes, D x With D z They represent the high-order difference operators in the x-direction and the z-direction respectively. The specific form is shown in formula (2). Δt represents the time step of the seismic wave.

[0060] The iterative form of the seismic wave stress component vector is:

[0061]

[0062] The finite difference iterative format for the electromagnetic wave vector is shown in equation (8):

[0063]

[0064] Δ tEM Represents one time step of an electromagnetic wave.

[0065] S4, N EM is the number of electromagnetic iteration updates, and the number of updates is calculated as follows:

[0066]

[0067] In formula (9), Δt is a time step of the seismic wave. After one seismic wave iteration through formula (5) and formula (7), formula (8) is used to perform N EM Electromagnetic wave iterations.

[0068] Example

[0069] See also Figure 1 , a time-domain numerical simulation method of a two-dimensional SHTE mode seismoelectric wave field based on high-order difference of the present invention is used for actual simulation, including:

[0070] 1) Set the calculation area to x: -3km~3km, z=-3km~3km, where the x-axis is horizontal and the z-axis is vertical. The nodes in the calculation area are evenly distributed, the node spacing is 10m, and the total number of nodes is 90601; the four boundaries of the calculation area adopt PML boundary conditions, and the artificial source is set at (0m, 0m);

[0071] 2) Set the porous medium parameters in the calculation area and set the layered interface at z = -200m. The seismic-electric coupling coefficient of the upper porous medium is 17.25*10 -10sC / kg, the conductivity of the underground medium is 0.001S / m, and the dielectric constant of the medium is 7.8*10 -11 , the magnetic permeability of the medium is 4π*10 -7 , the equivalent density of porous media is 2647kg / m 3 , the pore fluid density is 1000kg / m 3 , porosity is 0.15, additional mass density is 300kg / m 3 , the viscosity coefficient of the pore fluid is 0.001Pa.s, and the permeability of the porous medium is 0.1*10 -12 m 2 , the shear modulus is 9.6 GPa. The conductivity of the lower porous medium is set to 0.01 S / m, the fluid viscosity is set to 0.1 Pa.s, and the other parameters are the same as the upper layer.

[0072] 3) Set the differential accuracy to 8th order and use formula (3) to solve the spatial high-order differential coefficients;

[0073] 4) The loading source is Ricker wavelet, and the source form is shown in formula (10);

[0074]

[0075] f m is the main frequency of the Ricker wavelet.

[0076] 5) Set the seismic wave simulation time step to 5*10 -4 s, use the finite difference iteration format of the seismic wave vector to update the seismic wave variable once;

[0077] 6) Set the electromagnetic wave simulation time step to 1*10 -8 s, use the finite difference iteration format of the electromagnetic wave vector to update the electromagnetic wave variable once;

[0078] 7) Determine whether the preset electromagnetic time node 5*10 is reached -4 s, if not completed, repeat step 6), if completed, go to step 8);

[0079] 8) Determine whether the preset earthquake time node 1s is reached. If not, repeat step 6). If completed, output the time domain response results of the seismic wave and electromagnetic wave variables;

[0080] See also Figure 2 The comparison between the numerical simulation results of the two-dimensional SHTE mode under the porous medium layered model and the analytical solution shows that the numerical simulation results can effectively simulate the co-seismic electric field of the source electromagnetic wave, interface electromagnetic wave and SH wave;

[0081] Figure 3It is a snapshot of the wave field of solid phase particle velocity in the porous medium layered model, which shows the reflection and transmission phenomenon of seismic wave particle velocity at the interface;

[0082] Figure 4 It is a wavelength snapshot of the electric field component Ey under the porous medium layered model, from which the co-seismic electric field induced by the reflected and transmitted waves of seismic waves can be seen.

[0083] The above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. Any modifications, equivalent substitutions and improvements made within the spirit and principles of the present invention should be included in the protection scope of the present invention.

Claims

1. A time-domain numerical simulation method of two-dimensional SHTE mode seismoelectric wave field based on high-order differences, characterized in that: The method includes: S1. By neglecting the feedback term L(ω)E of the electromagnetic wave field to the seismic wave field in the Pride equation of the two-dimensional SHTE mode, the time domain control equation of the two-dimensional SHTE mode wave is derived; S2. The seismic wave calculation area is divided by staggered grids, and the eighth-order precision differential spatial discretization of the seismic wave variables based on the time-domain finite difference method is used to obtain the discrete format of the time-domain partial derivatives; the electromagnetic wave calculation area is divided by Yee grids, and the spatial differential discretization of the electromagnetic wave variables based on the spatial high-order precision differential method is used to obtain the discrete format of the spatial partial derivatives; S3, substituting the discrete format into the time domain control equation of the two-dimensional SHTE mode wave to obtain the finite difference iteration format of the seismic wave vector and the finite difference iteration format of the electromagnetic wave vector, setting the PML boundary for the seismic wave area and the electromagnetic wave simulation area, and loading the seismic source; S4, iterating once using the finite difference iterative format of the seismic wave vector and then iterating multiple times using the finite difference iterative format of the electromagnetic wave vector; S5. Repeat S4 until the number of seismic wave iterations ends, and display the calculation results.

2. The method for time-domain numerical simulation of seismoelectric wave field in two-dimensional SHTE mode based on high-order difference according to claim 1 is characterized in that: In S1, the time domain control equation of the two-dimensional SHTE mode wave is expressed as: In formula (1), F represents the loaded source, E y represents the y-component electric field strength, H x ,H z They represent the magnetic field strength of the x component and the z component respectively, L is the seismic-electric coupling coefficient, σ(t) is the conductivity of the underground medium, ε is the dielectric constant of the medium, μ is the magnetic permeability of the medium, and v y represents the solid particle velocity in the y direction, q y represents the relative particle velocity of the solid flow in the y direction, τ xy With τ zy is the component of the stress tensor in the rectangular coordinate system, ρ is the equivalent density of the porous medium, ρ f is the pore fluid density, ρ m is the added mass density, η is the viscosity coefficient of the pore fluid, κ0 is the permeability of the porous medium, and G represents the shear modulus.

3. The method for time-domain numerical simulation of two-dimensional SHTE mode seismoelectric wave field based on high-order difference according to claim 2 is characterized in that: In S2, the eighth-order precision difference space discretization form of the seismic wave variable is: v y is the velocity component of the solid phase particles of the seismic wave, i and j represent the grid numbers in the x and z directions respectively, n represents the number of the selected Taylor expansion node, Δx is the side length of the unit grid in the x direction, is the differential coefficient of the spatial derivative with eighth-order accuracy, and the solution is as follows: The spatial difference discrete form of electromagnetic wave variables is:

4. The method for time-domain numerical simulation of seismoelectric wave fields in two-dimensional SHTE mode based on high-order differences according to claim 3 is characterized by: The finite difference iteration format of the seismic wave velocity vector in S3 is: Load the source into F in equation (5) y Term, k represents the number of seismic wave time steps, seismic wave velocity vector Defined at the half-time node, the vector consisting of stress components Defined at an integer number of time nodes, D x With D z They represent the high-order difference operators in the x-direction and the z-direction respectively. The specific form is shown in formula (2). Δt represents the time step of the seismic wave. The iterative form of the seismic wave stress component vector is: The finite difference iterative format for the electromagnetic wave vector is shown in equation (8): Δt EM Represents one time step of an electromagnetic wave.

5. The method for time-domain numerical simulation of seismoelectric wave fields in two-dimensional SHTE mode based on high-order differences according to claim 4 is characterized in that: S4, N EM is the number of electromagnetic iteration updates, and the number of updates is calculated as follows: In formula (9), Δt is a time step of the seismic wave. After one seismic wave iteration through formula (5) and formula (7), formula (8) is used to perform N EM Electromagnetic wave iterations.

Images