CN107193043A - A kind of subsurface structure imaging method of relief surface - Google Patents
A kind of subsurface structure imaging method of relief surface Download PDFInfo
- Publication number
- CN107193043A CN107193043A CN201710337714.5A CN201710337714A CN107193043A CN 107193043 A CN107193043 A CN 107193043A CN 201710337714 A CN201710337714 A CN 201710337714A CN 107193043 A CN107193043 A CN 107193043A
- Authority
- CN
- China
- Prior art keywords
- data
- mrow
- gradient
- msub
- model
- 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.)
- Granted
Links
- 238000003384 imaging method Methods 0.000 title claims abstract description 49
- 238000005070 sampling Methods 0.000 claims abstract description 40
- 238000013508 migration Methods 0.000 claims abstract description 37
- 230000005012 migration Effects 0.000 claims abstract description 37
- 238000002310 reflectometry Methods 0.000 claims abstract description 37
- 238000000034 method Methods 0.000 claims abstract description 19
- 238000012952 Resampling Methods 0.000 claims abstract description 14
- 238000012986 modification Methods 0.000 claims description 23
- 230000004048 modification Effects 0.000 claims description 23
- 238000004364 calculation method Methods 0.000 claims description 12
- 230000017105 transposition Effects 0.000 claims description 3
- 229910017435 S2 In Inorganic materials 0.000 claims 1
- 238000000205 computational method Methods 0.000 claims 1
- 230000006378 damage Effects 0.000 abstract description 4
- 238000011478 gradient descent method Methods 0.000 description 5
- 238000004088 simulation Methods 0.000 description 5
- 238000002939 conjugate gradient method Methods 0.000 description 4
- 238000010276 construction Methods 0.000 description 4
- 230000021615 conjugation Effects 0.000 description 3
- 239000011159 matrix material Substances 0.000 description 3
- 238000007689 inspection Methods 0.000 description 2
- 238000010606 normalization Methods 0.000 description 2
- 238000012360 testing method Methods 0.000 description 2
- 238000004458 analytical method Methods 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 235000013399 edible fruits Nutrition 0.000 description 1
- 230000000694 effects Effects 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 230000036039 immunity Effects 0.000 description 1
- 238000012804 iterative process Methods 0.000 description 1
- 238000012545 processing Methods 0.000 description 1
- 230000000630 rising effect Effects 0.000 description 1
- 150000003839 salts Chemical class 0.000 description 1
- 239000004576 sand Substances 0.000 description 1
- 238000002945 steepest descent method Methods 0.000 description 1
- 238000003325 tomography Methods 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
- G01V1/282—Application of seismic models, synthetic seismograms
-
- 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
- G01V1/30—Analysis
- G01V1/306—Analysis for determining physical properties of the subsurface, e.g. impedance, porosity or attenuation profiles
-
- 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
- G01V1/34—Displaying seismic recordings or visualisation of seismic data or attributes
- G01V1/345—Visualisation of seismic data or attributes, e.g. in 3D cubes
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V20/00—Geomodelling in general
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V2210/00—Details of seismic processing or analysis
- G01V2210/60—Analysis
- G01V2210/62—Physical property of subsurface
- G01V2210/624—Reservoir parameters
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V2210/00—Details of seismic processing or analysis
- G01V2210/60—Analysis
- G01V2210/64—Geostructures, e.g. in 3D data cubes
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V2210/00—Details of seismic processing or analysis
- G01V2210/60—Analysis
- G01V2210/66—Subsurface modeling
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V2210/00—Details of seismic processing or analysis
- G01V2210/70—Other details related to processing
- G01V2210/74—Visualisation of seismic data
Landscapes
- Physics & Mathematics (AREA)
- Life Sciences & Earth Sciences (AREA)
- Engineering & Computer Science (AREA)
- Remote Sensing (AREA)
- General Life Sciences & Earth Sciences (AREA)
- General Physics & Mathematics (AREA)
- Geophysics (AREA)
- Acoustics & Sound (AREA)
- Environmental & Geological Engineering (AREA)
- Geology (AREA)
- Investigating Or Analysing Materials By Optical Means (AREA)
Abstract
The invention discloses a kind of subsurface structure imaging method of relief surface, comprise the following steps：S1, input observation data, migration velocity field, while pre-setting the initial value of reflectivity model；Observation data are carried out the sampling based on coefficient correlation or resampling, obtain offset data by S2；S3, data residual error is obtained according to current reflectivity model and offset data；S4, according to data residual error, determines gradient, and reflectivity model iteration is carried out based on gradient, obtains final reflectivity model.The inventive method is by introducing the sampling policy based on coefficient correlation, avoid stochastical sampling, offset data, which can be effectively prevented from, while offset data amount is reduced produces great variety with iterations, and then weaken replacing offset data to the destruction of gradient direction conjugacy and the jitter phenomenon of convergence curve, improve the computational efficiency and final image quality during iteration.
Description
Technical field
The present invention relates to a kind of subsurface structure imaging method of relief surface, belong to oil gas physical prospecting engineering field.
Background technology
With going deep into for China's exploration and development, exploration key area is gradually transferred to western exploratory area.The exploration in western exploratory area
There is following main feature in exploitation：(1) acutely, near surface construction is complicated for surface relief；(2) Background Construction is complicated, such as inverse to cover
Push away and cover and imaging etc. under salt；(3) reservoir buries relatively deep, and signal is weaker.Therefore, how real the Important Problems of western exploratory area exploration are
Efficient, high accuracy under the conditions of existing relief surface protect width imaging.
Current conventional inversion imaging method is solved by iterative algorithms such as steepest descent method, conjugate gradient methods, although can
To obtain high-quality imaging results, but amount of calculation is excessively huge.If carrying out stochastical sampling to observation data, updating every time
A calculating section data, can reduce amount of calculation during model.However, error function when stochastic sampling strategy causes inversion imaging
All there is randomness with gradient.It is adjacent when subsurface structure complexity or violent surface relief due to observation data stochastical sampling
Correlation is poor between offset data in iteration twice, and acutely shake often occurs in the convergence curve of stochastic gradient descent method, is carrying out
Image quality is poor when underground structure is imaged, and convergence rate is slower, and computational efficiency is low.
The content of the invention
It is an object of the present invention to provide a kind of subsurface structure imaging method of relief surface, it can solve current skill
Computational efficiency when being imaged under the conditions of problem present in art, raising relief surface, while not influenceing its image quality.
In order to solve the above technical problems, the present invention is adopted the following technical scheme that：A kind of subsurface structure of relief surface into
Image space method, comprises the following steps：
S1, input observation data, migration velocity field, while pre-setting the initial value of reflectivity model；
Observation data are carried out the sampling based on coefficient correlation or resampling, obtain offset data by S2；
S3, data residual error is obtained according to current reflectivity model and offset data；
S4, according to data residual error, determines stochastic gradient and/or random conjugate gradient, based on stochastic gradient and/or common at random
Yoke gradient carries out reflectivity model iteration, obtains final reflectivity model.
Compared with prior art, with the method for the present invention when being imaged under the conditions of relief surface, offset data amount is being reduced
While can effectively reduce change of the offset data with iterations, and then weaken and change offset data gradient direction is conjugated
The destruction of property and the jitter phenomenon of convergence curve, improve the convergence stability and its image quality during iteration.In addition, of the invention
Relief surface reverse-time migration will be introduced in Least squares inversion thought, can be pressed into as noise, balanced imaging amplitude, improve into
As resolution ratio, near surface construction and infrastructure can be in high precision imaged.
Brief description of the drawings
Fig. 1 is the FB(flow block) of one embodiment of the present of invention；
Fig. 2 is the migration velocity field of one embodiment of the present of invention；
Fig. 3 is the relief surface elevation information of one embodiment of the present of invention；
Fig. 4 is the observation geological data of one embodiment of the present of invention；
Fig. 5 is carried out to observation data based on obtaining after coefficient correlation big gun data sampling for one embodiment of the present of invention
Offset data；
Fig. 6 is the relief surface reverse-time migration result of one embodiment of the present of invention；
Fig. 7 is the relief surface least square reverse-time migration result based on stochastic gradient descent method；
Fig. 8 is the relief surface least square reverse-time migration result of one embodiment of the present of invention；
Fig. 9 is that the single track amplitude of one embodiment of the invention and other two kinds of imaging method results is contrasted；
Figure 10 is the normalization data residual error convergence curve of one embodiment of the invention and other two kinds of imaging method results；
Figure 11 is to contrast one embodiment of the invention and the calculating time of other two kinds of imaging methods；
Figure 12 is one embodiment of the present of invention schematic flow sheet.
The present invention is further illustrated with reference to the accompanying drawings and detailed description.
Embodiment
Embodiments of the invention 1：The subsurface structure imaging method of a kind of relief surface, as shown in figure 1, including following step
Suddenly：
S1, input observation data, migration velocity field, while pre-setting the initial value of reflectivity model；Here sight
Survey data and include but is not limited to relief surface elevation, observation system, field inspection record (i.e. big gun data) and can be with default inclined
Shifting parameter etc., wherein offset parameter include the horizontal sampled point n of migration velocity fieldxAnd longitudinal sampled point nz, spatial sampling interval dx,
Time sampling interval dt, time sampling points nt, dominant frequency f0, big gun number N in the offset data that obtains after samplingsAnd model modification
Direction controlling factor o and iteration ends threshold value.Common, the initial value of reflectivity model is arranged to 0, migration velocity field
It is to be obtained after carrying out conventional velocity analysis to observation data.
Observation data are carried out the sampling based on coefficient correlation or resampling, obtain offset data by S2；Here correlation
Coefficient calculation method is：
Wherein, dj(x,z),For the single-shot data and its average in last iteration hour offset data, di(x,z),It is surplus
Single-shot data and its average of the remainder in.In formula (1), first to di(x,t),dj(x, t) carries out removing average value processing, subtracts
Influence of the weak Non-zero Mean noise to coefficient correlation, therefore the formula has certain noise immunity.Then, d is calculated respectivelyi(x,
T) with the coefficient correlation of each single-shot data in last iteration hour offset data.Finally, all coefficient correlation summations are obtained
di(x, t) and offset data coefficient correlation.From formula (1)-(3), coefficient correlation can reflect a certain single-shot data with it is upper
The correlation of secondary iteration hour offset data, using coefficient correlation as referring to and then develop based on coefficient correlation such as down-sampling plan
Slightly：
It is exactly in the 1st iteration, according to fixed interval sampling N when sampling firstsIt is individual uniform along offset distance direction
The observation data of distribution are more than big gun interval (knot of being sampled when the sampling interval is equal to big gun interval as offset data, wherein sampling interval
Fruit is whole data volume；Sampled result is only part big gun data when only the sampling interval is more than big gun interval).When iterations is more than
1 and when meeting resampling condition, the coefficient correlation per big gun data with last iteration hour offset data in remaining data is calculated,
Choose the maximum N of coefficient correlationsIndividual single-shot data as current iteration offset data；When the number of single-shot data in remaining data
Amount is less than NsWhen, not enough part is randomly selected from the data calculated., directly will be upper when being unsatisfactory for resampling condition
An iteration hour offset data as current iteration offset data.It should be noted that remaining data here refers to observation
Remove in data after the offset data sampled during former all iteration, remaining data.Resampling condition is：When iteration time
When number k divided by model modification direction controlling factor o remainder is 1, resampling is carried out to observation data, offset data is changed；It is no
Then, without to observation data carry out resampling, directly using offset data during last iteration as current iteration offset data.
In with machine gun data sampling, offset data changes at random with iterations, when offset numbers in adjacent iteration twice
According to correlation it is poor when, its convergence curve often produces acutely shake, and convergence rate is slower, and computational efficiency is low.And based on phase relation
Several big gun data sampling strategies ensure that the difference of offset data in adjacent iteration twice is minimum, that is, the sampling policy is not
Influenceed by subsurface structure or surface conditions, can effectively reduce change of the offset data with iterations, and then weaken replacing
Offset data improves convergence stability during iteration to the destruction of gradient direction conjugacy and the jitter phenomenon of convergence curve,
And then improve computational efficiency.Meanwhile, resampling decision condition is added in the big gun data sampling strategy based on coefficient correlation,
When being unsatisfactory for resampling condition, offset data is constant, so ensures that the conjugacy of gradient direction is set up completely, improves iteration
When convergence efficiency and image quality.
S3, data residual error is obtained according to current reflectivity model and offset data；Specifically, for reflectance factor
Model, model data P (x), the calculation formula used in inverse time inverse migration simulation model data are obtained by the simulation of inverse time inverse migration
For：
P (x)=ω2∫m(x')f(ω)G0(x'；xs)G0(x；x')dx' (4)
Wherein, P (x) represents the model data that inverse time inverse migration is obtained, and m (x') represents reflectivity model, and f (ω) is shake
Source, ω is frequency, G0(x'；xs),G0(x；X') it is the corresponding Green's function of background velocity, xsFor shot position, x, x' is respectively
Underground any point position.
Formula (4) can be solved by two step forward simulations：
Wherein, P0(x') it is in migration velocity v0(x) the background wave field tried to achieve under, by its product with reflectivity model
Excited again as focus, you can try to achieve the model data of inverse time inverse migration simulation.Enter traveling-wave field according to formula (5) and formula (6) and prolong
When opening up, often it is required for applying free boundary condition to wave field by a continuation：
PΠ(x, z)=0 (7)
Wherein, Π represents the position of relief surface, i.e., the wave field on relief surface is entered as into 0.
Based on described above, formula (4) can be simplified shown as with matrix operator form：
D'=Lm (8)
Wherein, m is the matrix form of reflectivity model, and d' is the matrix form of model data, and L calculates for inverse time inverse migration
Positive operator in son, the i.e. reverse-time migration of relief surface least square.That is it is reverse-biased by doing the inverse time to reflectivity model
Movement, which is calculated, can obtain model data P (x), in the methods of the invention, only need to calculate offset data correspondence in each iteration
Model data, be expressed as：
d's=Lsm (9)
Wherein, d's,LsRespectively the corresponding model data of offset data and inverse time inverse migration operator, obtain model data
Afterwards, subtract each other with current offset data, obtain data residual errorWithRespectively kth time changes
For when offset data and its corresponding inverse time inverse migration operator, m(k-1)The reflectance factor mould obtained for -1 iterative calculation of kth
Type.
S4, according to data residual error, determines model modification gradient, and carrying out reflectivity model based on model modification gradient changes
In generation, final reflectivity model is obtained, gradient here includes stochastic gradient and/or random conjugation, specifically, passes through
The calculation formula of data residual computations stochastic gradient is：
Wherein,ForConjugate transposition, i.e. reverse-time migration operator.
The calculation formula of conjugate gradient modifying factor is at random：
Then the calculation formula of random conjugate gradient is：
Then, in the iteration of each step, stochastic gradient or random conjugate gradient will be selected based on decision condition, as
Model modification gradient carries out reflectivity model iteration, obtains final reflectivity model, specifically includes following steps：
S41, whether the remainder for judging iterations k divided by model modification direction controlling factor o is 1, if so, then select with
Machine gradient is model modification gradient, is used as the direction of reflectivity model iteration；If it is not, then selecting random conjugate gradient to be model
Gradient is updated, the direction of reflectivity model iteration is used as.
S42, when described model modification gradient is less than threshold value, into step S44, threshold value here is can be according to warp
Test what the requirement adjustment with available accuracy was set, when described model modification gradient is more than threshold value, according to model modification gradient
Calculate and update step-length, the calculation formula for updating step-length is：
Wherein,For model modification gradient, including stochastic gradient and/or random conjugate gradient, which is selected in specific iteration
Individual parameter is determined in step S41.
S43, based on step-length is updated, iteration updates reflectivity model, and the more new formula of reflectivity model is：
And return to step S2,
S44, judges the reflectivity model now obtained as final reflectivity model, by final reflectance factor mould
Type obtains subsurface structure imaging results.Using the method for the present invention, by introducing the sampling policy based on coefficient correlation, it is to avoid
Stochastical sampling, when being imaged under conditions of relief surface, offset data can be effectively prevented from while offset data amount is reduced
Great variety is produced with iterations, and then weakens and changes offset data to the destruction of gradient direction conjugacy and convergence curve
Jitter phenomenon, improves computational efficiency during iteration.In addition, the present invention will in Least squares inversion thought introduce relief surface it is inverse
Hour offset, can be pressed into as noise, balanced imaging amplitude, improve imaging resolution, high near surface construction and infrastructure
Precision is imaged；Using mixing gradient in iterative process, used random conjugate gradient to solve in a part of iteration, relative to
The convergence of machine gradient is more stablized, and imaging effect is more preferable.
Embodiment 2
Method in order to further illustrate the present invention, illustrates the side of the present invention by taking Canadian overthrust fault model as an example
Method.Idiographic flow is as shown in figure 12.Input observation data, including migration velocity field (as shown in Figure 2), relief surface elevation are (such as
Shown in Fig. 3), observation system, field inspection data (as shown in Figure 4), the threshold value of iteration ends and default bias parameter, and set
The value of initial reflection Modulus Model is 0, and observation system is distributed as：278 big guns are uniformly distributed with 30 meters of intervals in relief surface, per big gun
All it is 556 geophone station receptions, geophone station is at intervals of 15 meters, and offset parameter is as follows：Migration velocity field transverse direction sampled point be 556 and
Longitudinal sampled point is 250, and spatial sampling interval is 15 meters, and time sampling interval is 0.5 millisecond, and time sampling points are 4000
Individual, dominant frequency is 25 hertz, and the big gun number of offset data is 40 after sampling, and the model modification direction controlling factor is 5；Observation data are entered
Big gun data sampling of the row based on coefficient correlation, the observation data after sampling are as shown in Figure 5；For current reflectance model, lead to
Inverse time inverse migration simulation model data are crossed, are subtracted each other with offset data, data residual error is calculated；According to data residual computations stochastic gradient
Descent direction；Judge whether to need amendment stochastic gradient direction, if meeting Rule of judgment, stochastic gradient is modified to random conjugation
Gradient, as model modification direction, otherwise, regard stochastic gradient as model modification direction；Whether judgment models more new direction is full
Sufficient threshold condition, if being unsatisfactory for threshold condition, updates step-length according to model modification direction calculating and updates reflectivity model, then
It is secondary that the big gun data sampling based on coefficient correlation, inverse time inverse migration, computation model more new direction are carried out to observation data, until model
More new direction meets threshold condition, if meeting threshold condition, exports final reflectivity model, i.e., final imaging results
(as shown in Figure 8).
From imaging results, using the method for the present invention, relative to relief surface reverse-time migration result (such as Fig. 6 institutes
Show), imaging noise has been suppressed, the imaging energy in deep is compensate for, near surface structure imaging and medium and deep tomography become apparent from, carried
High imaging resolution；Relative to the relief surface least square reverse-time migration imaging results based on stochastic gradient descent method (such as
Shown in Fig. 7), imaging noise is weaker, and imaging amplitude is more balanced, and the imaging that near surface is constructed is become apparent from.
Fig. 9 is from the relief surface least square reverse-time migration imaging results based on conjugate gradient method, based on stochastic gradient
Distance=1000 meters in the relief surface least square reverse-time migration imaging results of descent method and the imaging results of the present invention
Locate the single track amplitude curve extracted, imaging results amplitude and the relief surface least square based on conjugate gradient method of the invention is inverse
Hour offset imaging results amplitude is close, relative to the relief surface least square reverse-time migration imaging based on stochastic gradient descent method
As a result amplitude, which has, preferably protects width.
Figure 10 is the normalization data residual error convergence curve of three kinds of imaging methods, and data residual error of the invention is with being based on conjugation
The data residual error of the relief surface least square reverse-time migration method of gradient method can stable convergence, and based under stochastic gradient
The data residual error of the relief surface least square reverse-time migration imaging method of drop method has obvious fluctuation, and convergent speed is slow
In the inventive method.
Figure 11 contrasts for the run time of several method, CPU models Intel (R) Xeon (R) CPU used in imaging test
E5-2650v2@2.60GHz.The time that calculates of imaging method of the present invention is far smaller than the relief surface based on conjugate gradient method most
A young waiter in a wineshop or an inn multiplies the calculating time of reverse-time migration method, but the two image quality is quite, relative to rising based on stochastic gradient descent method
Table least square reverse-time migration method of the throwing oneself on the ground calculating time is close, but image quality is higher.Therefore, the inventive method is with respect to showing
There is technical method, while taking into account higher computational efficiency and image quality, the exploration to western exploratory area is significant.
Claims (10)
1. the subsurface structure imaging method of a kind of relief surface, it is characterised in that comprise the following steps：
S1, input observation data, migration velocity field, while pre-setting the initial value of reflectivity model；
Observation data are carried out the sampling based on coefficient correlation or resampling, obtain offset data by S2；
S3, data residual error is obtained according to current reflectivity model and offset data；
S4, gradient is determined according to data residual error, and reflectivity model iteration is carried out based on gradient, obtains final reflectance factor mould
Type.
2. a kind of subsurface structure imaging method of relief surface according to claim 1, it is characterised in that the step S2
Middle coefficient correlation, its computational methods is：
<mrow>
<msub>
<mi>R</mi>
<mi>i</mi>
</msub>
<mo>=</mo>
<munderover>
<mo>&Sigma;</mo>
<mrow>
<mi>j</mi>
<mo>=</mo>
<mn>1</mn>
</mrow>
<msub>
<mi>N</mi>
<mi>s</mi>
</msub>
</munderover>
<mfrac>
<mrow>
<munder>
<mo>&Sigma;</mo>
<mrow>
<mi>x</mi>
<mo>,</mo>
<mi>t</mi>
</mrow>
</munder>
<mo>&lsqb;</mo>
<mrow>
<mo>(</mo>
<msub>
<mi>d</mi>
<mi>i</mi>
</msub>
<mo>(</mo>
<mrow>
<mi>x</mi>
<mo>,</mo>
<mi>t</mi>
</mrow>
<mo>)</mo>
<mo>-</mo>
<mover>
<msub>
<mi>d</mi>
<mi>i</mi>
</msub>
<mo>&OverBar;</mo>
</mover>
<mo>)</mo>
</mrow>
<mrow>
<mo>(</mo>
<msub>
<mi>d</mi>
<mi>j</mi>
</msub>
<mo>(</mo>
<mrow>
<mi>x</mi>
<mo>,</mo>
<mi>t</mi>
</mrow>
<mo>)</mo>
<mo>-</mo>
<mover>
<msub>
<mi>d</mi>
<mi>j</mi>
</msub>
<mo>&OverBar;</mo>
</mover>
<mo>)</mo>
</mrow>
<mo>&rsqb;</mo>
</mrow>
<msqrt>
<mrow>
<munder>
<mo>&Sigma;</mo>
<mrow>
<mi>x</mi>
<mo>,</mo>
<mi>t</mi>
</mrow>
</munder>
<msup>
<mrow>
<mo>(</mo>
<msub>
<mi>d</mi>
<mi>i</mi>
</msub>
<mo>(</mo>
<mrow>
<mi>x</mi>
<mo>,</mo>
<mi>t</mi>
</mrow>
<mo>)</mo>
<mo>-</mo>
<mover>
<msub>
<mi>d</mi>
<mi>i</mi>
</msub>
<mo>&OverBar;</mo>
</mover>
<mo>)</mo>
</mrow>
<mn>2</mn>
</msup>
<munder>
<mo>&Sigma;</mo>
<mrow>
<mi>x</mi>
<mo>,</mo>
<mi>t</mi>
</mrow>
</munder>
<msup>
<mrow>
<mo>(</mo>
<msub>
<mi>d</mi>
<mi>j</mi>
</msub>
<mo>(</mo>
<mrow>
<mi>x</mi>
<mo>,</mo>
<mi>t</mi>
</mrow>
<mo>)</mo>
<mo>-</mo>
<mover>
<msub>
<mi>d</mi>
<mi>j</mi>
</msub>
<mo>&OverBar;</mo>
</mover>
<mo>)</mo>
</mrow>
<mn>2</mn>
</msup>
</mrow>
</msqrt>
</mfrac>
</mrow>
For the single-shot data and its average in last iteration hour offset data, NsFor single-shot data in offset data
Quantity,For the single-shot data and its average in remaining data, wherein
nxFor hollow sampling number of single-shot data, ntCounted for time sampling, i, j represents big gun number.
3. a kind of subsurface structure imaging method of relief surface according to claim 2, it is characterised in that the step S2
In sampling based on coefficient correlation or resampling are carried out to observation data, specifically include：
When iterations is 1, according to the fixed sampling interval, N is selectedsIt is individual along the equally distributed observation data in offset distance direction
As offset data, and the sampling interval is more than big gun interval；
When iterations is more than 1 and meets decision condition, resampling is carried out.
4. a kind of subsurface structure imaging method of relief surface according to claim 3, it is characterised in that the judgement bar
Part, is specifically included：
Whether the remainder for judging iterations k divided by model modification direction controlling factor o is 1, if so, being weighed to observation data
New sampling, otherwise, not resampling.
5. the subsurface structure imaging method of a kind of relief surface according to claim 4, it is characterised in that described to adopt again
Sample, specific method includes：
Coefficient correlation in calculating remaining data per big gun data with last iteration hour offset data, chooses coefficient correlation maximum
NsIndividual single-shot data as current iteration offset data；
As the lazy weight N of single-shot data in remaining datasWhen, not enough part is randomly selected from the data calculated, this
In NsFor the big gun number of offset data after sampling.
6. the subsurface structure imaging method of a kind of relief surface according to claim 1, it is characterised in that in the step S4
Gradient is determined by data residual error, specifically included, stochastic gradient is determined according to data residual error：
WhereinRepresent stochastic gradient,Respectively kth time iteration when offset data and its corresponding inverse time inverse migration operator,
K is iterations；ForConjugate transposition, i.e. reverse-time migration operator；m(k-1)The reflection obtained for -1 iterative calculation of kth
Modulus Model.
7. a kind of subsurface structure imaging method of relief surface according to claim 6, it is characterised in that the step S4
In gradient is determined by data residual error, in addition to random conjugate gradient is determined by data residual error：Wherein
Random conjugate gradient modifying factorWhereinRepresent stochastic gradient,Random conjugate gradient is represented, k is represented
Iterations.
8. a kind of subsurface structure imaging method of relief surface according to claim 7, it is characterised in that the step S4
In reflectivity model iteration is carried out based on gradient, obtain final reflectivity model, specifically include following steps：
S41, based on decision condition, one in selection stochastic gradient or random conjugate gradient is model modification gradient, as anti-
Penetrate the direction of Modulus Model iteration；
S42, when described model modification gradient is less than threshold value, into step S44, when described model modification gradient is more than threshold
During value, step-length is updated according to model modification gradient calculation；
S43, based on step-length is updated, iteration updates reflectivity model, and return to step S2；
S44, judges the reflectivity model now obtained as final reflectivity model.
9. a kind of subsurface structure imaging method of relief surface according to claim 8, it is characterised in that the step
Decision condition in S41, is specifically included：Whether the remainder for judging iterations k divided by model modification direction controlling factor o is 1,
If so, then selection stochastic gradient is model modification gradient, if it is not, then selecting random conjugate gradient to be model modification gradient.
10. a kind of subsurface structure imaging method of relief surface according to claim 9, it is characterised in that the step
In S43 based on step-length is updated, iteration updates reflectivity model, and its circular is：
Wherein m(k)The reflectivity model obtained for kth time iterative calculation,α(k)Represent and update step-length, its
InRepresentative model updates gradient, including stochastic gradient or random conjugate gradient, []*Represent conjugate transposition,Represent kth
The corresponding inverse time inverse migration operator of offset data during secondary iteration.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201710337714.5A CN107193043B (en) | 2017-05-15 | 2017-05-15 | A kind of subsurface structure imaging method of relief surface |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201710337714.5A CN107193043B (en) | 2017-05-15 | 2017-05-15 | A kind of subsurface structure imaging method of relief surface |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN107193043A true CN107193043A (en) | 2017-09-22 |
| CN107193043B CN107193043B (en) | 2019-03-29 |
Family
ID=59873489
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201710337714.5A Expired - Fee Related CN107193043B (en) | 2017-05-15 | 2017-05-15 | A kind of subsurface structure imaging method of relief surface |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN107193043B (en) |
Cited By (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN108562937A (en) * | 2018-03-15 | 2018-09-21 | 东北石油大学 | A kind of seismic imaging method |
| CN108845355A (en) * | 2018-09-26 | 2018-11-20 | 中国矿业大学(北京) | Seismic migration imaging method and device |
| CN110531429A (en) * | 2019-08-02 | 2019-12-03 | 中国科学院电子学研究所 | A kind of time-domain electromagnetic data object inversion method based on supervision descent method |
Citations (8)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101082676A (en) * | 2007-07-11 | 2007-12-05 | 成都理工大学 | Earthquake post-stack forward method of undulating surface |
| CN101840000A (en) * | 2009-03-20 | 2010-09-22 | 中国石油天然气集团公司 | Plane wave pre-stack depth migration method |
| CN102116869A (en) * | 2011-02-12 | 2011-07-06 | 中国石油大学(华东) | High-precision prestack domain least square migration seismic imaging technology |
| US20120051176A1 (en) * | 2010-08-31 | 2012-03-01 | Chevron U.S.A. Inc. | Reverse time migration back-scattering noise removal using decomposed wavefield directivity |
| WO2014028030A1 (en) * | 2012-08-17 | 2014-02-20 | Landmark Graphics Corporation | Systems and methods for imaging seismic data |
| US20150355355A1 (en) * | 2012-11-13 | 2015-12-10 | Total E&P Usa, Inc. | Process for Creating Image Gathers |
| CN105301641A (en) * | 2015-10-30 | 2016-02-03 | 中国石油天然气集团公司 | Azimuthal anisotropy speed inversion method and apparatus |
| CN106033124A (en) * | 2016-06-29 | 2016-10-19 | 中国石油化工股份有限公司 | Multi-seismic resource sticky sound least square reverse time migration method based on stochastic optimization |
-
2017
- 2017-05-15 CN CN201710337714.5A patent/CN107193043B/en not_active Expired - Fee Related
Patent Citations (8)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101082676A (en) * | 2007-07-11 | 2007-12-05 | 成都理工大学 | Earthquake post-stack forward method of undulating surface |
| CN101840000A (en) * | 2009-03-20 | 2010-09-22 | 中国石油天然气集团公司 | Plane wave pre-stack depth migration method |
| US20120051176A1 (en) * | 2010-08-31 | 2012-03-01 | Chevron U.S.A. Inc. | Reverse time migration back-scattering noise removal using decomposed wavefield directivity |
| CN102116869A (en) * | 2011-02-12 | 2011-07-06 | 中国石油大学(华东) | High-precision prestack domain least square migration seismic imaging technology |
| WO2014028030A1 (en) * | 2012-08-17 | 2014-02-20 | Landmark Graphics Corporation | Systems and methods for imaging seismic data |
| US20150355355A1 (en) * | 2012-11-13 | 2015-12-10 | Total E&P Usa, Inc. | Process for Creating Image Gathers |
| CN105301641A (en) * | 2015-10-30 | 2016-02-03 | 中国石油天然气集团公司 | Azimuthal anisotropy speed inversion method and apparatus |
| CN106033124A (en) * | 2016-06-29 | 2016-10-19 | 中国石油化工股份有限公司 | Multi-seismic resource sticky sound least square reverse time migration method based on stochastic optimization |
Non-Patent Citations (1)
| Title |
|---|
| 李振春 等: ""地震偏移成像技术研究现状与发展趋势"", 《石油地球物理勘探》 * |
Cited By (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN108562937A (en) * | 2018-03-15 | 2018-09-21 | 东北石油大学 | A kind of seismic imaging method |
| CN108845355A (en) * | 2018-09-26 | 2018-11-20 | 中国矿业大学(北京) | Seismic migration imaging method and device |
| CN110531429A (en) * | 2019-08-02 | 2019-12-03 | 中国科学院电子学研究所 | A kind of time-domain electromagnetic data object inversion method based on supervision descent method |
Also Published As
| Publication number | Publication date |
|---|---|
| CN107193043B (en) | 2019-03-29 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CA2384810C (en) | Determining optimal well locations from a 3d reservoir model | |
| AU2015101990B4 (en) | Conditioning of object or event based reservoir models using local multiple-point statistics simulations | |
| CN105388518B (en) | A kind of centroid frequency and earthquake inversion of quality factor method in the united well of Frequency spectrum ratio | |
| CN107462924B (en) | A kind of absolute wave impedance inversion method independent of well-log information | |
| NO326598B1 (en) | Three-dimensional geological modeling | |
| CN106324662B (en) | A kind of full waveform inversion method and system for destination layer | |
| CN107193043A (en) | A kind of subsurface structure imaging method of relief surface | |
| RU2012152063A (en) | METHOD OF Q TOMOGRAPHY | |
| US6999879B2 (en) | Method for controlling seismic coverage using decision theory | |
| US20160003010A1 (en) | Method for operating a substerranean formation from which a fluid is produced | |
| Song et al. | Pattern search algorithms for nonlinear inversion of high-frequency Rayleigh-wave dispersion curves | |
| CA2908934C (en) | Gridless simulation of a fluvio-deltaic environment | |
| WO2017011469A1 (en) | Ensemble based decision making | |
| CN106033124A (en) | Multi-seismic resource sticky sound least square reverse time migration method based on stochastic optimization | |
| CN107894618A (en) | A kind of full waveform inversion gradient preprocess method based on model smoothing algorithm | |
| CN110031898A (en) | Data optimization methods and Kichhoff integral pre-stack depth migration method | |
| CN103543478A (en) | Geologic morphological interpolation KM (Kriging and Multiple-point geostatistics) method | |
| CN107024717A (en) | A kind of improved adaptive GA-IAGA for earthquake data before superposition parametric inversion | |
| Song et al. | Insights into performance of pattern search algorithms for high-frequency surface wave analysis | |
| Hicks Jr et al. | Identifying and quantifying significant uncertainties in basin modeling | |
| NO20131328A1 (en) | Computer estimation method and method for oil exploration and recovery using such a method | |
| Song et al. | Shuffled complex evolution approach for effective and efficient surface wave analysis | |
| CN105572732A (en) | Fracture-developed zone detection method based on gradual increase of attribute change rate | |
| US20170248718A1 (en) | Method of characterising a subsurface volume | |
| CN110515127A (en) | A kind of earthquake quality factor determines method, apparatus, equipment, medium |
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 | ||
| CB03 | Change of inventor or designer information |
Inventor after: Li Chuang Inventor after: Huang Jianping Inventor after: Li Zhenchun Inventor before: Li Chuang Inventor before: Huang Pingjian Inventor before: Li Zhenchun |
|
| CB03 | Change of inventor or designer information | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant | ||
| CF01 | Termination of patent right due to non-payment of annual fee |
Granted publication date: 20190329 Termination date: 20200515 |
|
| CF01 | Termination of patent right due to non-payment of annual fee |