CN104007465B - Pure P wave data based on elastic inversion portrays reservoir development Details Methods - Google Patents

Abstract

The present invention relates to a kind of pure P wave data based on elastic inversion and portray reservoir development Details Methods.Utilize prestack Bao Fu road collection to obtain pure P wave data, effectively eliminate the impact of AVO effect on data resolution, and avoid traditional inversion method neutron deficiency estimation and low frequency model error on the impact of result.Model and actual data application confirm: the pure compressional wave data that prestack inversion obtains is higher than original full STACK DATA resolution, eliminates tuning effect; Compared with original full STACK DATA, the seismic wavelet extracted from pure compressional wave data is more stable within the scope of effective band, inverting section precision is higher, the reservoir distribution feature of prediction tallies with the actual situation more.Advantage of the present invention is to can be applicable to real data, obtains the pure P wave data being more suitable for interpretive analysis.

Description

Pure P wave data based on elastic inversion portrays reservoir development Details Methods

Technical field

The invention belongs to petroleum exploration field, relate to one and be adapted to remove AVO(Amplitudeversusoffset) effect and the pure P wave data based on elastic inversion that improves data resolution characteristic portray reservoir development Details Methods.

Background technology

Seism processing focuses on and obtains high resolving power, high s/n ratio and high-resolution result.In order to improve data signal to noise ratio (S/N ratio), main by the CRP(Commonreflectionpoint after skew at present) road collection stacks up (Mu Yongguang etc., 2006), and present earthquake industry also mainly to make an explanation analysis based on full superposition of data.The CRP road of a certain reflection spot is made to integrate as amplitude as S (θ), wherein θ=1,2 ..., n represents the incident compressional angle degree of not people having a common goal.Because seismic trace has carried out normal-moveout correction, need to determine a standard track y, make the difference of standard track and each recording channel minimum.Utilize least square principle, calculate the error sum of squares of any seismic trace and standard track:

$Q=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{[S\left({\mathrm{\θ}}_i\right)-y]}^2---\left(1\right)$

Wherein i is collection sequence number.Extreme value is asked to the above-mentioned multivariate function, then has:

$\begin{array}{c}\frac{\∂Q}{\∂y}=\frac{\∂}{\∂y}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{[S\left({\mathrm{\θ}}_i\right)-y]}^2\\ =-2\underset{i=1}{\overset{n}{\mathrm{\Σ}}}[S\left({\mathrm{\θ}}_i\right)-y]=0\end{array}---\left(2\right)$

Have: $\underset{i=1}{\overset{n}{\mathrm{\Σ}}}S\left({\mathrm{\θ}}_i\right)-\mathrm{ny}=0---\left(3\right)$

Then standard track, the amplitude namely after level superposition is: $y=\frac{1}{n}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}S\left({\mathrm{\θ}}_i\right)---\left(4\right)$

Each trace weighting coefficient that above formula participates in superposition is identical, and between actual each road, quality often has difference.When differing greatly between road in CRP gather, equal weight superposition is difficult to obtain ideal effect, has self-adaptation stacking method so follow-up:

$y_2=\frac{1}{n}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{k}_{i}S\left({\mathrm{\θ}}_i\right)---\left(5\right)$

Wherein, k is each road weight factor, can be calculated obtain by the degree of correlation of itself and standard track.

The original intention of level superposition is the signal to noise ratio (S/N ratio) in order to improve data, but when dynamic comparatively speed is forbidden, complex structure, in the situations such as offset distance is too large, normal moveout correction usually exists residual move out time, and thus superposition can reduce the resolution (Mu Yongguang etc., 2006) of seismic signal.In addition, such average effect can reduce resolution and the accuracy of the incident compressional wave data with reflecting of true vertical.In fact, namely sidewalk collection is completely straight, and the AVO effect that lithology or fluid difference produce also can make superposition reduction data resolution characteristic and order of accuarcy.So the data of superposition being used as the analysis that makes an explanation of compressional wave Normal incidence reflectance data will inevitably produce error.

