A forward modeling method for seismic electromagnetic field induced by piezoelectric effect

Abstract

The invention discloses a forward modeling method for generating seismic electromagnetic fields induced by piezoelectric effect, and relates to the technical field of geophysics. Based on the frequency domain piezoelectric effect seismic electromagnetic coupling control equation, the equation group is decoupled into PSVTM mode equations and SHTE mode equations by means of coordinate transformation. Then, the two groups of equations are respectively solved in a three-dimensional horizontal layered model by means of a global matrix method to obtain the seismic electromagnetic wave field in the frequency-wave number domain, and finally, the solution in the time-space domain is obtained by using Hankel transformation and Fourier transformation. By using the invention, the seismic electromagnetic field induced by the piezoelectric effect can be quickly calculated.

Description

Technical Field

[0001] The invention relates to the field of geophysical technology, and in particular to a forward modeling method for generating seismic electromagnetic fields induced by piezoelectric effect. Background Art

[0002] When a force in a specific direction is applied to a quartz crystal (α-quartz, trigonal system 32 point group), its internal lattice structure will be deformed, causing the center of positive and negative charges to shift, forming microscopic electric dipoles. The directions of the electric dipoles generated inside a pure quartz crystal are basically the same. They will induce bound charges on the surface of the quartz, forming a macroscopic polarized electric field. This is the piezoelectric effect of the quartz crystal. When an earthquake occurs, due to the piezoelectric effect of the quartz crystal in the rock, the changing stress field in the crust will produce a changing macroscopic polarized electric field, thereby stimulating electromagnetic signals.

[0003] In the current forward modeling of seismic electromagnetic fields, the wavenumber integration method can be used to calculate the earthquake and electromagnetic fields induced by the piezoelectric effect in three-dimensional horizontal layered media in a rectangular coordinate system. The amount of data that needs to be processed is large. This method involves complex integral operations, especially when dealing with three-dimensional space problems. The amount of calculation is huge, resulting in slow calculation speed and low efficiency. Summary of the invention

[0004] Based on this, it is necessary to provide a forward modeling method for generating seismic electromagnetic fields induced by piezoelectric effect in response to the problems mentioned in the above background technology.

[0005] The object of the present invention can be achieved by the following technical solution: A forward modeling method for generating seismic electromagnetic field induced by piezoelectric effect, comprising the following steps:

[0006] 101: Based on the frequency domain piezoelectric effect seismic electromagnetic coupling equation, the control equation is transformed into the surface harmonic coordinate system by coordinate transformation. According to the symmetry, the three-dimensional equation group can be decoupled into two sets of independent two-dimensional equation groups, namely the PSVTM mode equation and the SHTE mode equation.

[0007] 102: By using the global matrix method to solve the PSVTM mode equation and the SHTE mode equation in the three-dimensional horizontal layered model, the displacement, stress, electric field and magnetic field in the frequency-wavenumber domain are obtained;

[0008] 103: The displacement, stress, electric field and magnetic field in the frequency-wavenumber domain are transformed into the corresponding solutions in the time-space domain through the Hankel transform of the wavenumber k and the inverse Fourier transform of the frequency.

[0009] In some embodiments, a set of orthogonal and complete coordinate systems is introduced to transform the piezoelectric effect seismic electromagnetic coupling equation into a surface harmonic coordinate system, specifically:

[0010] Introduce the surface harmonic coordinate vector:

[0011]

[0012]

[0013]

[0014] in, is a harmonic function. m (kr) is the m-order first-kind Bessel function, k represents the wave number, and i is the imaginary unit. m is related to the symmetry of the earthquake source: for a point force source, m can be 0, ±1; for a point source of the seismic moment tensor, m can be 0, ±1, ±2. r , e θ , e z are the basis vectors in cylindrical coordinates. is called the basis vector of the surface harmonic function. For any vector ξ(r,θ,z) in the cylindrical coordinate system, the surface harmonic coordinate expansion can be performed.

