A forward modeling method and device for seismic wave data in a viscoelastic medium

Abstract

An embodiment of this specification discloses a forward simulation method and device for seismic wave data in viscoelastic media. By obtaining an initial parameter field; meshing the initial parameter field to generate a discrete initial parameter field; and using a preset viscoelastic medium wavefield propagation operator to calculate velocity components V x and V z and stress components, the viscoelastic medium wavefield propagation operator includes a naturally decoupled amplitude attenuation term and a phase dispersion term, which can achieve separate simulation or simultaneous simulation of amplitude attenuation and phase dispersion. Because when performing attenuation compensation reverse time migration imaging, only the sign of the amplitude attenuation term needs to be changed while keeping the sign of the phase dispersion unchanged. Therefore, compared with the traditional viscoelastic seismic wave forward simulation in which the amplitude attenuation and phase dispersion based on the standard linear solid model are coupled together, it can efficiently simulate seismic wave data with amplitude attenuation and phase dispersion characteristics in viscoelastic media.

Description

Technical Field

[0001] This specification relates to the field of petroleum geophysical exploration, and particularly to a method and device for forward modeling of seismic wave data in viscoelastic media. Background Art

[0002] In seismic exploration and seismology, the properties and occurrence states of underground media can be obtained through seismic wave forward modeling, migration imaging, and velocity inversion techniques.

[0003] Generally, when seismic waves propagate in fluid-containing media, they exhibit characteristics of amplitude attenuation and phase dispersion. If the influence of these viscosities of underground media on the propagation of seismic waves is ignored, it will cause deviations in the results of velocity inversion and migration imaging. In the past, the amplitude attenuation term and phase dispersion term of the traditional viscoelastic wave equation were coupled together, which was not conducive to the implementation of velocity inversion and migration imaging algorithms.

[0004] Therefore, to accurately obtain the true properties of underground rocks and the true shape of structural interfaces, developing a viscoelastic wave equation that can efficiently simulate and decouple the amplitude attenuation and phase dispersion characteristics is the key to solving complex underground media imaging. Summary of the Invention

[0005] Embodiments of this specification provide a method and device for forward modeling of seismic wave data in viscoelastic media to solve the following technical problems: A forward modeling scheme for seismic wave data that can efficiently simulate the amplitude attenuation and phase dispersion characteristics in viscoelastic media is required.

[0006] To solve the above technical problems, one or more embodiments of this specification are implemented as follows:

[0007] In a first aspect, embodiments of this specification provide a method for forward modeling of seismic wave data in viscoelastic media, including:

[0008] Obtain an initial parameter field, where the initial parameter field includes v p 、v s 、Q p 、Q s and ρ, where v p represents the longitudinal wave velocity of the medium, v s represents the shear wave velocity of the medium, Q p is the P-wave quality factor, and Q s is the S-wave quality factor;

[0009] Perform grid discretization on the initial parameter field to generate a discrete initial parameter field;

[0010] Use a preset wave field propagation operator for viscoelastic media to calculate velocity components V x and V zAnd stress components, the viscoelastic medium wave field propagation operator is in the form of:

[0011]

[0012] Wherein,

[0013] γ p = arctan(1 / Q p ) / π, γ s = arctan(1 / Q s ) / π;

[0014] P is the discrete form of the velocity component V x and V x is the velocity field in the time domain representing the horizontal direction, G is the discrete form of the velocity component V z and V z is the velocity field in the time domain in the vertical direction, U is the discrete form of the stress component σ xx and V is the discrete form of the stress component σ xz and W is the discrete form of the stress component σ zz where i and j respectively represent the horizontal and vertical spatial discrete indices, m represents the time discrete index, Δx represents the transverse grid spacing, Δz represents the longitudinal grid spacing, Δt is the time sampling interval, represents the regular grid difference coefficient, N represents the difference order, FFT represents the fast Fourier transform, and FFT -1 represents the inverse transform of the fast Fourier transform.

[0015] In a second aspect, one or more embodiments of this specification provide an electronic device, including:

[0016] At least one processor; and,

[0017] A memory communicatively connected to the at least one processor; wherein,

[0018] The memory stores instructions executable by the at least one processor, and the instructions are executed by the at least one processor so that the at least one processor can execute the method as described in the first aspect.

