Numerical simulation method for fluid viscoelastic rock seismic dispersion attenuation wave field

Abstract

The invention belongs to the field of seismic exploration of oil and gas resources, and provides a numerical simulation method for a seismic dispersion attenuation wave field containing fluid viscoelastic rock. According to the method, a reflection coefficient calculation method for accurately describing the energy conversion relation of seismic waves at the viscoelastic rock interface is established, the reflection time difference caused by seismic wave speed dispersion can be reflected, and the dispersion attenuation characteristics of the seismic waves propagating in viscoelastic stratum rocks are completely described; the method can be used for numerical simulation of seismic wave field propagation characteristics of various different complex rocks, and an intuitive and efficient implementation way is provided for researching seismic wave propagation rules of complex viscoelastic media.

Description

Technical Field

[0001] The present invention belongs to the fields of conventional and unconventional oil and gas resource seismic exploration, carbon storage rock physics and seismic wave propagation simulation, and realizes a numerical simulation method for seismic dispersion attenuation wave field of viscoelastic rock containing fluid. Background Art

[0002] Studying the propagation characteristics of seismic wave fields in viscoelastic reservoir rocks is the basis for using the attribute information in seismic waves to detect underground strata structures, achieve reservoir parameter characterization and oil and gas detection, etc. Numerical simulation of seismic wave fields based on wave equations can be used to describe the propagation characteristics of seismic waves in viscoelastic media, but this type of method often faces huge computational overhead in complex rock physics models, which limits its application.

[0003] The numerical simulation method of seismic waves based on the convolution model has fast calculation speed and is easy to implement. It has been widely used in actual seismic exploration, such as well seismic calibration of synthetic seismic records and seismic inversion. At present, the convolution model has been extended to the viscoelastic framework so that it can partially describe the dispersion attenuation characteristics of seismic waves. The dispersion attenuation of seismic waves is manifested in that the seismic reflection amplitude decreases with the increase of propagation distance, and the reflection coefficient value and reflection time at the same reflection interface change with different frequencies. However, the convolution model algorithm under the existing viscoelastic framework often ignores the change of seismic reflection time on the specific formation interface with frequency, and fails to fully and completely describe the dispersion attenuation characteristics of the seismic wave field in the viscoelastic medium. The method of the patent of this invention is to solve this problem. Summary of the invention

[0004] The purpose of the present invention is to provide a numerical simulation method for the seismic dispersion attenuation wave field of viscoelastic rocks containing fluids, which is characterized by establishing a reflection coefficient calculation method that accurately describes the energy conversion relationship of seismic waves at the viscoelastic rock interface, can reflect the reflection time difference caused by the dispersion of seismic wave velocity, and fully describe the dispersion attenuation characteristics of seismic waves propagating in viscoelastic formation rocks. The method mainly includes the following main steps:

[0005] The rock physical parameters required for modeling fluid-containing fractured porous rocks are obtained through logging and core data, including: rock grain bulk modulus and density, dry rock bulk modulus and shear modulus, porosity, and pore fluid bulk modulus, density, viscosity, occurrence radius and saturation;

[0006] Given a set containing The frequency vector of the frequencies , set the index number ;

[0007] From the frequency vector Extract the index number The corresponding frequency value ;

[0008] According to the steps The rock physical parameters are calculated by using the dynamic rock physical model method to calculate the density of the rock. ,frequency The corresponding longitudinal wave phase velocity , longitudinal wave attenuation coefficient , shear wave phase velocity , shear wave attenuation coefficient , is the two-way travel time of the seismic wave;

[0009] Using the viscoelastic rock seismic wave reflection theory, the frequency The reflection coefficient The calculation method is as follows:

[0010]

[0011] in, is the incident angle of the seismic wave, represents the seismic reference reflection coefficient calculated by the present invention, It means taking the difference between the parameters of the upper and lower strata;

[0012] The frequency is calculated by the following formula: Correspondingly, the unsteady reflection coefficient that fully describes the dispersion and attenuation effect of seismic waves is :

[0013]

