CN104977609A - Prestack longitudinal wave and transverse wave combined inversion method based on rapid simulated annealing - Google Patents

Abstract

The invention discloses a prestack longitudinal wave and transverse wave combined inversion method based on rapid simulated annealing by use of information of longitudinal waves and transverse waves. By use of conversion of longitudinal wave and transverse wave data, the speed of transverse waves is obtained through calculation, a natural acquisition presetack reflection data body is drawn to form angle gathers of different angles according to incident angles, and PP and PS wave prestack combined inversion is carried out. According to the invention, the inversion longitudinal wave impedance match is good, through observation of the error analysis of a data synthetic angle gather and an inversion data synthetic angle gather, it is indicated that the two angle gathers match well, a prestack PP and PS wave combined random inversion algorithm has quite high precision, and it is demonstrated that inversed longitudinal wave and transverse wave impedance profiles and Vp/Vs profiles have quite good corresponding relations with earthquake data and logging data.

Description

A kind of prestack based on fast simulated annealing is ripple joint inversion method in length and breadth

Technical field

The present invention is applicable to seismic interpretation multi-components prestack inversion field.Along with the sternness increasingly of whole world exploration form, the precision management of oil reservoir and fine description and non-constructilvity hydrocarbon-bearing pool difficulty increasing, very large challenge is brought to seismic interpretation, therefore, utilize prestack information, especially make full use of the convection cell inverting of shear wave data and identification seems particularly important.

Background technology

The joint inversion of multi-components seismic reflection data is subject to PP (compressional wave incident, longitudinal wave reflection) and the incidence of PS(compressional wave, transverse wave reflection always) the coupling difficult problems of ripple lineups perplexed.Although in refutation process, can complete joint inversion by coupling PP and PS ripple lineups, its amount of calculation is large and precision is low.Conventional PS ripple disposal route is ripple downgoing P-wave (compressional wave) and the up shear wave of reflection being assumed to a kind of simple propagation, and transformed wave is considered to neither with p wave interval velocity, the virtual or effective ripple neither propagated with shear wave velocity.PP ripple is processed into picture in PP travel-time territory, and PS ripple is processed into picture in PS travel-time territory.Even in same depth point, the data such as PP and the PS wave velocity obtained, Dao Ji, superposition and migrated section all have the different travel-times.These data make the horizon calibration in joint inversion and ripple coupling is very difficult in length and breadth.

Agullo described a kind of three steps invers$ method in 2004: from the inverting of PP and PS ripple, try to achieve S-wave impedance respectively; Carry out being correlated with to the low-frequency component of these two S-wave impedances and obtain P-S wave velocity ratio γ, and matching treatment is carried out to PS ripple; The joint inversion of PP and PS ripple.Current linear or non-linear inversion is all complete in this way substantially.But horizon calibration problem important between PP and PS ripple is also unresolved.Have multiple method to demarcate the layer position of PP and PS ripple at present, the waveform relationship between large multi-method supposition PP and PS ripple is identical or have identical feature, but, this hypothesis and the larger contradiction of physical presence are especially more obvious in reservoir position.

Dariu described the simulated annealing of the global optimization of a kind of multi-components AVO inverting in 2003.This global optimization approach, by mating the gross energy of real data and generated data, asks for P-S wave velocity ratio, although effective to the experimental verification of simulation and real data, but will rely on known actual data scaling layer position before inverting.

Multi-component seismic technology has the superiority simultaneously utilizing PP and PS ripple information to carry out reservoir prediction, can identify lithology, fluid and crack, but converted shear wave velocity analysis and the joint inversion of PP and PS ripple have difficulties.

Summary of the invention

The object of the invention is to provide a kind of PP and PS ripple imaging problem solving same Time-Scale Domain,

Improve inversion accuracy, fully utilize the ripple joint inversion method in length and breadth of the prestack based on fast simulated annealing of the information of ripple in length and breadth.

The present invention is realized by following steps:

1) field earthquake is blown out and is collected compressional wave and converted shear wave data, utilizes converted shear wave data, calculates shear wave velocity by following formula;

Calculate shear wave velocity Vs formula:

$t_{\mathrm{ps}}=\sqrt{t_{p0}^2+{\left(\frac{x-x_c}{V_p}\right)}^2}+\sqrt{t_{s0}^2+{\left(\frac{x_c}{V_s}\right)}^2}---\left(1\right)$

