CN106125135A - Gas-bearing sandstone reservoir seismic response method for numerical simulation based on petrophysical model - Google Patents

Abstract

The present invention provides a kind of gas-bearing sandstone reservoir seismic response method for numerical simulation based on petrophysical model, based on spherical patchy saturation, rock physics and the seismic dynamics feature of conventional gas sand can preferably be described, apply viscid disperse fluctuation Equation Theory, introduce decay and energy loss that seismic wave causes because of the existence of fluid when propagating.Frequency-wavenumber domain viscid disperse wave equation is utilized to carry out forward simulation, effectively the geophysical character of petrophysical model can be showed, the viscid feature containing fluid media (medium) brought in conjunction with viscid dispersion theory and dispersive test, it is possible to more intuitively and exactly provide the seismic response features of conventional gas sand.

Description

Gas-bearing sandstone reservoir seismic response method for numerical simulation based on petrophysical model

Technical field

The present invention relates to geological technique field, particularly relate to a kind of gas-bearing sandstone reservoir earthquake based on petrophysical model Response method for numerical simulation.

Background technology

The attenuation of seismic wave feature analysis being currently based on petrophysical model rests on reflection coefficient discussion more, and traditional In Wave equation forward modeling, mostly the foundation of model is based on to the estimation of seimic wave velocity in actual formation, therefore owing to increasing The subjectivity of human intervention becomes inaccurate, and the geological model set up based on petrophysical model more conforms to reality, therefore Petrophysical model and the combination of Wave equation forward modeling are the development trends that method for numerical simulation is inevitable.

Patchy saturation is that one is full of fluid, and the petrophysical model of part gassiness partially aqueous, and it can be used Simulate the petrophysics property of conventional gas sand；Viscid disperse wave equation is mainly used to the simulation earthquake containing fluid media (medium) Response characteristic.

Summary of the invention

It is an object of the invention to solve the defect that above-mentioned prior art exists, it is provided that a kind of based on petrophysical model Gas-bearing sandstone reservoir seismic response method for numerical simulation, it is possible to effectively the geophysical character of petrophysical model is shown Come, more intuitively and exactly provide the seismic response features of conventional gas sand.

Not up to above-mentioned purpose, the present invention uses following scheme:

Gas-bearing sandstone reservoir seismic response method for numerical simulation based on petrophysical model, petrophysical model is space The spherical patchy saturation of periodic arrangement；The method comprises the following steps:

According to the effective bulk modulus of pore-fluid in the case of Wood law calculating low frequency, thus utilize Gassmann equation Obtain rock equivalent volume modulus；First calculate the bulk modulus of different fluid, recycling according to Gassmann equation under high frequency Hill Theoretical Calculation rock volume modulus；Utilize Johnson theory to provide dynamic volume modulus, and then calculate P wave complex-velocity Degree；Obtain the dispersion relation of P-wave phase velocity and quality factor according to Carcoine attenuation theory, and be introduced into frequency wave Number field obtains complex wave number, is combined with frequency and becomes the initial wave field that reflection coefficient formula obtains with wavelet convolution, substitutes into two dimension viscid Disperse wave equation carries out numerical simulation.

Further, process is implemented as follows:

Arranging spherical patchy saturation based on space periodic, ball interior is full of gas, and outside is full of water, when seismic wave passes When being multicast to speckle saturated media, owing to the existence of gas and water has resulted in different pore pressures, utilize Biot theoretical explanation Corresponding pore pressure homeostasis process, and provide corresponding seismic wave diffusion length:

$L=\sqrt{F_i/\mathrm{\ω}}---\left(1\right)$

Wherein ω is angular frequency, F_iFor diffusion coefficient, i=g, w；Wherein g represents gas, and w represents water:

Wherein, κ is the permeability of patchy saturation, and η is the coefficient of viscosity of pore-fluid,M_c、K_avAccording to containing fluid Some physical attributes of porous media are given:

$M_c\left(K_f\right)=K_G\left(K_f\right)+\frac{4}{3}\mathrm{\μ}---\left(5\right)$

In formula, K_m、K_s、K_fBeing respectively dry rock, solid particle and the bulk modulus of pore-fluid, μ is to be situated between containing fluid porous The modulus of shearing of matter, theoretical according to Hill, it is assumed that the modulus of shearing of remaining each blending ingredients is the most identical, and φ is porosity, K_G(K_f) Then calculate according to corresponding fluid modulus and Gassmann equation:

