A forward modeling method for seismoelectric wave fields in VTI porous media

Abstract

The present invention discloses a forward modeling method for the seismoelectric wave field in VTI porous media. A set of surface harmonic coordinate basis vectors are introduced to expand the frequency domain control equation in surface harmonic coordinates and decouple it into the PSVTM mode and the SHTE mode. The frequency-wavenumber domain seismoelectric wave field in the horizontal layered model is then solved by the global matrix method. Finally, the frequency and wavenumber are integrated using the Hankel transform and the fast Fourier transform to obtain the time-space domain seismoelectric wave field. Compared with the forward modeling method based on the isotropic porous media theory, the calculation results of the present invention are closer to the observed values. It has high computational efficiency, provides a basis for studying the contribution of elastic waves to electromagnetic waves in VTI porous media, and provides a forward modeling operator for seismic electromagnetic wave inversion. The present invention based on VTI media is significantly superior to other methods based on isotropic conditions, with higher sensitivity and higher accuracy.

Description

Technical Field

[0001] The present invention belongs to the field of geophysics, and in particular relates to a forward modeling method for seismoelectric wave fields in VTI porous media. Background Art

[0002] Current models for underground exploration are mostly based on the assumption that the subsurface medium is isotropic. However, the actual subsurface medium is complex and generally exhibits anisotropy in both mechanical and electrical properties. For example, aligned fractures or periodically distributed thin layers can cause rocks to exhibit transverse isotropy (TI), resulting in seismic wave velocity anisotropy. The presence of aligned fractures also causes anisotropy in electrical conductivity and permeability. Therefore, accounting for anisotropy is essential.

[0003] VTI media is a typical anisotropic medium that can more accurately characterize underground media conditions. Differences in conductivity and permeability often occur at the interface of underground media. Traditional methods based on isotropic media assume that the conductivity and permeability values ​​are the same in all directions, which can lead to differences between simulation results and experimental results, or even incorrect conclusions. Based on the seismoelectric effect in VTI media, considering conductivity and permeability as two values ​​in the vertical and horizontal directions can more accurately reflect the conditions at the underground medium interface, making the conclusions more consistent with the actual underground media conditions.

[0004] Traditional three-dimensional simulations are computationally intensive and difficult to implement. However, the global matrix method can rapidly calculate the seismoelectric wavefield response in horizontally layered VTI porous media and process three-dimensional wavefields in one-dimensional inhomogeneous media. Summary of the Invention

[0005] In view of the shortcomings of the existing technology, the purpose of the present invention is to provide a forward modeling method for the seismoelectric wave field in VTI porous media to simulate the three-dimensional seismoelectric wave field in horizontally layered VTI porous media.

[0006] The purpose of the present invention can be achieved through the following technical solutions:

[0007] A forward modeling method for seismoelectric wave fields in VTI porous media, the specific process includes:

[0008] S1. Introducing surface harmonic basis vectors, the frequency domain control equations are expanded in surface harmonic coordinates, resulting in two sets of decoupled control equations. One set is only related to qP waves, qSV waves, and TM waves, and the other set is only related to SH waves and TE waves, which are called PSVTM mode and SHTE mode, respectively. The introduction of cylindrical harmonic coordinates simplifies the three-dimensional problem into a two-dimensional problem, and at the same time transforms the vector from the space domain to the wavenumber domain.

[0009] S2. Solve the frequency-wavenumber domain wave field in the horizontal layered model using a global matrix method. Given the source and medium parameters, all boundary conditions and source contributions are used to obtain a system of linear equations with wave amplitude as the unknown in each layer of the medium. Solve this system of equations to obtain wave fields such as solid phase displacement, seepage displacement, electric field, and magnetic field.

[0010] S3. Use Hankel transform and fast Fourier transform to transform the frequency-wavenumber domain results into the time-space domain wave field.

[0011] Furthermore, a set of surface harmonic coordinate basis vectors is introduced:

[0012]

[0013]

[0014]

[0015] In the above formula, e r represents the unit vector in the radial direction in the cylindrical coordinate system, e θ represents the unit vector in the tangential direction in the cylindrical coordinate system, e z represents the unit vector in the vertical direction in the cylindrical coordinate system;

[0016] Y m (r,θ)=J m (kr)e imθ (m=0,±1,±2,...), J m (kr) is the Bessel function of the first kind, order m, and k is the horizontal wave number;

[0017] ξ is the solid phase displacement vector u, the osmotic displacement w, the body force densities F and f acting on the entire medium and fluid phase, the electric field vector E, the magnetic field vector H, and can also be the stress vector at any position on the horizontal plane: τ = τ rz e r +τ θz e θ +τ zz e z ;

[0018] m depends on the symmetry of the source in the circumferential direction. For explosive point sources or vertical point force sources with axial symmetry, m = 0; for horizontal point force sources, m = ± 1; for seismic moment tensors, m ≤ 2;

[0019] and represents the surface harmonic coordinate component of vector ξ in the frequency-wavenumber domain.

[0020] Furthermore, the global matrix is ​​processed to avoid numerical overflow.

[0021] For PSVTM mode, the depth (z n <z<z n-1 ) at W n The solution has the following form:

[0022]

[0023] And for SHTE mode:

[0024]

[0025] represents the amplitude of each downgoing wave at the upper interface of the nth layer of porous medium;

[0026] and represents the amplitude of each upgoing wave at the lower interface of the nth layer of porous medium;

[0027] W n represents the amplitude vector in the nth layer of medium;

[0028] q Pf is the vertical slowness of the fast longitudinal wave, q Ps is the vertical slowness of the slow longitudinal wave, q SV is the vertical slowness of the SV wave, q TM is the vertical slowness of the TM wave, q SH is the vertical slowness of the SH wave, q TE is the vertical slowness of the TE wave;

[0029] z n-1 represents the depth of the upper interface of the nth layer of VTI porous medium, z n Indicates the depth of the lower interface of the nth layer of VTI porous medium.

[0030] Furthermore, the fast Hankel transform and fast Fourier transform are used to integrate the frequency and wave number to obtain the time-space domain wave field, using the formula:

[0031]

[0032]

[0033]

[0034] Transform the frequency-wavenumber domain vector to the frequency-space domain;

[0035] Then by the formula

[0036]

[0037] Transform the frequency-space domain vector to the time-space domain.

[0038] Beneficial effects of the present invention:

[0039] 1. The forward modeling method for the seismoelectric wave field in VTI porous media proposed in this paper introduces a set of surface-harmonic coordinate basis vectors, expands the frequency-domain governing equations in surface-harmonic coordinates, and decouples them into the PSVTM model and the SHTE model. The global matrix method is then used to solve the frequency-wavenumber domain seismoelectric wave field in the horizontal layered model. Finally, the frequency and wavenumber are integrated using the Hankel transform and the fast Fourier transform to obtain the time-space domain seismoelectric wave field. Compared with the forward modeling method based on isotropic porous media theory, the calculation results of this invention are closer to the observed values.

[0040] 2. The forward modeling method for the seismoelectric wave field in VTI porous media proposed in this invention has high computational efficiency, provides a basis for studying the contribution of elastic waves to electromagnetic waves in VTI porous media, and provides a forward modeling operator for seismic electromagnetic wave inversion. BRIEF DESCRIPTION OF THE DRAWINGS

[0041] In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the following briefly introduces the drawings required for use in the embodiments or the description of the prior art. Obviously, for ordinary technicians in this field, other drawings can be obtained based on these drawings without any creative work.

[0042] Figure 1 1 is a flow chart of a forward modeling method according to an embodiment of the present invention;

[0043] Figure 2 This is a schematic diagram of a horizontal layered model according to Example 1 of the present invention;