[0015] Among them, ξ S,m (z,k),ξ T,m (z,k) and ξ R,m (z, k) is the expansion coefficient of the surface harmonic coordinates, which is a function of m and k and can be obtained by the following integral:

[0016]

[0017]

[0018]

[0019] [] * Indicates conjugation.

[0020] Conversely, if the plane harmonic coordinate components of ξ are known, its cylindrical coordinate components can also be obtained by the following expression:

[0021]

[0022]

[0023]

[0024] Among them, only one wavenumber integration is required to transform ξ from the wavenumber domain to the space domain.

[0025] In some embodiments, the specific process of decoupling the three-dimensional equation system into two sets of two-dimensional equation systems according to symmetry is:

[0026] The displacement vector u(r,θ,z), body force vector F(r,θ,z), horizontal plane stress vector Q(r,θ,z), electric field E(r,θ,z) and magnetic field H(r,θ,z) in the cylindrical coordinate system are respectively expressed as

[0027]

[0028] Their surface harmonic coordinate expansions can be obtained respectively. Substituting the expansions of u, F, Q, E, H and the polarization vector P into the piezoelectric effect seismic electromagnetic coupling control equation, two sets of independent first-order linear differential equations are obtained through mathematical derivation, which are the PSVTM mode equations:

[0029]

[0030] And the SHTE mode equation:

[0031]

[0032] In some embodiments, the specific process of deriving the surface harmonic coordinate expansion of the piezoelectric modulus is:

[0033] The governing equations of piezoelectric effect-induced seismic electromagnetics in the frequency domain are:

[0034]

[0035]

[0036]

[0037]

[0038] P = dτ′

[0039] in, is the Hamiltonian operator, τ is the stress tensor, u is the displacement vector, and f is the body density vector; is the volume strain, I is the unit matrix, ρ is the medium density, t is the time, ω is the angular frequency; λ and G are the Lamé coefficients, σ, ε and μ represent the conductivity, dielectric constant and magnetic permeability of the rock medium respectively, J f represents the applied current density, P is the electric polarization intensity vector of the medium; the stress vector τ′ is expressed as τ′=(τ rr ,τ θθ ,τ zz ,τ θz ,τ rz ,τ rθ ) T , the piezoelectric modulus d is selected to have '∞m' symmetry, expressed as

[0040]

[0041] Where d0 = 5 × 10 -16 C / N, in which case the quartz rock has a vertical axis of symmetry;

[0042] The cylindrical coordinate components of the electric polarization vector P can be obtained:

[0043] P r =-2d0τ rz

[0044] P θ =-2d0τ θz

[0045] P z =-d0τ rr -d0τ θθ +2d0τ zz

[0046] Further, the polarization vector P is expressed by displacement u and stress Q,

[0047] P r =-2d0Q r

[0048] P θ =-2d0Q θ

[0049]

[0050] Finally, the surface harmonic coordinate expansion of displacement u and stress Q is substituted into it to obtain the surface harmonic expansion of the piezoelectric effect electric polarization vector P:

[0051]

[0052]

[0053]

[0054] In some embodiments, the global matrix method is used to solve the PSVTM model equation and the SHTE model equation in the three-dimensional horizontal layered model, and the specific solution process is:

[0055] Using the formula

[0056]

[0057] And the formula

[0058]

[0059] The amplitudes of earthquake and electromagnetic waves are solved in the PSVTM model equation and the SHTE model equation respectively, where represents the eigenvector of the downlink wave in the nth layer, represents the eigenvector of the upgoing wave in the nth layer, P n represents the propagation matrix of waves in the nth layer, and They represent the amplitude coefficients of the downgoing wave and the upgoing wave in the nth layer respectively, S represents the discontinuity vector of the source; the superscript V represents the PSVTM mode, and H represents the SHTE mode; the linear equations can be solved to obtain the amplitude vectors of seismic and electromagnetic waves.