In formula:

T _psfor downgoing P-wave and up shear wave whilst on tour sum,

T _p0and t _s0be respectively P ripple and shear wave outward journey time,

Vp is P ripple,

X is geophone offset,

X _cfor transfer point is to the distance of shot point.

Shear wave velocity described in step 1) comprises stack velocity or migration before stack speed, and stack velocity or migration before stack speed are obtained by stack velocity spectrometry.

2) the angle gathers naturally gathering prestack reflection data body and to be taken into by incident angle different angles;

Step 2) described in angle gathers in PP ripple angle gathers be angle gathers PP ripple prestack reflection data body incident angle being taken into different angles, PS ripple angle gathers is angle gathers PS ripple prestack reflection data body incident angle being taken into different angles.

3) joint inversion of PP and PS ripple prestack is carried out;

Joint inversion described in step 3) is:

(1) PP and PS ripple is expressed as incident angle with minor function:

$\begin{array}{c}{R}_{\mathrm{PP}}\left(\mathrm{\θ}\right)\≈\left(\frac{1+{\tan }^2\mathrm{\θ}}{2}\right)\frac{\mathrm{\ΔI}}{I}-4\frac{{\mathrm{\β}}^2}{{\mathrm{\α}}^2}{\sin }^2\mathrm{\θ}\frac{\mathrm{\ΔJ}}{J}\\ -(\frac{1}{2}{\tan }^2\mathrm{\θ}-2\frac{{\mathrm{\β}}^2}{{\mathrm{\α}}^2}{\sin }^2\mathrm{\θ})\frac{\mathrm{\Δ\ρ}}{\mathrm{\ρ}}\end{array}$

In formula: R _pP(θ) be longitudinal wave reflection, R _ps(θ, φ) is converted shear wave reflection,

α, β, ρ are respectively by the compressional wave at interface and shear wave average velocity and average density,

Δ α, Δ β, Δ ρ are respectively the variable quantity of compressional wave by interface and shear wave velocity and density,

θ is average reflection by the compressional wave at interface and angle of transmission,

for average reflection and the angle of transmission of the transformed wave by interface,

I, J are respectively compressional wave and S-wave impedance,

Δ I, Δ J is respectively the variable quantity of p-wave impedance and S-wave impedance;

That is: compressional wave and converted shear wave data are that the formula being function with offset distance (angle) is expressed.

(2) establishing target function:

$\begin{array}{c}E={2\mathrm{\ω}}_p\frac{{\left|\right|d_p^{\mathrm{obs}}-d_p^{\mathrm{pre}}\left|\right|}_2^2}{{\left|\right|d_p^{\mathrm{obs}}-d_p^{\mathrm{pre}}\left|\right|}_2^2+{\left|\right|d_p^{\mathrm{obs}}-d_p^{\mathrm{pre}}\left|\right|}_2^2}+{2\mathrm{\ω}}_{\mathrm{ps}}\frac{{\left|\right|d_{\mathrm{ps}}^{\mathrm{obs}}-d_{\mathrm{ps}}^{\mathrm{pre}}\left|\right|}_2^2}{{\left|\right|d_{\mathrm{ps}}^{\mathrm{obs}}-d_{\mathrm{ps}}^{\mathrm{pre}}\left|\right|}_2^2+{\left|\right|d_{\mathrm{ps}}^{\mathrm{obs}}{-d}_{\mathrm{ps}}^{\mathrm{pre}}\left|\right|}_2^2}\\ +{\mathrm{\μ}}_1{\left|\right|m_p^{\mathrm{pri}}-m_p^{\mathrm{new}}\left|\right|}_l+{\mathrm{\μ}}_2{\left|\right|m_{\mathrm{ps}}^{\mathrm{pri}}-m_{\mathrm{ps}}^{\mathrm{new}}\left|\right|}_l\end{array}---\left(3\right)$

In formula: E objective function, l is l mould, be respectively PP ripple, the observation of PS ripple and predicted data, ( entirety is a parameter, and obs is the abbreviation of observe-observation, and pre is the abbreviation of predictive---prediction, usage of trade literary style) be respectively PP ripple, the iteration of PS ripple and preset model, ω _p, ω _ps, μ ₁, μ ₂for coefficient, ω _p, ω _psbe respectively the weight (example 1,1) of compressional wave and transformed wave, μ ₁, μ ₂the Error weight (between 0.001-1) being respectively compressional wave and transformed wave adjusts according to real data.