When frequency of seismic wave is of a sufficiently low, have sufficient time between fluid to reach pore pressure balance, based on Wood Law can obtain the effective bulk modulus of pore-fluid:

$K_R={\left(\frac{f_i}{K_i}\right)}^{-1}---\left(7\right)$

Wherein f_i(i=g, w) is the volume fraction of medium Fluid in Pore, replaces with gas and water saturation degree here, K_i(i =g, w) is corresponding gas and the bulk modulus of water；Then under low frequency background, according to Gassmann equation calculate rock etc. Effect bulk modulus:

Relative, when frequency of seismic wave is the highest, between fluid, the not free pore pressure that reaches balances, now in hole The pressure produced is uneven, if assuming, pore pressure is different but is constant, first can calculate according to Gassmann equation The bulk modulus of different fluid, the rock equivalent volume modulus further according under Hill Theoretical Calculation high frequency:

$K_H=\underset{i=g,w}{\Σ}\frac{M_c\left(K_i\right)}{S_i}-\frac{4}{3}\mathrm{\μ}---\left(9\right)$

According to the rock equivalent volume modulus under high and low frequency, when trying to achieve corresponding intermediate frequency containing gas and water rock etc. Effect bulk modulus:

$K_D=K_H-\frac{K_H-K_L}{1-\mathrm{\θ}+\mathrm{\θ}\sqrt{1+100i\mathrm{\ω}\mathrm{\σ}/{\mathrm{\θ}}^2}}---\left(10\right)$

In formula, i is imaginary unit, θ and σ can obtain according to the physical property etc. of dry rock, fluid respectively, provides phase here The computing formula answered:

$\mathrm{\θ}=\frac{K_H-K_L}{2K_L}\·\frac{\mathrm{\σ}}{T}---\left(11\right)$

$\mathrm{\σ}=\frac{{(K_H-K_L)}^2}{K_H^2{G}^2}---\left(12\right)$

$G={\[\frac{\underset{i=g,w}{\Σ}(Z_i+Q_i)M_c\left(K_i\right)}{\underset{i=g,w}{\Σ}{\mathrm{\φS}}_i{K}_G\left(K_i\right)M_c\left(K_i\right)}\]}^2\·\frac{S}{V}\sqrt{F}---\left(13\right)$

Z_i=φ²K_av(K_i) (14)

$F={\left(\frac{{\mathrm{\κK}}_H}{\underset{i=g,w}{\Σ}{\mathrm{\η}}_i\sqrt{D_i}}\right)}^2---\left(16\right)$

In above formula, D_i(i=g w) is the diffusion coefficient of corresponding fluid, can be obtained by the Biot Theoretical Calculation mentioned before Arriving, in space periodic arrangement patchy saturation, S/V, T are represented by:

$\frac{S}{V}=3\frac{R_g^2}{R_w^2}---\left(17\right)$

$T=\frac{1}{\mathrm{\κ}}\frac{K_L{\mathrm{\φ}}^2}{30R_w^3}\left\{\begin{array}{c}\[3{\mathrm{\η}}_w{g}_w^2+5({\mathrm{\η}}_g-{\mathrm{\η}}_w)g_g{g}_w-3{\mathrm{\η}}_g{g}_g^2\]R_g^5\\ +5g_w\[3{\mathrm{\η}}_w{g}_w-(2{\mathrm{\η}}_w+{\mathrm{\η}}_g)g_g\]R_g^2{R}_w^3\\ -15{\mathrm{\η}}_w{g}_w(g_w-g_g)R_g^3{R}_w^2-3{\mathrm{\η}}_w{g}_w^2{R}_w^5\end{array}\right\}---\left(18\right)$

Wherein:

It is calculated the complex velocity that P wave is propagated in gassiness, aqueous reservoir by above-mentioned formula:

$V_R=\sqrt{\frac{K_D+\frac{4}{3}\mathrm{\μ}}{{\mathrm{\ρ}}_b}}---\left(20\right)$

Wherein, ρ_bFor the body density of medium,ρ in formula_s、ρ_iIt is respectively solid particle With the density of fluid, S_iSaturation for corresponding pore-fluid；According to Carcoine attenuation theory, complex velocity is utilized to provide ground The seismic wave equivalent phase velocity in attenuation medium and quality factor:

$V_P={\[\mathrm{Re}\left(\frac{1}{V_R}\right)\]}^{-1}---\left(21\right)$

$Q=\frac{\mathrm{Re}\left(V_R^2\right)}{\mathrm{Im}\left(V_R^2\right)}---\left(22\right)$

Obtain seismic wave phase velocity and the dispersion relation of quality factor in speckle saturated media according to this, utilize above-mentioned frequency dispersion to close System, in conjunction with the viscid-disperse wave equation of two dimension, carries out numerical simulation.

Further, viscid-disperse Wave equation forward modeling to implement step as follows:

First according to the viscid-disperse wave equation of two dimension:

$\frac{{\∂}^{2}u}{\∂t^2}+\mathrm{\γ}\frac{\∂u}{\∂t}-\mathrm{\η}(\frac{{\∂}^{3}u}{\∂x^2\∂t}+\frac{{\∂}^{3}u}{\∂z^2\∂t})-v^2(\frac{{\∂}^{2}u}{\∂x^2}+\frac{{\∂}^{2}u}{\∂z^2})=0---\left(23\right)$

In above formula, u is particle displacement, and γ is the dispersion coefficient containing fluid media (medium), and η is the coefficient of viscosity containing fluid media (medium), v For seimic wave velocity in medium, provide the solution of the simple harmonic wave form of equation:

$u=e^{i(k_{z}z+k_{x}x)}{e}^{i\mathrm{\ω}t}---\left(24\right)$

In formula, ω is circular frequency, k_x、k_zBe respectively x, z direction wave number, unit be 1/m, i be imaginary unit.By above formula Substitution wave equation has:

Then have:

$k_z^2=\frac{{\mathrm{\ω}}^2(v^2-\mathrm{\γ}\mathrm{\η})}{v^4+{\mathrm{\ω}}^2{\mathrm{\η}}^2}{k}_x^2+i\frac{\mathrm{\ω}({\mathrm{\γv}}^2+{\mathrm{\η\ω}}^2)}{v^4+{\mathrm{\ω}}^2{\mathrm{\η}}^2}---\left(26\right)$

Make k_z=k+i α (27)

Both members square obtainsObtain k and the α expression formula about viscid dispersion coefficient；

Frequency wavenumber domain wave field extrapolation formula according to Gazdag etc.:

$P(z_0+dz)=P\left(z_0\right)e^{{\mathrm{ik}}_{z}dz}=P\left(z_0\right)e^{i(k+i\mathrm{\α})dz}---\left(28\right)$

Wherein, e^i(k+iα)dzFor phase shift factor, utilize above formula carry out prestack viscid-disperse wave equation forward numerical simulation；

According to the dispersion relation formula (21) of seismic wave phase velocity in patchy saturation, substituted into the formula calculating wave number Formula (26-28) tries to achieve phase shift factor, and the seismic response features to gas sand that designs a model carries out accurate analysis.

The present invention is based on spherical patchy saturation, it is possible to preferably describe rock physics and the earthquake of conventional gas sand Dynamic characteristic, apply viscid disperse fluctuate Equation Theory, introduce seismic wave propagate time because of the existence of fluid declining of causing Subtract and energy loss.Frequency wavenumber domain viscid disperse wave equation is utilized to carry out forward simulation, it is possible to effectively by rock physics The geophysical character of model shows, and the viscid feature containing fluid media (medium) brought in conjunction with viscid dispersion theory and disperse are special Levy, it is possible to more intuitively and exactly provide the seismic response features of conventional gas sand.

Accompanying drawing explanation

Fig. 1 geological model

Seismic wave phase velocity and the dispersion relation curve of inverse quality factor under Fig. 2 difference gas saturation

Three class sandstone reservoir seismic response under Fig. 3 difference gas saturation

Fig. 4 difference gas saturation interface echo frequency spectrum at present

Detailed description of the invention

For making the object, technical solutions and advantages of the present invention clearer, below technical scheme in the present invention carry out clearly Chu, it is fully described by, it is clear that described embodiment is a part of embodiment of the present invention rather than whole embodiments.Based on Embodiment in the present invention, it is every other that those of ordinary skill in the art are obtained under not making creative work premise Embodiment, broadly falls into the scope of protection of the invention.