[0044] Figure 3 1 is a comparison diagram of the wave field calculation results of the isotropic model and the VTI porous medium model in Example 1 of the present invention;

[0045] Figure 4 This is a schematic diagram of a two-layer model of Example 2 of the present invention;

[0046] Figure 5 This is a diagram of seismoelectric wave field signals received by the two-layer model observation points in Example 2 of the present invention;

[0047] Figure 6 Schematic diagram of the effect of elastic parameters on the seismoelectric signal of a VTI medium according to Example 2 of the present invention;

[0048] Figure 7 Schematic diagram of the effect of permeability on the seismoelectric signal of a VTI medium according to Example 2 of the present invention;

[0049] Figure 8This is a schematic diagram of the effect of conductivity on the seismoelectric signal of the VTI medium in Example 2 of the present invention. DETAILED DESCRIPTION

[0050] The following will clearly and completely describe the technical solutions in the embodiments of the present invention in conjunction with the accompanying drawings. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without making any creative efforts shall fall within the scope of protection of the present invention.

[0051] This application is suitable for underground exploration using the global matrix method, and can efficiently realize the numerical simulation of three-dimensional seismoelectric wave fields in horizontally layered VTI porous media.

[0052] Example 1

[0053] like Figure 1 As shown in FIG, a forward modeling method for the seismoelectric wave field in VTI porous media includes the following steps:

[0054] S1. Introducing surface harmonic basis vectors, the frequency domain control equations are expanded in surface harmonic coordinates, resulting in two sets of decoupled control equations. One set is only related to qP waves, qSV waves, and TM waves, and the other set is only related to SH waves and TE waves, which are called PSVTM mode and SHTE mode, respectively. The introduction of cylindrical harmonic coordinates simplifies the three-dimensional problem into a two-dimensional problem, and at the same time transforms the vector from the space domain to the wavenumber domain.

[0055] S2. Solve the frequency-wavenumber domain wave field in the horizontal layered model using a global matrix method. Given the source and medium parameters, all boundary conditions and source contributions are used to obtain a system of linear equations with wave amplitude as the unknown in each layer of the medium. Solve this system of equations to obtain wave fields such as solid phase displacement, seepage displacement, electric field, and magnetic field.

[0056] S3. Use Hankel transform and fast Fourier transform to transform the frequency-wavenumber domain results into the time-space domain wave field.

[0057] Specifically, a horizontal layering model is established. The Earth's interior is heterogeneous, but overall, the heterogeneity in the vertical direction is generally greater than in the horizontal direction. Actual stratigraphic structures often exhibit a structural pattern similar to horizontal layering. Therefore, the Earth's medium can be considered a horizontal layering model consisting of multiple parallel VTI layers. This is a commonly used approximate geological model that is homogeneous and isotropic only in the horizontal plane. Figure 2A schematic diagram of the horizontal layered model is given, where the z axis of the coordinate system is consistent with the symmetry axis of each layer. The model contains N+1 layers of porous media, the bottom layer of which is a uniform semi-infinite space, the free surface at z=z0, and the air above. The earthquake source is placed at a depth of z=z s The axis of the s ). The seismoelectric wave field in each layer of the medium in the model obeys the VTI porous medium governing equation.

[0058] The governing equations of seismic electromagnetic wave propagation in VTI porous media are converted into cylindrical coordinates. Next, a set of surface harmonic coordinate basis vectors are introduced:

[0059]

[0060]

[0061]

[0062] In the above formula, e r represents the unit vector in the radial direction in the cylindrical coordinate system, e θ represents the unit vector in the tangential direction in the cylindrical coordinate system, e z Represents the unit vector in the vertical direction in the cylindrical coordinate system.

[0063] Y m (r,θ)=J m (kr)e imθ (m=0,±1,±2,...), and J m (kr) is the first kind m-order Bessel function, where k is the horizontal wave number. In the frequency domain, any vector in the cylindrical coordinate system It can be expanded into the following form:

[0064]

[0065] In the above formula, ξ can be the solid phase displacement vector u, the osmotic displacement w, the body force density F and f acting on the entire medium and fluid phase, the electric field vector E, the magnetic field vector H, or the stress vector at any position on the horizontal plane: τ = τ rz e r +τ θz e θ +τ zz e z .

[0066] It can be the solid phase displacement vector u, seepage displacement and Body force vector F, or stress vector t = σ on the horizontal plane rz er +σ θz e θ +σ zz e z The summation range m in the formula depends on the symmetry of the source in the circumferential direction. For point sources with axial symmetry such as explosion point sources and vertical point force sources, m = 0; for horizontal point force sources, m = ±1; for seismic moment tensor, m≤2. and represents the surface harmonic coordinate component of vector ξ in the frequency-wavenumber domain.

[0067] For vector Its cylindrical coordinate components Harmonic coordinate components The following conversion relationships exist:

[0068]

[0069]

[0070]

[0071]

[0072]

[0073]

[0074] In the above equation, the symbol * represents the complex conjugate. Expanding all vectors in the governing equations using surface harmonic vectors yields two sets of decoupled governing equations: one set related only to qP waves, qSV waves, and TM waves, and the other set related only to SH waves and TE waves, referred to as the PSVTM and SHTE modes, respectively. The introduction of surface harmonic vectors decomposes the original three-dimensional problem in cylindrical coordinates into two two-dimensional problems, with the wave fields in both modes conforming to plane wave theory. It also transforms the vectors from the spatial domain to the wavenumber domain. The governing equations for each mode can be expressed as follows:

[0075]

[0076] In the above formula, N is the displacement-stress vector,

[0077] The superscript "V" indicates the PSVTM mode, and "H" indicates the SHTE mode. A is an 8th-order matrix in the PSVTM mode, and A is a 4th-order matrix in the SHTE mode.

[0078] F is a vector related to the source, called the body force vector:

[0079]

[0080]

[0081] Ignoring the body force vector, the governing equations degenerate into:

[0082]

[0083] Let iωΛ and represent the eigenvalue matrix and eigenvector matrix respectively, we can get: Λ is expressed as:

[0084]

[0085]

[0086] In the above formula, q Pf is the vertical slowness of the fast longitudinal wave, q Ps is the vertical slowness of the slow longitudinal wave, q SV is the vertical slowness of the SV wave, q TM is the vertical slowness of the TM wave, q SH is the vertical slowness of the SH wave, q TE is the vertical slowness of the TE wave. When the matrix A is determined, Λ and It can be obtained by eigenvalue decomposition and introducing linear transformation: The control equation can be obtained as follows: W is the amplitude vector, whose components represent the amplitudes of the upgoing and downgoing waves in each layer of the medium.

[0087] Use W n Represents the amplitude vector in the nth layer of the medium.

[0088] Process the global matrix to avoid numerical overflow,

[0089] For PSVTM mode, the depth (z n <z<z n-1 ) at W n The solution has the following form:

[0090] And for SHTE mode:

[0091]

[0092] In the above formula Represents the amplitude of each downgoing wave at the upper interface of the nth layer of porous medium.

[0093] z n-1 represents the depth of the upper interface of the nth layer of VTI porous medium, z n Indicates the depth of the lower interface of the nth layer of VTI porous medium.

[0094] and represents the amplitude of each upgoing wave at the lower interface of the nth layer of porous media. These amplitude coefficients can be solved using source expressions and boundary conditions. When the surface of a porous medium has an open pore boundary condition, it means that the pore fluid can flow freely across the interface. This means that the pore fluid pressure at the interface is zero. In this case, the following boundary conditions exist at the surface:

[0095] For PSVTM mode, and For SH mode,

[0096] Joint relation The boundary conditions at the free surface can be rewritten as: For PSVTM mode:

[0097]

