CN106950598B - A kind of migration velocity field method for evaluating reliability - Google Patents

A kind of migration velocity field method for evaluating reliability Download PDF

Info

Publication number
CN106950598B
CN106950598B CN201710187173.2A CN201710187173A CN106950598B CN 106950598 B CN106950598 B CN 106950598B CN 201710187173 A CN201710187173 A CN 201710187173A CN 106950598 B CN106950598 B CN 106950598B
Authority
CN
China
Prior art keywords
velocity field
field
scale
migration velocity
topological
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.)
Active
Application number
CN201710187173.2A
Other languages
Chinese (zh)
Other versions
CN106950598A (en
Inventor
符力耘
吴玉
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Institute of Geology and Geophysics of CAS
Original Assignee
Institute of Geology and Geophysics of CAS
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by Institute of Geology and Geophysics of CAS filed Critical Institute of Geology and Geophysics of CAS
Priority to CN201710187173.2A priority Critical patent/CN106950598B/en
Publication of CN106950598A publication Critical patent/CN106950598A/en
Application granted granted Critical
Publication of CN106950598B publication Critical patent/CN106950598B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • GPHYSICS
    • G01MEASURING; TESTING
    • G01VGEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/28Processing seismic data, e.g. for interpretation or for event detection
    • G01V1/30Analysis

Landscapes

  • Engineering & Computer Science (AREA)
  • Remote Sensing (AREA)
  • Physics & Mathematics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Acoustics & Sound (AREA)
  • Environmental & Geological Engineering (AREA)
  • Geology (AREA)
  • General Life Sciences & Earth Sciences (AREA)
  • General Physics & Mathematics (AREA)
  • Geophysics (AREA)
  • Geophysics And Detection Of Objects (AREA)

Abstract

The present invention relates to a kind of migration velocity field method for evaluating reliability, comprising the following steps: 1) by seimic wave propagation operator by wave field from earth's surface continuation to underground;2) continuation wave field is imaged by the way that operator is imaged;3) approximate processing is carried out to propogator matrix, so that propogator matrix only retains the main information of velocity field;4) calculate approximate processing after propogator matrix topological scale, by topological scale portray in migration velocity field structure;5) migration velocity field of different resolution scale is calculated, for different migration velocity fields, repeat implementation steps 1)~4), the topological scale curve of each migration velocity field is obtained, the degree of closeness of the topological scale curve of different migration velocity fields is to characterize the degree of closeness of migration result.

Description