The invention provides a kind of gas-bearing sandstone reservoir seismic response method for numerical simulation based on petrophysical model, should Model is that space periodic arranges spherical patchy saturation, and the method comprises the following steps: calculate low frequency feelings according to Wood law The effective bulk modulus of pore-fluid under condition, thus utilize Gassmann equation to obtain rock equivalent volume modulus；Head under high frequency First calculate the bulk modulus of different fluid according to Gassmann equation, recycle Hill Theoretical Calculation rock volume modulus；Utilize Johnson theory provides dynamic volume modulus, and then calculates P wave complex velocity；Obtain ground according to Carcoine attenuation theory Shake P phase velocity of wave and the dispersion relation of quality factor, and be introduced into frequency-wavenumber domain and obtain complex wave number, it is combined with frequency and becomes anti- Penetrate the initial wave field that coefficient formula obtains with wavelet convolution, substitute into two-dimentional viscid disperse wave equation and can carry out numerical simulation. The present invention gives the gas-bearing sandstone reservoir method for numerical simulation of a kind of reality of more fitting.Below for implementing process:

Arranging spherical patchy saturation based on space periodic, ball interior is full of gas, and outside is full of water.When seismic wave passes When being multicast to speckle saturated media, owing to the existence of gas and water has resulted in different pore pressures, utilize Biot theoretical explanation Corresponding pore pressure homeostasis process, and provide corresponding seismic wave diffusion length:

$L=\sqrt{F_i/\mathrm{\ω}}---\left(1\right)$

Wherein ω is angular frequency, F_i(i=g, w；Wherein g represents gas, and w represents water) be diffusion coefficient:

Wherein, κ is the permeability of patchy saturation, and η is the coefficient of viscosity of pore-fluid,M_c、K_avThen can be according to containing Some physical attributes of fluid porous media are given:

$M_c\left(K_f\right)=K_G\left(K_f\right)+\frac{4}{3}\mathrm{\μ}---\left(5\right)$

In formula, K_m、K_s、K_fBeing respectively dry rock, solid particle and the bulk modulus of pore-fluid, μ is to be situated between containing fluid porous The modulus of shearing of matter, theoretical according to Hill, it will be assumed that the modulus of shearing of remaining each blending ingredients is the most identical, and φ is porosity, K_G(K_f) then can calculate according to corresponding fluid modulus and Gassmann equation:

When frequency of seismic wave is of a sufficiently low, have sufficient time between fluid to reach pore pressure balance, based on Wood Law can obtain the effective bulk modulus of pore-fluid:

$K_R={\left(\frac{f_i}{K_i}\right)}^{-1}---\left(7\right)$

Wherein f_i(i=g, w) is the volume fraction of medium Fluid in Pore, replaces with gas and water saturation degree here, K_i(i =g, w) is corresponding gas and the bulk modulus of water.Then, under low frequency background, rock can be calculated according to Gassmann equation Equivalent volume modulus:

Relative, when frequency of seismic wave is the highest, between fluid, the not free pore pressure that reaches balances, now in hole The pressure produced is uneven, if assuming, pore pressure is different but is constant, first can calculate according to Gassmann equation The bulk modulus of different fluid, the rock equivalent volume modulus further according under Hill Theoretical Calculation high frequency:

$K_H=\underset{i=g,w}{\Σ}\frac{M_c\left(K_i\right)}{S_i}-\frac{4}{3}\mathrm{\μ}---\left(9\right)$

According to the rock equivalent volume modulus under high and low frequency, can be in the hope of the rock Han gas and water during corresponding intermediate frequency Equivalent volume modulus:

$K_D=K_H-\frac{K_H-K_L}{1-\mathrm{\θ}+\mathrm{\θ}\sqrt{1+100i\mathrm{\ω}\mathrm{\σ}/{\mathrm{\θ}}^2}}---\left(10\right)$

In formula, i is imaginary unit, θ and σ can obtain according to the physical property etc. of dry rock, fluid respectively, provides phase here The computing formula answered:

$\mathrm{\θ}=\frac{K_H-K_L}{2K_L}\·\frac{\mathrm{\σ}}{T}---\left(11\right)$

$\mathrm{\σ}=\frac{{(K_H-K_L)}^2}{K_H^2{G}^2}---\left(12\right)$

