CN104977609A - Prestack longitudinal wave and transverse wave combined inversion method based on rapid simulated annealing - Google Patents
Prestack longitudinal wave and transverse wave combined inversion method based on rapid simulated annealing Download PDFInfo
- Publication number
- CN104977609A CN104977609A CN201410145440.6A CN201410145440A CN104977609A CN 104977609 A CN104977609 A CN 104977609A CN 201410145440 A CN201410145440 A CN 201410145440A CN 104977609 A CN104977609 A CN 104977609A
- Authority
- CN
- China
- Prior art keywords
- wave
- ripple
- data
- prestack
- angle
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Pending
Links
Landscapes
- Geophysics And Detection Of Objects (AREA)
- Measurement Of Mechanical Vibrations Or Ultrasonic Waves (AREA)
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
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:
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:
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:
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, μ
1, μ
2for coefficient, ω
p, ω
psbe respectively the weight (example 1,1) of compressional wave and transformed wave, μ
1, μ
2the 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:
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:
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, μ
1, μ
2for coefficient, ω
p, ω
psbe respectively the weight (example 1,1) of compressional wave and transformed wave, μ
1, μ
2the 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
0, 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:
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, μ
1, μ
2for coefficient, ω
p, ω
psbe respectively the weight of compressional wave and transformed wave, μ
1, μ
2be 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
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.
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201410145440.6A CN104977609A (en) | 2014-04-11 | 2014-04-11 | Prestack longitudinal wave and transverse wave combined inversion method based on rapid simulated annealing |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201410145440.6A CN104977609A (en) | 2014-04-11 | 2014-04-11 | Prestack longitudinal wave and transverse wave combined inversion method based on rapid simulated annealing |
Publications (1)
Publication Number | Publication Date |
---|---|
CN104977609A true CN104977609A (en) | 2015-10-14 |
Family
ID=54274284
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201410145440.6A Pending CN104977609A (en) | 2014-04-11 | 2014-04-11 | Prestack longitudinal wave and transverse wave combined inversion method based on rapid simulated annealing |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN104977609A (en) |
Cited By (5)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN105807316A (en) * | 2016-04-25 | 2016-07-27 | 吉林大学 | Surface observation microseism speed model correcting method based on amplitude stack |
CN108333626A (en) * | 2018-01-10 | 2018-07-27 | 成都理工大学 | A kind of genetic algorithm Optimum Impedance Inversion Method based on best retention strategy |
CN108363099A (en) * | 2018-02-11 | 2018-08-03 | 中国石油化工股份有限公司 | The longitudinal and shear wave correlation technique analyzed based on simulated annealing and multiple dimensioned ridge ripple |
CN110780351A (en) * | 2018-07-31 | 2020-02-11 | 中国石油化工股份有限公司 | Longitudinal wave and converted wave prestack joint inversion method and system |
CN112180440A (en) * | 2019-07-03 | 2021-01-05 | 中国石油天然气集团有限公司 | AVO characteristic analysis-based prestack stochastic inversion method and system |
Citations (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US20030062902A1 (en) * | 2001-08-27 | 2003-04-03 | Stolarczyk Larry G. | Method and system for radio-imaging underground geologic structures |
CN102062873A (en) * | 2009-11-13 | 2011-05-18 | 中国石油化工股份有限公司 | Method for matching longitudinal and transverse waves |
CN102870009A (en) * | 2010-03-12 | 2013-01-09 | 地球物理维里达斯集团(美国)有限公司 | Methods and systems for performing azimuthal simultaneous elatic inversion |
-
2014
- 2014-04-11 CN CN201410145440.6A patent/CN104977609A/en active Pending
Patent Citations (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US20030062902A1 (en) * | 2001-08-27 | 2003-04-03 | Stolarczyk Larry G. | Method and system for radio-imaging underground geologic structures |
CN102062873A (en) * | 2009-11-13 | 2011-05-18 | 中国石油化工股份有限公司 | Method for matching longitudinal and transverse waves |
CN102870009A (en) * | 2010-03-12 | 2013-01-09 | 地球物理维里达斯集团(美国)有限公司 | Methods and systems for performing azimuthal simultaneous elatic inversion |
Non-Patent Citations (3)
Title |
---|
ZHIWEN DENG 等: "PP & PS wave image and joint stochastic prestack inversion in the same PP time scale domain", 《INTERNATIONAL PETROLEUM TECHNOLOGY CONFERENCE》 * |
司马立强: "《测井地质应用技术》", 31 October 2002, 石油工业出版社 * |
邢春颖 等: "三维三分量地震资料试处理", 《石油地球物理勘探》 * |
Cited By (8)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN105807316A (en) * | 2016-04-25 | 2016-07-27 | 吉林大学 | Surface observation microseism speed model correcting method based on amplitude stack |
CN108333626A (en) * | 2018-01-10 | 2018-07-27 | 成都理工大学 | A kind of genetic algorithm Optimum Impedance Inversion Method based on best retention strategy |
CN108333626B (en) * | 2018-01-10 | 2019-06-18 | 成都理工大学 | A kind of genetic algorithm Optimum Impedance Inversion Method based on best retention strategy |
CN108363099A (en) * | 2018-02-11 | 2018-08-03 | 中国石油化工股份有限公司 | The longitudinal and shear wave correlation technique analyzed based on simulated annealing and multiple dimensioned ridge ripple |
CN110780351A (en) * | 2018-07-31 | 2020-02-11 | 中国石油化工股份有限公司 | Longitudinal wave and converted wave prestack joint inversion method and system |
CN110780351B (en) * | 2018-07-31 | 2021-11-23 | 中国石油化工股份有限公司 | Longitudinal wave and converted wave prestack joint inversion method and system |
CN112180440A (en) * | 2019-07-03 | 2021-01-05 | 中国石油天然气集团有限公司 | AVO characteristic analysis-based prestack stochastic inversion method and system |
CN112180440B (en) * | 2019-07-03 | 2023-05-26 | 中国石油天然气集团有限公司 | Pre-stack random inversion method and system based on AVO feature analysis |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
EP3028071B1 (en) | Method and device for the generation and application of anisotropic elastic parameters in horizontal transverse isotropic (hti) media | |
US6826484B2 (en) | 3D prestack time migration method | |
KR101548976B1 (en) | Estimation of soil properties using waveforms of seismic surface waves | |
US6839658B2 (en) | Seismic processing with general non-hyperbolic travel-time corrections | |
US7315783B2 (en) | Traveltime calculation in three dimensional transversely isotropic (3D TI) media by the fast marching method | |
CN103995288B (en) | Gauss beam prestack depth migration method and device | |
CA2899796C (en) | Seismic velocity model updating and imaging with elastic wave imaging | |
CN104407378B (en) | Anisotropy parameter inversion method and device | |
US20100118652A1 (en) | Determination of depth moveout and of residual radii of curvature in the common angle domain | |
Igel et al. | SH‐wave propagation in the whole mantle using high‐order finite differences | |
CN106556861B (en) | A kind of azimuthal AVO inversion method based on Omnibearing earthquake auto data | |
AU2007229430A1 (en) | Method for multi-azimuth prestack time migration for general heterogeneous, anisotropic media | |
Witten et al. | Extended wave-equation imaging conditions for passive seismic data | |
CN104977609A (en) | Prestack longitudinal wave and transverse wave combined inversion method based on rapid simulated annealing | |
Górszczyk et al. | Graph‐space optimal transport concept for time‐domain full‐waveform inversion of ocean‐bottom seismometer data: Nankai trough velocity structure reconstructed from a 1D model | |
EP3067718A1 (en) | Boundary layer tomography method and device | |
Gong et al. | Combined migration velocity model-building and its application in tunnel seismic prediction | |
CN103472482B (en) | Multi-wave seismic data time domain matching process and system based on gene order system | |
Zhu et al. | Amplitude and phase versus angle for elastic wide-angle reflections in the τ‐p domain | |
US20150149093A1 (en) | Method and system for efficient extrapolation of a combined source-and-receiver wavefield | |
Ge et al. | Estimation of weak anisotropy parameters and fracture density of asymmetric fractures in monoclinic medium | |
Kim et al. | Prestack elastic generalized-screen migration for multicomponent data | |
CN105988134A (en) | Multi-wave synthesis record calibration method | |
Wang et al. | 3D isotropic elastic reverse time migration using signed magnitudes of elastic data components | |
Lee | Efficient ray tracing algorithms based on wavefront construction and model based interpolation method |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
C06 | Publication | ||
PB01 | Publication | ||
C10 | Entry into substantive examination | ||
SE01 | Entry into force of request for substantive examination | ||
RJ01 | Rejection of invention patent application after publication | ||
RJ01 | Rejection of invention patent application after publication |
Application publication date: 20151014 |