After relative amplitude preserved processing, still there is certain amplitude variations between CRP road collection upper track, i.e. AVO feature in seismic data; Rule and the underground medium character of these changes are closely related, define the inversion method predictably descending dielectric attribute according to it.Existing statement amplitude mainly all comes from theory---the Knot-Zoeppritz equation of plane wave incidence and transmission along with the theoretical foundation that offset distance changes.The reflection coefficient formula of Application comparison Aki & Richards compressional wave incidence is widely (Aki & Richards, 1980):

$R_{\mathrm{pp}}\left(\mathrm{\θ}\right)\≈\frac{1\mathrm{\Δ\ρ}}{2\mathrm{\ρ}}(1-\frac{{4\mathrm{\β}}^2}{{\mathrm{\α}}^2}{\sin }^2\mathrm{\θ})+\frac{1}{{2\cos }^2\mathrm{\θ}}\frac{\mathrm{\Δ\α}}{\mathrm{\α}}-\frac{4{\mathrm{\β}}^2}{{\mathrm{\α}}^2}\frac{\mathrm{\Δ\β}}{\mathrm{\β}}{\sin }^2\mathrm{\θ}---\left(6\right)$

Wherein, the p-and s-wave velocity of interface upper and lower medium and density are respectively α ₁, α ₂, β ₁, β ₂, ρ ₁, ρ ₂, θ ₁, θ ₂reflection angle and angle of transmission respectively.Δα=α ₂-α ₁，Δβ=β ₂-β ₁，Δρ=ρ ₂-ρ ₁，α=（α ₁+α ₂）/2，β=（β ₁+β ₂）/2，θ=（θ ₁+θ ₂）/2，ρ=（ρ ₁+ρ ₂）/2。

Abbreviation is carried out to (1) formula, can following formula be obtained:

$R_{\mathrm{pp}}\left(\mathrm{\θ}\right)\″\frac{1}{2}(1+{\tan }^2\mathrm{\θ})(\frac{\mathrm{\Δ\α}}{\mathrm{\α}}+\frac{\mathrm{\Δ\ρ}}{\mathrm{\ρ}})-4\frac{{\mathrm{\β}}^2}{{\mathrm{\α}}^2}(\frac{\mathrm{\Δ\β}}{\mathrm{\β}}+\frac{\mathrm{\Δ\ρ}}{\mathrm{\ρ}}){\sin }^2\mathrm{\θ}-(\frac{1}{2}{\tan }^2\mathrm{\θ}-2\frac{{\mathrm{\β}}^2}{{\mathrm{\α}}^2}{\sin }^2\mathrm{\θ})\frac{\mathrm{\Δ\ρ}}{\mathrm{\ρ}}---\left(7\right)$

Usually, compressional wave normal-incidence reflection coefficient is:

$R_p=\frac{{\mathrm{\α}}_2{\mathrm{\ρ}}_2-{\mathrm{\α}}_1{\mathrm{\ρ}}_1}{{\mathrm{\α}}_2{\mathrm{\ρ}}_2+{\mathrm{\α}}_1{\mathrm{\ρ}}_1}=\frac{{\mathrm{\α}}_2{\mathrm{\ρ}}_2-{\mathrm{\α}}_2{\mathrm{\ρ}}_1+{\mathrm{\α}}_2{\mathrm{\ρ}}_1-{\mathrm{\α}}_1{\mathrm{\ρ}}_1}{{\mathrm{\α}}_2{\mathrm{\ρ}}_2+{\mathrm{\α}}_1{\mathrm{\ρ}}_1}=\frac{{\mathrm{\α}}_2\mathrm{\Δ\ρ}+{\mathrm{\Δ\α\ρ}}_1}{{\mathrm{\α}}_2{\mathrm{\ρ}}_2+{\mathrm{\α}}_1{\mathrm{\ρ}}_1}---\left(8\right)$

Because α ₁=α ₂-Δ α, α=(α ₁+ α ₂)/2, obviously, α ₂=2 α-α ₁=2 α-α ₂+ Δ α, then α ₂=α+Δ α/2,

Similarly, ρ ₁=ρ-Δ ρ/2,

So, α ₂ρ ₂+ α ₁ρ ₁=2 α ρ+Δ α ρ ≈ 2 α ρ.

Like this,