$G={\[\frac{\underset{i=g,w}{\Σ}(Z_i+Q_i)M_c\left(K_i\right)}{\underset{i=g,w}{\Σ}{\mathrm{\φS}}_i{K}_G\left(K_i\right)M_c\left(K_i\right)}\]}^2\·\frac{S}{V}\sqrt{F}---\left(13\right)$

Z_i=φ²K_av(K_i) (14)

$F={\left(\frac{{\mathrm{\κK}}_H}{\underset{i=g,w}{\Σ}{\mathrm{\η}}_i\sqrt{D_i}}\right)}^2---\left(16\right)$

In above formula, D_i(i=g w) is the diffusion coefficient of corresponding fluid, can be obtained by the Biot Theoretical Calculation mentioned before Arrive.In space periodic arrangement patchy saturation, S/V, T are represented by:

$\frac{S}{V}=3\frac{R_g^2}{R_w^2}---\left(17\right)$

$T=\frac{1}{\mathrm{\κ}}\frac{K_L{\mathrm{\φ}}^2}{30R_w^3}\left\{\begin{array}{c}\[3{\mathrm{\η}}_w{g}_w^2+5({\mathrm{\η}}_g-{\mathrm{\η}}_w)g_g{g}_w-3{\mathrm{\η}}_g{g}_g^2\]R_g^5\\ +5g_w\[3{\mathrm{\η}}_w{g}_w-(2{\mathrm{\η}}_w+{\mathrm{\η}}_g)g_g\]R_g^2{R}_w^3\\ -15{\mathrm{\η}}_w{g}_w(g_w-g_g)R_g^3{R}_w^2-3{\mathrm{\η}}_w{g}_w^2{R}_w^5\end{array}\right\}---\left(18\right)$

Wherein:

It is calculated the complex velocity that P wave is propagated in gassiness, aqueous reservoir by above-mentioned formula:

$V_R=\sqrt{\frac{K_D+\frac{4}{3}\mathrm{\μ}}{{\mathrm{\ρ}}_b}}---\left(20\right)$

Wherein, ρ_bFor the body density of medium,ρ in formula_s、ρ_iIt is respectively solid particle With the density of fluid, S_iSaturation for corresponding pore-fluid.According to Carcoine attenuation theory, complex velocity is utilized to provide ground The seismic wave equivalent phase velocity in attenuation medium and quality factor:

$V_P={\[\mathrm{Re}\left(\frac{1}{V_R}\right)\]}^{-1}---\left(21\right)$

$Q=\frac{\mathrm{Re}\left(V_R^2\right)}{\mathrm{Im}\left(V_R^2\right)}---\left(22\right)$

The most i.e. can get seismic wave phase velocity and the dispersion relation of quality factor in speckle saturated media.

Utilize above-mentioned dispersion relation, in conjunction with the viscid-disperse wave equation of two dimension, numerical simulation can be carried out.To be situated between below Continue viscid-disperse Wave equation forward modeling implements process.

First according to the viscid-disperse wave equation of two dimension:

$\frac{{\∂}^{2}u}{\∂t^2}+\mathrm{\γ}\frac{\∂u}{\∂t}-\mathrm{\η}(\frac{{\∂}^{3}u}{\∂x^2\∂t}+\frac{{\∂}^{3}u}{\∂z^2\∂t})-v^2(\frac{{\∂}^{2}u}{\∂x^2}+\frac{{\∂}^{2}u}{\∂z^2})=0---\left(23\right)$

In above formula, u is particle displacement, and γ is the dispersion coefficient containing fluid media (medium), and η is the coefficient of viscosity containing fluid media (medium), v For seimic wave velocity in medium.Provide the solution of the simple harmonic wave form of equation:

$u=e^{i(k_{z}z+k_{x}x)}{e}^{i\mathrm{\ω}t}---\left(24\right)$

In formula, ω is circular frequency, k_x、k_zBe respectively x, z direction wave number, unit be 1/m, i be imaginary unit.By above formula Substitution wave equation has:

Then have:

Make k_z=k+i α (27)

Both members square obtainsStudy according to forefathers, high-order α can be saved²?.Then K and the α expression formula about viscid dispersion coefficient can be obtained.

Frequency wavenumber domain wave field extrapolation formula according to Gazdag etc. can have:

$P(z_0+dz)=P\left(z_0\right)e^{{\mathrm{ik}}_{z}dz}=P\left(z_0\right)e^{i(k+i\mathrm{\α})dz}---\left(28\right)$

