A forward modeling method for seismic wavefields excited by a point source in horizontally layered dual-porosity media

Abstract

The present invention discloses a forward modeling method for seismic wave fields excited by a point source in a horizontally layered dual-porosity medium. The method converts the governing equations of seismic wave propagation in the time-domain dual-porosity medium into the frequency domain through Fourier transform. Then, a set of surface-harmonic coordinate basis vectors are introduced to expand the frequency-domain governing equations in the surface-harmonic coordinates and decouple them into PSV modes and SH modes. The frequency-wavenumber domain wave field in the horizontally layered model is solved through a global matrix method. Finally, the frequency and wavenumber are integrated using fast Hankel transform and fast Fourier transform to obtain the time-space domain wave field. Compared with the forward modeling method based on the theory of single-porosity media, the calculated results in the low-frequency band are closer to the observed values. The three-dimensional seismic wave field in the dual-porosity medium is obtained through an analytical algorithm with high computational efficiency. The method provides a basis for studying the influencing factors of P-wave and S-wave attenuation in the dual-porosity medium and provides a forward modeling operator for seismic wave inversion.

Description

Technical Field

[0001] The invention belongs to the field of geophysical technology, and in particular relates to a forward modeling method for a seismic wave field excited by a point source in a horizontally layered dual-porosity medium. Background Art

[0002] Traditional elastic wave theory assumes that the medium is a single-phase elastic solid. However, natural seismology and exploration seismology study seismic waves propagating through underground rocks. Underground rocks are often porous, with fluids coexisting in their pores. These fluids can move relative to the solid skeleton, making them multiphase rather than single-phase elastic solids. The presence of pore fluids affects the velocity and attenuation of seismic waves. Therefore, it is necessary to consider the presence of rock pores.

[0003] The single-porosity theory, however, ignores the mesoscale heterogeneity of rock, resulting in predicted P-wave attenuation in the seismic exploration frequency band (1-100 Hz) that is far less than the observed value. The dual-porosity theory, on the other hand, accounts for the energy attenuation caused by mesoscale local flow between high-pressure and low-pressure areas of rock during seismic wave propagation, resulting in predicted P-wave attenuation that is closer to the observed results.

[0004] Currently, the application of dual-porosity media theory is mostly in two dimensions. Two-dimensional forward modeling can only handle P-waves and SV-waves, not SH-waves. Furthermore, seismic wave propagation in two-dimensional geometries differs significantly from that in three-dimensional geometries, and two-dimensional simulations cannot provide the absolute amplitude of seismic signals. Therefore, developing forward modeling methods for three-dimensional dual-porosity media is highly valuable. Furthermore, dual-porosity media can better simulate underground oil and gas reservoirs, and the rich information contained in their full seismic waveforms can have important applications in acoustic logging and other fields.

[0005] Traditional three-dimensional simulations are computationally intensive and difficult to implement. However, the global matrix method can rapidly calculate the seismic wave response in horizontally layered dual-porosity media and process three-dimensional wave fields in one-dimensional inhomogeneous media. Summary of the Invention

[0006] In view of the shortcomings of the prior art, the present invention aims to provide a forward modeling method for seismic wave fields excited by a point source in horizontally layered dual-porosity media, so as to simulate the three-dimensional seismic wave fields in horizontally layered dual-porosity media.

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

[0008] A forward modeling method for seismic wave fields excited by a point source in a horizontally layered dual-porosity medium. The specific process includes:

[0009] S1. Transform the governing equations of seismic wave propagation in dual-porosity media in the time domain into the frequency domain by Fourier transform.

[0010] S2. Introducing a set of surface-harmonic coordinate basis vectors, the frequency-domain governing equations are expanded in surface-harmonic coordinates, resulting in two sets of decoupled governing equations. One set is related only to P-waves and SV-waves, and the other is related only to SH-waves, referred to as PSV and SH modes, respectively. The introduction of surface-harmonic coordinates decomposes the original three-dimensional problem into a two-dimensional plane wave problem and transforms the vectors from the spatial domain to the wavenumber domain.

[0011] S3. 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 upgoing and downgoing wave amplitudes in each layer of porous medium as unknown parameters. Solve this system of equations to obtain wave fields such as solid phase displacement and seepage displacement.