[0098]

[0099]

[0100]

[0101] For SHTE mode:

[0102]

[0103]

[0104] J + (1,1:2)=D(2,3:4)

[0105]

[0106]

[0107] At the interface between the nth layer and the n+1th layer, z=z n At this point, the displacement-stress vector is continuous, that is, N n (z n )=N n+1 (z n ). It can be further expressed as:

[0108] So there is

[0109] In the above formula yes The block matrix form of .

[0110] For PSVTM mode, where h n =z n -z n-1 Indicates the thickness of the nth layer of porous medium.

[0111] where h n+1 =z n+1 -z n . is an identity matrix. For the SHTE mode,

[0112] The contribution of the earthquake source is as follows: Assuming that at depth z = z s There is a virtual interface at which the vector N is discontinuous and the relationship S=N s+1 (z s )-N s (z s )=F1+AF2. Where S is called the discontinuity vector, and F1 and F2 can be obtained by the following formula:

[0113] The seismic moment tensor source can be expressed in terms of equivalent body density as:

[0114] When m=0

[0115]

[0116]

[0117] When m=±1

[0118]

[0119]

[0120] When m=±2

[0121]

[0122]

[0123] Given the source and medium parameters, using all boundary conditions and source contributions, we can obtain a and is a system of linear equations with unknown parameters

[0124]

[0125] Solving this set of equations, we can obtain the amplitudes of the upgoing and downgoing waves in each layer of porous media, and then use the formula The frequency-wavenumber domain wavefield in each layer is calculated.

[0126] In order to obtain the wave field in the time-space domain, the frequency and wavenumber are integrated using the fast Hankel transform and the fast Fourier transform.

[0127] Using the formula

[0128]

[0129]

[0130]

[0131] Transform the frequency-wavenumber domain vector to the frequency-space domain, ξ z Can be applied to u z , w z , τ zz ,-P,E z and H z ,ξ r Can be applied to u r ,w r ,τ rz ,E r and -H r ,ξ θ Can be applied to u θ ,w θ ,τ θz ,E θ and -H θ .

[0132] Then by the formula

[0133]

[0134] The frequency-space domain vector is transformed into the time-space domain, where ξ can be a vector such as u, w, E, H, etc.

[0135] To find the components in a rectangular coordinate system, use the following coordinate transformation:

[0136]

[0137] In order to verify the correctness of the method in this invention, a half-space model is designed with the following model parameters: c 11 =62.43GPa,c 13 =14.68GPa,c 33 =62.43GPa,c 44 =23.85GPa,c 66 =23.85GPa、φ=0.15、κ H0 =0.1Darcy,κV0 =0.1Darcy,α H∞ =3,α V∞ =3,ρ=2650kg / m 3 , ρ f =1000kg / m 3 , eta = 0.001Pa·s, C0 = 0.001mol·L -1 , K s =35.7GPa, K f =2.25GPa,Λ H =12.65×10 -6 m、Λ V =12.65×10 -6 m. This is equivalent to the isotropic medium case, and the VTI porous medium is degenerated into an isotropic medium. A double-couple source is used to excite seismic waves, and its seismic moment component is M xz =M zx =1.5×10 18 In order to avoid the influence of surface reflection waves, the source is placed at a depth z far away from the surface. s = 10km. The position of the observation point is (x = 50km, y = 0, z = 50km). Figure 3 The solid line in the middle is the result of the analytical solution of the isotropic model, and the dotted line is the result of the degenerate VTI porous medium model (the two lines overlap in the figure). It can be seen that the results of the two models are very consistent, verifying the effectiveness of this application.

[0138] Example 2

[0139] The sensitivity of VTI porous media is often stronger than that of isotropic media. When there is a discontinuous interface in the subsurface, more accurate conclusions can be obtained based on the VTI porous media. To further verify the effectiveness of the method, the seismic response of VTI porous media is compared with that of isotropic models. A two-layer model is used. Figure 4 , the first layer is isotropic medium, and the second layer is divided into isotropic and VTI conditions for comparison.