Described iteration and preset be calculative variable, be rate pattern or impedance model.

(3) iterative computation obtains p-wave impedance I, S-wave impedance J and the density p of each common midpoint or common transfer point (CMP (CCP)) position.

Described iterative computation is following formula:

M _i+1=M _i+ΔM _i（4）

In formula:

M _ithe p-wave impedance obtained for wanting inverting, S-wave impedance, density, at every turn to model vector M _iupgrade: if objective function E reaches minimum, then stop iteration, then obtain final model vector solution.

Described M _imodel obtains initial value by well-log information.

Present invention achieves PP and the PS wave datum process problem of same Time-Scale Domain, solve the difficult problem that time calibrating is inaccurate, through test, found by the contrast of observation with inverting p-wave impedance, both mate well.The error analysis that observation data synthesis angle gathers and invert data synthesize angle gathers shows, two angle gathers coincide good; Found by the contrast of observation with inverting S-wave impedance, both mate well, and the error analysis that observation data synthesis angle gathers and invert data synthesize angle gathers shows, two angle gathers coincide good.Illustrate that prestack PP & PS ripple associating stochastic inverse algorithm has very high precision.Tested by real data, the P-wave And S impedance profile of inverting and the display of Vp/Vs section have good corresponding relation with seismic data and well-log information.

Accompanying drawing explanation

Fig. 1: PP (left side) & PS (right side) pre-stack time migration sectional view;

Fig. 2: PP angle gathers (left side), PS angle gathers (right side);

Fig. 3: compressional wave recording error is analyzed;

Fig. 4: transformed wave recording error is analyzed;

Fig. 5: observation and inverting p-wave impedance contrast;

Fig. 6: observation and inverting S-wave impedance contrast.

Embodiment

The present invention is described in detail below in conjunction with the drawings and specific embodiments.

1) field earthquake is blown out and is collected compressional wave and converted shear wave data, utilizes converted shear wave data, calculates shear wave velocity by following formula;

Calculate shear wave velocity Vs formula:

$t_{\mathrm{ps}}=\sqrt{t_{p0}^2+{\left(\frac{x-x_c}{V_p}\right)}^2}+\sqrt{t_{s0}^2+{\left(\frac{x_c}{V_s}\right)}^2}---\left(1\right)$

In formula:

T _psfor downgoing P-wave and up shear wave whilst on tour sum,

T _p0and t _s0be respectively P ripple and shear wave outward journey time,

Vp is P ripple,

X is geophone offset,

Xc is the distance of transfer point to shot point.

Shear wave velocity described in step 1) comprises stack velocity or migration before stack speed, and stack velocity or migration before stack speed are obtained by stack velocity spectrometry.

Fig. 1 is PP & PS ripple Kirchhoff pre-stack time migration section, and same depth PP & PS ripple lineups are consistent for reflection interval.

2) the angle gathers naturally gathering prestack reflection data body and to be taken into by incident angle different angles.As shown in Figure 2, the left side is PP angle gathers, and the right side is PS angle gathers, and Taoist monastic name 1-8 position is PP ripple, Taoist monastic name 11-18 is PS ripple, angular range: 5-40 °, is spaced apart 5 °;

Step 2) described in angle gathers in PP ripple angle gathers be angle gathers PP ripple migration before stack reflectance data body incident angle being taken into different angles, PS ripple angle gathers is angle gathers PS ripple migration before stack reflectance data body incident angle being taken into different angles.

3) joint inversion of PP and PS ripple prestack is carried out;

Joint inversion described in step 3) is:

First, in wave-field simulation, compressional wave and transformed wave Aki & Richards (1980) approximate expression of Zoeppritz equation are the function of incident angle.

(1) PP and PS wave table is shown:

$\begin{array}{c}{R}_{\mathrm{PP}}\left(\mathrm{\θ}\right)\≈\left(\frac{1+{\tan }^2\mathrm{\θ}}{2}\right)\frac{\mathrm{\ΔI}}{I}-4\frac{{\mathrm{\β}}^2}{{\mathrm{\α}}^2}{\sin }^2\mathrm{\θ}\frac{\mathrm{\ΔJ}}{J}\\ -(\frac{1}{2}{\tan }^2\mathrm{\θ}-2\frac{{\mathrm{\β}}^2}{{\mathrm{\α}}^2}{\sin }^2\mathrm{\θ})\frac{\mathrm{\Δ\ρ}}{\mathrm{\ρ}}\end{array}$

In formula: R _pP(θ) be longitudinal wave reflection, Rps (θ, φ) is converted shear wave reflection,

α, β, ρ are respectively by the compressional wave at interface and shear wave average velocity and average density,

Δ α, Δ β, Δ ρ are respectively the variable quantity of compressional wave by interface and shear wave velocity and density,

θ is average reflection by the compressional wave at interface and angle of transmission,

for average reflection and the angle of transmission of the transformed wave by interface,

I, J are respectively compressional wave and S-wave impedance,

Δ I, Δ J is respectively the variable quantity of p-wave impedance and S-wave impedance;

Namely compressional wave and converted shear wave data are that the formula being function with offset distance (angle) is expressed.

(2) establishing target function:

In formula: E objective function, l is l mould, be respectively PP ripple, the observation of PS ripple and predicted data, ( entirety is a parameter, and obs is the abbreviation of observe-observation, and pre is the abbreviation of predictive---prediction, usage of trade literary style) be respectively PP ripple, the iteration of PS ripple and preset model, ω _p, ω _ps, μ ₁, μ ₂for coefficient, ω _p, ω _psbe respectively the weight (example 1,1) of compressional wave and transformed wave, μ ₁, μ ₂the Error weight (between 0.001-1) being respectively compressional wave and transformed wave adjusts according to real data.The first two fraction is that data is coincide, and latter two fraction and model are identical, after two parameters relevant with data signal to noise ratio (S/N ratio), due to, we expect that identical data is more some more, affect less, therefore by model, signal to noise ratio (S/N ratio) is higher, and data reliability is high, and latter two parameter value should be less.

Described iteration and preset be calculative variable, be rate pattern or impedance model.

(3) iterative computation obtains p-wave impedance I, S-wave impedance J and the density p of each common midpoint or common transfer point (CMP (CCP)) position.

Described iterative computation is following formula:

M _i+1=M _i+ΔM _i（4）

In formula: M _ithe p-wave impedance obtained for wanting inverting, S-wave impedance, density,

Described M _imodel obtains initial value by well-log information.That is survey line well-log information was utilized to set up initial model M ₀, then, at every turn to model vector M _iupgrade: finally make objective function E reach minimum, stop iteration.

It is the analysis of compressional wave recording error shown in Fig. 3, left-hand component is recorded as and utilizes observation data to synthesize angle gathers, center section is recorded as the angle gathers obtaining parameter calculating according to inverting, and right-hand component is the difference that observation data synthesis angle gathers and invert data synthesize angle gathers.

It is the analysis of transformed wave recording error shown in Fig. 4, left-hand component is recorded as and utilizes observation data to synthesize angle gathers, center section is recorded as the angle gathers obtaining parameter calculating according to inverting, right-hand component is the difference that observation data synthesis angle gathers and invert data synthesize angle gathers, and these two angle gathers coincide good.Utilize the method inverting to make this error of objective function E be very little as seen, observation data and generated data are closely.

Obtain final model vector solution, be illustrated in figure 5 observation and contrast with inverting p-wave impedance, left-hand component is observation p-wave impedance, and center section is inverting p-wave impedance, and right-hand component is for stacking contrast; Figure 6 shows that observation and inverting S-wave impedance contrast, left-hand component is observation S-wave impedance, and center section is inverting S-wave impedance, and right-hand component is for stacking contrast.The impedance matching that the impedance that inverting obtains obtains with observation is good.Illustrate that the prestack PP & PS ripple associating stochastic inverse algorithm based on fast simulated annealing has very high precision.

Claims (6)

1., based on a prestack ripple joint inversion method in length and breadth for fast simulated annealing, feature is realized by following steps:

1) field earthquake is blown out and is collected compressional wave and converted shear wave data, utilizes converted shear wave data to calculate shear wave velocity Vs;

2) the angle gathers naturally gathering prestack reflection data body and to be taken into by incident angle different angles;

3) joint inversion of PP and PS ripple prestack is carried out;

Joint inversion described in step 3) is:

(1) PP and PS ripple is expressed as incident angle with minor function:

In formula: R _pP(θ) be longitudinal wave reflection, R _ps(θ, φ) is converted shear wave reflection,

α, β, ρ are respectively by the compressional wave at interface and shear wave average velocity and average density, and Δ α, Δ β, Δ ρ are respectively the variable quantity of compressional wave by interface and shear wave velocity and density,

θ is average reflection by the compressional wave at interface and angle of transmission,

for average reflection and the angle of transmission of the transformed wave by interface,

I, J are respectively compressional wave and S-wave impedance,

Δ I, Δ J is respectively the variable quantity of p-wave impedance and S-wave impedance;

(2) following objective function is built:

$\begin{array}{c}E={2\mathrm{\ω}}_p\frac{{\left|\right|d_p^{\mathrm{obs}}-d_p^{\mathrm{pre}}\left|\right|}_2^2}{{\left|\right|d_p^{\mathrm{obs}}-d_p^{\mathrm{pre}}\left|\right|}_2^2+{\left|\right|d_p^{\mathrm{obs}}-d_p^{\mathrm{pre}}\left|\right|}_2^2}+{2\mathrm{\ω}}_{\mathrm{ps}}\frac{{\left|\right|d_{\mathrm{ps}}^{\mathrm{obs}}-d_{\mathrm{ps}}^{\mathrm{pre}}\left|\right|}_2^2}{{\left|\right|d_{\mathrm{ps}}^{\mathrm{obs}}-d_{\mathrm{ps}}^{\mathrm{pre}}\left|\right|}_2^2+{\left|\right|d_{\mathrm{ps}}^{\mathrm{obs}}{-d}_{\mathrm{ps}}^{\mathrm{pre}}\left|\right|}_2^2}\\ +{\mathrm{\μ}}_1{\left|\right|m_p^{\mathrm{pri}}-m_p^{\mathrm{new}}\left|\right|}_l+{\mathrm{\μ}}_2{\left|\right|m_{\mathrm{ps}}^{\mathrm{pri}}-m_{\mathrm{ps}}^{\mathrm{new}}\left|\right|}_l\end{array}---\left(3\right)$

In formula: E objective function, || || _lfor l mould, be respectively PP ripple, the observation of PS ripple and predicted data, be respectively PP ripple, the iteration of PS ripple and preset model, ω _p, ω _ps, μ ₁, μ ₂for coefficient, ω _p, ω _psbe respectively the weight of compressional wave and transformed wave, μ ₁, μ ₂be respectively the Error weight of compressional wave and transformed wave, adjust according to real data;

(3) iterative computation obtains p-wave impedance I, S-wave impedance J and the density p of each common midpoint or common transfer point (CMP (CCP)) position.

2. method according to claim 1, feature is that the calculating shear wave velocity Vs described in step 1) adopts following formula:

$t_{\mathrm{ps}}=\sqrt{t_{p0}^2+{\left(\frac{x-x_c}{V_p}\right)}^2}+\sqrt{t_{s0}^2+{\left(\frac{x_c}{V_s}\right)}^2}---\left(1\right)$ In formula:

T _psfor downgoing P-wave and up shear wave whilst on tour sum,

T _p0and t _s0be respectively P ripple and shear wave outward journey time,

Vp is P ripple,

X is geophone offset,

X _cfor transfer point is to the distance of shot point.

3. method according to claim 1, feature is that the shear wave velocity Vs described in step 1) comprises stack velocity or migration before stack speed, and stack velocity or migration before stack speed are obtained by stack velocity spectrometry.

4. method according to claim 1, feature is step 2) described in angle gathers in PP ripple angle gathers be angle gathers PP ripple prestack reflection data body incident angle being taken into different angles, PS ripple angle gathers is angle gathers PS ripple prestack reflection data body incident angle being taken into different angles.

5. method according to claim 1, feature is iteration described in step 3) and preset model is calculative variable, is rate pattern or impedance model.

6. method according to claim 1, feature is the iterative computation described in step 3) is following formula:

M _i+1=M _i+ΔM _i（4）

In formula: M _ithe p-wave impedance obtained for wanting inverting, S-wave impedance, density, at every turn to model vector M _iupgrade: if objective function E reaches minimum, then stop iteration, then obtain final model vector solution.