A kind of migration velocity field method for evaluating reliability
Technical field
The present invention relates to a kind of migration velocity field method for evaluating reliability, belong to seismic exploration technique field.
Background technique
With gradually going deep into for seismic prospecting, seismic migration imaging is gradually transitioned into Depth Domain from time-domain.Speed and depth The uncertainty of degree makes depth migration very sensitive to velocity field.In order to obtain high-precision migration imaging as a result, it is desirable to right The reliability of migration velocity field is analyzed and evaluated.
The quality that migration velocity field reliability mainly passes through migration result at present is evaluated, and can be given birth to by migration before stack At various common imaging gathers, mainly offset away from domain common imaging gather (ODCIG) and angle domain common image gathers (ADCIG), according to the relationship between common imaging gather and migration velocity, when image gather is evened up or is focused, It is considered that migration velocity is accurately and reliably.But this method calculation amount is very big, needs to carry out pre-stack depth migration generation Common imaging gather.
Summary of the invention
It can be before migration imaging to the reliability of velocity field in view of the above-mentioned problems, the object of the present invention is to provide one kind Carry out the migration velocity field method for evaluating reliability of Fast Evaluation.
To achieve the above object, the invention adopts the following technical scheme: a kind of migration velocity field method for evaluating reliability, packet Include following steps: 1) by seimic wave propagation operator by wave field from earth's surface continuation to underground;2) by imaging operator to continuation wave Field is imaged;3) approximate processing is carried out to propogator matrix, so that propogator matrix only retains the main information of velocity field;4) it calculates The topological scale of propogator matrix after approximate processing, by topological scale portray in migration velocity field structure;5) it calculates not Implementation steps 1 are repeated for different migration velocity fields with the migration velocity field of resolution-scale)~4), obtain each offset speed The topological scale curve of field is spent, the degree of closeness of the topological scale curve of different migration velocity fields is to characterize approaching for migration result Degree.
As follows to the process of wave field progress continuation in the step 1): the earthquake record of earth's surface can be expressed as p0(x, z= 0, ω), wave field downward continuation Δ z can be indicated are as follows:
p1(x, z=Δ z, ω)=A0p0(x, z=0, ω) (1)
In formula, A0Indicate the propogator matrix of first layer;p0Indicate surface seismic data;X indicates lateral distance;Z indicates deep Degree;ω indicates frequency;Δ z indicates step size;p1Indicate the seismic data of underground first layer;(1) formula is applied multiple times by wave field From earth's surface downward continuation n-layer:
pn(x, z=n Δ z, ω)=An-1…Ai…A0p0(x, z=0, ω) (2)
In formula, AiIndicate the propogator matrix of i+1 layer, it is related with i-th layer of speed;pnIndicate the earthquake of underground n-th layer Data.
The process that continuation wave field is imaged in the step 2) is as follows:
Imaging operator R is acted on into continuation wave field (2) formula, it may be assumed that
In=Rpn (3)
In formula: InIndicate that reflectivity is imaged in n-th layer.
Propogator matrix A=An-1...A0, to the approximate processing process of propogator matrix in the step 3) are as follows:
Adiag=diag (A) (5)
In formula, diag indicates to take the diagonal element of matrix.
In the step 4), calculates topological scale and uses elimination trend Fluctuation Method, calculating process is as follows:
To a discrete series Xn, n=1,2 ..., N, mean value isCorresponding cumulative departure sequence be y (m), m=1, 2 ..., M,
Sequences y (m) is divided into the M/l minizone that length is l, removes its background trend between each cell, then The fluctuation root mean square F (l) of cumulative departure sequences y (m) is calculated,
Wherein ylIndicate the background trend of first of minizone.It is full between F (l) and l for topological scale invariance medium The following relationship of foot,
F(l)∝lα, (8)
Wherein α indicates topological scale;Under log-log coordinate, topological scale is calculated by fit slope;For propagating square Battle array Adiag, each column regard a discrete series as, a topological scale can be calculated using the above method, final to obtain One topological scale curve.
The invention adopts the above technical scheme, which has the following advantages: the present invention can be (folded to different disposal method Acceleration analysis, migration velocity modeling, full waveform inversion etc.) obtain rate pattern carry out Analysis of Topological Structure, offset at Reliability evaluation is carried out to various rate patterns as before, topological scale curve differs smaller its migration result of rate pattern and gets over It is close.
Detailed description of the invention
Fig. 1 is three kinds of more representational complex fault block rate patterns, wherein figure (a) is original offset velocity field, figure It (b) is longitudinal smoothed offset velocity field, figure (c) is lateral smoothed offset velocity field;
Fig. 2 is topological scale curve corresponding to three kinds of migration velocity fields;
Fig. 3 is three kinds of complex fault block imaging sections, wherein figure (a) is original offset velocity field imaging section, and figure (b) is Longitudinal smoothed offset velocity field imaging section, figure (c) is lateral smoothed offset velocity field imaging section;
Fig. 4 is the migration velocity field changed original offset velocity field after shallow-layer velocity amplitude (increasing 300m/s)
Fig. 5 is that original offset velocity field and topological scale corresponding to the migration velocity field after change shallow-layer velocity amplitude are bent Line;
Fig. 6 is the migration velocity field imaging section changed after shallow-layer velocity amplitude.
Specific embodiment
The present invention is described in detail below with reference to the accompanying drawings and embodiments.
The invention proposes a kind of migration velocity field method for evaluating reliability, comprising the following steps:
1) pass through seimic wave propagation operator by wave field from earth's surface continuation to underground, detailed process is as follows:
The earthquake record of earth's surface can be expressed as p0(x, z=0, ω) can indicate wave field downward continuation Δ z are as follows:
p1(x, z=Δ z, ω)=A0p0(x, z=0, ω) (1)
In formula, A0Indicate the propogator matrix of first layer;p0Indicate surface seismic data;X indicates lateral distance;Z indicates deep Degree;ω indicates frequency;Δ z indicates step size;p1Indicate the seismic data of underground first layer.
(1) formula is applied multiple times can be by wave field from earth's surface downward continuation n-layer:
pn(x, z=n Δ z, ω)=An-1…Ai…A0p0(x, z=0, ω) (2)
In formula, AiIndicate the propogator matrix of i+1 layer, it is related with i-th layer of speed;pnIndicate the earthquake of underground n-th layer Data.
2) continuation wave field is imaged by the way that operator is imaged, detailed process is as follows:
Imaging operator R is acted on into continuation wave field (2) formula, it may be assumed that
In=Rpn (3)
In formula: InIndicate that reflectivity is imaged in n-th layer.
3) approximate processing is carried out to propogator matrix, so that propogator matrix only retains the main information of velocity field, detailed process It is as follows:
Formula (1) (2) (3) above-mentioned establishes the connection between migration velocity field and imaging results, can by above-mentioned formula To find out, velocity field is that final imaging results are influenced by propogator matrix, between lower surface analysis propogator matrix and velocity field Relationship, for n-th layer medium, wave field is from pn(x, ω) continuation is to pn+1(x, ω) can be indicated are as follows:
In formula: k indicates wave number;C (x) indicates n-th layer wavefield velocity;J is imaginary unit.
By (4) formula it is found that matrix AnThe corresponding convolution operator of every a line, frequency domain response is
In order to analyze the influence covered medium and propagate wave field, need to estimate A=An-1...A0.For 2 dimension media, A matrix Size is nx×nx, wherein nxFor lateral offset sampling number, directly calculating An-1...A0It is very time-consuming, and need to compare Big memory.Therefore, approximation can be carried out to propogator matrix, so that propogator matrix retains the main information of velocity field, tool Body way is the diagonal element that can only retain matrix A:
Adiag=diag (A) (5)
In formula, diag indicates to take the diagonal element of matrix;
4) propogator matrix A is calculateddiagTopological scale, by topological scale portray in migration velocity field knot Structure.
Topological scale measurement is zoomed in or out to research object, its form, complexity, scrambling etc. The degree of variation.Topological scale invariance refers to a regional area optional to velocity field, no matter its amplification, diminution or deformation, The various characteristics such as complexity, scrambling do not change.One velocity field kept with topological structure must meet Scale invariance, i.e. migration result are identical.It carries out smooth to migration velocity field or changes an interval velocity, will lead to velocity field Different variations occurs for form, complexity, scrambling, can be with this variation of quantificational expression by topological scale.
Elimination trend Fluctuation Method can effective calculated curve topological scale, calculating process is as follows: giving a discrete series Xn, n=1,2 ..., N, mean value isCorresponding cumulative departure sequence be y (m), m=1,2 ..., M,
Sequences y (m) is divided into the M/l minizone that length is l, removes its background trend between each cell, then The fluctuation root mean square F (l) of cumulative departure sequences y (m) is calculated,
Wherein ylIndicate the background trend of first of minizone.It is full between F (l) and l for topological scale invariance medium The following relationship of foot,
F(l)∝lα, (8)
Wherein α indicates topological scale.Under log-log coordinate, topological scale can be calculated by fit slope.
For propogator matrix Adiag, each column are considered as a discrete series, can calculate using the above method To a topological scale, a topological scale curve may finally be obtained.
5) different resolutions can be obtained by existing methods such as stack velocity analysis, migration velocity modeling, full waveform inversions The migration velocity field of rate scale repeats implementation steps 1 for different migration velocity fields)~4) obtain each migration velocity field Topological scale curve, topological scale curve are used to quantitatively portray the topological structure of velocity field, the topology mark of different migration velocity fields Line of writing music is closer, show in migration velocity field topological structure it is closer, corresponding migration result is also closer.
Technical effect of the invention is illustrated with a specific embodiment below:
1) three kinds of migration velocity fields are inputted
As shown in Figure 1, three kinds of more representational velocity fields of this example selective analysis: (a) original offset velocity field is made For a standard of other two kinds of velocity field reliability evaluations, (b) longitudinal smoothed offset velocity field, (c) lateral smoothed offset speed Spend field.Smoothed out velocity field is regarded as the macroscopical migration velocity field established by general velocity analysis, two kinds of smooth manners Its topological scale of the velocity field of foundation is different, and corresponding imaging precision is also different.
2) three kinds of velocity field topology scale curves are calculated
For the immanent structure feature of quantitative analysis migration velocity field, we calculate the corresponding propagation square of three kinds of velocity fields Battle array and its topological scale curve (as shown in Figure 2).As can be seen that it is longitudinal to velocity field progress smooth, original speed is maintained substantially The topological scale features for spending field illustrate that longitudinal smooth velocity field is compared with raw velocity field, maintain similar immanent structure;And The laterally smooth basic large scale variation for keeping raw velocity field topology scale curve, but destroy its small dimensional variation feature, i.e., Laterally smooth more larger than longitudinal smooth topological scale destruction to raw velocity field, corresponding imaging precision is lower.
3) in order to verify above-mentioned conclusion, suitable frequency (dominant frequency: 30Hz) is chosen, to the two-dimension earthquake number of three kinds of velocity fields According to post-stack migration imaging is carried out, as a result as shown in Figure 3.As it can be seen that since longitudinal smooth velocity field maintains opening up for raw velocity field Feature is flutterred, migration result and original image section are almost the same, deviation very little;And laterally smooth velocity field destroys original speed The small dimensional variation feature of field topology scale, migration result and original image section slightly deviation are spent, breakpoint is mainly reflected in It is unclean with the convergence of small fault block, there is diffracted wave by a small margin.
4) imaging is seriously affected in order to further protrude destruction raw velocity field topology scale, we change shallow-layer One layer of velocity amplitude (increasing 300m/s), to test migration result to the sensibility of a certain interval velocity, as a result as shown in Figure 4.Change Velocity field after becoming an interval velocity is regarded as the very big migration velocity field of complex region error established by general velocity analysis, The serious topological scale for destroying true velocity field.Fig. 5 is the corresponding propogator matrix of velocity field and its topological scale curve of calculating, It indicates that significant changes have occurred in the immanent structure of migration velocity field, serious imaging is caused to misplace.To raw velocity field and change Velocity field after shallow-layer velocity amplitude carries out migration imaging after two-dimension earthquake stacked data, as a result as shown in Figure 6.As it can be seen that due to changing One interval velocity is larger to the immanent structure feature destruction of raw velocity field, and many place wave fields are not restrained, and stratum deformation is serious, Also great changes will take place for the corresponding depth in stratum after imaging.
Therefore, the topological characteristic for changing raw velocity field, is affected to final off-set construction.This case migration imaging example Illustrate the validity of this patent migration velocity method for evaluating reliability.
The various embodiments described above are merely to illustrate the present invention, and wherein the implementation steps etc. of method may be changed, All equivalents and improvement carried out based on the technical solution of the present invention, should not exclude in protection scope of the present invention Except.