$R_p=\frac{{\mathrm{\α}}_2\mathrm{\Δ\ρ}+{\mathrm{\Δ\α\ρ}}_1}{{\mathrm{\α}}_2{\mathrm{\ρ}}_2+{\mathrm{\α}}_1{\mathrm{\ρ}}_1}\≈\frac{(\mathrm{\α}+\mathrm{\Δ\α}/2)\mathrm{\Δ\ρ}}{2\mathrm{\α\ρ}}+\frac{(\mathrm{\ρ}-\mathrm{\Δ\ρ}/2)\mathrm{\Δ\α}}{2\mathrm{\α\ρ}}\≈\frac{1}{2}(\frac{\mathrm{\Δ\ρ}}{\mathrm{\ρ}}+\frac{\mathrm{\Δ\α}}{\mathrm{\α}})---\left(9\right)$

Similarly, shear wave normal-incidence reflection coefficient is: $R_s\≈\frac{1}{2}(\frac{\mathrm{\Δ\ρ}}{\mathrm{\ρ}}+\frac{\mathrm{\Δ\β}}{\mathrm{\β}})---\left(10\right)$

(9) and (10) are brought into (7) formula, (Gidlowetal., 1992 can be obtained; Sun, 1999):

$R_{\mathrm{pp}}\left(\mathrm{\θ}\right)\≈\frac{1}{2}(1+{\tan }^2\mathrm{\θ})R_p-8\frac{{\mathrm{\β}}^2{\sin }^2\mathrm{\θ}}{{\mathrm{\α}}^2}{R}_s-(\frac{1}{2}{\tan }^2\mathrm{\θ}-2\frac{{\mathrm{\β}}^2}{{\mathrm{\α}}^2}{\sin }^2\mathrm{\θ})\frac{\mathrm{\Δ\ρ}}{\mathrm{\ρ}}---\left(11\right)$

Formula (7) or (11) are mainly used in prestack inversion and obtain the elastic parameters such as p-and s-wave velocity, in length and breadth wave reflection coefficient, density for data interpretation, but its result is usually controlled by the quality (Zhangetal. of wavelet estimation and initial low frequency model, 2011), must the participating in directly of log data and layer position information, there is larger personal error in net result.

Summary of the invention

The present invention proposes a kind of pure P wave data based on elastic inversion and portrays reservoir development Details Methods, and the method can remove AVO(Amplitudeversusoffset) effect and improve data resolution characteristic.By method provided by the invention, the pure P wave data higher than original full STACK DATA resolution can be obtained, clearly portray reservoir development details.

The specific embodiments that the present invention realizes above-mentioned purpose is as follows: a kind of pure P wave data based on elastic inversion portrays reservoir development Details Methods, comprises the steps:

Step 1: carry out relative amplitude preserved processing to seismic data, extracts high resolving power, high s/n ratio, high-fidelity common reflection point (CRP) road collection;

Step 2: based on high resolving power, high s/n ratio, high-fidelity common reflection point (CRP) road collects acquisition point angle superposition of data body;

Step 3: the pure P wave data acquiring method bipartition angle degree superposition of data set up based on the present invention calculates, and calculates the pure P wave data on each time-sampling point on each CRP point;

Step 4: export pure P wave data body;

Step 5: carry out data interpretation based on the pure P wave data calculated, and carry out geological Significance analysis.

The beneficial effect that compared with the prior art the present invention reaches is:

Direct pure P wave data acquiring method is established based on elastic inversion, prestack Bao Fu road collection is utilized to obtain pure P wave data, effectively eliminate the impact of AVO effect on data resolution, and avoid traditional inversion method neutron deficiency estimation and low frequency model error on the impact of result.Model and actual data application confirm: the pure compressional wave data that prestack inversion obtains is higher than original full STACK DATA resolution, eliminates tuning effect; Compared with original full STACK DATA, the seismic wavelet extracted from pure compressional wave data is more stable within the scope of effective band, inverting section precision is higher, the reservoir distribution feature of prediction tallies with the actual situation more.Advantage of the present invention is to can be applicable to real data, obtains the pure P wave data being more suitable for interpretive analysis.

Accompanying drawing explanation

Fig. 1 is different close contrast figure (a-d, first to fourth class AVO) of four class AVO exceptions.