[0019] One or more embodiments of this specification adopting the above at least one technical solution can achieve the following beneficial effects: By obtaining an initial parameter field, the initial parameter field includes v p , v s , Q p , Q s and ρ, wherein v p represents the longitudinal wave velocity of the medium, v s represents the shear wave velocity of the medium, Qp is the P-wave quality factor, Q s is the S-wave quality factor; the initial parameter field is meshed to generate a discrete initial parameter field; a preset viscoelastic medium wave field propagation operator is used to calculate the velocity components V x and V z and stress components. The viscoelastic medium wave field propagation operator contains a naturally decoupled amplitude attenuation term and phase dispersion term, which can realize separate simulation or simultaneous simulation of amplitude attenuation and phase dispersion. Because when doing attenuation compensation reverse time migration imaging, only the sign of the amplitude attenuation term needs to be changed while keeping the sign of the phase dispersion unchanged. Therefore, compared with the traditional viscoelastic seismic wave forward simulation in which the amplitude attenuation and phase dispersion are coupled based on the standard linear solid model, it can efficiently simulate seismic wave data with amplitude attenuation and phase dispersion characteristics in the viscoelastic medium. BRIEF DESCRIPTION OF THE DRAWINGS

[0020] In order to more clearly illustrate the technical solutions in the embodiments of this specification or the prior art, the following will briefly introduce the drawings required for use in the description of the embodiments or the prior art. Obviously, the drawings in the following description are only some embodiments recorded in this specification. For those of ordinary skill in the art, other drawings can be obtained based on these drawings without creative efforts.

[0021] Figure 1 is a schematic flow chart of a method for forward simulation of seismic wave data in a viscoelastic medium provided by an embodiment of this specification;

[0022] Figure 2 is a schematic diagram of decoupled amplitude attenuation and phase dispersion characteristic wave field simulation in a viscoelastic medium provided by an embodiment of this specification, where (a) is the horizontal component and (b) is the vertical component.

[0023] Figure 3 is a schematic diagram of the relationship between the viscoelastic parameters and spatial positions of a geological seismic model provided by an embodiment of this specification; (a) is the P-wave velocity, (b) is the P-wave quality factor, and (c) is the density.

[0024] Figure 4 is a schematic diagram of the comparison between the wave field simulated by the method provided by an embodiment of this specification and the wave field simulated by the elastic wave simulation scheme;

[0025] Figure 5 is a schematic diagram of the longitudinal single-trace comparison between the wave field simulated by the method provided by an embodiment of this specification and the wave field simulated by the elastic wave equation. DETAILED DESCRIPTION OF THE EMBODIMENTS

[0026] Embodiments of this specification provide a method and device for forward simulation of seismic wave data in a viscoelastic medium.

[0027] To enable those skilled in the art to better understand the technical solutions in this specification, the technical solutions in the embodiments of this specification will be clearly and completely described below in conjunction with the accompanying drawings in the embodiments of this specification. Obviously, the described embodiments are only a part of the embodiments of this application, rather than all of the embodiments. Based on the embodiments of this specification, all other embodiments obtained by those of ordinary skill in the art without creative efforts shall fall within the protection scope of this application.

[0028] First, it is necessary to explain the viscoelastic wave propagation operator used in the simulation process of this application.

[0029] The known complex value expressions of the longitudinal wave velocity and the transverse wave velocity in the underground viscoelastic medium are:

[0030]

[0031] It is also known that the wavenumber-frequency domain viscoelastic wave equation in seismology can be expressed as:

[0032]

[0033] Among them, z and x respectively represent the longitudinal and transverse spatial coordinates, ω is the angular frequency, ω0 is the reference angular frequency, v p and v s are the longitudinal wave and transverse wave velocities of the medium, γ p = arctan(1 / Q p ) / π, γ s = arctan(1 / Q s ) / π is a dimensionless quantity, and its value is between 0 and 1. Q p and Q s are the quality factors of the P wave and the S wave respectively, representing the amount of energy attenuation. The larger Q p and Q s , the weaker the attenuation. i represents the imaginary unit, and respectively represent the frequency domain wave fields in the horizontal and vertical directions, s x and s z are the frequency domain source terms in the horizontal and vertical directions, ρ represents the medium density, k x and k zrespectively represent the wavenumbers in the horizontal and vertical directions. Solving formula (2) in the frequency domain requires a huge amount of computational effort and memory. Therefore, we consider transforming it into the time domain and then solving it in the time domain. After the traditional transformation into the time domain, it is a time-fractional wave equation, and its numerical solution requires storing the wave field values at all past times, which also requires a huge amount of memory. To obtain the space-fractional viscoelastic wave equation and reduce the computational effort and memory, we propose the following approximation:

[0034]

[0035] where the real part terms in equations (5) and (6) above are used to simulate the phase dispersion characteristics, and the imaginary part terms simulate the amplitude attenuation characteristics, and k represents the spatial wavenumber.

[0036] Substitute equations (5) and (6) into equations (3) and (4), and transform them into the time-space domain, the first-order velocity-stress form viscoelastic wave equation in the time domain can be obtained as

[0037]

[0038]

[0039] where V x and V z respectively represent the velocity fields in the horizontal and vertical directions in the time domain, σ xx , σ xz , σ zz represent the stress wave fields, represents the Laplace operator, f x and f z respectively represent the source terms in the horizontal and vertical directions in the time domain, The terms related to the coefficients α p and α s can describe the phase dispersion characteristics during the propagation of seismic waves and are called dispersion terms; the terms related to the coefficients β p and β s can describe the amplitude attenuation characteristics during the propagation of seismic waves and are called attenuation terms. Therefore, in the wave equations (7) and (8), the amplitude attenuation term and the phase dispersion term are decoupled. Therefore, when assuming that the coefficient of the phase dispersion term is 0, the amplitude attenuation characteristics can be simulated. When assuming that the coefficient of the amplitude attenuation term is 0, the phase dispersion characteristics can be simulated.

[0040] Since there is a spatially variable fractional Laplace operator in equations (9), (10), and (11) above. Therefore, when using the pseudospectral method for numerical solution, it is necessary to decouple the wavenumber and the fractional order to avoid generating numerical simulation noise. The decoupling of the wavenumber and the fractional order can use the following relationship

[0041] k2γ(x) ≈1 + 2ln(k)γ(x) + 2(ln(k)) 2 γ 2 (x), (k > 0), (12)

[0042] k γ(x) ≈1 + ln(k)γ(x) + (ln(k)) 2 γ 2 (x) / 2, (k > 0). (13)

[0043] Substituting (12) and (13) into (9) - (11) gives

[0044]

[0045]

[0046] In this case, the finite difference - pseudo - spectral hybrid numerical simulation method is adopted (i.e., formulas (7) and (8) use finite differences, and formulas (14) - (15) use the pseudo - spectral method) for numerical solution. Therefore, the visco - elastic medium wave - field propagation operator can be expressed as:

[0047]

[0048] where P, G, U, V, W respectively represent the numerical simulation representation forms of velocity components V x and V z and stress components σ xx , σ xz , σ zz , i and j respectively represent the horizontal and vertical spatial discretization indices, m represents the time discretization index, Δx represents the horizontal grid spacing, Δz represents the vertical grid spacing, Δt is the time sampling interval, represents the regular grid difference coefficient, N represents the difference order. When N = 5, it is the 10 - order spatial difference accuracy. FFT represents the fast Fourier transform, and FFT -1 represents the inverse transform of the fast Fourier transform.

[0049] The foregoing part explains and illustrates the wave - field propagation operator adopted in the embodiments of this specification. The specific usage is as Figure 1 shown. In the first aspect, the embodiments of this specification provide a forward seismic data simulation method in a visco - elastic medium, Figure 1 which is a schematic flow diagram of a forward seismic data simulation method in a visco - elastic medium provided by the embodiments of this specification, including:

[0050] S101, obtaining an initial parameter field, where the initial parameter field includes v p , v s , Q p , Qs and ρ, where v p represents the longitudinal wave velocity of the medium, v s represents the shear wave velocity of the medium, Q p is the P-wave quality factor, Q s is the S-wave quality factor;

[0051] Specifically, data such as field geological reconnaissance, logging, geophysical exploration, and empirical models can be used to draw the geological and geophysical models of the study area, and the initial parameter field can be obtained by filling the velocity values into them.

[0052] S103, mesh the initial parameter field to generate a discrete initial parameter field;