Claims (1)

1. a kind of migration velocity field method for evaluating reliability, comprising the following steps:
1) by seimic wave propagation operator by wave field from earth's surface continuation to underground;
2) continuation wave field is imaged by the way that operator is imaged;
3) approximate processing is carried out to propogator matrix, so that propogator matrix only retains the main information of velocity field;
4) calculate approximate processing after propogator matrix topological scale, by topological scale portray in migration velocity field knot Structure;
5) migration velocity field for calculating different resolution scale repeats implementation steps 1 for different migration velocity fields)~4), The topological scale curve of each migration velocity field is obtained, the degree of closeness of the topological scale curve of different migration velocity fields characterizes inclined Move the degree of closeness of result;
It is as follows to the process of wave field progress continuation in the step 1):
The earthquake record of earth's surface can be expressed as p0(x, z=0, ω) can indicate wave field downward continuation Δ z are as follows:
p1(x, z=Δ z, ω)=A0p0(x, z=0, ω) (1)
In formula, A0Indicate the propogator matrix of first layer;p0Indicate surface seismic data;X indicates lateral distance;Z indicates depth;ω Indicate frequency;Δ z indicates step size;p1Indicate the seismic data of underground first layer;
(1) formula is applied multiple times by wave field from earth's surface downward continuation n-layer:
pn(x, z=n Δ z, ω)=An-1…Ai…A0p0(x, z=0, ω) (2)
In formula, AiIndicate the propogator matrix of i+1 layer, it is related with i-th layer of speed;nIndicate the seismic data of underground n-th layer;
The process that continuation wave field is imaged in the step 2) is as follows:
Imaging operator R is acted on into continuation wave field (2) formula, it may be assumed that
In=Rpn (3)
In formula: InIndicate that reflectivity is imaged in n-th layer;
Propogator matrix A=An-1...A0,
To the approximate processing process of propogator matrix in the step 3) are as follows:
Adiag=diag (A) (5)
In formula, diag indicates to take the diagonal element of matrix
In the step 4), calculates topological scale and uses elimination trend Fluctuation Method, calculating process is as follows:
To a discrete series Xn, n=1,2 ..., N, mean value isCorresponding cumulative departure sequence be y (m), m=1,2 ..., M,
Sequences y (m) is divided into the M/ minizone that length is l, its background trend is removed between each cell, then calculates tired The fluctuation root mean square F (l) of product biased sequence y (m),
Wherein ylThe background trend for indicating first of minizone meets as follows for topological scale invariance medium, between F (l) and l Relationship,
F(l)∝lα, (8)
Wherein α indicates topological scale;Under log-log coordinate, topological scale is calculated by fit slope;
For propogator matrix Adiag, each column regard a discrete series as, a topology can be calculated using the above method Scale, it is final to obtain a topological scale curve.
CN201710187173.2A 2017-03-27 2017-03-27 A kind of migration velocity field method for evaluating reliability Active CN106950598B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201710187173.2A CN106950598B (en) 2017-03-27 2017-03-27 A kind of migration velocity field method for evaluating reliability

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201710187173.2A CN106950598B (en) 2017-03-27 2017-03-27 A kind of migration velocity field method for evaluating reliability