[0012] S4. Use fast Hankel transform and fast Fourier transform to integrate the frequency and wave number, and finally obtain the time-space domain wave field.

[0013] Furthermore, the method of converting the time domain dual-porosity medium control equations into the frequency domain using Fourier transform includes:

[0014] Introducing the Fourier transform pair,

[0015] The time domain control equations are converted to the frequency domain through forward transformation, and A(r,θ,z,t) represents the time domain variables such as the solid phase displacement u r 、u θ 、u z ; represents the corresponding frequency domain variable.

[0016] Furthermore, the introduction of the surface harmonic coordinate system includes:

[0017] Introducing a set of surface harmonic coordinate basis vectors

[0018]

[0019]

[0020]

[0021] 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;

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

[0023]

[0024] in is the solid phase displacement vector u, seepage displacement and Body force vector F, or stress vector t = σ on the horizontal plane rz e r +σ θ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 components of the vector in the frequency-wavenumber domain.

[0025] Furthermore, the global matrix is ​​processed to avoid numerical overflow, and the amplitude vector in the nth layer of porous media in the PSV mode and SH mode is calculated respectively:

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

[0027]

[0028] For SH mode:

[0029]

[0030] and 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;

[0031] q Pf is the vertical slowness of the fast longitudinal wave, q Ps1 is the vertical slowness of the first slow longitudinal wave, q Ps2is the vertical slowness of the second slow longitudinal wave, q SV is the vertical slowness of the SV wave, q SH is the vertical slowness of the SH wave;

[0032] z n-1 represents the depth of the upper interface of the nth layer of porous media, z n Indicates the depth of the lower interface of the nth layer of porous media.

[0033] Furthermore, the frequency and wavenumber are integrated using the fast Hankel transform and the fast Fourier transform:

[0034] For u z , σ zz ,

[0035]

[0036] For u r , σ rz

[0037]

[0038] For u θ , σ θz

[0039]

[0040] Beneficial effects of the present invention:

[0041] 1. The forward modeling method for point-source excitation of seismic wave fields in horizontally layered dual-porosity media proposed in this invention converts the governing equations for seismic wave propagation in the time domain to the frequency domain through Fourier transform. A set of surface-harmonic coordinate basis vectors is then introduced to expand the frequency-domain governing equations in surface-harmonic coordinates and decouple them into PSV and SH modes. The frequency-wavenumber domain wavefield in the horizontally layered model is then solved using a global matrix method. Finally, the frequency and wavenumber are integrated using fast Hankel transform and fast Fourier transform to obtain the time-space domain wavefield. Compared with forward modeling methods based on single-porosity media theory, the calculated results in the low-frequency band are closer to the observed values.

[0042] 2. The forward modeling method for point source excitation seismic wave fields in horizontally layered dual-porosity media proposed in the present invention obtains the three-dimensional seismic wave field in dual-porosity media through an analytical algorithm. It has high computational efficiency, provides a basis for studying the influencing factors of longitudinal and shear wave attenuation in dual-porosity media, and provides a forward modeling operator for seismic wave inversion. BRIEF DESCRIPTION OF THE DRAWINGS

[0043] 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.

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

[0045] Figure 2 is a schematic diagram of a horizontal layered model according to an embodiment of the present invention;

[0046] Figure 3 3. This is a comparison diagram of wave field calculation results of the Biot poroelastic model and the dual-porosity medium model according to an embodiment of the present invention;

[0047] Figure 4 3 is a comparison diagram of wave fields at three observation points in a single-porous medium and a double-porous medium according to an embodiment of the present invention. DETAILED DESCRIPTION

[0048] 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.

[0049] The present application is preferably applicable to underground exploration using a global matrix method, and is characterized in that it can efficiently realize the numerical simulation of three-dimensional seismic wave fields in horizontally layered dual-porosity media.

[0050] like Figure 1 As shown in Figure 2, the specific process of the forward modeling method for the seismic wave field excited by a point source in a horizontally layered dual-porosity medium includes:

[0051] S1. Transform the governing equations of seismic wave propagation in dual-porosity media in the time domain into the frequency domain by Fourier transform.