The P wave data that Fig. 2 is the road collection of four class exceptions, vertical incidence data, superposition of data, inverting obtain and error comparison diagram (a-d, first to fourth class AVO; E1 and E2 is respectively superposition and inversion error).

Fig. 3 is the prestack inversion method flow diagram based on seismic constraint modeling.

Fig. 4 is the full superposition of data (a) in wellblock in the middle ancient times 8 and pure P-wave section (b) comparison diagram.

Fig. 5 is many wells wavelet comparison diagram (a, the wavelet form that full superposition and pure compressional wave two overlap data extraction; B, spectral amplitude; C, phase spectrum).

Fig. 6 be the middle ancient times 8 well full superposition and pure compressional wave two overlap data and respective inversion result figure (a, full superposition of data; B, full superposition p-wave impedance inversion result; C, pure P wave data; D, pure compressional wave p-wave impedance inversion result).

Fig. 7 is that wellblock hawk in the middle ancient times 8 mountain group entirely superposes and to overlap p-wave impedance prediction of result reservoir distribution figure (a, full stack result with pure compressional wave two; B, pure compressional wave result; Structure map is group top, hawk mountain).

Embodiment

Below in conjunction with example and accompanying drawing, the specific embodiment of the invention is described.

Pure P wave data based on elastic inversion portrays a reservoir development Details Methods, comprises the steps:

Step 1: carry out relative amplitude preserved processing to seismic data, extracts high resolving power, high s/n ratio, high-fidelity common reflection point (CRP) road collection;

Step 2: based on high resolving power, high s/n ratio, high-fidelity common reflection point (CRP) road collects selects high-quality part to obtain point angle superposition of data body;

Step 3: the pure P wave data acquiring method bipartition angle degree superposition of data set up based on the present invention calculates, and calculates the pure P wave data on each time-sampling point on each CRP point;

Step 4: export pure P wave data body;

Step 5: carry out data interpretation based on the pure P wave data calculated, and carry out geological Significance analysis.

The ultimate principle of the pure P wave data acquiring method based on elastic inversion that embodiment provides is as follows:

Prestack CRP road collection S (θ) is first directly considered as R by the method principle _pp(θ), new pure P wave data computing formula is set up:

$S_{\mathrm{pp}}\left(\mathrm{\θ}\right)\≈\frac{1}{2}(1+{\tan }^2\mathrm{\θ})S_p-8\frac{{\mathrm{\β}}^2{\sin }^2\mathrm{\θ}}{{\mathrm{\α}}^2}{S}_s-(\frac{1}{2}{\tan }^2\mathrm{\θ}-2\frac{{\mathrm{\β}}^2}{{\mathrm{\α}}^2}{\sin }^2\mathrm{\θ})S_D---\left(12\right)$

Wherein, S _pp(θ) be pre-stack seismic road collection, S _p=R _p* W _t, S _s=R _p* W _t, S _d=R _d* W _tbe respectively pure P wave data, shear wave data and density data.W _tfor seismic wavelet, R _d=Δ ρ/ρ.Although the quality of prestack inversion depends on the quality of previous earthquake data, inverting itself can not change the impact of data quality on inversion result, so road collection S (θ) of prestack is directly considered as R _pp(θ) set up formula and do not produce personal error.

Can find out, formula (12) is Nonlinear System of Equations, and a given known mudstone relation (beta/alpha), can try to achieve compressional wave S linearly _p, shear wave S _swith density S _dparameter.But shear wave S _swith density S _dthe levels of precision of parameter is directly related to the accuracy of initial mudstone relation (beta/alpha), but this is most unknown in practical situations both.Accordingly, compressional wave S _pparameter and mudstone relation (beta/alpha) have nothing to do, very stable.Therefore, based on mudstone relation (beta/alpha) hypothesis near 0.5, non-linear formula (12) is rewritten as linear formula further:

$S_{\mathrm{pp}}\left(\mathrm{\θ}\right)\≈\frac{1}{2}(1+{\tan }^2\mathrm{\θ})S_p-{2\sin }^2\mathrm{\θ}{S}_s-\frac{1}{2}({\tan }^2\mathrm{\θ}-{\sin }^2\mathrm{\θ})S_D---\left(12\right)$