Publications (2)

Publication Number Publication Date
CN106950598A CN106950598A (en) 2017-07-14
CN106950598B true CN106950598B (en) 2019-01-29

Family

ID=59473792

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201710187173.2A Active CN106950598B (en) 2017-03-27 2017-03-27 A kind of migration velocity field method for evaluating reliability

Country Status (1)

Country Link
CN (1) CN106950598B (en)

Families Citing this family (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN108680968B (en) * 2018-07-24 2020-01-07 中国石油天然气集团有限公司 Evaluation method and device for seismic exploration data acquisition observation system in complex structural area
CN109188513B (en) * 2018-09-30 2020-03-10 中国石油天然气股份有限公司 Method and device for generating depth domain data volume and storage medium

Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101082676A (en) * 2007-07-11 2007-12-05 成都理工大学 Earthquake post-stack forward method of undulating surface
CN101276001A (en) * 2008-04-25 2008-10-01 符力耘 Underground non-uniform medium seismic investigation complexity quantitative evaluating method
CN101609167A (en) * 2009-07-17 2009-12-23 中国石化集团胜利石油管理局 Crosshole seismic wave equation pre stack depth migration formation method based on relief surface
CN102890290A (en) * 2012-09-25 2013-01-23 中国石油天然气股份有限公司 Prestack depth migration method under condition of undulating surface
CN104570124A (en) * 2013-10-29 2015-04-29 中国石油化工股份有限公司 Continuation imaging method suitable for cross-well seismic large-angle reflection conditions
CN104678440A (en) * 2015-02-15 2015-06-03 山东科技大学 Well-constrained two-dimensional seismic variable velocity field nonlinear error correction method
CN105807315A (en) * 2016-03-14 2016-07-27 中国石油大学(华东) Elastic vector reverse time migration imaging method

Patent Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101082676A (en) * 2007-07-11 2007-12-05 成都理工大学 Earthquake post-stack forward method of undulating surface
CN101276001A (en) * 2008-04-25 2008-10-01 符力耘 Underground non-uniform medium seismic investigation complexity quantitative evaluating method
CN101609167A (en) * 2009-07-17 2009-12-23 中国石化集团胜利石油管理局 Crosshole seismic wave equation pre stack depth migration formation method based on relief surface
CN102890290A (en) * 2012-09-25 2013-01-23 中国石油天然气股份有限公司 Prestack depth migration method under condition of undulating surface
CN104570124A (en) * 2013-10-29 2015-04-29 中国石油化工股份有限公司 Continuation imaging method suitable for cross-well seismic large-angle reflection conditions
CN104678440A (en) * 2015-02-15 2015-06-03 山东科技大学 Well-constrained two-dimensional seismic variable velocity field nonlinear error correction method
CN105807315A (en) * 2016-03-14 2016-07-27 中国石油大学(华东) Elastic vector reverse time migration imaging method

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
New seismic attribute:Fractal scaling exponent based on gray detrended fluctuation analysis;Huang Ya-Ping et al;《APPLIED GEOPHYSICS》;20160930;第12卷(第3期);343-352 *
改进的消除波动趋势分析法;庞宇磊 等;《科学技术与工程》;20100131;第10卷(第3期);634-637,642 *

Also Published As

Publication number Publication date
CN106950598A (en) 2017-07-14

Similar Documents

Publication Publication Date Title
CN103293552B (en) A kind of inversion method of Prestack seismic data and system
CN105277978B (en) A kind of method and device for determining near-surface velocity model
CN108254780A (en) A kind of microseism positioning and anisotropic velocity structure tomographic imaging method
CN106526674A (en) Three-dimensional full waveform inversion energy weighted gradient preprocessing method
CN107329171A (en) Depth Domain reservoir seismic inversion method and device
CN109459787B (en) coal mine underground structure imaging method and system based on seismic channel wave full-waveform inversion
CN106950598B (en) A kind of migration velocity field method for evaluating reliability
CN106249297A (en) Fracturing microseism seismic source location method and system based on Signal estimation
CN103604836B (en) A kind of method and apparatus of measuring gas hydrates reservoir saturation degree
Chen et al. A compact program for 3D passive seismic source‐location imaging
Guidio et al. Passive seismic inversion of SH wave input motions in a truncated domain
CN109490947B (en) High-temperature medium seismic wave propagation simulation method
CN113158315A (en) Rock-soil body parameter three-dimensional non-stationary condition random field modeling method based on static cone penetration data
Putz-Perrier et al. Spatial distribution of brittle strain in layered sequences
CN112100906A (en) Data-driven large-scale density modeling method, computing device and storage medium
CN105093300A (en) Geologic body boundary identification method and apparatus
CN113219531A (en) Method and device for identifying gas-water distribution of tight sandstone
CN103901472B (en) Frequency domain forward modeling method and device
CN112257241B (en) Triangular net Fresnel time difference tomography inversion method
CN114200538B (en) Ground stress direction determining method, device, storage medium and computer equipment
CN113552624B (en) Porosity prediction method and device
CN114624779A (en) Pre-stack multi-parameter inversion method for balanced model constraint
US20200233111A1 (en) Estimation of Reservoir Flow Properties From Seismic Data
Rosenblad et al. Potential phase unwrapping errors associated with SASW measurements at soft-over-stiff sites
CN105259577A (en) Method and device for determining angle information of formation boundary

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant