CN108845355A - Seismic migration imaging method and device - Google Patents
Seismic migration imaging method and device Download PDFInfo
- Publication number
- CN108845355A CN108845355A CN201811128030.5A CN201811128030A CN108845355A CN 108845355 A CN108845355 A CN 108845355A CN 201811128030 A CN201811128030 A CN 201811128030A CN 108845355 A CN108845355 A CN 108845355A
- Authority
- CN
- China
- Prior art keywords
- wave field
- indicate
- model
- current
- field residual
- 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
- 238000003384 imaging method Methods 0.000 title claims abstract description 92
- 230000005012 migration Effects 0.000 title claims abstract description 89
- 238000013508 migration Methods 0.000 title claims abstract description 89
- 238000000034 method Methods 0.000 claims abstract description 49
- 238000002310 reflectometry Methods 0.000 claims abstract description 42
- 238000002939 conjugate gradient method Methods 0.000 claims abstract description 16
- 230000006870 function Effects 0.000 claims description 90
- 238000004088 simulation Methods 0.000 claims description 3
- 230000001419 dependent effect Effects 0.000 claims 1
- 238000007689 inspection Methods 0.000 claims 1
- 230000008569 process Effects 0.000 abstract description 12
- 230000000694 effects Effects 0.000 abstract description 9
- 238000004364 calculation method Methods 0.000 abstract description 7
- 238000004422 calculation algorithm Methods 0.000 abstract description 3
- 238000010586 diagram Methods 0.000 description 8
- 238000005516 engineering process Methods 0.000 description 3
- 238000004458 analytical method Methods 0.000 description 2
- 238000009434 installation Methods 0.000 description 2
- 230000009286 beneficial effect Effects 0.000 description 1
- 238000004590 computer program Methods 0.000 description 1
- 238000001514 detection method Methods 0.000 description 1
- 235000013399 edible fruits Nutrition 0.000 description 1
- 238000001727 in vivo Methods 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 230000009467 reduction Effects 0.000 description 1
- 230000017105 transposition Effects 0.000 description 1
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
- G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
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 provides a kind of seismic migration imaging method and devices, including：Wave field residual error objective function is established based on seismic observation data and analogue data to be solved；Current wave field residual values are determined in conjunction with current reflectance model and wave field residual error objective function；Offset method based on time shift image-forming condition deviates current wave field residual values, obtains target function gradient；Updated reflectivity model is determined using conjugate gradient method based on target function gradient；It is updated iteration using updated reflectivity model as current reflectance model, until reaching default the number of iterations, obtains target reflection factor model, and then obtain migration imaging result.The process employs least-squares algorithms, effectively increase the precision of seismic migration imaging, imaging effect is good, time shift image-forming condition is used during migration imaging, calculation amount is small, the efficiency for improving migration imaging alleviates traditional migration imaging inefficiency, and the technical problem that imaging effect is bad.
Description
Technical field
The present invention relates to the technical fields of seismic prospecting, more particularly, to a kind of seismic migration imaging method and device.
Background technique
Method of seismic prospecting is the important means that the mankind obtain underground space information, wherein seismic migration imaging is earthquake
The committed step of Data processing.Offset can make back wave accurately playback, diffracted wave convergence, to intuitively show subterranean
The true form made.
The seismic imaging that migration and imaging techniques traditional at present generally use under single-point scattering approximation is theoretical, is mainly base
In the matching analysis of source wavefield and geophone station wave field, forward and reverse is carried out respectively to source wavefield and geophone station wave field first
Continuation reapplies suitable image-forming condition and reflecting surface is imaged.Suitable image-forming condition for migration velocity analysis and at
As effect has great influence.It is to be protected in image-forming condition using a kind of wide image-forming condition at present that sky, which moves image-forming condition,
Offset distance information is stayed, realizes the angle domain imaging of wave equation migration.But sky moves image-forming condition and needs huge cross-correlation calculation
Amount, computational efficiency is lower, is unfavorable for the processing of real data；And traditional migration imaging is easy to produce noise and offset is false
As.
To sum up, traditional migration imaging inefficiency, and imaging effect is bad.
Summary of the invention
In view of this, the purpose of the present invention is to provide a kind of seismic migration imaging method and device, it is traditional to alleviate
Migration imaging inefficiency, and the technical problem that imaging effect is bad.
In a first aspect, the embodiment of the invention provides a kind of seismic migration imaging methods, including：
Seismic observation data and analogue data to be solved based on pending area establish wave field residual error objective function；
Current wave field residual values are determined in conjunction with current reflectance model and the wave field residual error objective function, wherein institute
Current reflectance model is stated for determining the analogue data to be solved；
Offset method based on time shift image-forming condition deviates the current wave field residual values, obtains target with determination
Functional gradient；
Updated reflectivity model is determined using conjugate gradient method based on the target function gradient；
It is updated iteration using the updated reflectivity model as the current reflectance model, until reaching
To default the number of iterations, target reflection factor model is obtained, and then obtains migration imaging result.
With reference to first aspect, the embodiment of the invention provides the first possible embodiments of first aspect, wherein base
Establishing wave field residual error objective function in the seismic observation data of pending area and analogue data to be solved includes：
Obtain the seismic observation data and migration velocity body of the pending area；
The wave field residual error target between the seismic observation data and the analogue data to be solved is established by L2 norm
Function.
With reference to first aspect, the embodiment of the invention provides second of possible embodiments of first aspect, wherein knot
It closes current reflectance model and the wave field residual error objective function determines that current wave field residual values include：
Formula d is calculated according to analogue datamod=LmkCalculate the analogue data to be solved, wherein dmodIndicate it is described to
Analogue data is solved, L indicates the seismic wave field forward operator that velocity information is depended under Born approximate condition, mkWork as described in expression
Front-reflection Modulus Model；
According to the calculating formula of the wave field residual error objective functionWork as prewave described in calculating
Field residual values, wherein E (mk) indicate the current wave field residual values, dmodIndicate the analogue data to be solved, dobsIndicate institute
State seismic observation data.
With reference to first aspect, the embodiment of the invention provides the third possible embodiments of first aspect, wherein base
The current wave field residual values are deviated in the offset method of time shift image-forming condition, target function gradient packet is obtained with determination
It includes：
Geophone station wave field is subjected to backward extension by the migration velocity body, and passes through the migration velocity body for focus
Wave field carries out positive continuation；
Based on time shift image-forming condition to the geophone station wave field after backward extension and the source wavefield after positive continuation carry out at
Picture obtains the migrated section about time shift amount；
Using the migrated section of zero moment as the target function gradient in the migrated section about time shift amount.
With reference to first aspect, the embodiment of the invention provides the 4th kind of possible embodiments of first aspect, wherein base
The geophone station wave field after the backward extension and the source wavefield after the positive continuation are imaged in time shift image-forming condition
Including：
According to the imaging formula R (m based on time shift image-forming conditionk, t, τ) and=Pr(mk,t+τ)*Ps(mk, t- τ) and it is imaged,
Obtain the migrated section about time shift amount, wherein R (m, t, τ) indicates that the migrated section about time shift amount, * indicate
Time-domain cross-correlation, t indicate the time, and τ indicates the time shift amount in time shift image-forming condition, PrDetection after indicating the backward extension
Point wave field, PsSource wavefield after indicating the positive continuation, mkIndicate the current reflectance model.
With reference to first aspect, the embodiment of the invention provides the 5th kind of possible embodiments of first aspect, wherein base
Determine that updated reflectivity model includes using conjugate gradient method in the target function gradient：
It is determined based on the target function gradient and updates step-length and more new direction；
Formula m is updated according to reflectivity modelk+1=mk+αkdkCalculate the updated reflectivity model, wherein
mk+1Indicate the updated reflectivity model, mkIndicate current reflectance model, αkIndicate the update step-length, dkTable
Show the more new direction.
With reference to first aspect, the embodiment of the invention provides the 6th kind of possible embodiments of first aspect, wherein base
Update step-length is determined in the target function gradient and more new direction includes：
According to more new direction formulaDetermine the more new direction, whereindkIndicate the more new direction, dk-1Table
Show a more new direction, g (mk) indicate in the current reflectance model mkUnder target function gradient, g (mk-1) indicate
Upper reflectivity model mk-1Under target function gradient；
The update step-length is determined according to step-length condition is updated, wherein the update step-length condition is：When meeting E (mk+αkdk)≤E(mk)+αkcg(mk)dkWhen, then αk=αk-1, otherwise, αk=λ αk-1, αkIndicate the update step-length, αk-1Indicate upper one
Update step-length, E (mk+αkdk) indicate that in more new direction and update step-length be respectively dkAnd αkWhen, updated reflectivity model
Under wave field residual values, g (mk) indicate in the current reflectance model mkUnder target function gradient, 0<c<1,0<λ<1.
Second aspect, the embodiment of the invention also provides a kind of seismic migration imaging devices, including：
Establish module, for based on pending area seismic observation data and analogue data to be solved establish wave field residual error
Objective function；
First determining module, for working as prewave in conjunction with current reflectance model and wave field residual error objective function determination
Field residual values, wherein the current reflectance model is for determining the analogue data to be solved；
Offset module deviates the current wave field residual values for the offset method based on time shift image-forming condition,
Target function gradient is obtained with determination；
Second determining module, for determining that updated reflection is using conjugate gradient method based on the target function gradient
Exponential model；
Update iteration module, for using the updated reflectivity model as the current reflectance model into
Row updates iteration, until reaching default the number of iterations, obtains target reflection factor model, and then obtain migration imaging result.
In conjunction with second aspect, the embodiment of the invention provides the first possible embodiments of second aspect, wherein institute
It states and establishes module and include：
Acquiring unit, for obtaining the seismic observation data and migration velocity body of the pending area；
Unit is established, for establishing between the seismic observation data and the analogue data to be solved by L2 norm
Wave field residual error objective function.
In conjunction with second aspect, the embodiment of the invention provides second of possible embodiments of second aspect, wherein institute
Stating the first determining module includes：
First computing unit, for calculating formula d according to analogue datamod=LmkThe analogue data to be solved is calculated,
Wherein, dmodIndicate the analogue data to be solved, L indicates depending on the seismic wave field of velocity information just under Born approximate condition
Calculation, mkIndicate the current reflectance model；
Second computing unit, for the calculating formula according to the wave field residual error objective function
Calculate the current wave field residual values, wherein E (mk) indicate the current wave field residual values, dmodIndicate the simulation to be solved
Data, dobsIndicate the seismic observation data.
The embodiment of the present invention brings following beneficial effect：
In the present embodiment, first to establish wave field residual for the seismic observation data based on pending area and analogue data to be solved
Poor objective function, and then current reflectance model and wave field residual error objective function is combined to determine current wave field residual values, then
Offset method based on time shift image-forming condition deviates current wave field residual values, obtains target function gradient with determination, into
One step determines updated reflectivity model using conjugate gradient method based on target function gradient, finally by updated reflection
Modulus Model is updated iteration as current reflectance model, until reaching default the number of iterations, obtains target reflection system
Exponential model, and then obtain migration imaging result.As can be seen from the above description, seismic migration imaging method of the invention uses most
Small square law effectively increases the precision of seismic migration imaging, and imaging effect is good, in addition, when using during migration imaging
Image-forming condition is moved, calculation amount is small, improves the efficiency of migration imaging, alleviate traditional migration imaging inefficiency, and at
As ineffective technical problem.
Other features and advantages of the present invention will illustrate in the following description, also, partly become from specification
It obtains it is clear that understand through the implementation of the invention.The objectives and other advantages of the invention are in specification, claims
And specifically noted structure is achieved and obtained in attached drawing.
To enable the above objects, features and advantages of the present invention to be clearer and more comprehensible, preferred embodiment is cited below particularly, and cooperate
Appended attached drawing, is described in detail below.
Detailed description of the invention
It, below will be to specific in order to illustrate more clearly of the specific embodiment of the invention or technical solution in the prior art
Embodiment or attached drawing needed to be used in the description of the prior art be briefly described, it should be apparent that, it is described below
Attached drawing is some embodiments of the present invention, for those of ordinary skill in the art, before not making the creative labor
It puts, is also possible to obtain other drawings based on these drawings.
Fig. 1 is a kind of flow chart of seismic migration imaging method provided in an embodiment of the present invention；
Fig. 2 is that the seismic observation data provided in an embodiment of the present invention based on pending area and analogue data to be solved are built
The method flow diagram of vertical wave field residual error objective function；
Fig. 3 is that the offset method provided in an embodiment of the present invention based on time shift image-forming condition carries out current wave field residual values
The method flow diagram of offset；
Fig. 4 is the schematic diagram of common imaging gather provided in an embodiment of the present invention；
Fig. 5 determines updated reflection using conjugate gradient method based on target function gradient to be provided in an embodiment of the present invention
The method flow diagram of Modulus Model；
Fig. 6 (a) is the schematic diagram for the migration imaging result that seismic migration imaging method provided in an embodiment of the present invention obtains；
Fig. 6 (b) is the schematic diagram for the migration imaging result that conventional migration technique imaging method provided in an embodiment of the present invention obtains；
Fig. 7 is a kind of schematic diagram of seismic migration imaging device provided in an embodiment of the present invention.
Specific embodiment
In order to make the object, technical scheme and advantages of the embodiment of the invention clearer, below in conjunction with attached drawing to the present invention
Technical solution be clearly and completely described, it is clear that described embodiments are some of the embodiments of the present invention, rather than
Whole embodiments.Based on the embodiments of the present invention, those of ordinary skill in the art are not making creative work premise
Under every other embodiment obtained, shall fall within the protection scope of the present invention.
For convenient for understanding the present embodiment, first to a kind of seismic migration imaging side disclosed in the embodiment of the present invention
Method describes in detail.
Embodiment one：
According to embodiments of the present invention, a kind of embodiment of seismic migration imaging method is provided, it should be noted that attached
The step of process of figure illustrates can execute in a computer system such as a set of computer executable instructions, though also,
So logical order is shown in flow charts, but in some cases, it can be to be different from shown by sequence execution herein
Or the step of description.
Fig. 1 is a kind of flow chart of seismic migration imaging method according to an embodiment of the present invention, as shown in Figure 1, this method
Include the following steps：
Step S102, seismic observation data and analogue data to be solved based on pending area establish wave field residual error target
Function；
In embodiments of the present invention, pending area is the region of field acquisition.Seismic observation data is to pending district
Domain carries out the data that seismic wave collects, which is given data.Analogue data to be solved is related to reflectivity model,
After obtaining reflectivity model, corresponding analogue data to be solved is assured that, is hereinafter described in detail again.
Step S104 determines current wave field residual values in conjunction with current reflectance model and wave field residual error objective function,
In, current reflectance model is for determining analogue data to be solved；
It, can be according to current anti-if obtaining current reflectance model after obtaining wave field residual error objective function
It penetrates Modulus Model and determines analogue data to be solved, then analogue data to be solved is substituted into wave field residual error objective function, it can be really
Settled preceding wave field residual values.
Step S106, the offset method based on time shift image-forming condition deviate current wave field residual values, with determining
To target function gradient；
After obtaining current wave field residual values, the offset method of time shift image-forming condition is based further on to current wave field residual error
Value is deviated, and is obtained target function gradient, is hereinafter described in detail again.
Step S108 determines updated reflectivity model using conjugate gradient method based on target function gradient；
After obtaining target function gradient, be based further on target function gradient using conjugate gradient method determine it is updated
Reflectivity model.
Updated reflectivity model is updated iteration by step S110, until
Reach default the number of iterations, obtain target reflection factor model, and then obtains migration imaging result.
Specifically, reach the target reflection factor model obtained after default the number of iterations be migration imaging as a result, according to
Target reflection factor model is drawn, and the image of final migration imaging result is obtained.
In the present embodiment, first to establish wave field residual for the seismic observation data based on pending area and analogue data to be solved
Poor objective function, and then current reflectance model and wave field residual error objective function is combined to determine current wave field residual values, then
Offset method based on time shift image-forming condition deviates current wave field residual values, obtains target function gradient with determination, into
One step determines updated reflectivity model using conjugate gradient method based on target function gradient, finally by updated reflection
Modulus Model is updated iteration as current reflectance model, until reaching default the number of iterations, obtains target reflection system
Exponential model, and then obtain migration imaging result.As can be seen from the above description, seismic migration imaging method of the invention uses most
Small square law effectively increases the precision of seismic migration imaging, and imaging effect is good, in addition, when using during migration imaging
Image-forming condition is moved, calculation amount is small, improves the efficiency of migration imaging, alleviate traditional migration imaging inefficiency, and at
As ineffective technical problem.
Above content has carried out brief introduction to seismic migration imaging method of the invention, below to the tool being directed to
Hold in vivo and describes in detail.
In an optional embodiment of the invention, with reference to Fig. 2, seismic observation data based on pending area and to
Solution analogue data is established wave field residual error objective function and is included the following steps：
Step S201 obtains the seismic observation data and migration velocity body of pending area；
Specifically, migration velocity body is handled seismic observation data, and in embodiments of the present invention, offset
Body of velocity is Depth Domain body of velocity.It all include SEGY standard trace header letter in obtained seismic observation data and migration velocity body
Breath.
Step S202 establishes the wave field residual error target between seismic observation data and analogue data to be solved by L2 norm
Function.
After obtaining seismic observation data, can be established by L2 norm seismic observation data and analogue data to be solved it
Between wave field residual error objective function.
Specifically, the wave field residual error objective function established is：
Wherein, E indicates wave field residual error, dmodIndicate analogue data to be solved, dobsIndicate seismic observation data, | | | |2Table
Show L2 norm.In fact, the process is the process so that the seismic observation data of analogue data approaching to reality to be solved.
The process for establishing wave field residual error objective function is described in detail in above content, below to the current wave field of determination
The process of residual values is described in detail.
It is true in conjunction with current reflectance model and wave field residual error objective function in an optional embodiment of the invention
Settled preceding wave field residual values include the following steps：
(1) formula d is calculated according to analogue datamod=LmkCalculate analogue data to be solved, wherein dmodIndicate to be solved
Analogue data, L indicate the seismic wave field forward operator that velocity information is depended under Born approximate condition, mkIndicate current reflective system
Exponential model；
When calculating first time, current reflectance model is initial reflection Modulus Model, the initial reflection Modulus Model
For preset reflectivity model.
(2) according to the calculating formula of wave field residual error objective functionCalculate current wave field residual error
Value, wherein E (mk) indicate current wave field residual values, dmodIndicate analogue data to be solved, dobsIndicate seismic observation data.
As described in (1), analogue data d to be solved is being obtainedmodAfterwards, by analogue data d to be solvedmodSubstitute into wave field residual error mesh
In the calculating formula of scalar functions, it will be able to current wave field residual values be calculated.
The process of the current wave field residual values of determination is described in detail in above content, terraced to objective function is determined below
The process of degree describes in detail.
In an optional embodiment of the invention, with reference to Fig. 3, the offset method based on time shift image-forming condition is to current
Wave field residual values are deviated, to determine that obtaining target function gradient includes the following steps：
Geophone station wave field is carried out backward extension by migration velocity body, and will be shaken by migration velocity body by step S301
Source wave field carries out positive continuation；
Step S302, based on time shift image-forming condition to the geophone station wave field after backward extension and the focus wave after positive continuation
Field is imaged, and the migrated section about time shift amount is obtained；
Specifically, according to the imaging formula R (m based on time shift image-forming conditionk, t, τ) and=Pr(mk,t+τ)*Ps(mk, t- τ) into
Row imaging, obtains the migrated section about time shift amount, wherein R (m, t, τ) indicates the migrated section about time shift amount, when * is indicated
Between domain cross-correlation, t indicate the time, τ indicate time shift image-forming condition in time shift amount, PrGeophone station wave after indicating backward extension
, PsSource wavefield after indicating positive continuation, mkIndicate current reflectance model.
Step S303, using the migrated section of zero moment as target function gradient in the migrated section about time shift amount.
Specifically, g (mk)=LT(dmod-dobs)=R (m, τ, t=0), wherein g (mk) indicate zero moment migrated section,
LTIndicate the transposition of forward operator.
In fact, R obtained in above-mentioned steps S402 (m, t, τ) is 3D data volume, Z axis is depth, and X-axis is time t,
Y-axis is τ, and the progress common imaging gather superposition (as shown in Figure 4) of section when t=0 is taken (when superposition, it is lateral to carry out same depth
Superposition is built up together), migrated section has just been obtained after superposition.Common imaging gather superposition when why taking t=0, be because
When speed is accurate, it can be focused at 0 moment, that is, most strong in 0 moment energy.
The process for determining objective function is described in detail in above content, below to the updated reflection coefficient of determination
The process of model describes in detail.
It is true using conjugate gradient method based on target function gradient with reference to Fig. 5 in an optional embodiment of the invention
Fixed updated reflectivity model includes the following steps：
Step S501 is determined based on target function gradient and is updated step-length and more new direction；
Specifically, in conjugate gradient method, the more new formula of reflectivity model is：
mk+1=mk+αkdk, wherein mk+1Indicate updated reflectivity model, mkIndicate current reflectance model, αk
It indicates to update step-length, dkIndicate more new direction.
Specifically：
(1) according to more new direction formulaDetermine more new direction, whereindkIndicate more new direction, dk-1In expression
One more new direction, g (mk) indicate in current reflectance model mkUnder target function gradient, g (mk-1) indicate in a upper reflection
Modulus Model mk-1Under target function gradient；
(2) according to updating, step-length condition is determining to update step-length, wherein updating step-length condition is：When meeting E (mk+αkdk)≤
E(mk)+αkcg(mk)dkWhen, then αk=αk-1, otherwise, αk=λ αk-1, αkIt indicates to update step-length, αk-1Indicate that upper one updates step-length, E
(mk+αkdk) indicate that in more new direction and update step-length be respectively dkAnd αkWhen, the wave field under updated reflectivity model is residual
Difference, g (mk) indicate in current reflectance model mkUnder target function gradient, 0<c<1,0<λ<1.
It is, more new direction is：Wherein, βkIt can be calculated by following formula
It obtains：
Wherein, (g (mk),g(mk)-g(mk-1)) indicate vector g (mk) and to
Measure g (mk)-g(mk-1) inner product, similarly, this is no longer going to repeat them for the meaning of other parameters.
It updates step-length to be calculated using searching algorithm, sets a larger initial step length αk, reduction factor λ (0<λ<1), than
Example factor c (0<c<1), if it meets condition：E(mk+αkdk)≤E(mk)+αkcg(mk)dk, then αk=αk-1, otherwise αk=λ αk-1,
The meaning of each parameter can the content with reference to documented by other positions in the present invention.
Step S502 updates formula m according to reflectivity modelk+1=mk+αkdkUpdated reflectivity model is calculated,
Wherein, mk+1Indicate updated reflectivity model, mkIndicate current reflectance model, αkIt indicates to update step-length, dkIt indicates
More new direction.
After obtaining updated reflectivity model, using updated reflectivity model as current reflectance mould
Type is updated iteration, until reaching default the number of iterations, obtains target reflection factor model, and then obtain migration imaging knot
Fruit.As shown in Fig. 6 (a), the migration imaging that seismic migration imaging method as according to the present invention obtains is as a result, Fig. 6 (b) is normal
The migration imaging result that rule offset imaging method obtains.By comparison it is found that the resolution for the migration imaging result that Fig. 6 (a) is obtained
Rate significantly improves, and becomes apparent especially for portraying for infrastructure.
The invention proposes a kind of seismic migration imaging method and apparatus, this method constructs Green's letter with time shift image-forming condition
Number carries out wave-field simulation using the linear Born approximate form of seismic wave field, passes through the wave field to analogue data and observation data
Residual error carries out least-squares iteration, obtains least-squares migration imaging section.The present invention is based on time shift image-forming conditions, utilize minimum
Two, which multiply migration technology, carries out earthquake data offset imaging, and since time shift image-forming condition calculates, cost is relatively low, and then alleviates existing
It is calculated as this higher technical problem in technology, improves the efficiency of migration imaging, while least-squares migration algorithm, it can be effective
Improve the precision of seismic migration imaging.
Embodiment two：
The embodiment of the invention also provides a kind of seismic migration imaging device, which is mainly used for holding
Seismic migration imaging method provided by row above content of the embodiment of the present invention, it is inclined to earthquake provided in an embodiment of the present invention below
It moves imaging device and does specific introduction.
Fig. 7 is a kind of schematic diagram of seismic migration imaging device according to an embodiment of the present invention, as shown in fig. 7, the earthquake
Migration imaging device mainly includes：Establish module 11, the first determining module 12, offset module 13, the second determining module 14 and more
New iteration module 15, wherein：
Establish module, for based on pending area seismic observation data and analogue data to be solved establish wave field residual error
Objective function；
First determining module, for combining current reflectance model and wave field residual error objective function to determine that current wave field is residual
Difference, wherein current reflectance model is for determining analogue data to be solved；
Offset module deviates current wave field residual values for the offset method based on time shift image-forming condition, with true
Surely target function gradient is obtained；
Second determining module, for determining updated reflection coefficient mould using conjugate gradient method based on target function gradient
Type；
Iteration module is updated, for being updated updated reflectivity model as current reflectance model repeatedly
In generation, obtains target reflection factor model, and then obtain migration imaging result until reaching default the number of iterations.
In the present embodiment, first to establish wave field residual for the seismic observation data based on pending area and analogue data to be solved
Poor objective function, and then current reflectance model and wave field residual error objective function is combined to determine current wave field residual values, then
Offset method based on time shift image-forming condition deviates current wave field residual values, obtains target function gradient with determination, into
One step determines updated reflectivity model using conjugate gradient method based on target function gradient, finally by updated reflection
Modulus Model is updated iteration as current reflectance model, until reaching default the number of iterations, obtains target reflection system
Exponential model, and then obtain migration imaging result.As can be seen from the above description, seismic migration imaging device of the invention uses most
Small square law effectively increases the precision of seismic migration imaging, and imaging effect is good, in addition, when using during migration imaging
Image-forming condition is moved, calculation amount is small, improves the efficiency of migration imaging, alleviate traditional migration imaging inefficiency, and at
As ineffective technical problem.
Optionally, establishing module includes：
Acquiring unit, for obtaining the seismic observation data and migration velocity body of pending area；
Unit is established, for establishing the wave field residual error between seismic observation data and analogue data to be solved by L2 norm
Objective function.
Optionally, the first determining module includes：
First computing unit, for calculating formula d according to analogue datamod=LmkCalculate analogue data to be solved, wherein
dmodIndicate analogue data to be solved, L indicates the seismic wave field forward operator that velocity information is depended under Born approximate condition, mk
Indicate current reflectance model；
Second computing unit, for the calculating formula according to wave field residual error objective functionMeter
Calculate current wave field residual values, wherein E (mk) indicate current wave field residual values, dmodIndicate analogue data to be solved, dobsIndicate ground
Shake observation data.
Optionally, offset module includes：
Continuation unit for geophone station wave field to be carried out backward extension by migration velocity body, and passes through migration velocity body
Source wavefield is subjected to positive continuation；
Imaging unit, for based on time shift image-forming condition to the geophone station wave field after backward extension and the shake after positive continuation
Source wave field is imaged, and the migrated section about time shift amount is obtained；
Setup unit, for terraced using the migrated section of zero moment as objective function in the migrated section about time shift amount
Degree.
Optionally, imaging unit is also used to：
According to the imaging formula R (m based on time shift image-forming conditionk, t, τ) and=Pr(mk,t+τ)*Ps(mk, t- τ) and it is imaged,
Obtain the migrated section about time shift amount, wherein R (m, t, τ) indicates that the migrated section about time shift amount, * indicate that time-domain is mutual
Correlation, t indicate the time, and τ indicates the time shift amount in time shift image-forming condition, PrGeophone station wave field after indicating backward extension, PsIt indicates
Source wavefield after positive continuation, mkIndicate current reflectance model.
Optionally, the second determining module includes：
Determination unit updates step-length and more new direction for determining based on target function gradient；
Third computing unit, for updating formula m according to reflectivity modelk+1=mk+αkdkCalculate updated reflection
Modulus Model, wherein mk+1Indicate updated reflectivity model, mkIndicate current reflectance model, αkIt indicates to update step
It is long, dkIndicate more new direction.
Optionally it is determined that unit is also used to：
According to more new direction formulaDetermine more new direction, whereindkIndicate more new direction, dk-1In expression
One more new direction, g (mk) indicate in current reflectance model mkUnder target function gradient, g (mk-1) indicate in a upper reflection
Modulus Model mk-1Under target function gradient；
Step-length is updated according to updating step-length condition and determining, wherein updating step-length condition is：When meeting E (mk+αkdk)≤E
(mk)+αkcg(mk)dkWhen, then αk=αk-1, otherwise, αk=λ αk-1, αkIt indicates to update step-length, αk-1Indicate that upper one updates step-length, E
(mk+αkdk) indicate that in more new direction and update step-length be respectively dkAnd αkWhen, the wave field under updated reflectivity model is residual
Difference, g (mk) indicate in current reflectance model mkUnder target function gradient, 0<c<1,0<λ<1.
The technical effect and preceding method embodiment phase of device provided by the embodiment of the present invention, realization principle and generation
Together, to briefly describe, Installation practice part does not refer to place, can refer to corresponding contents in preceding method embodiment.
The computer program product of seismic migration imaging method and device provided by the embodiment of the present invention, including store
The computer readable storage medium of program code, the instruction that said program code includes can be used for executing in previous methods embodiment
The method, specific implementation can be found in embodiment of the method, and details are not described herein.
It is apparent to those skilled in the art that for convenience and simplicity of description, the system of foregoing description
It with the specific work process of device, can refer to corresponding processes in the foregoing method embodiment, details are not described herein.
In addition, in the description of the embodiment of the present invention unless specifically defined or limited otherwise, term " installation ", " phase
Even ", " connection " shall be understood in a broad sense, for example, it may be being fixedly connected, may be a detachable connection, or be integrally connected；It can
To be mechanical connection, it is also possible to be electrically connected；It can be directly connected, can also can be indirectly connected through an intermediary
Connection inside two elements.For the ordinary skill in the art, above-mentioned term can be understood at this with concrete condition
Concrete meaning in invention.
It, can be with if the function is realized in the form of SFU software functional unit and when sold or used as an independent product
It is stored in a computer readable storage medium.Based on this understanding, technical solution of the present invention is substantially in other words
The part of the part that contributes to existing technology or the technical solution can be embodied in the form of software products, the meter
Calculation machine software product is stored in a storage medium, including some instructions are used so that a computer equipment (can be a
People's computer, server or network equipment etc.) it performs all or part of the steps of the method described in the various embodiments of the present invention.
And storage medium above-mentioned includes：USB flash disk, mobile hard disk, read-only memory (ROM, Read-Only Memory), arbitrary access are deposited
The various media that can store program code such as reservoir (RAM, Random Access Memory), magnetic or disk.
In the description of the present invention, it should be noted that term " center ", "upper", "lower", "left", "right", "vertical",
The orientation or positional relationship of the instructions such as "horizontal", "inner", "outside" be based on the orientation or positional relationship shown in the drawings, merely to
Convenient for description the present invention and simplify description, rather than the device or element of indication or suggestion meaning must have a particular orientation,
It is constructed and operated in a specific orientation, therefore is not considered as limiting the invention.In addition, term " first ", " second ",
" third " is used for descriptive purposes only and cannot be understood as indicating or suggesting relative importance.
Finally it should be noted that：Embodiment described above, only a specific embodiment of the invention, to illustrate the present invention
Technical solution, rather than its limitations, scope of protection of the present invention is not limited thereto, although with reference to the foregoing embodiments to this hair
It is bright to be described in detail, those skilled in the art should understand that：Anyone skilled in the art
In the technical scope disclosed by the present invention, it can still modify to technical solution documented by previous embodiment or can be light
It is readily conceivable that variation or equivalent replacement of some of the technical features；And these modifications, variation or replacement, do not make
The essence of corresponding technical solution is detached from the spirit and scope of technical solution of the embodiment of the present invention, should all cover in protection of the invention
Within the scope of.Therefore, the protection scope of the present invention shall be subject to the protection scope of the claims.
Claims (10)
1. a kind of seismic migration imaging method, which is characterized in that including：
Seismic observation data and analogue data to be solved based on pending area establish wave field residual error objective function；
Current wave field residual values are determined in conjunction with current reflectance model and the wave field residual error objective function, wherein described to work as
Front-reflection Modulus Model is for determining the analogue data to be solved；
Offset method based on time shift image-forming condition deviates the current wave field residual values, obtains objective function with determination
Gradient；
Updated reflectivity model is determined using conjugate gradient method based on the target function gradient；
It is updated iteration using the updated reflectivity model as the current reflectance model, until reaching pre-
If the number of iterations, target reflection factor model is obtained, and then obtains migration imaging result.
2. the method according to claim 1, wherein seismic observation data based on pending area and to be solved
Analogue data establishes wave field residual error objective function：
Obtain the seismic observation data and migration velocity body of the pending area；
The wave field residual error objective function between the seismic observation data and the analogue data to be solved is established by L2 norm.
3. the method according to claim 1, wherein in conjunction with current reflectance model and the wave field residual error mesh
Scalar functions determine that current wave field residual values include：
Formula d is calculated according to analogue datamod=LmkCalculate the analogue data to be solved, wherein dmodIndicate described to be solved
Analogue data, L indicate the seismic wave field forward operator that velocity information is depended under Born approximate condition, mkIndicate described current anti-
Penetrate Modulus Model；
According to the calculating formula of the wave field residual error objective functionCalculate the current wave field residual error
Value, wherein E (mk) indicate the current wave field residual values, dmodIndicate the analogue data to be solved, dobsIndicate the earthquake
Observe data.
4. according to the method described in claim 2, it is characterized in that, the offset method based on time shift image-forming condition is to described current
Wave field residual values are deviated, to determine that obtaining target function gradient includes：
Geophone station wave field is subjected to backward extension by the migration velocity body, and passes through the migration velocity body for source wavefield
Carry out positive continuation；
The geophone station wave field after backward extension and the source wavefield after positive continuation are imaged based on time shift image-forming condition, obtained
To the migrated section about time shift amount；
Using the migrated section of zero moment as the target function gradient in the migrated section about time shift amount.
5. according to the method described in claim 4, it is characterized in that, based on time shift image-forming condition to the inspection after the backward extension
Source wavefield after wave point wave field and the positive continuation carries out imaging and includes：
According to the imaging formula R (m based on time shift image-forming conditionk, t, τ) and=Pr(mk,t+τ)*Ps(mk, t- τ) and it is imaged, it obtains
The migrated section about time shift amount, wherein R (m, t, τ) indicates that the migrated section about time shift amount, * indicate the time
Domain cross-correlation, t indicate the time, and τ indicates the time shift amount in time shift image-forming condition, PrGeophone station wave after indicating the backward extension
, PsSource wavefield after indicating the positive continuation, mkIndicate the current reflectance model.
6. the method according to claim 1, wherein true using conjugate gradient method based on the target function gradient
Determining updated reflectivity model includes：
It is determined based on the target function gradient and updates step-length and more new direction；
Formula m is updated according to reflectivity modelk+1=mk+αkdkCalculate the updated reflectivity model, wherein mk+1
Indicate the updated reflectivity model, mkIndicate current reflectance model, αkIndicate the update step-length, dkIt indicates
The more new direction.
7. according to the method described in claim 6, it is characterized in that, based on the determining update step-length of the target function gradient and more
New direction includes：
According to more new direction formulaDetermine the more new direction, whereindkIndicate the more new direction, dk-1Table
Show a more new direction, g (mk) indicate in the current reflectance model mkUnder target function gradient, g (mk-1) indicate
Upper reflectivity model mk-1Under target function gradient；
The update step-length is determined according to step-length condition is updated, wherein the update step-length condition is：When meeting E (mk+αkdk)
≤E(mk)+αkcg(mk)dkWhen, then αk=αk-1, otherwise, αk=λ αk-1, αkIndicate the update step-length, αk-1Indicate that upper one updates
Step-length, E (mk+αkdk) indicate that in more new direction and update step-length be respectively dkAnd αkWhen, under updated reflectivity model
Wave field residual values, g (mk) indicate in the current reflectance model mkUnder target function gradient, 0<c<1,0<λ<1.
8. a kind of seismic migration imaging device, which is characterized in that including：
Establish module, for based on pending area seismic observation data and analogue data to be solved establish wave field residual error target
Function；
First determining module, for determining that current wave field is residual in conjunction with current reflectance model and the wave field residual error objective function
Difference, wherein the current reflectance model is for determining the analogue data to be solved；
Offset module deviates the current wave field residual values for the offset method based on time shift image-forming condition, with true
Surely target function gradient is obtained；
Second determining module, for determining updated reflection coefficient mould using conjugate gradient method based on the target function gradient
Type；
Iteration module is updated, for carrying out the updated reflectivity model as the current reflectance model more
New iteration obtains target reflection factor model, and then obtain migration imaging result until reaching default the number of iterations.
9. device according to claim 8, which is characterized in that the module of establishing includes：
Acquiring unit, for obtaining the seismic observation data and migration velocity body of the pending area；
Unit is established, for establishing the wave field between the seismic observation data and the analogue data to be solved by L2 norm
Residual error objective function.
10. device according to claim 8, which is characterized in that first determining module includes：
First computing unit, for calculating formula d according to analogue datamod=LmkCalculate the analogue data to be solved, wherein
dmodIndicate the analogue data to be solved, L indicates that the seismic wave field under Born approximate condition dependent on velocity information just calculates
Son, mkIndicate the current reflectance model；
Second computing unit, for the calculating formula according to the wave field residual error objective functionMeter
Calculate the current wave field residual values, wherein E (mk) indicate the current wave field residual values, dmodIndicate the simulation number to be solved
According to dobsIndicate the seismic observation data.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201811128030.5A CN108845355A (en) | 2018-09-26 | 2018-09-26 | Seismic migration imaging method and device |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201811128030.5A CN108845355A (en) | 2018-09-26 | 2018-09-26 | Seismic migration imaging method and device |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN108845355A true CN108845355A (en) | 2018-11-20 |
Family
ID=64188069
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201811128030.5A Pending CN108845355A (en) | 2018-09-26 | 2018-09-26 | Seismic migration imaging method and device |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN108845355A (en) |
Cited By (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN111624658A (en) * | 2020-05-29 | 2020-09-04 | 中国石油天然气集团有限公司 | Depth domain imaging simulation method and system |
| CN112698280A (en) * | 2020-12-09 | 2021-04-23 | 南京长峰航天电子科技有限公司 | Bistatic SAR real-time echo simulation method based on DSP and FPGA architecture |
| CN112773396A (en) * | 2021-01-13 | 2021-05-11 | 佟小龙 | Medical imaging method based on full waveform inversion, computer equipment and storage medium |
| CN112946744A (en) * | 2019-12-11 | 2021-06-11 | 中国石油天然气集团有限公司 | Least square offset imaging method and system based on dynamic time difference warping |
| CN113866825A (en) * | 2021-08-23 | 2021-12-31 | 中国石油大学(华东) | Angle domain least square reflectivity inversion method based on coherent superposition |
| CN112698280B (en) * | 2020-12-09 | 2024-05-31 | 南京长峰航天电子科技有限公司 | Double-base SAR real-time echo simulation method based on DSP and FPGA architecture |
Citations (12)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN103675895A (en) * | 2012-08-30 | 2014-03-26 | 中国石油化工股份有限公司 | Method for utilizing GPU (Graphic Processing Unit) to increase computing efficiency of wave field continuation |
| CN105487112A (en) * | 2014-09-18 | 2016-04-13 | 中国石油化工股份有限公司 | Method for constructing stratum reflection coefficient |
| CN106033124A (en) * | 2016-06-29 | 2016-10-19 | 中国石油化工股份有限公司 | Multi-seismic resource sticky sound least square reverse time migration method based on stochastic optimization |
| CN107193043A (en) * | 2017-05-15 | 2017-09-22 | 中国石油大学(华东) | A kind of subsurface structure imaging method of relief surface |
| CN107229071A (en) * | 2017-05-25 | 2017-10-03 | 中国石油大学(华东) | A kind of subsurface structure inversion imaging method |
| CN107356972A (en) * | 2017-06-28 | 2017-11-17 | 中国石油大学(华东) | A kind of imaging method of anisotropic medium |
| CN107589443A (en) * | 2017-08-16 | 2018-01-16 | 东北石油大学 | Method and system based on elastic wave least square reverse-time migration imaging |
| CN107783190A (en) * | 2017-10-18 | 2018-03-09 | 中国石油大学(北京) | A kind of least square reverse-time migration gradient updating method |
| CN107870363A (en) * | 2016-09-27 | 2018-04-03 | 中国石油化工股份有限公司 | Least-squares migration imaging optimization method and system |
| CN108241173A (en) * | 2017-12-28 | 2018-07-03 | 中国石油大学(华东) | A kind of seismic data offset imaging method and system |
| CN108333628A (en) * | 2018-01-17 | 2018-07-27 | 中国石油大学(华东) | Elastic wave least square reverse-time migration method based on regularization constraint |
| CN108802813A (en) * | 2018-06-13 | 2018-11-13 | 中国石油大学(华东) | A kind of multi-component seismic data offset imaging method and system |
-
2018
- 2018-09-26 CN CN201811128030.5A patent/CN108845355A/en active Pending
Patent Citations (12)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN103675895A (en) * | 2012-08-30 | 2014-03-26 | 中国石油化工股份有限公司 | Method for utilizing GPU (Graphic Processing Unit) to increase computing efficiency of wave field continuation |
| CN105487112A (en) * | 2014-09-18 | 2016-04-13 | 中国石油化工股份有限公司 | Method for constructing stratum reflection coefficient |
| CN106033124A (en) * | 2016-06-29 | 2016-10-19 | 中国石油化工股份有限公司 | Multi-seismic resource sticky sound least square reverse time migration method based on stochastic optimization |
| CN107870363A (en) * | 2016-09-27 | 2018-04-03 | 中国石油化工股份有限公司 | Least-squares migration imaging optimization method and system |
| CN107193043A (en) * | 2017-05-15 | 2017-09-22 | 中国石油大学(华东) | A kind of subsurface structure imaging method of relief surface |
| CN107229071A (en) * | 2017-05-25 | 2017-10-03 | 中国石油大学(华东) | A kind of subsurface structure inversion imaging method |
| CN107356972A (en) * | 2017-06-28 | 2017-11-17 | 中国石油大学(华东) | A kind of imaging method of anisotropic medium |
| CN107589443A (en) * | 2017-08-16 | 2018-01-16 | 东北石油大学 | Method and system based on elastic wave least square reverse-time migration imaging |
| CN107783190A (en) * | 2017-10-18 | 2018-03-09 | 中国石油大学(北京) | A kind of least square reverse-time migration gradient updating method |
| CN108241173A (en) * | 2017-12-28 | 2018-07-03 | 中国石油大学(华东) | A kind of seismic data offset imaging method and system |
| CN108333628A (en) * | 2018-01-17 | 2018-07-27 | 中国石油大学(华东) | Elastic wave least square reverse-time migration method based on regularization constraint |
| CN108802813A (en) * | 2018-06-13 | 2018-11-13 | 中国石油大学(华东) | A kind of multi-component seismic data offset imaging method and system |
Non-Patent Citations (1)
| Title |
|---|
| 刘守伟 等: "时空移动成像条件及偏移速度分析", 《地球物理学报》 * |
Cited By (7)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN112946744A (en) * | 2019-12-11 | 2021-06-11 | 中国石油天然气集团有限公司 | Least square offset imaging method and system based on dynamic time difference warping |
| CN111624658A (en) * | 2020-05-29 | 2020-09-04 | 中国石油天然气集团有限公司 | Depth domain imaging simulation method and system |
| CN112698280A (en) * | 2020-12-09 | 2021-04-23 | 南京长峰航天电子科技有限公司 | Bistatic SAR real-time echo simulation method based on DSP and FPGA architecture |
| CN112698280B (en) * | 2020-12-09 | 2024-05-31 | 南京长峰航天电子科技有限公司 | Double-base SAR real-time echo simulation method based on DSP and FPGA architecture |
| CN112773396A (en) * | 2021-01-13 | 2021-05-11 | 佟小龙 | Medical imaging method based on full waveform inversion, computer equipment and storage medium |
| CN113866825A (en) * | 2021-08-23 | 2021-12-31 | 中国石油大学(华东) | Angle domain least square reflectivity inversion method based on coherent superposition |
| CN113866825B (en) * | 2021-08-23 | 2023-12-01 | 中国石油大学(华东) | Angular domain least square reflectivity inversion method based on coherent superposition |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN108845355A (en) | Seismic migration imaging method and device | |
| Bistacchi et al. | Photogrammetric digital outcrop reconstruction, visualization with textured surfaces, and three-dimensional structural analysis and modeling: Innovative methodologies applied to fault-related dolomitization (Vajont Limestone, Southern Alps, Italy) | |
| NO20161080A1 (en) | Seismic data migration system and method | |
| EA008733B1 (en) | Method for seismic imaging in geologically complex formations | |
| CN106772592B (en) | Diffracted wave focuses the analysis method and device of energy | |
| EA032008B1 (en) | Two stage seismic velocity model generation | |
| BR112017020991B1 (en) | SEISMIC DATA PROCESSING METHOD AND SYSTEM, AND COMPUTER READABLE RECORDING MEDIA | |
| CN110031898B (en) | Data optimization method and integral method prestack depth migration method | |
| BR102013014787A2 (en) | Method, system, and one or more computer readable storage media | |
| CN107894613A (en) | Elastic wave vector imaging method, device, storage medium and equipment | |
| CN107229071B (en) | A kind of subsurface structure inversion imaging method | |
| CN109633749A (en) | Non-linear Fresnel zone seismic traveltime tomography method based on scattering integral method | |
| CN106324665A (en) | Method and system of inverting fracture density | |
| EP3436850B1 (en) | Determining displacement between seismic images using optical flow | |
| CN106772593B (en) | The imaging method and device of diffracted wave | |
| CN109116413B (en) | Imaging domain stereo chromatography velocity inversion method | |
| CN109636912A (en) | Tetrahedron subdivision finite element interpolation method applied to three-dimensional sonar image reconstruction | |
| CN105137479B (en) | A kind of computational methods and device of bin degree of covering | |
| CN106842302B (en) | A kind of method and device of batch editor first arrival | |
| CN104199088B (en) | Incident angle gather extraction method and system | |
| CN108845350A (en) | The method and device of inverting two-dimension speed model | |
| CN105353406B (en) | A kind of method and apparatus for generating angle gathers | |
| CN104316961A (en) | Method for obtaining geological parameters of weathered layer | |
| Gong et al. | Combined migration velocity model-building and its application in tunnel seismic prediction | |
| CN105093300B (en) | A kind of boundary recognition of geological body method and device |
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 | ||
| RJ01 | Rejection of invention patent application after publication |
Application publication date: 20181120 |
|
| RJ01 | Rejection of invention patent application after publication |