[0140] The first layer still uses the above parameters, and the second layer parameters are set to: 11 =30.34GPa,c 13 =3.82GPa,c 33 =30.34GPa,c 44 =13.26GPa,c 66 =13.26GPa、φ=0.15、κ H0 =0.1Darcy,κ V0 =0.1Darcy,α H∞ =3,αV∞ =3,ρ=2650kg / m 3 , ρ f =1000kg / m 3 , eta = 0.001Pa·s, C0 = 0.001mol·L -1 , K s =37.9GPa, K f =2.25GPa,Λ H =12.65×10 -6 m、Λ V =12.65×10 -6 m.

[0141] Some of the second layer parameters can be adjusted to meet different VTI medium conditions. xx =M yy =M zz =2.54×10 7 The explosion source is Nm, the source time function is Ricker wavelet with a central frequency of 200 Hz, 30 observation points are symmetrically distributed on both sides of the source (-75m, -70m, ... -5m, 5m, ... 70m, 75m), and the burial depth is 0.5m.

[0142] Consider the impact of these types of situations:

[0143] ① Elastic modulus

[0144] Considering anisotropy parameters and Set separately

[0145] (a1)c 11 =0.5c 33 (ν=-0.25)

[0146] (a2)c 11 =c 33 (ν=0, isotropic)

[0147] (a3)c 11 =2c 33 (ν=0.5)

[0148] (a4)c 11 =3c 33 (ν=1)

[0149] (b1)c 66 =0.5c 44 (γ=-0.25)

[0150] (b2)c 66 =c 44 (γ=0i.e., isotropic)

[0151] (b3)c 66 =2c 44 (γ=0.5)

[0152] (b4)c 66 =3c 44 (γ=1)

[0153] ②. Penetration rate

[0154] Considering the horizontal permeability κ H0 and vertical permeability κ V0 , set them separately

[0155] (c1)κ H0 =0.01Darcy

[0156] (c2)κ H0 =1Darcy(κ V0 Keep 0.01Darcy)

[0157] (d1)κ V0 =0.01Darcy

[0158] (d2)κ V0 =1Darcy(κ H0 Keep 0.01Darcy)

[0159] ③ Conductivity

[0160] The curvature determines the conductivity, considering the horizontal curvature α H∞ and vertical curvature α V∞ , set them separately

[0161] (e1)α H∞ =2

[0162] (e2)α H∞ =4(α V∞ Keep 3)

[0163] (f1)α V∞ =2

[0164] (f2)α V∞ =4(α H∞ Keep 3)

[0165] Figure 5 The medium dd event is the interface electromagnetic signal reflected by the underground discontinuous interface, which is amplified 1000 times. Figure 6 , Figure 7 and Figure 8The three factors are shown below, each corresponding to the influence of three different VTI media. All three factors significantly affect the electromagnetic signal at the interface, indicating that the VTI-based method is significantly superior to the isotropic method, achieving higher sensitivity and accuracy.

[0166] Solving the three-dimensional wave equations of VTI porous media using a full three-dimensional numerical solver is computationally intensive and difficult to implement. However, this method directly solves the governing equations of VTI porous media in a horizontally layered model to obtain a three-dimensional seismoelectric wave field, which has high computational efficiency. The propagation of seismic waves in two-dimensional geometric structures is very different from that in three-dimensional geometric structures, and two-dimensional simulations cannot give the absolute amplitude of seismic signals. The present invention provides a basis for studying the contribution of elastic waves to electromagnetic waves in VTI porous media and provides a forward operator for seismic electromagnetic wave inversion.

[0167] Throughout this specification, references to terms such as "one embodiment," "example," or "specific example" indicate that the specific features, structures, materials, or characteristics described in conjunction with that embodiment or example are included in at least one embodiment or example of the present invention. In this specification, schematic representations of these terms do not necessarily refer to the same embodiment or example. Furthermore, the specific features, structures, materials, or characteristics described may be combined in any suitable manner in any one or more embodiments or examples.

[0168] The basic principles, main features, and advantages of the present invention are shown and described above. Those skilled in the art should understand that the present invention is not limited to the above embodiments. The above embodiments and descriptions are merely illustrative of the principles of the present invention. Various changes and modifications may be made to the present invention without departing from the spirit and scope of the present invention, and such changes and modifications fall within the scope of the invention as claimed.

Claims

1. A forward modeling method for seismoelectric wave fields based on VTI porous media, characterized by: The specific process includes: S1. Introducing surface harmonic basis vectors, the frequency domain control equations are expanded in surface harmonic coordinates, resulting in two sets of decoupled control equations. One set is only related to qP waves, qSV waves, and TM waves, and the other set is only related to SH waves and TE waves, which are called PSVTM mode and SHTE mode, respectively. The introduction of cylindrical harmonic coordinates simplifies the three-dimensional problem into a two-dimensional problem, and at the same time transforms the vector from the space domain to the wavenumber domain. S2. Solve the frequency-wavenumber domain wave field in the horizontal layered model using a global matrix method. Given the source and medium parameters, all boundary conditions and source contributions are used to obtain a system of linear equations with the wave amplitude as the unknown in each layer of the medium. Solve this system of equations to obtain the solid phase displacement, seepage displacement, electric field, and magnetic field wave fields. S3. Use Hankel transform and fast Fourier transform to transform the frequency-wavenumber domain results into the time-space domain wave field; Introduce a set of surface harmonic coordinate basis vectors: In the above formula, e r represents the unit vector in the radial direction in the cylindrical coordinate system, e θ represents the unit vector in the tangential direction in the cylindrical coordinate system, e z represents the unit vector in the vertical direction in the cylindrical coordinate system; Y m (r,θ)=J m (kr)e imθ J m (kr), m=±1,±2,...,J m (kr) is the Bessel function of the first kind, order m, and k is the horizontal wave number; ξ is the solid phase displacement vector u, the osmotic displacement w, the body force densities F and f acting on the entire medium and fluid phase, the electric field vector E, the magnetic field vector H, and can also be the stress vector at any position on the horizontal plane: τ = τ rz e r +τ θz e θ +τ zz e z ; m depends on the symmetry of the source in the circumferential direction. For explosive point sources or vertical point force sources with axial symmetry, m = 0; for horizontal point force sources, m = ± 1; for seismic moment tensors, m ≤ 2; and represents the surface harmonic coordinate component of vector ξ in the frequency-wavenumber domain; Process the global matrix to avoid numerical overflow, For PSVTM mode, depth z, Z n <Z<Z n-1 W n The solution has the following form: And for SHTE mode: represents the amplitude of each downgoing wave at the upper interface of the nth layer of porous medium; and represents the amplitude of each upgoing wave at the lower interface of the nth layer of porous medium; W n represents the amplitude vector in the nth layer of medium; q Pf is the vertical slowness of the fast longitudinal wave, q Ps is the vertical slowness of the slow longitudinal wave, q SV is the vertical slowness of the SV wave, q TM is the vertical slowness of the TM wave, q SH is the vertical slowness of the SH wave, q TE is the vertical slowness of the TE wave; z n-1 represents the depth of the upper interface of the nth layer of VTI porous medium, z n Indicates the depth of the lower interface of the nth layer of VTI porous medium.

2. The forward modeling method of seismoelectric wave field based on VTI porous media according to claim 1 is characterized in that: The time-space domain wave field is obtained by integrating the frequency and wave number using the fast Hankel transform and the fast Fourier transform, using the formula: Transform the frequency-wavenumber domain vector to the frequency-space domain; Then by the formula Transform the frequency-space domain vector to the time-space domain.

Images