S accurately can be tried to achieve by above formula when road collection number is greater than three _p.Now owing to not removing seismic wavelet, the S of gained _pbe real P wave data.Owing to eliminating the interference of AVO effect, and there is not any human factor, so the more traditional full superposition of data of its result is better.

Fig. 1 ~ Fig. 2 numerical model forward simulation information contrast changes the difference of pure compressional wave data that method asks for and full STACK DATA, the outstanding validity changing method.Introduce gas sand to be extremely analyzed by the four quasi-representative AVO that shale bed covers, its basic parameter is as shown in table 1.

The basic elastic parameter of table 1 four class AVO exception

Fig. 1 is the reflection coefficient of just drilling four kinds of gas sands in the table 1 of acquisition according to the various approximate formulas of Knot-Zoeppritz equation along with offset distance changes data.Obviously, each class gas sand all has certain AVO feature, namely each class is in reflection coefficient (can be considered amplitude) the numerical value difference all to some extent at different offset distance or incident angle place, and the poststack result adopting all average or average weighted modes to obtain is inevitable different from numerical value time vertical incidence (zero degree).

Select the Ricker wavelet of 50Hz, the above-mentioned four class gas sand road collection of 0 to 42 degree scopes are just being drilled respectively according to Aki & Richards formula, then adopt the mode of average superposition to obtain its corresponding superposition of data respectively, adopt formula of the present invention (12) to obtain its pure P wave data (S _p).

Fig. 2 is the Sp contrast that the reflectance data at four class gas sand zero incident angle places, stack result and the present invention obtain.Error_1 and Error_2 represents the reflectance data at zero incident angle place and the S of stack result and inverting respectively _pdifference.Can find out: prestack inversion obtain pure P wave data almost there is no error, and superposition result along with the different errors of AVO feature also different.Usually, in the conventional data ranges of incidence angles of 0-42 degree, stack result slightly can reduce than the zero degree incident angle of first kind AVO, seriously increases than the result of Equations of The Second Kind AVO, seriously increases than the result of the 3rd class, slightly reduces than the result of the 4th class.

Fig. 3 is the process flow diagram of the pure P wave data acquiring method based on elastic inversion:

Step 1: carry out relative amplitude preserved processing to seismic data, extracts high resolving power, high s/n ratio, high-fidelity common reflection point (CRP) road collection;

Step 2: based on resolution, high s/n ratio, high-fidelity common reflection point (CRP) road collects acquisition point angle superposition of data body;

Step 3: the pure P wave data acquiring method bipartition angle degree superposition of data set up based on the present invention calculates, and calculates the pure P wave data on each time-sampling point on each CRP point;

Step 4: export pure P wave data body;

Step 5: carry out data interpretation based on the pure P wave data calculated, and carry out geological Significance analysis.

Fig. 4 is the pure P-wave section comparison diagram for the original full superposition of data in wellblock in the middle ancient times 8 and prestack inversion obtain, and figure Green curve is interface, carbonatite top (well compactness top To3s).Can find out: after AVO tuning effect is removed, the pure P-wave section of right side inverting obviously improves than the full stacked section resolution in left side.First: the position, Silurian reflection horizon of such as figure upper yellow arrows position instruction be the very wide negative reflection of redness on full stacked section, and after removing AVO effect correspondence clearly two covers reflect; Secondly, the Sang Tamu group that on figure, black arrow is corresponding and Silurian unconformability contact relation on pure P-wave section clearly; Again, the reflection (red arrow) of on carbonatite interface and its underpart is also more clear on pure P-wave section.

Fig. 5 is the carbonatite target interval many wells average wavelet contrast utilizing full superposition and the pure compressional wave two that obtains to overlap data to extract, and wavelet is extracted according to its composite traces and the highest principle of the geological data degree of correlation.Can find out: the red seismic wavelet extracted from pure P wave data increases, as shown in the black arrow in Fig. 5 b than the blue seismic wavelet frequency extracted from full superposition of data; In addition, can find out that the former is more stable within the scope of the effective band of 0-40hz, close to zero phase from the phase spectrum of Fig. 5 c.In fact, although utilize prestack inversion to remove AVO effect thus the degree improving seismic data resolution does not have deconvolution in processing procedure obvious to the raising of dominant frequency, but the process later stage improve data quality and improve that reservoir prediction precision is most important, unheeded step again.