Wherein, e^i(k+iα)dzFor phase shift factor, utilize above formula can carry out prestack viscid-disperse Wave equation forward modeling numerical value Simulation.

According to the dispersion relation (formula 21) of seismic wave phase velocity in patchy saturation, substituted into above-mentioned calculating wave number Formula (formula 26-28) tries to achieve phase shift factor, designs a model and the seismic response features of gas sand can be carried out accurate point Analysis.

Accompanying drawing 1 is three layers of geological model of design, and it is shale with bottom at the middle and upper levels, and centre is gas-bearing sandstone reservoir.If Counting a undercompacted sandstone, its physical parameter such as permeability, porosity given, wherein the bulk modulus of rock particles is 38GPa, density is 2.65g/cm³, the bulk modulus of dry rock and modulus of shearing are respectively 1.56GPa and 1.10GPa, hole Degree is 33%, and in hole, the bulk modulus of institute's gassiness is 0.03GPa, and density is 0.15g/cm³, the bulk modulus of water is 2.42GPa, density is 1.0g/cm³。

Utilize described technology path, first available undercompaction sandstone reservoir based on patchy saturation as shown in Figure 2 Middle P-wave phase velocity and the dispersion relation curve of inverse quality factor.It can be seen that when gas saturation increases, be situated between In matter, fluid creates tremendous influence to the decay of seismic wave, causes its speed to there occurs the reduction that hundreds of meters is per second, and Low-frequency range there occurs obvious dispersion phenomenon.The dispersion curve of inverse quality factor then indicates, and gas saturation increases energy The peak value enough making inverse quality factor moves to high frequency direction, and during crest frequency, seimic wave velocity decay is the strongest.

Seismic wave phase velocity and the dispersion relation curve of inverse quality factor under Fig. 2 difference gas saturation

And then according to the viscid-disperse wave equation of two dimension, substitute into the velocity dispersion relation of patchy saturation, and based on figure Geological model shown in 1, can obtain the 128th road prestack numeric modeling result as shown in Figure 3, can be significantly from figure To due to the increase of air content, the decay of seismic wave is strengthened by fluid, causes interface under reservoir to reflect the now significantly time and prolongs Late and amplitude variations, wherein amplitude is changed by reservoir natural impedance is affected and is increased.

Extract the frequency spectrum of reflecting interface seismic wave under reservoir, as shown in Figure 4, it can be seen that the dominant frequency of seismic wave is along with storage In Ceng, the increase Han Fluid Volume is gradually moved to low frequency, reflects that seismic wave becomes strong in the decay of high band.

According to said method, based on spherical patchy saturation, utilize viscid-disperse Wave equation forward modeling, it is possible to effectively Ground reflects the decay and loss, petrophysical model and seismic wave that in classical gas sand, seismic wave occurs because of the existence of fluid The combination that dynamic equation is just being drilled can provide earthquake model forward modeling method more accurately, and this can be for being further appreciated by containing fluid layer Seismic response provide certain theories integration.

Last it is noted that above example is only in order to illustrate technical scheme, it is not intended to limit；Although With reference to previous embodiment, the present invention is described in detail, it will be understood by those within the art that: it still may be used So that the technical scheme described in foregoing embodiments to be modified, or wherein portion of techniques feature is carried out equivalent； And these amendment or replace, do not make appropriate technical solution essence depart from various embodiments of the present invention technical scheme spirit and Scope.

Claims (3)

1. gas-bearing sandstone reservoir seismic response method for numerical simulation based on petrophysical model, it is characterised in that: rock physics Model is that space periodic arranges spherical patchy saturation；

The method comprises the following steps:

According to the effective bulk modulus of pore-fluid in the case of Wood law calculating low frequency, thus Gassmann equation is utilized to obtain Rock equivalent volume modulus；First calculate the bulk modulus of different fluid under high frequency according to Gassmann equation, recycle Hill Theoretical Calculation rock volume modulus；Utilize Johnson theory to provide dynamic volume modulus, and then calculate P wave complex velocity； Obtain the dispersion relation of P-wave phase velocity and quality factor according to Carcoine attenuation theory, and be introduced into frequency wave number Territory obtains complex wave number, is combined with frequency and becomes the initial wave field that reflection coefficient formula obtains with wavelet convolution, substitutes into two dimension viscid more Scattered wave equation carries out numerical simulation.

