阅读说明：本技术 一种粘声各向异性介质中的正演成像方法 (Forward imaging method in visco-acoustic anisotropic medium ) 是由 慕鑫茹 黄建平 李振春 庄苏斌 彭炜颋 于 2021-09-29 设计创作，主要内容包括：本说明书实施例公开了一种粘声各向异性介质中的正演成像方法。发明通过波动方程系数描述了传播过程中的相位频散或者振幅衰减影响的性质,可很方便的实现振幅补偿逆时偏移成像。相比传统的波动方程,具有相同的数值模拟精度,且在极其复杂的介质中仍能稳定模拟波场传播。此外,新方程具有解耦的振幅衰减项与相位频散项,相比传统波动方程,在实现逆时偏移成像时更加稳定,可很好地应用于地震勘探。(The embodiment of the specification discloses a forward imaging method in a visco-acoustic anisotropic medium. The invention describes the property of phase dispersion or amplitude attenuation influence in the propagation process through the wave equation coefficient, and can conveniently realize amplitude compensation reverse time migration imaging. Compared with the traditional wave equation, the method has the same numerical simulation precision, and can still stably simulate wave field propagation in an extremely complex medium. In addition, the new equation has decoupled amplitude attenuation terms and phase frequency dispersion terms, and compared with the traditional wave equation, the new equation is more stable in reverse time migration imaging, and can be well applied to seismic exploration.)

1. A method of forward imaging in a visco-acoustic anisotropic medium, comprising:

acquiring an initial parameter field, wherein parameters in the initial parameter field comprise viscosity parameters and anisotropy parameters, and the viscosity parameters comprise v_p、Q、ε_QAnd delta_QWherein v is_pDenotes the propagation velocity of longitudinal waves in the medium, Q denotes a quality factor, ε_QAnd delta_QRepresenting pairs of attenuation anisotropy parameters including ε, δ, φ, wherein ε, δ characterize the velocity anisotropyThe strength phi is an anisotropic dip angle parameter

Mesh generation is carried out on the initial parameter field;

determining, for each grid point, the wave equation coefficient a related to the viscosity parameter and the anisotropy parameter₁、a₂、a₃、b₁、b₂And b₃Wherein, in the step (A),

wherein the content of the first and second substances,γ_ij＝arctan(1/Q_ij)/π；

calculating according to the wave equation coefficient and the initial parameter field by adopting a preset wave field propagation operator to generate a longitudinal wave field value at each grid point, wherein the wave field propagation operator is in the form of:

where p represents the longitudinal wavefield value, t represents the wavefield travel time,representing the laplacian, f representing the seismic source term, x being the abscissa, and z being the ordinate.

2. The method of claim 1, meshing the parameter field, comprising:

the distance delta d between the grids obtained by dividing meets the condition:wherein v is_{p max}Representing the maximum longitudinal wave velocity and at representing the time sampling interval.

3. A computer device comprising a memory, a processor and a computer program stored on the memory and executable on the processor, wherein the processor implements the method of any of claims 1 to 2 when executing the program.

Technical Field

The specification relates to the field of exploration geophysics, in particular to a forward imaging method in a visco-acoustic anisotropic medium.

Background

The subterranean medium typically exhibits viscosity due to the internal friction between the actual earth medium rocks and the fluid filling the rocks. Seismic waves exhibit amplitude attenuation and phase dispersion characteristics when propagating in a viscous medium. If these viscosity influences are not considered in seismic imaging, the amplitude of the obtained imaging result is weak, and the continuity of the seismic wave in-phase axis is poor. In addition, the continuously sedimentary formations generally exhibit anisotropic properties, and seismic waves exhibit a property of faster propagation speed in a certain direction when propagating in an anisotropic medium. If these anisotropic characteristics are not corrected during offset imaging, the offset imaging result may cause imaging position inaccuracy and energy dispersion, thereby generating imaging noise. In many formations where both properties are always present simultaneously, conventional viscoacoustic anisotropic wave equations exhibit instability when simulating wave propagation in complex media.

Therefore, an imaging method that can simultaneously consider viscosity and anisotropy characteristics and stably simulate seismic wave propagation in a complex medium is needed.

Disclosure of Invention

The invention aims to provide an imaging method capable of stably simulating seismic wave propagation in a complex medium.

In order to solve the technical problems, the invention adopts the following technical scheme:

a method of forward imaging in a visco-acoustic anisotropic medium, comprising:

acquiring an initial parameter field, wherein parameters in the initial parameter field comprise viscosity parameters and anisotropy parameters, and the viscosity parameters comprise v_p、Q、ε_QAnd delta_QWherein v is_pDenotes the propagation velocity of longitudinal waves in the medium, Q denotes a quality factor, ε_QAnd delta_QRepresents a pair of attenuation anisotropy parameters including epsilon, delta and phi, wherein epsilon and delta represent the strength of the velocity anisotropy, and phi is an anisotropic dip angle parameter

Mesh generation is carried out on the parameter field;

determining, for each grid point, the wave equation coefficient a related to the viscosity parameter and the anisotropy parameter₁、a₂、a₃、b₁、b₂And b₃Wherein, in the step (A),

wherein the content of the first and second substances,γ_ij＝arctan(1/Q_ij)/π；

calculating according to the wave equation coefficient and the initial parameter field by adopting a preset wave field propagation operator to generate a longitudinal wave field value at each grid point, wherein the wave field propagation operator is in the form of:

where p represents the longitudinal wavefield value, t represents the wavefield travel time,representing the laplacian, f representing the seismic source term, x being the abscissa, and z being the ordinate.

The embodiment of the specification adopts at least one technical scheme which can achieve the following beneficial effects:

compared with the prior art, the method describes the property of phase dispersion or amplitude attenuation influence in the propagation process through the wave equation coefficient, and can conveniently realize amplitude compensation reverse time migration imaging. Compared with the traditional wave equation, the method has the same numerical simulation precision, and can still stably simulate wave field propagation in an extremely complex medium. In addition, the new equation has decoupled amplitude attenuation terms and phase frequency dispersion terms, and compared with the traditional wave equation, the new equation is more stable in reverse time migration imaging, and can be well applied to seismic exploration.

Drawings

Fig. 1 is a schematic flow chart of a forward imaging method in a viscoelastic anisotropic medium according to an embodiment of the present disclosure;

FIG. 2 is a schematic diagram of an amplitude attenuation and phase dispersion wavefield simulation with decoupling in a viscoelastic-acoustic anisotropic medium provided by an embodiment of the present disclosure;

FIG. 3 is a schematic diagram of the relationship between the viscosity and anisotropy parameters and the position of a model provided in an embodiment of the present disclosure;

FIG. 4 is a diagram illustrating a comparison between a wave field snapshot obtained by a conventional method and a method provided by an embodiment of the present disclosure;

FIG. 5 is a schematic diagram of the relationship between the viscosity and anisotropy parameters and the position of another model provided in the embodiments of the present disclosure;

FIG. 6 is a diagram illustrating a comparison between a wave field snapshot obtained by the method provided by the embodiments of the present specification and a wave field snapshot obtained by a conventional method.

Detailed Description

In order to make the objects, technical solutions and advantages of the present application more apparent, the technical solutions of the present application will be described in detail and completely with reference to the following specific embodiments of the present application and the accompanying drawings. It should be apparent that the described embodiments are only some of the embodiments of the present application, and not all of the embodiments. All other embodiments obtained by a person skilled in the art based on the embodiments in the present specification without any inventive step are within the scope of the present application.

The wavefield propagation operators used in the simulation process of the present application will first be described.

The complex-valued exact phase velocity formula in known visco-acoustic VTI media (VTI stands for special anisotropic media without dip anisotropy effects) can be expressed as:

wherein, V_pRepresenting the complex-valued phase velocity of the qP wave, theta representing the phase angle, M_i,jRepresents a complex-valued stiffness coefficient, which can be expressed as:

where ω denotes the angular frequency, ω₀Denotes the reference angular frequency, i being the imaginary unit, γ_ij＝arctan(1/Q_ij) And/pi is a dimensionless quantity. For any positive value Q_ijIn the presence of 0 < gamma_ij< 0.5, quality factor Q_ijCan be expressed as:

Q₃₃＝Q_p, (5)

wherein epsilon_QAnd delta_QThomsen parameter pairs, C, representing attenuation correlations_ijIs the elastic stiffness coefficient, which can be expressed as:

wherein v is_pDenotes the longitudinal wave velocity along the direction of the symmetry axis, and ε and δ denote Thomsen parameters. Because the root number term exists in the formula (1), the derivation cannot be directly carried out by using the formula (1), and then an explicit wave equation expression is obtained. Applying a first order Taylor approximation to equation (1)The root term can be removed and Z represents a complex function, so that:

because of the existence of the relation sin theta ═ V_p(θ)k_xω and cos θ ═ V_p(θ)k_zω, where k_xAnd k_zFor the spatial wave numbers along the horizontal direction and the vertical direction, the viscous sound anisotropic wave equation in the spatial wave number domain can be obtained through a plurality of mathematical operations as follows:

wherein the content of the first and second substances,

equation (10) is a fluctuation equation of the spatial wave number domain of the viscoelastic sound VTI medium, and by introducing an anisotropic dip angle phi into equation (10), the fluctuation equation of the spatial wave number domain of the viscoelastic sound TTI medium (the TTI medium represents an anisotropic medium containing dip angle anisotropic influence) can be obtained as follows:

transforming equation (11) to the time-space domain can obtain the following time-domain sticky-sound TTI medium wave equation:

where p represents the longitudinal wavefield value, t represents the wavefield travel time,representing the laplacian, f the seismic source term. Note the sum coefficient a₁，a₂，a₃The related wave equation item can describe the viscosity of the seismic wavesInfluence of phase dispersion generated during propagation in a medium, and coefficient b₁，b₂，b₃The related wave equation terms can describe the amplitude attenuation effect generated in the process of propagating the seismic waves in the viscous medium. When we assume b₁，b₂，b₃When the values of (a) and (b) are all 0, the mode (12) describes the phase dispersion characteristic in the wave propagation process. When we assume gamma₁₁，γ₁₃，γ₃₃When the values of (a) and (b) are all 0, the mode (12) describes the amplitude attenuation characteristic in the wave propagation process. It is because of the nature of wave equation (12) that can separately describe the effects of phase dispersion or amplitude attenuation that can be conveniently implemented in amplitude compensated reverse time shift imaging.

In simulating the wave field propagation using equation (12), to facilitate numerical dispersion using a computer, we equate equation (12) to the following equation:

wherein the content of the first and second substances,

q₁＝FFT^-1(ln(k)FFT(p)),q₂＝FFT^-1((ln(k))²FFT(p)),

it should be noted that these q variables are merely introduced quantities for the calculation of the auxiliary equation, and have no practical physical significance.

At this time, the equation (13) is discretized by using a time second-order and space high-order difference format, so that the seismic wave propagation operator in the visco-acoustic anisotropic medium is obtained as follows:

wherein the content of the first and second substances,

wherein p is a seismic wave stress field value, q is an auxiliary seismic wave stress field value, i and j respectively represent the positions of transverse and longitudinal grid points, k represents time dispersion, Δ x represents the transverse spacing of a discrete grid, Δ z represents the longitudinal spacing of the discrete grid, Δ t represents the time sampling spacing of differential dispersion, a₀And a_nnDifferential coefficients representing finite difference dispersion, FFT representing fast Fourier transform operator, FFT^-1Inverse transform operator, k, representing a fast Fourier transform_xAnd k_zRepresenting wave numbers along the transverse direction and the longitudinal direction, respectively, and N represents the difference order of finite difference dispersion. m, N, u, mx, nz, mz, ux, uz represent intermediate variables in numerical calculation of a defined finite difference method, which have no practical physical meaning and represent spatial difference accuracy of order 10 when the value of N takes 5.

The foregoing section explains and explains the visco-acoustic anisotropic seismic wave field propagation operator adopted in the embodiment of the present specification, and a specific use manner is as shown in fig. 1, where fig. 1 is a schematic flow chart of a seismic wave data forward modeling method in a visco-acoustic anisotropic medium provided by the embodiment of the present specification, and the method includes:

s101, acquiring an initial parameter field, wherein parameters in the initial parameter field comprise viscosity parameters and anisotropic parameters. Viscosity parameters include v_pAnd Q, ε_Q，δ_QWherein v is_pRepresenting the propagation velocity of longitudinal waves in a medium, Q representing a quality factor, the size of Q determining the energy attenuation strength of seismic waves, epsilon_Q，δ_QRepresenting attenuation anisotropy parameter pairs for representing the strength of attenuation anisotropy, wherein the anisotropy parameters comprise epsilon, delta and phi, the epsilon and the delta are Thomsen anisotropy parameter pairs for representing the strength of velocity anisotropy, and the phi is an anisotropy dip angle parameter and describes the directional characteristic of the velocity anisotropy;

specifically, a geological geophysical model is drawn according to data obtained by methods of field geological investigation, mechanical drilling analysis, geophysical inversion, well logging detection and the like, and then the model is filled with viscosity and anisotropy parameter values to obtain a required parameter field for computer numerical simulation.

S103, performing mesh generation on the initial parameter field;

the model can be divided into regular grids or irregular grids according to actual needs, irregular grid subdivision is generally carried out on the condition that the undulating surface needs to be processed, and simple regular grid division is only required on the condition of the horizontal surface, namely, the grids are divided into rectangles. In actual industrial production, in order to reduce the calculation cost, the smaller the calculation amount, the smaller the memory requirement, and the better, that means the smaller the number of the mesh, the better. Meanwhile, in order to ensure stable simulation of the finite difference algorithm, the time sampling interval and the grid interval need to meet certain conditions. Therefore, in order to ensure efficient, accurate and stable finite difference numerical simulation, under the condition of satisfying the stability simulation and minimizing the numerical dispersion, the speed, the time sampling interval and the grid interval should satisfy the following relations:wherein v is_pmaxRepresents the maximum longitudinal wave velocity, Δ t represents the time sampling interval, and Δ d represents the divided grid spacing.

S105, determining wave equation coefficients related to the viscosity parameter and the anisotropy parameter for each grid point, wherein the coefficients comprise: a is₁，a₂，a₃，b₁，b₂，b₃And can be determined according to equation (10), and will not be described herein.

And S107, calculating according to the wave equation coefficient and the initial parameter field by adopting a preset wave field propagation operator to generate a longitudinal wave field value at each grid point.

According to a visco-acoustic anisotropic medium seismic wave field propagation operator given in a formula (12), the coefficient of a wave field propagation equation related to a space grid position and the viscosity and anisotropy parameter field are obtained through calculation, after a distributed shot point and a distributed wave detection point position are determined, finite difference numerical simulation is carried out through a computer, and therefore the seismic wave stress field value of the position of the wave detector is recorded.

The seismic wavefield propagation operator has been explained in detail above. In practical application, the presetting can be performed in the form of a functional module or an algorithm module. The shot and wave detector locations may be located at any position on the grid, typically on the first layer of the grid (i.e., the earth's surface).

As described above, in the numerical simulation process, equation (14) is solved by using a time second-order and space tenth-order difference cellular mode, equation (23) -equation (33) is calculated by using a pseudo-spectrum method, and numerical simulation based on the newly derived visco-acoustic anisotropic wave equation is realized, so that stress field data (namely, longitudinal wave field value) about seismic wave field propagation is obtained.

FIG. 2 is a comparison diagram of a wavefield simulation provided herein that simulates a wavefield without amplitude attenuation and phase dispersion effects, with phase dispersion effects only, with amplitude attenuation effects only, and with amplitude attenuation and phase dispersion effects simultaneously. In fig. 2, a is a simulated wave field containing neither amplitude attenuation nor phase dispersion, i.e. an anisotropic acoustic wave simulation result, and as shown in the figure, the wave field propagation speed contains directivity, and is faster along a certain direction. Part b of fig. 2 is the result of the seismic wavefield simulation with only the effect of dispersion, which is consistent in amplitude compared to part a of fig. 2, but the wavefront surface has a hysteresis effect. The section in figure 3 is the result of a seismic wavefield simulation with only the effect of amplitude attenuation, which is consistent with the wavefront surface but with amplitude attenuation compared to section a in figure 2. Part d of FIG. 2 is the result of a seismic wavefield simulation with both phase dispersion and amplitude attenuation, with both wavefront surface lag and amplitude energy reduction compared to part a of FIG. 2. According to the analysis, the viscoacoustic anisotropic wave equation provided by the invention can simulate the viscosity and anisotropic influence in the underground wave field propagation process, and simultaneously can independently simulate the amplitude attenuation influence and the phase frequency dispersion influence caused by the viscosity, and the characteristic is convenient for realizing attenuation-compensated reverse time migration imaging. In the figure, the abscissa represents the length x, and the ordinate represents the depth z.

In order to verify the accuracy of the viscous acoustic anisotropic forward modeling method in the simulation of the seismic wave field in the complex model, a complex air chimney model is firstly constructed and applied to the scheme provided by the invention. The model has strong viscosity and anisotropy characteristics, and a parametric model is shown in fig. 3, when seismic waves pass through a region containing viscosity, the seismic waves have the characteristics of amplitude attenuation and phase dispersion, and when the seismic waves pass through the region containing anisotropy, the characteristic that the wave field propagation speed is faster along a certain direction due to the influence of anisotropy is generated. Fig. 3 is a schematic diagram of the relationship between the viscosity and anisotropy parameters and the position of a complex model provided in the embodiments of the present disclosure. In the figure, the abscissa represents the length x, and the ordinate represents the depth z.

The conventional visco-acoustic anisotropic wave equation is chosen here as a reference, which has a high accuracy of seismic wavefield simulation, but whose amplitude attenuation effects are coupled with phase dispersion effects. Fig. 4 is a schematic diagram of a comparison between the method provided in the embodiment of the present disclosure and a snapshot of a wave field obtained by a conventional method at a certain time, where the abscissa is length x and the ordinate is depth z. Wherein, a part a in fig. 4 is a wave field snapshot simulated by a conventional method, b part in fig. 4 is a wave field snapshot simulated by an embodiment of the present specification, and c part in fig. 4 is a wave field superposition result obtained by superposing b part in fig. 4, which is represented by a red curve, on a part a in fig. 4, which is represented by a black curve. The black curve and the red curve are well matched, so that the accuracy of the visco-acoustic anisotropic medium forward modeling method provided by the patent is proved.

In order to verify the stability of the sticky sound anisotropic medium forward modeling method in the seismic wave field modeling in an extremely complex model, a parameter model containing viscosity and anisotropy is firstly established according to a BP2007 model and applied to the scheme provided by the invention. The model contains an anisotropic tilt model of a sharp tilt change. When the traditional visco-acoustic anisotropic seismic wave equation is used for simulating the propagation of seismic waves in the model containing severe inclination angle changes, high-frequency noise which increases exponentially with time is generated, and finally instability is generated. Fig. 5 is a schematic diagram of the relationship between the viscosity and anisotropy parameters and the position of the second model provided in the embodiments of the present disclosure. In the figure, the abscissa represents the length x, and the ordinate represents the depth z.

The conventional visco-acoustic anisotropic wave equation is still selected as a reference, and is easy to generate instability when being used for simulating a wave field in a model containing complex anisotropic dip angles. Fig. 6 is a schematic diagram of a comparison between the method provided in the embodiment of the present disclosure and a snapshot of a wave field obtained by a conventional method at a certain time, where the abscissa is length x and the ordinate is depth z. Wherein, part a in fig. 6 is a wave field snapshot simulated by the conventional method, and part b in fig. 6 is a wave field snapshot simulated by the embodiment of the present specification. The part a in fig. 6 contains high frequency instability as shown by an arrow, and the part b in fig. 6 does not generate high frequency instability, so that the stability of the method for simulating the visco-acoustic anisotropic medium forward simulation provided by the patent in numerical simulation in an extremely complex model is proved.

Correspondingly, the embodiment of the application also provides computer equipment, which comprises a memory, a processor and a computer program which is stored on the memory and can run on the processor, wherein when the processor executes the program, the forward modeling method for the acoustic wave seismic data in the viscous medium is realized.

The embodiments in the present specification are described in a progressive manner, and the same and similar parts among the embodiments are referred to each other, and each embodiment focuses on the differences from the other embodiments. Especially, as for the device, apparatus and medium type embodiments, since they are basically similar to the method embodiments, the description is simple, and the related points may refer to part of the description of the method embodiments, which is not repeated here.

Correspondingly, the embodiment of the application also provides computer equipment, which comprises a memory, a processor and a computer program which is stored on the memory and can run on the processor, wherein when the processor executes the program, the forward modeling method for the acoustic wave seismic data in the viscous medium is realized.

The embodiments in the present specification are described in a progressive manner, and the same and similar parts among the embodiments are referred to each other, and each embodiment focuses on the differences from the other embodiments. Especially, as for the device, apparatus and medium type embodiments, since they are basically similar to the method embodiments, the description is simple, and the related points may refer to part of the description of the method embodiments, which is not repeated here.

The embodiments in the present specification are described in a progressive manner, and the same and similar parts among the embodiments are referred to each other, and each embodiment focuses on the differences from the other embodiments. Especially, as for the device, apparatus and medium type embodiments, since they are basically similar to the method embodiments, the description is simple, and the related points may refer to part of the description of the method embodiments, which is not repeated here.

The foregoing description has been directed to specific embodiments of this disclosure. Other embodiments are within the scope of the following claims. In some cases, the actions or steps or modules recited in the claims may be performed in a different order than in the embodiments and still achieve desirable results. In addition, the processes depicted in the accompanying figures do not necessarily require the particular order shown, or sequential order, to achieve desirable results. In some embodiments, multitasking and parallel processing may also be possible or may be advantageous.