[0060] In some embodiments, the solution in the time-space domain is obtained using Hankel transform and inverse Fourier transform, and the specific process is as follows:

[0061] Using the formula

[0062]

[0063]

[0064]

[0065] Transform the displacement u, stress Q, electric field E, and magnetic field H in the wave number domain to the space domain, i.e., cylindrical coordinates;

[0066] Then use the inverse Fourier transform formula

[0067]

[0068] Transform the displacement u, stress Q, electric field E, and magnetic field H in the frequency domain to the time domain.

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

[0070] 1. The present invention proposes that the piezoelectric effect of rocks induces earthquakes to generate electromagnetic signals, gives the piezoelectric effect seismic electromagnetic coupling control equation, introduces a set of orthogonal and complete surface harmonic coordinate basis vectors, transforms the piezoelectric effect seismic electromagnetic control equation into the surface harmonic coordinate system, and obtains two sets of independent two-dimensional equations, namely the PSVTM mode equation and the SHTE mode equation. The global matrix method is used to solve the two sets of equations respectively, and the seismic wave field and the electromagnetic wave field can be calculated simultaneously;

[0071] 2. The present invention analytically derives analytical expressions of earthquake and electromagnetic wave characteristic vectors by using plane wave theory, and applies them to the process of solving differential equations using the global matrix method, greatly improving the accuracy and efficiency of the three-dimensional problem in simulation calculations. BRIEF DESCRIPTION OF THE DRAWINGS

[0072] In order to more clearly illustrate the technical solutions in the embodiments of the present application or the conventional technology, the drawings required for use in the embodiments or the conventional technology descriptions are briefly introduced below. Obviously, the drawings described below are only some embodiments of the present application. For ordinary technicians in this field, other drawings can be obtained based on these drawings without creative work.

[0073] Figure 1 A flow chart of a forward modeling method for generating seismic electromagnetic fields induced by piezoelectric effect provided by the present invention;

[0074] Figure 2 It is a schematic diagram of a three-dimensional horizontal layered model of the present invention;

[0075] Figure 3 The seismic and electromagnetic responses induced by the piezoelectric effect of the full-space model of the present invention;

[0076] Figure 4 It is the early electromagnetic wave radiation pattern induced by the piezoelectric effect in the full-space model of the present invention;

[0077] Figure 5 It is a schematic diagram of a half-space model of the present invention;

[0078] Figure 6 The seismic and electromagnetic responses induced by the piezoelectric effect in the half-space model of the present invention;

[0079] Figure 7 This is the influence of the electrical conductivity of the present invention on the electromagnetic response of the piezoelectric effect. DETAILED DESCRIPTION

[0080] In order to make the above-mentioned objects, features and advantages of the present invention more obvious and easy to understand, the specific embodiments of the present invention are described in detail below in conjunction with the accompanying drawings. In the following description, many specific details are set forth to facilitate a full understanding of the present invention. However, the present invention can be implemented in many other ways different from those described herein, and those skilled in the art can make similar improvements without violating the connotation of the present invention, so the present invention is not limited by the specific embodiments disclosed below.

[0081] like Figure 1 FIG. 1 is a schematic diagram of a forward modeling method for generating seismic electromagnetic fields induced by piezoelectric effect, the method comprising:

[0082] Establishment of a three-dimensional horizontal layered medium model. Underground rocks are roughly distributed in layers. In the process of propagating from the earthquake source to the observation point, seismic waves will be reflected and transmitted many times at the interface between different underground rock layers and the free surface. The data actually recorded by the receiving instrument will inevitably contain this part of the reflection and transmission signals. Therefore, in order to be able to explain the actual observation data and understand seismic electromagnetic phenomena, it is reasonable and necessary to approximate the underground rock medium as a layered structure.

[0083] Figure 2 A schematic diagram of the horizontal layered model is given, assuming that each layer of rock in the model is an isotropic uniform elastic medium, ignoring the influence of other structures such as pores, fluids, and cracks. The parameters of rocks in different layers are different, including elastic parameters (density, elastic modulus, etc.), and electrical parameters (conductivity, dielectric constant, etc.). A rectangular coordinate system is used, with the origin placed on the surface, the x-axis representing the north direction; the y-axis representing the east direction; and the z-axis representing the bottom. In order to facilitate calculations using symmetry, the source is placed on the z-axis. Observation points are generally placed near the surface. Above the free surface is the air layer, Figure 2 The first layer is the layer below the ground, and so on. The last layer is a semi-infinite space, which is the Nth layer. For convenience, such a model with N+1 layers is still called an N-layer model, which will not be explained later. For the i-th layer, z i is the depth coordinate of its lower interface.

[0084] In order to use the coordinate transformation method to quickly solve the three-dimensional problem, a set of orthogonal and complete coordinate systems are introduced to transform the piezoelectric effect seismic electromagnetic coupling equation into the surface harmonic coordinate system. This step also transforms the equation in the frequency-space domain into the frequency-wavenumber domain. Introduce the surface harmonic coordinate vector:

[0085]

[0086]

[0087]

[0088] in, is a harmonic function. m (kr) is the m-order first-kind Bessel function, k represents the wave number, and i is the imaginary unit. m is related to the symmetry of the earthquake source: for a point force source, m can be 0, ±1; for a point source of the seismic moment tensor, m can be 0, ±1, ±2. r , e θ , e z are the basis vectors in cylindrical coordinates. is called the basis vector of the surface harmonic function. For any vector ξ(r,θ,z) in the cylindrical coordinate system, the surface harmonic coordinate expansion can be performed.

[0089] Among them, ξ S,m (z,k),ξ T,m (z,k) and ξ R,m (z,k) is the coefficient of the surface harmonic coordinate expansion, which is a function of m and k and can be obtained by the following integral:

[0090]

[0091]

[0092]

[0093] [] * Indicates conjugation.

[0094] Conversely, if the plane harmonic coordinate components of ξ are known, its cylindrical coordinate components can also be obtained by the following expression:

[0095]

[0096]

[0097]

[0098] Among them, ξ only needs to be transformed from the wavenumber domain to the space domain by one wavenumber integration. The displacement vector u(r,θ,z), body force vector F(r,θ,z), horizontal plane stress vector Q(r,θ,z), electric field E(r,θ,z) and magnetic field H(r,θ,z) in the cylindrical coordinate system are respectively expressed by the formula Their surface harmonic coordinate expansions can be obtained respectively. Substituting the expansions of u, F, Q, E, H and the polarization vector P into the piezoelectric effect seismic electromagnetic coupling control equation, two sets of independent first-order linear differential equations are obtained through mathematical derivation, which are the PSVTM mode equations:

[0099]

[0100] And the SHTE mode equation:

[0101]

[0102] Define the eigenvalues ​​and eigenvectors of the coefficient matrix A of the differential equation system in the mth level: A m D m =D m Λ m . Where Λ m Yes A mThe diagonal matrix formed by the eigenvalues ​​of m Yes A m The transformation matrix composed of the eigenvectors of m Can be written as:

[0103]

[0104]

[0105] in, and They represent the eigenvectors of the descending P, SV, SH, TM and TE waves respectively; and They represent the eigenvectors of the upgoing P, SV, SH, TM and TE waves respectively. They represent the basic solutions of the displacement-stress-electromagnetic response excited by the downgoing and upgoing waves. In order to improve the accuracy and speed of the calculation, the analytical expressions of these eigenvectors are derived.

[0106] A m The matrix Λ consisting of the eigenvalues ​​of m In PSVTM and SHTE modes, they can be expressed as:

[0107]

[0108]

[0109] in, and They represent P, S and the wave number of electromagnetic waves propagating in the vertical direction in the mth layer, that is, In the formula There are two roots, take the one with a positive imaginary part. x is the horizontal wave number. s j is the slowness of different waves. Rewrite the displacement-stress-electromagnetic vector b(z) as b m (z) = D m (z)W m (z), where W m (z) is the amplitude vector of seismic and electromagnetic waves at depth z in the mth layer. m (z) is the equation for the unknown quantity:

[0110] At this time, the coefficient matrix of the equation is a diagonal matrix, and its solution is easy to obtain, namely W m (z) = θα m Where Θ is the wave propagation matrix, which can be written as:

[0111]

[0112]

[0113] α m is the amplitude vector of the downgoing wave and the upgoing wave at the interface in the mth layer, which can also be expressed as:

[0114]

[0115]

[0116] α m The wave amplitude W at any depth can be obtained by solving the linear equations by loading source terms and boundary conditions. m (z), and finally obtain the displacement-stress-electromagnetic vector b m (z).

[0117] For seismic waves and electromagnetic waves propagating underground, displacement, stress and electromagnetic field at the interface between any two layers z = z m The continuity condition must be satisfied, that is, b m (z m )=b m+1 (z m ). At the free interface z = z0, the stress in the horizontal plane must be zero and the horizontal component of the electromagnetic field must be continuous, that is,

[0118]

[0119]

[0120] In addition, the contribution of the source needs to be considered. The source contribution S can be expressed as the source depth z = z s It is represented by a displacement-stress-electromagnetic discontinuity vector at , S = b s (z s+ )-b s (z s- ). After derivation, the rectangular coordinate expression of the seismic moment tensor source term can also be obtained, PSVTM mode:

[0121]

[0122] And SHTE mode:

[0123]

[0124] With the given source terms and boundary conditions, a large linear system of equations can be assembled:

[0125]

[0126] Solve the linear equations of the PSVTM mode and the SHTE mode respectively to obtain the amplitude coefficient vector α of the seismic and electromagnetic waves in the stratum. Then, the unknown vector b(z) of displacement-stress-electric field-magnetic field in the frequency-wavenumber domain is obtained by the formula W(z)=Θα and b(z)=D(z)W(z). Finally, the inverse Fourier transform and Hankel transform are used to obtain their results in the time-space domain.

[0127] Considering the full-space model and ignoring the existence of interfaces such as the ground surface, the piezoelectric electromagnetic response of body waves directly excited by the earthquake source is studied. The medium parameters used are: V P =6500km / s, V S =3900km / s,ρ=2800kg / m 3 ,σ=0.001S / m,ε=6.2×10 -11 F / m. Simulate the seismic electromagnetic response excited by the fault slip source, using a double couple source, and the seismic moment component is M xz =M zx =1.531×10 18 Nm, relative to an earthquake of magnitude Mw 6.1. The focal depth is 10km, the origin is placed at the focal point, and the coordinates of the observation point relative to the focal point are (50km, 70km, -10km). Figure 3 The numerical simulation results are given, and it can be seen that earthquakes can stimulate two types of electromagnetic signals. One is the co-seismic electromagnetic signal (black line) that arrives at the observation point at the same time as the seismic wave, and the other is the early electromagnetic signal (blue line) that arrives earlier than the seismic wave. The propagation speed of the early electromagnetic wave is faster, and the signal can be observed almost at the moment the earthquake occurs, so it is helpful for early warning of earthquakes.

[0128] Figure 3 The signal of early electromagnetic waves is weak, which may be because the observation point is far away from the epicenter, resulting in a large attenuation of early electromagnetic waves. Use multiple receiving points to observe early electromagnetic waves simultaneously, and obtain the radiation pattern in the horizontal receiving plane 10km above the epicenter, such as Figure 4 The results show that different components of the early electromagnetic field have different radiation patterns, and the early electromagnetic wave signals are stronger near the epicenter, where the electric field can reach 1μV / m and the magnetic field can reach 0.001nT, both of which can be observed. This proves that the early electromagnetic waves of earthquakes caused by the piezoelectric effect may help us with earthquake early warning.