Fig. 6 be the middle ancient times 8 well two cover data and the p-wave impedance Comparative result of corresponding inverting.Visible pure P-wave section and inversion result thereof are all slightly high than full stack result resolution, and the waveform character of beading reflection is also more perfect.Especially, its reflection of black arrow upper 3.95 second position process removal AVO effect is schemed more clear, the solution cavity response that correspondence one is little.

Fig. 7 is wellblock hawk in the middle ancient times 8 mountain group carbonate reservoir distribution plan of the p-wave impedance inversion result prediction according to full superposition of data and pure P wave data, it is group top, hawk mountain that figure constructs map, overall by southwest directional structure vectorical structure reduction northeastward, grow west and south TZ23 structural belt, TZ10 structural belt and northern palaeohigh etc. three structure high-order bits therebetween respectively.Visible: the reservoir predicted according to full superposition of data is mainly distributed in northern protuberance and tower near No. 10 zone of fracture, and TZ23 structural belt reservoir distribution is less, and comprises the wavelet tuning effect illusion at northern low position; And be mainly distributed on three structural belts according to the reservoir that pure compressional wave is predicted, be that some favourable beading distributions have appearred in south the most significantly, and the tuning interference of the wavelet eliminating northeast, with brill well logging information and geological condition more identical.So far proved that the pure P wave data that AVO inverting obtains is higher than traditional full stack result resolution, the reservoir distribution according to its prediction is more reasonable.

Claims (1)

1. the pure P wave data based on elastic inversion portrays a reservoir development Details Methods, it is characterized in that, described method comprises the steps:

Step 1: carry out relative amplitude preserved processing to seismic data, extracts common reflection point (CRP) the road collection of high resolving power, high s/n ratio, high-fidelity;

Step 2: obtain a point angle superposition of data based on high resolving power, high s/n ratio, high-fidelity reflection spot (CRP) road collection;

Step 3: bipartition angle degree superposition of data calculates, calculates the pure P wave data on each time-sampling point on each road collection point;

Step 4: export pure P wave data body;

Step 5: carry out data interpretation based on the pure P wave data calculated, and carry out geological Significance analysis;

Wherein: described bipartition angle degree superposition of data body calculates, and the step calculating the pure P wave data on each road collection point on each time-sampling point is:

The first step: prestack CRP road collection S (θ) is directly considered as R _pp(θ), new pure P wave data nonlinear computation formula is set up:

$S_{pp}\left(\mathrm{\θ}\right)\≈\frac{1}{2}(1+{\tan }^2\mathrm{\θ})S_p-8\frac{{\mathrm{\β}}^2{\sin }^2\mathrm{\θ}}{{\mathrm{\α}}^2}{S}_s-(\frac{1}{2}{\tan }^2\mathrm{\θ}-2\frac{{\mathrm{\β}}^2}{{\mathrm{\α}}^2}{\sin }^2\mathrm{\θ})S_D$

Wherein, S _pp(θ) be pre-stack seismic road collection; S _p=R _p* W _t, S _s=R _s* W _t, S _d=R _d* W _tbe respectively pure P wave data, shear wave data and density data; W _tfor seismic wavelet, R _d=Δ ρ/ρ, θ are incident compressional angle degree, and α is velocity of longitudinal wave, and β is shear wave velocity, and ρ is Media density, R _ppfor the reflection coefficient of compressional wave incidence, R _pcompressional wave normal-incidence reflection coefficient, R _sfor shear wave normal-incidence reflection coefficient;

Second step: based on mudstone relation (beta/alpha) hypothesis near 0.5, is rewritten as linear formula further by non-linear formula (12):

$S_{pp}\left(\mathrm{\θ}\right)\≈\frac{1}{2}(1+{\tan }^2\mathrm{\θ})S_p-2{\sin }^2{\mathrm{\θS}}_s-\frac{1}{2}({\tan }^2\mathrm{\θ}-{\sin }^2\mathrm{\θ})S_D$

P wave data S is accurately obtained by described linear formula when road collection number is greater than three _p.