[0052] S2. Introducing a set of surface-harmonic coordinate basis vectors, the frequency-domain governing equations are expanded in surface-harmonic coordinates, resulting in two sets of decoupled governing equations. One set is related only to P-waves and SV-waves, and the other is related only to SH-waves, referred to as PSV and SH modes, respectively. The introduction of surface-harmonic coordinates decomposes the original three-dimensional problem into a two-dimensional plane wave problem and also transforms the vectors from the spatial domain to the wavenumber domain.

[0053] S3. 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 upgoing and downgoing wave amplitudes in each layer of porous medium as unknown parameters. Solve this system of equations to obtain wave fields such as solid phase displacement and seepage displacement.

[0054] S4. Use fast Hankel transform and fast Fourier transform to integrate the frequency and wave number, and finally obtain the time-space domain wave field.

[0055] Specifically, we first establish a horizontal layering model. The Earth's interior is heterogeneous, but overall, the heterogeneity is generally greater in the vertical direction 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 horizontally layered medium model consisting of multiple parallel layers. This is a commonly used approximate geological model that is homogeneous and isotropic horizontally and piecewise uniform vertically.

[0056] Horizontal layered models such as Figure 2 As shown, 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, of which the bottom layer is a uniform semi-infinite space, z=z0 is the free surface, and above it is air. The earthquake source is placed at a depth of z=z s On the z-axis, its coordinates are (0,0,z s ). The seismic waves in each layer of the medium in the model obey the governing equations of dual-porosity media.

[0057] The governing equations for seismic wave propagation in dual-porosity media are converted into cylindrical coordinates, and the body force term F is considered in the momentum conservation equation. r , F θ , F z And by introducing the following Fourier transform pair:

[0058]

[0059]

[0060] The time domain governing equations are transformed into the frequency domain, where A(r,θ,z,t) represents the time domain variables such as the solid phase displacement u r 、u θ 、u z wait, represents the corresponding frequency domain variable.

[0061] Then introduce a set of surface harmonic coordinate basis vectors:

[0062]

[0063]

[0064]

[0065] where 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.

[0066] 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:

[0067]

[0068] in 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 e r +σ θ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 components of the vector in the frequency-wavenumber domain.

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

[0070]

[0071]

[0072]

[0073]

[0074]

[0075]

[0076] The symbol * denotes complex conjugation. Expanding all vectors in the governing equations using surface harmonic vectors yields two decoupled sets of governing equations: one related only to P- and SV-waves, and the other only to SH-waves, referred to as the PSV and SH modes, respectively. The introduction of surface harmonic vectors decomposes the original three-dimensional cylindrical coordinate problem into a two-dimensional plane wave problem, with the wave fields in both modes conforming to plane wave theory. It also transforms vectors from the spatial domain to the wavenumber domain.

[0077] The governing equations for each set of modes can be expressed as follows:

[0078]

[0079] where B is the displacement-stress vector, The superscript "V" indicates the PSV mode, and "H" indicates the SH mode. In the PSV and SH modes, the matrices are 8th and 2nd order, respectively. F is a vector related to the source, called the body force vector:

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

[0081]

[0082] Let iωΛ represent the eigenvalue matrix of A and D represent the eigenvector matrix of A, then we can get: A=DiωΛD -1 Λ is expressed as:

[0083] where q Pf is the vertical slowness of the fast longitudinal wave, q Ps1 is the vertical slowness of the first slow longitudinal wave, q Ps2 is the vertical slowness of the second slow longitudinal wave, q SV is the vertical slowness of the SV wave, q SH is the vertical slowness of the SH wave.

[0084] Once the matrix A is determined, Λ and D can be obtained through eigenvalue decomposition. Introducing the linear transformation: Β = DW, the governing equations 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. n Represents the amplitude vector in the nth layer of the medium.

[0085] Process the global matrix to avoid numerical overflow. For PSV mode, the depth z (z n <z<zn-1 ) at W n The solution has the following form:

[0086]

[0087] For SH mode:

[0088]

[0089] in and 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 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:

[0090] For PSV mode, For SH mode, in, represents the stress component, and represent the pore pressure in the background phase and the embedded body, respectively.