[0053] The grid is divided into small rectangles or irregular shapes such as triangles. Irregular shapes can simulate undulating interfaces, while regular grids have better adaptability to continuous interfaces with slow slope changes. In practical applications, regular rectangular grids are mostly used because they are simple to mesh, have a small computational amount, and can meet most production requirements. Whether using the finite difference method or the pseudo-spectral method for numerical simulation, the selection of the grid spacing and the time sampling interval must satisfy the limitations of the stability conditions.

[0054] where, Δd max is the maximum value of Δx and Δz. V max is the maximum longitudinal wave velocity value.

[0055] At the same time, to avoid numerical dispersion, in areas with small velocities, small grid spacings and time sampling intervals need to be selected. Therefore, a variable grid simulation method is also considered to improve the numerical simulation calculation efficiency.

[0056] S105, using a preset viscoelastic medium wave field propagation operator, calculate the velocity components V x and V z and the stress components according to the discrete initial parameter field. The form of the viscoelastic medium wave field propagation operator is as shown in equation (17):

[0057]

[0058] The wave field propagation operator has been described in detail above. In practical applications, it can be programmed and set in the form of functional modules or algorithm modules. The positions of the shot point and the geophone can be set at any position in the grid to achieve different observation purposes. In seismic exploration, generally, both the shot point and the geophone are set on the first layer of the grid (i.e., the surface). In seismology, the shot point is generally set deep underground, and the geophones are arranged on the surface.

[0059] As described above, in the numerical simulation process, for the difference calculation part, a second-order time and tenth-order space difference scheme is applied. For the pseudo-spectral method calculation part, a second-order difference is used for time, and the space is solved in the wavenumber domain (theoretically, the space can reach infinite-order difference accuracy). Then, based on the numerical simulation of the newly derived viscoelastic medium wave equation, data on the stress field component or velocity field component of the wave field propagation in the viscoelastic medium can be obtained.

[0060] Since in Equation (17), it contains a naturally decoupled amplitude attenuation term (the term related to coefficients β p and β s ), and a phase dispersion term (the term related to coefficients α p and α s ), therefore, the calculated velocity components V x and V z and the stress components are also naturally decoupled in terms of amplitude attenuation and dispersion.

[0061] Figure 2 This is a schematic comparison diagram of the wave field simulation results of the horizontal component ( Figure 2 (a)) and vertical component ( Figure 2 (b)) provided in this specification, which do not contain phase dispersion and do not contain amplitude attenuation (i.e., elastic waves), only contain amplitude attenuation, only contain phase dispersion, and contain both amplitude attenuation and phase dispersion. Among them, Quadrant Ⅰ is the wave field simulation result without attenuation and without dispersion (i.e., the elastic wave simulation result), Quadrant Ⅱ is the wave field simulation result only containing phase dispersion, Quadrant Ⅲ is the wave field simulation result only containing attenuation, and Quadrant Ⅳ is the wave field simulation result containing both attenuation and dispersion. Figure 2 It shows that the viscoelastic wave equation proposed by the present invention can simulate the wave field with amplitude attenuation characteristics or the wave field with phase dispersion characteristics respectively, and can also simulate the wave field with both amplitude attenuation and phase dispersion characteristics simultaneously. In the figure, the abscissa is the length x, and the ordinate is the depth z.

[0062] For the solution provided in the embodiments of this specification, its equation has a naturally decoupled amplitude attenuation term and phase dispersion term, and can achieve separate simulation or simultaneous simulation of amplitude attenuation and phase dispersion. Because when performing attenuation compensation reverse time migration imaging, we only need to change the sign of the amplitude attenuation term while keeping the sign of the phase dispersion unchanged. Therefore, compared with the traditional viscoelastic seismic wave forward simulation in which the amplitude attenuation and phase dispersion are coupled based on the standard linear solid model, it can efficiently simulate the seismic wave data with amplitude attenuation and phase dispersion characteristics in the viscoelastic medium.

[0063] Numerical simulation tests and practical applications also prove that the viscoelastic medium seismic wave forward simulation proposed by this solution can be well applied to seismic exploration and seismology research. At the same time, compared with the traditional time fractional wave equation, the space fractional wave equation provided in this embodiment does not need to save the wave field at all times during numerical solution, and the calculation is simple and the amount of calculation is small.

[0064] To verify the effectiveness and robustness of the viscoelastic wavefield forward simulation method of this patent in simulating an actual medium model, we constructed a complex gas cloud model and applied the solution proposed in the present invention to seismic wave simulation therein. The velocity parameters of this model are as Figure 3 shown. There is a gas reservoir in the middle of the model, and this gas reservoir has a relatively small quality factor value. When seismic waves pass through it, strong amplitude attenuation and phase dispersion characteristics will occur.

[0065] Figure 3 It is a schematic diagram showing the relationship between the velocity parameter values and spatial positions in a viscoelastic medium provided in the embodiments of this specification.

[0066] Among them, the shear wave velocity can be given according to the relationship v s = v p / 1.732. The shear wave quality factor can be given according to the relationship Q s = Q p / 1.1. In the figure, the abscissa is the length x, and the ordinate is the depth z.

[0067] Here, the results of elastic wave equation simulation are selected for reference. Figure 4 It is a schematic diagram comparing the horizontal component (i.e., the (a) part and (c) part in Figure 4 ) and the vertical component (i.e., the (b) part and (d) part in Figure 4 ) of a certain shot record obtained by the method provided in the embodiments of this specification and the elastic medium wave equation. In the figure, the abscissa is the length x, and the ordinate is the time t. Figure 4 The (a) part and (c) part in Figure 4 are the wavefields simulated by this method, and the (b) part and (d) part in Figure 5 are the wavefields simulated by elastic waves. Figure 4 It is a schematic diagram comparing the waveforms of the horizontal component (a) and the vertical component (b) of a certain vertical single-trace record in

[0068] . It can be seen that compared with the elastic wavefield, the wavefield simulated by the solution proposed in this patent has a reduced amplitude and a phase dispersion. Through the test of this gas cloud model, the stability and effectiveness of the seismic wave forward simulation method in viscoelastic media proposed in this patent are proved.

[0069] Each embodiment in this specification is described in a progressive manner. For the same or similar parts among the embodiments, reference can be made to each other. Each embodiment focuses on the differences from other embodiments. In particular, for the embodiments of devices, equipment, and media, since they are basically similar to the method embodiments, the description is relatively simple. For the relevant parts, reference can be made to the partial description of the method embodiments, and details will not be repeated here.

Claims

1. A forward simulation method for seismic wave data in a viscoelastic medium, comprising: Obtain an initial parameter field, where the initial parameter field includes v p , v s , Q p , Q s and ρ, where v p represents the longitudinal wave velocity of the medium, v s represents the shear wave velocity of the medium, Q p is the P-wave quality factor, Q s is the S-wave quality factor; Meshing the initial parameter field to generate a discrete initial parameter field; Using a preset viscoelastic medium wavefield propagation operator, the velocity components V x , V z and stress components are calculated according to the discrete initial parameter field, and the form of the viscoelastic medium wavefield propagation operator is: Among them, γ p = arctan(1 / Q p ) / π, γ s = arctan(1 / Q s ) / π; P is the discrete form of the velocity component V x , where V x is the velocity field in the time domain representing the horizontal direction, G is the discrete form of the velocity component V z , and V z is the velocity field in the time domain representing the vertical direction, U is the discrete form of the stress component σ xx , V is the discrete form of the stress component σ xz , W is the discrete form of the stress component σ zz , i and j respectively represent the horizontal and vertical spatial discrete indices, m represents the time discrete index, Δx represents the lateral grid spacing, Δz represents the longitudinal grid spacing, Δt is the time sampling interval, represents the regular grid difference coefficient, N represents the difference order, FFT represents the fast Fourier transform, and FFT -1 represents the inverse transform of the fast Fourier transform.

2. The method according to claim 1, wherein Meshing the initial parameter field includes: Meshing the initial parameter field with irregular-shaped meshes for simulating an undulating interface; Alternatively, meshing the initial parameter field with regular-shaped meshes for simulating a continuous interface with a slowly changing slope.

3. The method according to claim 1, wherein Meshing the initial parameter field includes: Determine the grid spacings Δx, Δz and the time sampling interval Δt of the grid according to the stability conditions, where the stability conditions include: where Δd max is the maximum value of Δx and Δz, and V max is the maximum P-wave velocity value.

4. An electronic device, comprising: At least one processor; And A memory communicatively connected to the at least one processor; wherein The memory stores instructions executable by the at least one processor, and the instructions are executed by the at least one processor to enable the at least one processor to execute the method according to any one of claims 1 to 3.

Images