[0129] Using the same medium parameters, the seismic electromagnetic response generated by the earthquake point source in the half-space model is simulated. The calculation model is as follows Figure 5 As shown in the figure, the origin of the coordinates is on the surface, the earthquake source is 10km underground, and the coordinates of the receiving point are (50km, 70km, 0.1m). Figure 6The simulated seismic displacement field, electric field and magnetic field are given. The results show that there are also two kinds of electromagnetic responses, namely co-seismic electromagnetic signals and early electromagnetic signals. However, compared with the full-space model, the electromagnetic response obtained by the half-space model has a more complex waveform, because the influence of the free surface and air layer on seismic waves and electromagnetic waves is considered. The seismic electromagnetic response induced by the piezoelectric effect is not only related to the model adopted, but also to the parameters of the medium.

[0130] Consider the influence of the conductivity of the medium on the seismic electromagnetic response induced by the piezoelectric effect. Figure 5 In the half-space model, the rock conductivity σ is set to σ=0.0001S / m, σ=0.001S / m, σ=0.01S / m, σ=0.1S / m and σ=1S / m. The corresponding seismic electromagnetic response is calculated. Figure 7 The curves of the maximum amplitude of the electromagnetic field with different conductivity are given. The results show that the co-seismic electric field decreases with the increase of conductivity, and the co-seismic magnetic field hardly changes with the change of conductivity.

[0131] The method of the present invention provides a forward modeling method for generating an electromagnetic field of an earthquake induced by a piezoelectric effect. Solving a three-dimensional forward modeling problem directly in a rectangular coordinate system requires a huge amount of calculation and has a very low calculation efficiency. By using a coordinate transformation method to transform the control equation to a harmonic coordinate system for solution, a set of three-dimensional equations can be decoupled into two sets of independent two-dimensional equations, greatly reducing the amount of calculation. At the same time, the equation is transformed from the frequency-space domain to the frequency-wavenumber domain, and then the global matrix method is used in the wavenumber domain to quickly solve the equation. Finally, the Hankel transform and the inverse Fourier transform are used to obtain a solution in the time-space domain. The piezoelectric effect can induce earthquakes to generate co-seismic electromagnetic signals and early electromagnetic signals, the latter of which is helpful for early warning of earthquakes.

[0132] The above-mentioned embodiments only express several implementation methods of the present invention, and the descriptions thereof are relatively specific and detailed, but they cannot be understood as limiting the scope of the invention patent. It should be pointed out that, for ordinary technicians in this field, several variations and improvements can be made without departing from the concept of the present invention, and these all belong to the protection scope of the present invention. Therefore, the protection scope of the patent of the present invention shall be subject to the attached claims.

Claims