Gas-bearing sandstone reservoir seismic response method for numerical simulation based on petrophysical model the most according to claim 1, It is characterized in that implementing process as follows:

Arranging spherical patchy saturation based on space periodic, ball interior is full of gas, and outside is full of water, when seimic wave propagation arrives During speckle saturated media, owing to the existence of gas and water has resulted in different pore pressures, utilize Biot theoretical explanation corresponding Pore pressure homeostasis process, and provide corresponding seismic wave diffusion length:

$L=\sqrt{F_i/\mathrm{\ω}}---\left(1\right)$

Wherein ω is angular frequency, F_iFor diffusion coefficient, i=g, w；Wherein g represents gas, and w represents water:

Wherein, κ is the permeability of patchy saturation, and η is the coefficient of viscosity of pore-fluid,M_c、K_avAccording to containing fluid porous Some physical attributes of medium are given:

$M_c\left(K_f\right)=K_G\left(K_f\right)+\frac{4}{3}\mathrm{\μ}---\left(5\right)$

In formula, K_m、K_s、K_fBeing respectively dry rock, solid particle and the bulk modulus of pore-fluid, μ is Saturate porous medium Modulus of shearing, theoretical according to Hill, it is assumed that the modulus of shearing of remaining each blending ingredients is the most identical, and φ is porosity, K_G(K_f) then root Calculate according to corresponding fluid modulus and Gassmann equation:

When frequency of seismic wave is of a sufficiently low, have sufficient time between fluid to reach pore pressure balance, based on Wood law Can obtain the effective bulk modulus of pore-fluid:

$K_R={\left(\frac{f_i}{K_i}\right)}^{-1}---\left(7\right)$

Wherein f_i(i=g, w) is the volume fraction of medium Fluid in Pore, replaces with gas and water saturation degree here, K_i(i=g, W) it is corresponding gas and the bulk modulus of water；Then, under low frequency background, the equivalent of rock is calculated according to Gassmann equation Product module amount:

Relative, when frequency of seismic wave is the highest, the not free pore pressure balance that reaches, now generation in hole between fluid Pressure be uneven, if assuming, pore pressure is different but is constant, can first calculate difference according to Gassmann equation The bulk modulus of fluid, the rock equivalent volume modulus further according under Hill Theoretical Calculation high frequency:

$K_H=\underset{i=g,w}{\Σ}\frac{M_c\left(K_i\right)}{S_i}-\frac{4}{3}\mathrm{\μ}---\left(9\right)$

According to the rock equivalent volume modulus under high and low frequency, when trying to achieve corresponding intermediate frequency, contain the equivalent of gas and water rock Product module amount:

$K_D=K_H-\frac{K_H-K_L}{1-\mathrm{\θ}+\mathrm{\θ}\sqrt{1+100i\mathrm{\ω}\mathrm{\σ}/{\mathrm{\θ}}^2}}---\left(10\right)$

In formula, i is imaginary unit, θ and σ can obtain according to the physical property etc. of dry rock, fluid respectively, is given corresponding here Computing formula:

$\mathrm{\θ}=\frac{K_H-K_L}{2K_L}\·\frac{\mathrm{\σ}}{T}---\left(11\right)$

$\mathrm{\σ}=\frac{{(K_H-K_L)}^2}{K_H^2{G}^2}---\left(12\right)$

$G={\[\frac{\underset{i=g,w}{\Σ}(Z_i+Q_i)M_c\left(K_i\right)}{\underset{i=g,w}{\Σ}{\mathrm{\φS}}_i{K}_G\left(K_i\right)M_c\left(K_i\right)}\]}^2\·\frac{S}{V}\sqrt{F}---\left(13\right)$

Z_i=φ²K_av(K_i) (14)

$F={\left(\frac{{\mathrm{\κK}}_H}{\underset{i=g,w}{\Σ}{\mathrm{\η}}_i\sqrt{D_i}}\right)}^2---\left(16\right)$

In above formula, D_i(i=g, w) is the diffusion coefficient of corresponding fluid, can be obtained by the Biot Theoretical Calculation mentioned before, In space periodic arrangement patchy saturation, S/V, T are represented by:

$\frac{S}{V}=3\frac{R_g^2}{R_w^2}---\left(17\right)$