[0014] The frequency is calculated by the following formula: The corresponding seismic wavelet components :

[0015]

[0016] in, is the main frequency of the seismic wavelet;

[0017] Calculate the frequency according to the following formula Corresponding seismic wave field :

[0018]

[0019] in, Represents the convolution operation;

[0020] make , repeat steps to , until After that, the final seismic dispersion attenuation wave field data is calculated according to the following formula. :

[0021] . BRIEF DESCRIPTION OF THE DRAWINGS

[0022] Figure 1 The four-layer geological model designed for the present invention has the first, second and fourth layers being completely elastic mudstone layers with thicknesses of 86m, 152m and 200m respectively, and the third layer being a viscoelastic gas-bearing sandstone layer with a thickness of 162m. Figure 2 1 is a non-steady-state reflection coefficient for completely describing the dispersion attenuation effect of seismic waves at different frequencies according to an embodiment of the present invention. The horizontal axis is frequency in Hertz (Hz), and the vertical axis is time in milliseconds (ms). Figure 3 It is the seismic dispersion attenuation wave field record at different incident angles. The horizontal axis is time in milliseconds (ms), and the vertical axis is amplitude, which is dimensionless. DETAILED DESCRIPTION

[0023] In order to make the technical solution of the present invention clearer and the advantages clearer, the specific implementation methods of the technical solution in the present invention will be clearly and completely described below in conjunction with the diagrams. It should be understood that the specific implementation methods described here are only used to explain the present invention and are not used to limit the present invention.

[0024] Build as Figure 1 The geological model shown in the figure, based on the well logging and core data, obtains the P-wave and S-wave velocities and density of each layer of mudstone in the surrounding rock (Table 1); obtains the rock grain bulk modulus and density, dry rock bulk modulus, dry rock shear modulus and porosity of the gas-bearing sandstone layer (Table 2); obtains the bulk modulus, density, viscosity, occurrence radius and saturation of the pore fluid (Table 3);

[0025] Table 1 Elastic parameters of mudstone

[0026] Physical parameters Mudstone 1 Mudstone 2 Mudstone 3 P-wave velocity (m / s) 1710 2035 2035 Shear wave velocity (m / s) 1038 1237 1237 <![CDATA[Density (g / cm 3 )]]> 2.08 2.109 2.109

[0027] Table 2 Physical parameters of gas-bearing sandstone rock skeleton

[0028] Physical parameters Numeric Rock particle bulk modulus (MPa) 37 <![CDATA[Rock particle density (g / cm 3 )]]> 2.65 Dry rock bulk modulus (MPa) 4.8 Dry rock shear modulus (MPa) 5.7 Porosity (%) 30

[0029] Table 3 Physical parameters of pore fluid in gas-bearing sandstone

[0030] Physical parameters water gas Bulk modulus (GPa) 2.25 0.012 <![CDATA[Density (g / cm 3 ).]]> 1.04 0.078 Viscosity (Pa • s) 0.03 0.0015 Radius (m) 0.4 0.29 Saturation (%) 60 40

[0031] Given a frequency vector containing N=125 frequencies , set the index number ;

[0032] From the frequency vector Extract the index number The corresponding frequency value ;

[0033] According to the steps The rock physical parameters are calculated by using the dynamic rock physical model method to calculate the density of the rock. ,frequency The corresponding longitudinal wave phase velocity , longitudinal wave attenuation coefficient , shear wave phase velocity , shear wave attenuation coefficient , is the two-way travel time of the seismic wave;

[0034] Using the viscoelastic rock seismic wave reflection theory, the frequency The reflection coefficient The calculation method is as follows:

[0035]

[0036] in, is the incident angle of the seismic wave, represents the seismic reference reflection coefficient calculated by the present invention, It means taking the difference between the parameters of the upper and lower strata;

[0037] The frequency is calculated by the following formula: Correspondingly, the unsteady reflection coefficient that fully describes the dispersion and attenuation effect of seismic waves is ,like Figure 2 As shown:

[0038]

[0039] The frequency is calculated by the following formula: The corresponding seismic wavelet components :

[0040]

[0041] in, is the main frequency of the seismic wavelet;

[0042] Calculate the frequency according to the following formula Corresponding seismic wave field :

[0043]

[0044] in, Represents the convolution operation;

[0045] make , repeat steps to , until After that, the final seismic dispersion attenuation wave field data is calculated according to the following formula. ,like Figure 3 As shown:

[0046] .

[0047] Figure 2 This is the reflection coefficient of different frequencies when the seismic wave incident angle θ=0º of the embodiment of the present invention. It can be seen from the figure that at the same frequency, the reflection coefficient of the bottom interface of the gas-bearing sandstone layer (located between 350ms-400ms) has obvious amplitude attenuation compared with the reflection coefficient of the top interface of the gas-bearing sandstone layer (located at 250ms); at the bottom interface of the gas-bearing sandstone layer, the reflection coefficients of different frequencies, their time positions gradually decrease with the increase of frequency.

[0048] Figure 3 The figures are seismic dispersion attenuation wave field records at different incident angles (θ=0º, 15º, 30º) of the embodiment of the present invention. It can be seen from the figure that the seismic waves at the bottom of the gas-bearing sandstone layer (between 350ms and 400ms) have obvious amplitude attenuation and phase change, and the attenuation degree and phase change characteristics are different at different incident angles. The amplitude attenuation is caused by the attenuation of seismic waves of different frequencies, and the phase change is caused by the velocity dispersion of different frequencies.

[0049] The above embodiments show that the numerical simulation method proposed in the present invention can completely describe the dispersion attenuation characteristics of the viscoelastic rock seismic wave field.

[0050] The advantages and novelties of the present invention are: A reflection coefficient calculation method was established that can accurately describe the energy conversion relationship of seismic waves at the viscoelastic rock interface; It fully reflects the difference in seismic reflection time caused by velocity dispersion of seismic waves, and can fully describe the dispersion and attenuation characteristics of the seismic wave field of viscoelastic rocks; The method of the present invention can be used for numerical simulation of seismic wave field propagation characteristics of various complex rocks without solving complex wave equations, and provides an intuitive and efficient implementation approach for studying the seismic wave propagation laws of complex viscoelastic media.

[0051] The above embodiments are only used to illustrate the present invention, wherein the various implementation steps of the method can be changed, and any equivalent transformations and improvements based on the technical solution of the present invention should not be excluded from the protection scope of the present invention.

Claims

1. A numerical simulation method for seismic dispersion attenuation wave field of viscoelastic rock containing fluid, mainly comprising the following steps: The rock physical parameters required for modeling fluid-containing fractured porous rocks are obtained through logging and core data, including: rock grain bulk modulus and density, dry rock bulk modulus and shear modulus, porosity, and pore fluid bulk modulus, density, viscosity, occurrence radius and saturation; Given a set containing The frequency vector of the frequencies , set the index number ; From the frequency vector Extract the index number The corresponding frequency value ; According to the steps The rock physical parameters are calculated by using the dynamic rock physical model method to calculate the density of the rock. ,frequency The corresponding longitudinal wave phase velocity , longitudinal wave attenuation coefficient , shear wave phase velocity , shear wave attenuation coefficient , is the two-way travel time of the seismic wave; Using the viscoelastic rock seismic wave reflection theory, the frequency The reflection coefficient The calculation method is as follows: in, is the incident angle of the seismic wave, represents the seismic reference reflection coefficient calculated by the present invention, It means taking the difference between the parameters of the upper and lower strata; The frequency is calculated by the following formula: Correspondingly, the unsteady reflection coefficient that fully describes the dispersion and attenuation effect of seismic waves is : The frequency is calculated by the following formula: The corresponding seismic wavelet components : in, is the main frequency of the seismic wavelet; Calculate the frequency according to the following formula Corresponding seismic wave field : in, Represents the convolution operation; make , repeat steps to , until After that, the final seismic dispersion attenuation wave field data is calculated according to the following formula. : 。

Images