1. A forward modeling method for generating seismic electromagnetic field induced by piezoelectric effect, characterized in that: The following steps are involved: 101: Based on the seismic electromagnetic coupling equation caused by the piezoelectric effect in the frequency domain, the control equation is transformed into the surface harmonic coordinate system by the coordinate transformation method. According to the symmetry, the three-dimensional equation group can be decoupled into two sets of independent two-dimensional equation groups, namely the PSVTM mode equation and the SHTE mode equation; 102: Using the global matrix method, the PSVTM and SHTE mode equations are solved in a three-dimensional horizontal layered model to obtain displacement, stress, electric field, and magnetic field in the frequency-wavenumber domain; 103: The displacement, stress, electric field and magnetic field in the frequency-wavenumber domain are transformed into the corresponding solutions in the time-space domain by the Hankel transform of the wavenumber k and the inverse Fourier transform of the frequency; Based on the seismic electromagnetic coupling equation caused by the piezoelectric effect in the frequency domain, the specific process of converting the control equation into the surface harmonic coordinate system by coordinate transformation is as follows: 201: Introducing the surface harmonic coordinate basis vector: in, is the surface harmonic function, J m (kr) is the m-th order Bessel function of the first kind, k represents the wave number, i is the imaginary unit, e r , e θ , e z are the basis vectors of cylindrical coordinates; It is called the basis vector of the surface harmonic function; 202: For any vector ξ(r,θ,z) in the cylindrical coordinate system, perform surface harmonic coordinate expansion to obtain: Among them, ξ S,m (z,k),ξ T,m (z,k) and ξ R,m (z, k) is the expansion coefficient of the surface harmonic coordinates, which is a function of m and k and can be obtained by the following integral: [] * It means taking conjugate; 203: Conversely, if the plane harmonic coordinate components of ξ are known, its cylindrical coordinate components can be obtained by the following expression: Among them, ξ only needs to be transformed from the wavenumber domain to the space domain by one wavenumber integration, and the displacement vector u(r,θ,z), body force vector F(r,θ,z), horizontal plane stress vector Q(r,θ,z), electric field E(r,θ,z) and magnetic field H(r,θ,z) in the cylindrical coordinate system are used respectively: The corresponding surface harmonic coordinate expansion is derived; The specific process of deriving the surface harmonic coordinate expansion of the piezoelectric modulus is as follows: The control equations of piezoelectric effect-induced seismic electromagnetic in the frequency domain are: P = dτ′ in, is the Hamiltonian operator, τ is the stress tensor, u is the displacement vector, and f is the body density vector; is the volume strain, I is the unit matrix, ρ is the medium density, t is the time, ω is the angular frequency; λ and G are the Lamé coefficients, σ, ε and μ represent the conductivity, dielectric constant and magnetic permeability of the rock medium respectively, J f represents the applied current density, P is the electric polarization intensity vector of the medium; the stress vector τ′ is expressed as τ′=(τ rr ,τ θθ ,τ zz ,τ θz ,τ rz ,τ rθ ) T , the piezoelectric modulus d is selected to have '∞m' symmetry, expressed as: Where d0 = 5 × 10 -16 C / N, at this time the quartz rock has a vertical symmetry axis; thus the cylindrical coordinate components of the polarization vector P are obtained: P r =-2d0τ rz P θ =-2d0τ θz P z =-d0τ rr -d0τ θθ +2d0τ zz Further, the polarization vector P is expressed by displacement u and stress Q, P r =-2d0Q r P θ =-2d0Q θ Finally, the surface harmonic coordinate expansion of displacement u and stress Q is substituted into it to obtain the surface harmonic expansion of the piezoelectric effect electric polarization vector P: The global matrix method is used to solve the PSVTM mode equation and the SHTE mode equation in the three-dimensional horizontal layered model. The specific solution process is as follows: Using the formula And the formula The amplitudes of seismic waves and electromagnetic waves are solved in the PSVTM model equation and the SHTE model equation respectively, where: represents the eigenvector of the downlink wave in the nth layer, represents the eigenvector of the upgoing wave in the nth layer, P n represents the propagation matrix of waves in the nth layer, and denote the amplitude coefficients of the downgoing and upgoing waves in the nth layer, respectively; S denotes the source discontinuity vector; the superscript V denotes the PSVTM mode, and H denotes the SHTE mode; the amplitude vectors of the seismic and electromagnetic waves are obtained by solving the equations; The solution in the time-space domain is obtained using Hankel transform and inverse Fourier transform. The specific process is as follows: Using the formula Transform the displacement u, stress Q, electric field E, and magnetic field H in the wave number domain to the spatial domain, i.e., cylindrical coordinates; Then use the inverse Fourier transform formula Transform the displacement u, stress Q, electric field E, and magnetic field H in the frequency domain to the time domain.

Images