$T=\frac{1}{\mathrm{\κ}}\frac{K_L{\mathrm{\φ}}^2}{30R_w^3}\left\{\begin{array}{c}\[3{\mathrm{\η}}_w{g}_w^2+5({\mathrm{\η}}_g-{\mathrm{\η}}_w)g_g{g}_w-3{\mathrm{\η}}_g{g}_g^2\]R_g^5\\ +5g_w\[3{\mathrm{\η}}_w{g}_w-(2{\mathrm{\η}}_w+{\mathrm{\η}}_g)g_g\]R_g^2{R}_w^3\\ -15{\mathrm{\η}}_w{g}_w(g_w-g_g)R_g^3{R}_w^2-3{\mathrm{\η}}_w{g}_w^2{R}_w^5\end{array}\right\}---\left(18\right)$

Wherein:

It is calculated the complex velocity that P wave is propagated in gassiness, aqueous reservoir by above-mentioned formula:

$V_R=\sqrt{\frac{K_D+\frac{4}{3}\mathrm{\μ}}{{\mathrm{\ρ}}_b}}---\left(20\right)$

Wherein, ρ_bFor the body density of medium,ρ in formula_s、ρ_iIt is respectively solid particle and stream The density of body, S_iSaturation for corresponding pore-fluid；According to Carcoine attenuation theory, complex velocity is utilized to provide seismic wave Equivalent phase velocity in attenuation medium and quality factor:

$V_P={\[\mathrm{Re}\left(\frac{1}{V_R}\right)\]}^{-1}---\left(21\right)$

$Q=\frac{\mathrm{Re}\left(V_R^2\right)}{\mathrm{Im}\left(V_R^2\right)}---\left(22\right)$

Obtain seismic wave phase velocity and the dispersion relation of quality factor in speckle saturated media according to this, utilize above-mentioned dispersion relation, In conjunction with the viscid-disperse wave equation of two dimension, carry out numerical simulation.

Gas-bearing sandstone reservoir seismic response method for numerical simulation based on petrophysical model the most according to claim 2, It is characterized in that: viscid-disperse Wave equation forward modeling to implement step as follows:

First according to the viscid-disperse wave equation of two dimension:

$\frac{{\∂}^{2}u}{\∂t^2}+\mathrm{\γ}\frac{\∂u}{\∂t}-\mathrm{\η}(\frac{{\∂}^{3}u}{\∂x^2\∂t}+\frac{{\∂}^{3}u}{\∂z^2\∂t})-v^2(\frac{{\∂}^{2}u}{\∂x^2}+\frac{{\∂}^{2}u}{\∂z^2})=0---\left(23\right)$

In above formula, u is particle displacement, and γ is the dispersion coefficient containing fluid media (medium), and η is the coefficient of viscosity containing fluid media (medium), and v is for being situated between Seimic wave velocity in matter, provides the solution of the simple harmonic wave form of equation:

$u=e^{i(k_{z}z+k_{x}x)}{e}^{i\mathrm{\ω}t}---\left(24\right)$

In formula, ω is circular frequency, k_x、k_zBe respectively x, z direction wave number, unit be 1/m, i be imaginary unit.Above formula is substituted into Wave equation has:

Then have:

$k_z^2=\frac{{\mathrm{\ω}}^2(v^2-\mathrm{\γ}\mathrm{\η})}{v^4+{\mathrm{\ω}}^2{\mathrm{\η}}^2}-k_x^2+i\frac{\mathrm{\ω}({\mathrm{\γv}}^2+{\mathrm{\η\ω}}^2)}{v^4+{\mathrm{\ω}}^2{\mathrm{\η}}^2}---\left(26\right)$

Make k_z=k+i α (27)

Both members square obtainsObtain k and the α expression formula about viscid dispersion coefficient；

Frequency wavenumber domain wave field extrapolation formula according to Gazdag etc.:

$P(z_0+dz)=P\left(z_0\right)e^{{\mathrm{ik}}_{z}dz}=P\left(z_0\right)e^{i(k+i\mathrm{\α})dz}---\left(28\right)$

Wherein, e^i(k+iα)dzFor phase shift factor, utilize above formula carry out prestack viscid-disperse wave equation forward numerical simulation；

According to the dispersion relation formula (21) of seismic wave phase velocity in patchy saturation, substituted into the formula formula calculating wave number (26-28) trying to achieve phase shift factor, the seismic response features to gas sand that designs a model carries out accurate analysis.