[0091] Combining the relationship Β = DW, the boundary condition at the free surface can be rewritten as: For PSV mode, J + =D(3:6,1:4),J - =D(3:6,5:8); for SH mode, J + =D(2,1),J - =D(2,2).

[0092] 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, B n (z n )=Β n+1 (z n ). It can be further expressed as: D n W n (z n )=D n+1 W n+1 (z n ). So there is

[0093]

[0094] in It's D n For the PSV mode,

[0095] where h n =z n -z n-1 represents the thickness of the nth layer of porous medium, z n-1 represents the depth of the upper interface of the nth layer of porous media, z n Indicates the depth of the lower interface of the nth layer of porous media.

[0096] where h n+1 =z n+1 -z n I is an identity matrix. For the SH mode,

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

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

[0099] When m=0,

[0100]

[0101]

[0102]

[0103]

[0104]

[0105]

[0106] When m=±1,

[0107]

[0108]

[0109]

[0110]

[0111]

[0112]

[0113] When m=±2,

[0114]

[0115]

[0116]

[0117]

[0118]

[0119]

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

[0121]

[0122] By solving this set of equations, we can obtain the amplitudes of the upgoing and downgoing waves in each layer of porous media, and then calculate the frequency-wavenumber domain wave field in each layer using the formula Β=DW.

[0123] In order to obtain the time-space domain wave field, the frequency and wave number are integrated using the fast Hankel transform and the fast Fourier transform. z , σ zz , have

[0124]

[0125] For u r , σ rz ,have

[0126]

[0127] For u θ , σ θz ,have

[0128]

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

[0130]

[0131] In order to verify the correctness of the method in this invention, a half-space model is designed, and the model parameters are as follows: s =2650kg / m 3 , K s =38GPa, μ s =44GPa、φ 10 =0.1,φ 20 =0.3, ν1=1, ν2=0, c1=10, c2=200, c S =10, κ1=0.01Darcy, κ2=1Darcy,ρ s =1040kg / m 3 , K f =2.5GPa, η = 0.001Pa·s, R0 = 2.1cm. At this time, the embedded volume ratio is 0, and the dual-porosity medium is degenerated into a single-porosity medium. A double-couple source is used to excite seismic waves, and its seismic moment component is M xz =M zx =1.5×10 12 In order to avoid the influence of surface reflection waves, the source is placed at a depth z far away from the surface. s =1km. The location of the observation point is (x=120m, y=1m, z=1120m). Figure 3 The solid black line in the middle shows the results of the analytical solution of the Biot poroelastic model, while the dotted black line shows the results of the degenerate dual-porosity model. As can be seen, the results of the two models are very consistent, validating the effectiveness of the method described in this paper.

[0132] The longitudinal wave attenuation predicted by the single-porous medium theory is often smaller than the observed value, especially in the seismic wave frequency band (1-100Hz). The dual-porous medium theory predicts stronger attenuation in the low-frequency band, which can better explain the observation results. To further verify the effectiveness of the method described in this article, the seismic wave responses in the single-pore model and the dual-pore model are compared. A half-space model is still used, and the parameters of the single-porous medium are the same as above. The volume ratios of the background phase and the embedded body of the dual-porous medium are ν1=0.963 and ν2=0.037, respectively, and the other medium parameters are the same as those of the single-porous medium. A double-couple source with a depth of 1000km is used to excite seismic waves, and its seismic moment component is M xz =M zx =1.5×10 15The earthquake source time function is a Ricker wavelet with a center frequency of 30 Hz and a time delay of 0.03 s. The three observation points are located on the same straight line (x = 1 km, y = 0.5 km, z = 1001 km; x = 2 km, y = 1 km, z = 1002 km; x = 3 km, y = 1.5 km, z = 1003 km).

[0133] Figure 4 The black dotted line and the black solid line are the waveforms of the radial component of the solid-phase displacement in the single-porous medium and the dual-porous medium at three observation points, respectively. It can be seen that at the first receiving point, the amplitude of the longitudinal wave in the dual-porous medium is only one-tenth of that in the single-porous medium. At the third receiving point, the longitudinal wave is almost invisible. However, the amplitude of the shear wave in the two media is basically the same. This is because the mesoscopic local flow in the dual-porous medium affects the attenuation of seismic waves. The propagation of the longitudinal wave causes changes in the pore pressure, resulting in the generation of local flow, thereby dissipating the energy of the longitudinal wave. However, the shear wave does not cause changes in the pore pressure during propagation, and therefore does not cause the generation of local flow. This result is consistent with the laws of physics and proves the correctness of the method described in the present invention.

[0134] This application uses a full three-dimensional numerical solver to solve the three-dimensional wave equations in dual-porosity media, which is computationally intensive and difficult to implement. However, this method directly solves the governing equations for dual-porosity media in a horizontally layered model to obtain a three-dimensional seismic wave field, demonstrating high computational efficiency. Seismic wave propagation in two-dimensional geometries differs significantly from that in three-dimensional geometries, and two-dimensional simulations cannot provide the absolute amplitude of seismic signals. This provides a foundation for studying the factors influencing the attenuation of longitudinal and shear waves in dual-porosity media and offers a forward operator for seismic wave inversion.

[0135] 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.

[0136] 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 method for forward modeling of seismic wave fields excited by a point source in a horizontally layered dual-porosity medium, characterized by: The specific process includes: S1. Transform the governing equations of seismic wave propagation in dual-porosity media in the time domain into the frequency domain by Fourier transform. S2. Introducing a set of surface-harmonic coordinate basis vectors, the frequency-domain governing equations are expanded in surface-harmonic coordinates, resulting in two sets of decoupled governing equations. One set is related only to P-waves and SV-waves, and the other is related only to SH-waves, referred to as PSV and SH modes, respectively. The introduction of surface-harmonic coordinates decomposes the original three-dimensional problem into a two-dimensional plane wave problem and transforms the vectors from the spatial domain to the wavenumber domain. S3. 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 upgoing and downgoing wave amplitudes in each layer of porous medium as unknown parameters. Solve this system of equations to obtain the solid phase displacement and seepage displacement wave fields. S4. Use fast Hankel transform and fast Fourier transform to integrate the frequency and wave number, and finally obtain the time-space domain wave field; Introducing 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, ..., is the first kind of m-order Bessel function, k is the horizontal wave number, in the frequency domain, any vector in the cylindrical coordinate system Expands into the following form: in is the solid phase displacement vector u, seepage displacement and Body force vector F, or stress vector t = σ on the horizontal plane rz e r +σ θz e θ +σ zz e z The summation range m in the formula depends on the symmetry of the source in the annular direction. For explosive point sources and 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 the vector in the frequency-wavenumber domain; The global matrix is ​​processed to avoid numerical overflow, and the amplitude vector in the nth layer of porous media in the PSV mode and SH mode are calculated respectively: For PSV mode, depth z, Z n <Z<Z n-1 W n The solution has the following form: For SH mode: and represents the amplitude of each downgoing wave at the upper interface of the nth layer of porous media, and represents the amplitude of each upgoing wave at the lower interface of the nth layer of porous medium; q Pf is the vertical slowness of the fast longitudinal wave, q Ps1 is the vertical slowness of the first slow longitudinal wave, q Ps2 is the vertical slowness of the second slow longitudinal wave, q SV is the vertical slowness of the SV wave, q SH is the vertical slowness of the SH wave; z n-1 represents the depth of the upper interface of the nth layer of porous media, z n Indicates the depth of the lower interface of the nth layer of porous media.

2. The forward modeling method for point source excitation seismic wavefield in horizontally layered dual-porosity media according to claim 1, characterized in that: The Fourier transform is used to convert the governing equations of seismic wave propagation in dual-porosity media in the time domain to the frequency domain, including: Introducing the Fourier transform pair, The time domain control equations are converted to the frequency domain through forward transformation, and A(r,θ,z,t) represents the time domain variables such as the solid phase displacement u r 、u θ 、u z ; represents the corresponding frequency domain variable.

3. The forward modeling method for point source excitation seismic wavefield in horizontally layered dual-porosity media according to claim 1, characterized in that: Integrate over frequency and wavenumber using the fast Hankel transform and the fast Fourier transform: For u z , σ zz , For u r , σ rz For u θ , σ θz

Images