Two-dimensional numerical simulation method of seismic electric wave field based on high-order difference for SHTE mode

Abstract

This invention belongs to the field of geophysical exploration and presents a time-domain numerical simulation method for the seismic and electromagnetic wave fields of a two-dimensional SHTE model based on higher-order finite-difference methods. It employs staggered grids to partition the seismic wave computation region, and uses the finite-difference time-domain method to obtain the discrete format of the time-domain partial derivatives by spatially discretizing the seismic wave variables with eighth-order precision. After further region partitioning, the spatial partial derivatives of the electromagnetic wave variables are discretized using spatial difference methods to obtain the discrete format of the spatial partial derivatives. These discrete formats are then substituted into the time-domain governing equations of the two-dimensional SHTE model waves to obtain the finite-difference iterative formats for the seismic wave vector and the electromagnetic wave vector, respectively. The source is then loaded. The method iterates once using the finite-difference iterative format for the seismic wave vector, followed by multiple iterations of the finite-difference iterative format for the electromagnetic wave vector. This process is repeated until the required number of seismic wave iterations is reached. This invention enables time-domain numerical simulation of the seismic and electromagnetic wave fields of a two-dimensional SHTE model based on higher-order finite-difference methods.

Description

Technical Field

[0001] This invention belongs to the field of geophysical exploration, specifically a time-domain numerical simulation method for the seismic wave field of the two-dimensional SHTE mode based on high-order difference. Background Technology

[0002] The seismoelectric effect refers to the phenomenon in two-phase porous media where an electric double layer forms at the interface between the solid framework and the pore fluid. Seismic waves propagating through this porous medium induce the movement of charged ions in the pore fluid, thereby generating an electromagnetic wave field. Seismoelectric detection is a geophysical exploration method based on this effect. It not only possesses the high resolution advantage of traditional seismic exploration methods but also allows for the resolution of interfaces with electrical differences. Therefore, it can be used to detect complex oil and gas reservoirs that have indistinct mechanical properties but significant electrical differences compared to the surrounding strata.

[0003] Research on forward modeling of seismic wave fields can guide seismic detection methods, leading to more effective seismic detection. Currently, domestic and international research on forward modeling methods for seismic wave fields mainly focuses on analytical methods and numerical simulation methods. For analytical solutions, Pride and Haartsen (1996) solved the Pride equations in full space, providing Green's functions for solid displacements, fluid-solid displacements, and electric fields generated by point force and current sources. Gao and Hu (2010) derived the full-space Green's functions for displacements, electric fields, and magnetic fields generated by seismic dual-couple force sources. Haartsen and Pride (1997) developed an algorithm to simulate the seismic wave field generated by a point source in a layered medium. Compared to analytical methods, numerical simulation methods are more suitable for forward modeling complex underground models with inhomogeneous medium parameters and irregular interfaces. Haine and Pride (2006) used the finite-difference time-domain method to simulate the seismic response of non-homogeneous porous media, neglecting the displacement current and conduction current terms in Maxwell's equations when simulating the electromagnetic response, and using Poisson's equations for electromagnetic field simulation. However, the quasi-static approximation method cannot effectively simulate the accompanying electric field generated by shear waves, resulting in numerical errors. Therefore, time-domain numerical simulation based on Maxwell's full equations is required to ensure simulation accuracy. Gao et al. (2019) used the frequency-domain finite-difference method to realize the seismic electrical simulation of the two-dimensional SHTE model. However, the results calculated in the frequency domain require a time-frequency conversion method to transform the time-domain signal into the frequency domain. To ensure simulation accuracy, a wide frequency band needs to be selected for time-frequency conversion, which results in high computational cost and low efficiency. Using high-precision finite-difference in the direct time domain for seismoelectric wavefield simulation can efficiently and accurately obtain the seismoelectric response of complex underground porous media.

[0004] Chinese Patent Publication No. CN116184490A discloses a forward modeling method for the seismoelectric wave field of VTI porous media. The seismoelectric wave field of a horizontally layered model of VTI porous media is solved using the global matrix method.

[0005] Chinese Patent Publication No. CN117111174A discloses an interface detection method for seismic logging, which treats the electric field as a quasi-static field during the solution process, thereby calculating the electric field induced by acoustic waves.

[0006] In summary, current analytical methods primarily calculate full-space and layered models. For models with non-homogeneous media and irregular interfaces, numerical simulation methods are required. However, frequency-domain seismoelectric wavefield simulation based on Maxwell's full equations requires significant computational resources to ensure accurate time-frequency conversion. Furthermore, the use of quasi-static approximations in time-domain numerical simulations introduces numerical errors. Therefore, researching two-dimensional SHTE model seismoelectric wavefield simulation based on the full equations can overcome simulation errors caused by quasi-static approximations. Simultaneously, higher-order difference methods can reduce errors in numerical approximations, enabling high-precision simulation of the seismic electrical response of irregular anomalies in subsurface porous media, thus guiding instrument design and observation. Summary of the Invention

[0007] The technical problem this invention aims to solve is to provide a time-domain numerical simulation method for the two-dimensional SHTE mode seismic-electric wavefield based on higher-order differences. Addressing the numerical simulation errors caused by quasi-static approximation methods, and to achieve simulation of the seismic electrical response based on Maxwell's complete equations, the propagation equations of the two-dimensional SHTE mode seismic-electric wavefield are obtained as governing equations by neglecting the feedback term of electromagnetic waves to seismic waves. To improve the accuracy of the seismic electrical response simulation, the spatial higher-order difference coefficients are solved based on the spatial higher-order difference approximation method, giving the difference form of the propagation equations of the two-dimensional SHTE mode seismic-electric wavefield. Boundary conditions for seismic waves and electromagnetic waves are set separately, and after loading the seismic source, the governing equations for seismic waves and electromagnetic waves are iterated alternately, ultimately achieving high-precision simulation of the two-dimensional SHTE mode seismic-electric wavefield.

[0008] This invention is implemented as follows:

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

[0010] S1. By ignoring 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 computation region is partitioned using staggered grids, and the spatial partial derivative terms of the seismic wave variables are obtained by spatial discretization of the eighth-order precision difference based on the finite-difference time-domain method. The electromagnetic wave computation region is partitioned using Yee grids, and the spatial partial derivative terms of the electromagnetic wave variables are obtained by spatial discretization of the electromagnetic wave variables based on the high-order precision spatial difference method.

[0012] S3. Substitute the discrete scheme into the time-domain control equation of the two-dimensional SHTE mode wave to obtain the finite difference iterative scheme of the seismic wave vector and the finite difference iterative scheme of the electromagnetic wave vector respectively. Set the PML boundary for the seismic wave region and the electromagnetic wave simulation region, and load the source.

[0013] S4. After iterating once using the finite difference iterative scheme of the seismic wave vector, perform multiple iterations using the finite difference iterative scheme of the electromagnetic wave vector.

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

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

[0016]

[0017] In equation (1), F represents the loaded seismic source, and E y H represents the electric field intensity of the y-component. x H z Let x and z represent the magnetic field strengths of the x-component and z-component, respectively; L be the seismoelectric coupling coefficient; σ(t) be the electrical conductivity of the subsurface medium; ε be the permittivity of the medium; μ be the magnetic permeability of the medium; and v be the magnetic flux density. y q represents the solid particle velocity in the y-direction. y τ represents the solid-fluid relative particle velocity in the y-direction. xy With τ zy Let ρ be the component of the stress tensor in a rectangular coordinate system, and ρ be the equivalent density of the porous medium. f ρ is the pore fluid density. m Here, η represents the added mass density, κ0 represents the viscosity coefficient of the pore fluid, and G represents the shear modulus.

[0018] Furthermore, in S2, the eighth-order precision difference space discretization form of the seismic wave variables is as follows:

[0019]

[0020] v y Let represent the solid-phase particle velocity components of the seismic wave, i and j represent the grid indices in the x and z directions respectively, n represent the selected Taylor expansion node number, and Δx represent the side length of a unit grid in the x direction. The difference coefficients of the spatial derivative with eighth-order precision are calculated as follows:

[0021]

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

[0023]

[0024] Furthermore: the finite difference iterative scheme for the seismic wave velocity vector in S3 is:

[0025]

[0026] Load the seismic source into F in equation (5) y The term, k represents the seismic wave time step, and the seismic wave velocity vector. Defined at half a time node, the vector composed of stress components Defined at an integer number of time points, D x With D z These represent the higher-order difference operators in the x and z directions, respectively, and their specific forms are shown in Equation (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 scheme for the electromagnetic wave vector is shown in equation (8):

[0030]

[0031] Δ tEM It represents a time step of an electromagnetic wave.

[0032] Furthermore, in S4, N EM The number of electromagnetic iterations is given by equation (9).

[0033]

[0034] In equation (9), Δt is a time step of the seismic wave. After one seismic wave iteration using equations (5) and (7), equation (8) is used to perform N. EM Subsequent electromagnetic wave iteration.

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

[0036] This invention overcomes the simulation error of the shear wave-accompanying electric field caused by the use of quasi-static approximation in the time-domain seismoelectric field simulation of current research methods, and realizes the time-domain numerical simulation of the two-dimensional SHTE mode seismoelectric field based on the high-order difference method. Attached Figure Description

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

[0038] Figure 2 This is a comparison between the numerical simulation results and analytical solutions of the two-dimensional SHTE mode under the layered model of porous media provided in the embodiments of the present invention;

[0039] Figure 3 This is a snapshot of the wave field of solid-phase particle velocities under the porous medium layered model provided in the embodiments of the present invention.

[0040] Figure 4 This is a wavelength snapshot of the electric field component Ey under the porous medium layered model provided in the embodiments of the present invention. Detailed Implementation

[0041] To make the objectives, technical solutions, and advantages of this invention clearer, the invention will be further described in detail below with reference to embodiments. It should be understood that the specific embodiments described herein are merely illustrative and not intended to limit the invention.

[0042] A two-dimensional SHTE mode seismic electromagnetic wave field time-domain numerical simulation method based on high-order difference, the method includes...

[0043] S1. By ignoring 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 computation region is partitioned using staggered grids, and the spatial partial derivative terms of the seismic wave variables are obtained by spatial discretization of the eighth-order precision difference based on the finite-difference time-domain method. The electromagnetic wave computation region is partitioned using Yee grids, and the spatial partial derivative terms of the electromagnetic wave variables are obtained by spatial discretization of the electromagnetic wave variables based on the high-order precision spatial difference method.

[0045] S3. Substitute the discrete scheme into the time-domain control equation of the two-dimensional SHTE mode wave to obtain the finite difference iterative scheme of the seismic wave vector and the finite difference iterative scheme of the electromagnetic wave vector respectively. Set the PML boundary for the seismic wave region and the electromagnetic wave simulation region, and load the source.

[0046] S4. After iterating once using the finite difference iterative scheme of the seismic wave vector, perform multiple iterations using the finite difference iterative scheme of the electromagnetic wave vector.

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

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

[0049]

[0050] In equation (1), F represents the loaded seismic source, and E y H represents the electric field intensity of the y-component. x H z Let x and z represent the magnetic field strengths of the x-component and z-component, respectively; L be the seismoelectric coupling coefficient; σ(t) be the electrical conductivity of the subsurface medium; ε be the permittivity of the medium; μ be the magnetic permeability of the medium; and v be the magnetic flux density. y q represents the solid particle velocity in the y-direction. y τ represents the solid-fluid relative particle velocity in the y-direction. xy With τ zy Let ρ be the component of the stress tensor in a rectangular coordinate system, and ρ be the equivalent density of the porous medium. f For pore fluid density, Here, η represents the added mass density, κ0 represents the viscosity coefficient of the pore fluid, and G represents the shear modulus.

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

[0052]

[0053] v y Let represent the solid-phase particle velocity components of the seismic wave, i and j represent the grid indices in the x and z directions respectively, n represent the index of the Taylor expansion node, and Δx represent the side length of a unit grid in the x direction. The difference coefficients of the spatial derivative with eighth-order precision are calculated as follows:

[0054]

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

[0056]

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

[0058]

[0059] Load the seismic source into F in equation (5) y The term, k represents the seismic wave time step, and the seismic wave velocity vector. Defined at half a time node, the vector composed of stress components Defined at an integer number of time points, D x With D z These represent the higher-order difference operators in the x and z directions, respectively, and their specific forms are shown in Equation (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 scheme for the electromagnetic wave vector is shown in equation (8):

[0063]

[0064] Δ tEM It represents a time step of an electromagnetic wave.

[0065] S4, N EM The number of electromagnetic iterations is given by equation (9).

[0066]

[0067] In equation (9), Δt is a time step of the seismic wave. After one seismic wave iteration using equations (5) and (7), equation (8) is used to perform N. EM Subsequent electromagnetic wave iteration.

[0068] Example

[0069] See Figure 1 The present invention employs a two-dimensional SHTE mode seismoelectric wave field time-domain numerical simulation method based on high-order difference for practical simulation, including:

[0070] 1) Set the computational domain as x: -3km~3km, z=-3km~3km, where the x-axis is horizontal and the z-axis is vertically downward. Nodes are evenly distributed within the computational domain, with a node spacing of 10m, and a total of 90601 nodes. PML boundary conditions are used on the four boundaries of the computational domain, and the artificial source is set at (0m, 0m).

[0071] 2) Set the porosity parameters in the computational domain, and establish a layered interface at z = -200m. The seismoelectric coupling coefficient of the upper porosity medium is 17.25 * 10⁻⁶. -10 The conductivity of the underground medium is 0.001 S / m, and the dielectric constant of the medium is 7.8 × 10⁻⁶. -11 The permeability of the medium is 4π*10 -7 The equivalent density of the porous medium is 2647 kg / m³.3 The pore fluid density is 1000 kg / m³ 3 The porosity is 0.15, and the added mass density is 300 kg / m³. 3 The pore fluid viscosity coefficient is 0.001 Pa·s, and the permeability of the pore 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 coefficient is set to 0.1 Pa·s, and the other parameters are the same as those of the upper layer.

[0072] 3) Set the difference precision to 8th order and use equation (3) to solve for the higher-order difference coefficients in space;

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

[0074]

[0075] f m This is the dominant frequency of the Reichschild wavelet.

[0076] 5) Set the seismic wave simulation time step to 5*10. -4 s, using the finite difference iterative scheme of the seismic wave vector to perform a seismic wave variable update;

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

[0078] 7) Determine whether the preset electromagnetic time node 5*10 has been reached. -4 If not completed, repeat step 6); if completed, proceed to step 8.

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

[0080] See Figure 2 The numerical simulation results of the two-dimensional SHTE mode under the layered model of porous media are compared with the analytical solution. It can be found that the numerical simulation results can effectively simulate the co-seismic electric field of the source electromagnetic wave, the interface electromagnetic wave and the SH wave.

[0081] Figure 3 It is a snapshot of the wave field of solid particles velocity under a layered model of porous media, which shows the reflection and transmission phenomena of seismic wave particles at the interface.

[0082] Figure 4This is a snapshot of the wavelength of the electric field component Ey under a layered model of porous media, showing the co-seismic electric field induced by the reflected and transmitted waves of the seismic wave.

[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 within the protection scope of the present invention.

Claims

1. A two-dimensional SHTE mode seismoelectric wave field time-domain numerical simulation method based on high-order difference, characterized in that, The method includes: S1. By using the feedback term of the electromagnetic wave field to the seismic wave field in the Pride equation of the two-dimensional SHTE mode. Ignoring this, the time-domain control equations of the two-dimensional SHTE mode wave are derived. S2. The seismic wave computation region is partitioned using staggered grids, and the spatial partial derivative terms of the seismic wave variables are obtained by spatial discretization of the eighth-order precision difference based on the finite-difference time-domain method. The electromagnetic wave computation region is partitioned using Yee grids, and the spatial partial derivative terms of the electromagnetic wave variables are obtained by spatial discretization of the electromagnetic wave variables based on the high-order precision spatial difference method. S3. Substitute the discrete scheme into the time-domain control equation of the two-dimensional SHTE mode wave to obtain the finite difference iterative scheme of the seismic wave vector and the finite difference iterative scheme of the electromagnetic wave vector respectively. Set the PML boundary for the seismic wave region and the electromagnetic wave simulation region, and load the source. S4. After iterating once using the finite difference iterative scheme of the seismic wave vector, perform multiple iterations using the finite difference iterative scheme of the electromagnetic wave vector. S5. Repeat S4 until the number of seismic wave iterations ends, and then display the calculation results. In S1, the time-domain control equation of the two-dimensional SHTE mode wave is expressed as follows: （1）， In equation (1), Indicates the loading source, express Component electric field intensity , They represent Quantity and The magnetic field strength of the component, The electro-seismic coupling coefficient is... The electrical conductivity of the underground medium. The dielectric constant of the medium. The magnetic permeability of the medium, express Directional solid-phase particle velocity, express Directional solid-state fluid relative to particle velocity, and These are the components of the stress tensor in a rectangular coordinate system. The equivalent density of the porous medium, For pore fluid density, To add mass density, The viscosity coefficient of the pore fluid. The permeability of the porous medium, Indicates shear modulus; In S2, the eighth-order precision difference space discretization form of the seismic wave variables is as follows: （2）， , Represent , Grid number of direction, Indicates the index of the selected Taylor expansion node. for Side length of the unit grid in the direction, The difference coefficients of the spatial derivative with eighth-order precision are calculated as follows: （3）， The spatial difference discretization form of electromagnetic wave variables is: （4）。 2. The time-domain numerical simulation method for two-dimensional SHTE mode seismoelectric wave fields based on higher-order differences according to claim 1, characterized in that: The finite difference iteration scheme for the seismic wave velocity vector in S3 is: （5）， （6）， Load the seismic source into equation (5) item, Represents the time step of the seismic wave and the velocity vector of the seismic wave. Defined at half a time node, the vector composed of stress components Defined at an integer number of time points, and Represent direction and The higher-order difference operator for the direction is shown in equation (2). This represents the time step of the seismic waves; The iterative form of the seismic wave stress component vector is: （7）， The finite difference iterative scheme for the electromagnetic wave vector is shown in equation (8): （8）， It represents a time step of an electromagnetic wave.

3. The two-dimensional SHTE mode seismoelectric wave field time-domain numerical simulation method based on high-order difference as described in claim 2, is characterized in that, In S4, The number of electromagnetic iterations is given by equation (9). （9）， After one seismic wave iteration using equations (5) and (7), equation (8) is used for... Subsequent electromagnetic wave iteration.

Images