CN106950596A - A kind of finite difference contrast source full waveform inversion method based on wavelet iterative estimate - Google Patents
A kind of finite difference contrast source full waveform inversion method based on wavelet iterative estimate Download PDFInfo
- Publication number
- CN106950596A CN106950596A CN201710231308.0A CN201710231308A CN106950596A CN 106950596 A CN106950596 A CN 106950596A CN 201710231308 A CN201710231308 A CN 201710231308A CN 106950596 A CN106950596 A CN 106950596A
- Authority
- CN
- China
- Prior art keywords
- contrast
- wavelet
- wave field
- source
- contrast source
- 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
- 238000000034 method Methods 0.000 title claims abstract description 86
- 238000004364 calculation method Methods 0.000 claims abstract description 13
- 238000004088 simulation Methods 0.000 claims description 24
- 239000011159 matrix material Substances 0.000 claims description 12
- 230000017105 transposition Effects 0.000 claims description 2
- 230000008569 process Effects 0.000 description 7
- 238000010276 construction Methods 0.000 description 6
- 238000000354 decomposition reaction Methods 0.000 description 5
- 238000001228 spectrum Methods 0.000 description 5
- 230000036039 immunity Effects 0.000 description 4
- 230000008901 benefit Effects 0.000 description 3
- 238000005553 drilling Methods 0.000 description 3
- 238000000205 computational method Methods 0.000 description 2
- 238000007796 conventional method Methods 0.000 description 2
- 239000004576 sand Substances 0.000 description 2
- 238000012360 testing method Methods 0.000 description 2
- 238000012795 verification Methods 0.000 description 2
- 229910017435 S2 In Inorganic materials 0.000 description 1
- 238000004458 analytical method Methods 0.000 description 1
- 230000008859 change Effects 0.000 description 1
- 238000004891 communication Methods 0.000 description 1
- 238000002939 conjugate gradient method Methods 0.000 description 1
- 230000021615 conjugation Effects 0.000 description 1
- 238000009795 derivation Methods 0.000 description 1
- 235000013399 edible fruits Nutrition 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 238000000605 extraction Methods 0.000 description 1
- 238000009499 grossing Methods 0.000 description 1
- 238000007689 inspection Methods 0.000 description 1
- 238000012545 processing Methods 0.000 description 1
- 238000005070 sampling Methods 0.000 description 1
- 230000009466 transformation 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
- G01V1/282—Application of seismic models, synthetic seismograms
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)
- Radar Systems Or Details Thereof (AREA)
- Geophysics And Detection Of Objects (AREA)
Abstract
Source full waveform inversion method is contrasted the invention discloses a kind of finite difference based on wavelet iterative estimate, be the described method comprises the following steps：S1, inputs initial velocity model, builds finite difference operator, and calculate background wave field；S2, calculates approximate focus ratio, thus calculates and obtains scattered field, and obtains contrast source initial value, contrast initial value and total wave field initial value；S3, each big gun of calculating initially contrasts the Grad corresponding to source, builds conjugated gradient direction；S4, asks for the step-length that every big gun contrast source updates, and be iterated renewal to contrast source；S5, judges that the size of target function value now and allowable error determines whether iteration stops；S6, calculates according to contrast value now and obtains spread speed of the ripple in medium.The method of the present invention is effectively overcome when used wavelet differs larger with true wavelet, the difficult point of inverting acquired results substantial deviation true model, and amount of calculation is significantly reduced.
Description
Technical field
Source full waveform inversion method is contrasted the present invention relates to a kind of finite difference based on wavelet iterative estimate, belongs to frequency
Domain acoustic full waveform inversion technique field.
Background technology
Dynamic characteristic extraction rate parameter of the full waveform inversion (FWI) according to geological data, is chromatographed during compared to travelling
It is expected to obtain more accurate subsurface velocities structure Deng the method based on data motion feature.At present, come from the mode of realization
See and be broadly divided into two major classes, time-domain and frequency domain.Wherein, several frequencies are only needed just to obtain higher anti-of precision in frequency domain
Result is drilled, and is easily achieved Multi-Scale Calculation and is promoted rapidly.The topmost amount of calculation of frequency domain full waveform inversion is
Substantial amounts of forward simulation, one huge system of linear equations of demand solution, because direct method (as LU is decomposed) has calculating many big guns
Advantage, therefore be widely used, its main amount of calculation is the decomposition operation of matrix, and amount of calculation is larger.Simultaneously because lacking
Accurate source wavelet information, when simulating wavelet and true wavelet has big difference, its inversion result will be completely offset from truly
Model.
The content of the invention
It is an object of the present invention to provide a kind of finite difference contrast source full waveform inversion side based on wavelet iterative estimate
Method, it can solve problem present in current techniques, significantly reduce amount of calculation, while can not accurately be obtained in source wavelet
When, still correct rate pattern result can be obtained with inverting.
In order to solve the above technical problems, the present invention is adopted the following technical scheme that：A kind of having based on wavelet iterative estimate
Differential pair is limited than source full waveform inversion method, is comprised the following steps：
S1, builds finite difference operator, and calculate background wave field according to initial velocity model；
S2, initial approximation focus is calculated using the simulation wave field and true wavefield corresponding to arbitrary a certain big gun at receiving point
Than, scattered field is derived from, and according to scattering field computation contrast source initial value, as at the beginning of contrast source calculation of initial value obtains contrast
Initial value and total wave field initial value；
S3, the Grad according to corresponding to each big gun of inversion objective function calculating initially contrasts source, builds conjugate gradient side
To；
S4, asks for the step-length that every big gun contrast source updates, and be iterated renewal to contrast source；
S5, the contrast source value after being updated according to iteration updates total wave field, contrast and approximate focus ratio, and judge now
Target function value and allowable error relation, if target function value be more than allowable error, with the contrast source after renewal, contrast
Degree, total wave field and approximate focus are used for initial value and are transferred to step S3, otherwise into step S6；
S6, calculates according to contrast value now and obtains the spread speed of ripple in media as well.
In a kind of foregoing finite difference contrast source full waveform inversion method based on wavelet iterative estimate, the step S2
In approximate focus be than λ
R ∈ S, whereinAnd uref(r) represent that foregoing a certain big gun is corresponding respectively to receive
The column vector that simulation and observed wave field are constituted at point.
In a kind of foregoing finite difference contrast source full waveform inversion method based on wavelet iterative estimate, the step S2
The computing formula of middle scattered field is：r∈S。
In a kind of foregoing finite difference contrast source full waveform inversion method based on wavelet iterative estimate, the step S2
Middle use anti-pass broadcasting method obtains contrast source initial value, and its computational methods is：
In a kind of foregoing finite difference contrast source full waveform inversion method based on wavelet iterative estimate, the step S2
Middle contrast value computational methods are：
Wherein, Re represents to take real,
uj,nIt is total wave field of reconstruct
In a kind of foregoing finite difference contrast source full waveform inversion method based on wavelet iterative estimate, the step S3
In, inversion objective function Fn is collectively formed by data field and inverting domain：
Its Linear Operators LbAct on contrast source and produce and dissipate wave field；Projection operator MSAnd MD, wherein MSRepresent overall area
Data projection on T is on data field S, and MDThen by the data projection on T to inverting domain D；djRepresent the scattering received
, ηSWithIt is expressed as being expressed as follows respectivelyηSWithPlay a part of equalization data domain and inverting domain weight.
In a kind of foregoing finite difference contrast source full waveform inversion method based on wavelet iterative estimate, the step S4
The step-length of middle nth iterationIf nth iteration is total to
Yoke gradient direction is pj,n, contrast source wjMore new formula be:
In a kind of foregoing finite difference contrast source full waveform inversion method based on wavelet iterative estimate, the step S6
The middle spread speed for calculating ripple calculates formula and is：R ∈ D, v (r) therein and vb(r) divide
Wei not spread speed of the ripple in true model and background model.
Compared with prior art, the present invention has advantages below：
1st, when used wavelet differs larger with true wavelet, the inverting acquired results that current techniques are present are seriously inclined
From true model, and the efficiency of inverse process basically identical with using true wavelet can be obtained using the method for the present invention, i.e. this hair
Bright method is very low to the dependence of wavelet.
2nd, compared to conventional full wave shape inversion method, the object that the method for the present invention updates is the contrast source of every big gun, and carries on the back
Scape model keeps constant in single-frequency refutation process, therefore need to only carry out the decomposition of an impedance matrix for a frequency, from
And significantly reduce the amount of calculation of forward simulation in refutation process.
3. the method that the present invention is used has certain noise immunity；When signal to noise ratio is 30dB, construction occurs in that discontinuous
Phenomenon；When signal to noise ratio is 25dB, rate pattern is disturbed seriously, but not substantial deviation realistic model, and this also embodies
The method of the present invention to data containing random noise has preferable stability.
Brief description of the drawings
Fig. 1 is the flow chart of one embodiment of the present of invention；
Fig. 2 is the initial velocity model in the embodiment of the present invention, and wherein a is the true initial velocity obtained by common-shot record
Model, b is the initial velocity model smoothly obtained by true model；
Fig. 3 is embodiment of the present invention neutron waveform and amplitude spectrum, wherein (a) true wavelet；(b) wavelet is simulated；(c)
The amplitude spectrum of true wavelet；(d) amplitude spectrum of wavelet is simulated；
Fig. 4 is 24Hz single-shots true wavefield and reconstruct wave field real part in the embodiment of the present invention, wherein (a) true wavefield；
(b) the reconstruct wave field of KSS methods；(c) the reconstruct wave field of IES methods；
Fig. 5 is 24Hz contrast source real parts and contrast in the embodiment of the present invention, wherein the single-shot contrast of (a) KSS methods
Source；(b) the single-shot contrast source of IES methods；(c) contrast of KSS methods；(d) contrast of IES methods；
Fig. 6 is final inversion result in the embodiment of the present invention, wherein (a) KSS methods；(b) IES methods；(c) using simulation
Wavelet and it is left intact；
Fig. 7 is single track rate curve contrast (a) x=1.5km in the embodiment of the present invention；(b) x=3.0km；
Fig. 8 is each frequency subbands relative error broken line graph in the embodiment of the present invention；
Fig. 9 is the big gun record figure of the random noise containing 35dB in the embodiment of the present invention；
Figure 10 is the inversion result under different signal to noise ratios in the embodiment of the present invention, (a) 40dB；(b)35dB；(c)30dB；
(d)25dB；
Figure 11 is gradient initial model (a) and inversion result (b) in the embodiment of the present invention.
Figure 12 is a kind of flow chart of embodiment of the present invention.
The present invention is further illustrated with reference to the accompanying drawings and detailed description.
Embodiment
In a practical situation, the common big gun collection record (time-domain observes data) field acquisition obtained, is become using Fourier
Get the observed wave field of each frequency at detector position in return, the wavelet that dynamite source is excited is exactly true wavelet.And in inverting
A simulation wavelet is needed also exist for, for generating simulation wave field, but true wavelet is not known generally, so artificially set
Simulation wavelet is easily caused inverting failure when being differed greatly with true wavelet, the present invention solves this using the method for wavelet iteration
Problem.
A kind of embodiment of the present invention：A kind of finite difference contrast source full waveform inversion side based on wavelet iterative estimate
Method, as shown in figure 1, comprising the following steps：
S1, builds finite difference operator, and calculate background wave field according to initial velocity model；In the present embodiment, grid is big
Small is 101 × 251, and transverse and longitudinal spatial sampling interval is all 15m.Focus and wave detector are located at earth's surface, and the first bombard is located at 60m, big gun
Spacing is 75m, totally 49 big gun, and wave detector is set to earth's surface and received entirely, so as to obtain common big gun collection record, is obtained just by velocity analysis
Beginning rate pattern, as shown in Figure 2 a, then to common big gun collection record do the sight that Fourier transformation obtains each frequency at detector position
Survey wave field.Observed wave field is the true wavefield received at detector position.It can only be obtained in seismic prospecting at detector position
Wave field, and the true wavefields in the most places in underground are unknown.
In isotropic medium, the following Frequency-Space Domain Helmholtz equations of communication satisfaction of seismic wave
Wherein, uj(r) represent by true focus f1The total wave field produced,The locus of focus is represented, k (r) is space
Wave number at the r of position.Using the impedance matrix H of 9 difference scheme tectonic setting models of mixingb, and then calculate linear operatorkbThe wave number of background model is represented, k is thus calculatedb, and then obtain background wave fieldIt is wherein first anti-
With initial velocity model as background model when drilling, the rate pattern conduct obtained in inverting thereafter with previous frequency inversion
Background model.
S2, initial approximation focus is calculated using the simulation wave field and observed wave field corresponding to arbitrary a certain big gun at receiving point
Than, scattered field is derived from, and according to scattering field computation contrast source initial value, as at the beginning of contrast source calculation of initial value obtains contrast
Initial value and total wave field initial value, simulation wave field can be obtained by the simulation wavelet inputted.
In actual seismic exploration, true wavelet is generally difficult to accurate acquisition, true wavelet f1With simulation wavelet f2Not phase
Deng by f2Reasons for its use wave fieldDirectly difference can not be made with total wave field and obtain scattered field, approximate focus is defined for this and compares λ
For
In formula whereinThe column vector that simulation wave field is constituted at receiving point corresponding to the foregoing a certain big gun of expression, by just
Drill and obtain, uref(r) column vector then constituted for observed wave field, asterisk " * " represents conjugate transposition, r representation spaces position, and S is inspection
Ripple device region, it is initial approximate focus ratio that iteration calculates obtained λ before starting.
When background model is close to true model, λ ≈ f2/f1.So, the scattered field d now receivedjFor
, wherein j represents shot position, uj(r) observed wave field is represented,Represent background wave field.
When each frequency starts inverting, we it is also required to provide the initial contrast source of each big gun, that is, contrast source
Initial value.The initial value in contrast source is obtained using anti-pass broadcasting method herein：
Wherein, source w is contrastedj(r) by contrast χ (r) and total wave field uj(r) constitute
wj(r)=χ (r) uj(r) (5)
Contrast is expressed as follows
V (r) and v in formulab(r) it is respectively spread speed of the ripple in truth and background model.
Following contrast more new formula is recycled to obtain primary contrast's value (when taking n=0),
Wherein, Re represents to take real, uj,nIt is total wave field of reconstruct, when taking n=0, as total wave field initial value,
The initial contrast source and (8) then calculated by formula (4) calculates obtained total wave field equivalent to by simulation focus f2
Produce.Initial approximate focus is substantially correctly oriented than simply providing one for inverting, it is necessary to count again after model iteration
λ is calculated to continue the renewal of model, its method is used after model iterationRecalculate λ.
S3, under the framework of contrast source inversion method, is referred to as D, wave detector region is referred to as by the target area of exploration
Data field S, both constitutes overall area T.
The object function of inverting is collectively formed by data field and inverting domain
In above formula, linear operator LbAct on contrast source and produce and dissipate wave field；Projection operator MSAnd MD, wherein MSRepresenting will be total
Data projection on the T of region is on data field S, and MDThen by the data projection on T to inverting domain D；djRepresent that what is received dissipates
Penetrate field, ηSWithEffect with equalization data domain and inverting domain weight, it is ensured that the stable effective progress of inverting, they are represented respectively
It is as follows
It can be seen that contrast source inversion method is very different with conventional full wave shape inverting from object function, one is determined
, it is necessary to which the parameter updated is contrast source and contrast after individual frequency, the renewal of speed parameter is not related to.That is, working as
Only need to carry out its impedance matrix just drilling in all iteration under LU decomposition, current frequency after one background model of selection
The result decomposed for the first time can be used by calculating, so as to significantly reduce amount of calculation, improve inverting efficiency.
Contain two class variable (χ, w in object functionj), it is constant to assume initially that χ, to wjIt is updated.In order to simplify table
Show, the data field error for defining nth iteration is ρj,n, and the inverting domain error of nth iteration is rj,n, then object function pair
wjDerivation, can obtain each big gun wjGradient required for updating：
In above formula, upper line expression takes complex conjugate,Represent LbAdjoint operator, MS*Represent the matrix in data field
Projection is to overall area T, similarly, MD*Represent the matrix projection in inverting domain to overall area T；It is worth noting that, impedance matrix
Decomposition result may also be used for build operatorAnd need not decompose again.
Conjugate gradient method is used in the present embodiment, is tried to achieve after gradient direction, construction Polak-Ribiere directions are used as conjugation
Gradient direction Pj(i.e. CG directions Pj)。
S4, asks for the step-length that every big gun contrast source updates, and be iterated renewal to contrast source；
It is the step-length of nth iteration, orderThen
If the conjugated gradient direction of nth iteration is pj,n, then wjMore new formula it is as follows
S5, the contrast source value after being updated according to iteration updates total wave field by formula (8)；The Criterion of Selecting of contrast is to make
Object function totally reaches minimum, and only hasIn include contrast, by rightMoleculeMinimize, it is 0 to make it, the more new formula (7) of contrast is obtained, while according to public affairs
Formula (2) updates approximate focus ratio, updates after contrast and approximate focus ratio, judges target function value now and allowable error
Relation, if target function value is more than allowable error, is compared to the contrast source after renewal, contrast, total wave field and approximate focus
Step S3 is transferred to for initial value, otherwise into step S6, allowable error here can be according to different as an empirical value
Situation is set.
S6, calculates according to contrast value now and obtains the spread speed of ripple in media as well, computing formula is as follows：
V (r) and v in formulab(r) it is respectively spread speed in the spread speed and background model of ripple in media as well.
When used wavelet differs larger with true wavelet, the inverting acquired results substantial deviation that current techniques are present
True model.And direction of the method for the present invention using initial approximate focus than providing a renewal for inverting, make inverting court
The direction being in the main true to carry out, after model obtains renewal, the simulation wave field of generation is closer to truth, so as to approximately shake
Source is than being also more nearly real focus ratio.Utilize more accurately approximate focus ratio, the model that can obtain inverting
It is more accurate, it is achieved thereby that a kind of benign circulation.In the case that wavelet difference is larger, it is also possible to obtain correct
Inversion result, so as to reduce the dependence to wavelet.
Embodiment 2
As an alternative embodiment of the invention, inventor by the solution of the present invention by being applied to overthrust models
(i.e. Overthrust model), as shown in figure 12, further demonstrates the method for the present invention, initial smooth rate pattern (such as Fig. 2 b) by
Overthrust model smoothings, while being used to generate observation using the first derivative of Gaussian function as true wavelet (Fig. 3 a)
Wave field, Fig. 3 c are its corresponding amplitude spectrums.And the simulation wavelet that background wave field is generated in inverting is rake that dominant frequency is 10Hz
Ripple (Fig. 3 b), from its amplitude spectrum (Fig. 3 d) it can be seen that the frequency band of true wavelet is significantly wider than simulation wavelet, and both amplitudes
The contrast source inversion method that differs greatly is different from conventional full waveform inversion, and its object updated is the contrast source of every big gun, then
Scattered field, scattered field and the reconstruct wave field synthesized by ambient field and the more high then inverting knot of true wavefield degree of agreement are updated with this
Fruit is better.By inputting initial velocity model, finite difference operator is built, and calculate background wave field；Then approximate focus is calculated
Than thus calculating and obtaining scattered field, and obtain contrast source initial value, contrast initial value and total wave field initial value；Calculate each
Big gun initially contrasts the Grad corresponding to source, builds conjugated gradient direction；The step-length that every big gun contrast source updates further is asked for, and
Renewal is iterated to contrast source；Target function value and the size of allowable error according to judging now determine whether iteration stops
Only；Finally calculated when the iterations cease according to contrast value now and obtain spread speed of the ripple in medium.And by what is obtained
Inversion result and inversion method conventional at present compare, and detailed comparison scheme is as follows：
Three kinds of different inversion methods are applied to overthrust models：1) it is used for inverting, note using true wavelet
For KSS (know source signature)；2) from simulation wavelet, the finite difference pair based on wavelet iterative estimate is used
Than source full waveform inversion strategy, be designated as IES (iterative estimation ofsource signature), i.e., it is of the invention
Scheme；3) from simulation wavelet, and without any processing.Used frequency be all from 3Hz to 30Hz, using 3Hz as
Interval, totally 10 frequencies carry out inverting, the inversion result of low frequency as adjacent high frequency input.
25th big gun corresponding true wavefield real part when Fig. 4 a are 24Hz, Fig. 4 b and Fig. 4 c is that KSS and IES methods institute is right respectively
The real part for the reconstruct wave field answered, the reconstruct wave field for noting IES methods is to compare gained using the total wave field divided by approximate focus that update,
It is different when this point is from known wavelet.It is seen that the reconstruct wave field of two methods and true wavefield are closely, this explanation
Inverting must compare smoothly.The contrast source real part of KSS and the same big gun of IES methods when Fig. 5 (a-b) is 24Hz, its shape is not in
Rule, but remain to find out the vestige of wave field in Fig. 4.After the contrast source of all big guns is obtained with reconstruct wave field, Fig. 5 can be obtained
(c-d) contrast (Fig. 5 c shown in：KKS；Fig. 5 d：IES).
Fig. 6 (a-c) is the inversion result of three kinds of situations, and the inversion result when wavelet is correct is best (Fig. 6 a), Neng Gouzhun
The tectonic information really reflected in true model, this also embodies the validity of contrast source inversion method.Fig. 6 b are the anti-of IES methods
Result is drilled, it is basically identical with inversion result of wavelet when correct, it is simply slightly poor in deep-level location convergence situation.Due to true wavelet
Have big difference with simulation wavelet, when the refutation strategy without using wavelet iterative estimate, its inversion result (Fig. 6 c) is completely offset from
True model, so that it can also be seen that importance of the wavelet in full waveform inversion.
Each rate curve extracted at 1.5km and 3.0km from Fig. 6 (a-b), draws the single track speed as shown in Fig. 7 (a-b)
Spend curve comparison figure.The rate curve of two methods is essentially coincided, and is coincide preferably with true velocity, is especially become in construction
It is not very violent layer position to change.We define the relative error function of wavelet
Wherein, λtrueIt is the ratio for simulating wavelet and true wavelet.Fig. 8 is that each frequency is initial and final wavelet is relative by mistake
Poor broken line graph.From initial relative error, it can be seen that the approximate focus ratio as defined in formula (7) is relatively true focus
Ratio, this is also the premise that IES methods are capable of correct inverting.With the progress of iteration, model becomes closer to truly
Situation, so that wavelet is also updated, its relative error is less and less.
In addition, inventor has also done some other checking work：
1. the noise immunity of verification method, random white noise is added in being recorded using the AWGN functions in Matlab to big gun, point
Not Huo get signal to noise ratio be 40dB, 35dB, 30dB and 25dB synthetic seismic data, Fig. 9 be signal to noise ratio be 35dB certain single-shot remember
Record, its less lineups of medium and deep amplitude are fallen into oblivion by noise substantially.
The inversion result of IES methods using noisy data test above, wherein, inversion result during 40dB and 35dB
(Figure 10 a-b) in addition to the high speed position of fault convergence at middle part is slightly poor, model can restore truth substantially, embody
IES methods have good noise immunity.When signal to noise ratio is 30dB, due to the influence of noise, many places construction occurs in that discontinuous
Phenomenon (Figure 10 c), but all in all, inversion result is nevertheless not deviating from realistic model.When signal to noise ratio is 25dB, construction is disconnected
Discontinuously continuous phenomenon is serious, but relative to the inversion result in Fig. 6 c, Figure 10 d are not absorbed in local minimum seriously, embodies
The refutation strategy of wavelet iterative estimate is to random noise and insensitive.
2. verification method is to the adaptability of initial model.The initial model selected in above-mentioned test is smooth in Fig. 2 b
Model, full waveform inversion is stronger to the dependence of initial model, especially with wavelet iterative estimate refutation strategy when, initially
Model not only affects the nonlinear degree of inversion problem, also determine initial approximation focus than accuracy.We adopt below
Initial model (Figure 11 a), model top speed 2700m/s, bottom speed 6000m/ are used as with non-linear stronger gradient former
S, middle even variation.IES inversion methods result as shown in figure 11b, its efficiency of inverse process with it is basically identical in Fig. 6 b, this is not only anti-
Reflected contrast source inversion method it is relatively low to the dependence of initial model, the refutation strategy for also embodying wavelet iterative estimate exists
Still there is good stability when initial model is poor.
The present invention is on the basis of finite difference contrast source full waveform inversion, by the inverting plan for introducing wavelet iterative estimate
Slightly, it is proposed that a kind of iteration efficiency high and the frequency domain full waveform inversion method small to wavelet dependence, have the following advantages that：
(1) compared to conventional method, the object that this method updates is the contrast source of every big gun, and background model is protected in single-frequency refutation process
Hold constant, therefore need to only the decomposition of an impedance matrix be carried out for a frequency, so as to significantly reduce in refutation process just
Drill the amount of calculation of simulation；(2) when simulation wavelet differs larger with true wavelet used in routine techniques, conventional method inverting
Acquired results are by substantial deviation true model, and the refutation strategy of the present invention can be obtained with using true wavelet basically identical
Efficiency of inverse process.(3) found by the generated data inversion result of relatively more different signal to noise ratios：Signal to noise ratio be 40dB and 35dB when to anti-
Drilling result influences smaller, and demonstrating this method has certain noise immunity；And during 30dB, construction occurs in that discontinuous phenomenon；
When signal to noise ratio is 25dB, rate pattern is disturbed seriously, but not substantial deviation realistic model, and this also embodies this method
There is preferable stability to data containing random noise；(4) inversion result of gradient initial model and smooth initial model is basic
Unanimously, i.e. this method is relatively low to initial model dependence.
Claims (8)
1. a kind of finite difference contrast source full waveform inversion method based on wavelet iterative estimate, it is characterised in that including following
Step：
S1, builds finite difference operator, and calculate background wave field according to initial velocity model；
S2, initial approximate focus is calculated using the simulation wave field and observed wave field corresponding to arbitrary a certain big gun at receiving point
Than thus calculating and obtaining scattered field, and according to scattering field computation contrast source initial value, contrasted by contrasting source calculation of initial value
Spend initial value and total wave field initial value；
S3, according to inversion objective function, each big gun of calculating initially contrasts the Grad corresponding to source, builds conjugated gradient direction；
S4, asks for the step-length that every big gun contrast source updates, and be iterated renewal to contrast source；
S5, the contrast source value after being updated according to iteration updates total wave field, contrast and approximate focus ratio, and judge mesh now
The relation of offer of tender numerical value and allowable error, if target function value be more than allowable error, with the contrast source after renewal, contrast,
Total wave field and approximate focus are used for initial value and are transferred to step S3, otherwise into step S6；
S6, calculates according to contrast value now and obtains the spread speed of ripple in media as well.
2. a kind of finite difference contrast source full waveform inversion method based on wavelet iterative estimate according to claim 1,
Characterized in that, approximate focus is than λ in the step S2
R ∈ S, whereinSimulation wave field composition at receiving point corresponding to the foregoing a certain big gun of expression
Column vector, uref(r) column vector then constituted for the observed wave field of a certain big gun, asterisk " * " represents conjugate transposition, r representation spaces
Position, S is wave detector region.
3. a kind of finite difference contrast source full waveform inversion method based on wavelet iterative estimate according to claim 2,
Scattered field d is calculated with the following method characterized in that, being adopted in the step S2j(r)：R ∈ S, its
Middle j represents shot position, uj(r) observed wave field is represented,Represent background wave field.
4. a kind of finite difference contrast source full waveform inversion method based on wavelet iterative estimate according to claim 3,
Characterized in that, obtaining contrast source w using anti-pass broadcasting method in the step S2jInitial value, specific calculation is as follows：
Wherein,Represent linear operator LbAdjoint operator, MSRepresent on the matrix projection in overall area T to data field S, MS*
Represent the matrix projection in data field S to overall area T；LbBy background model impedance matrix HbInverse and background wave number kbJointly
Constitute,
5. a kind of finite difference contrast source full waveform inversion method based on wavelet iterative estimate according to claim 3,
Characterized in that, in the step S2, when taking n=0, calculating total wave field initial valueWith
And contrast initial valueWherein, Re represents to take real, MDRepresent in overall area T
Matrix projection on the D of inverting region, upper line represents complex conjugate, and n refers to the number of times of iteration.
6. a kind of finite difference contrast source full waveform inversion method based on wavelet iterative estimate according to claim 5,
Characterized in that, in the step S3, inversion objective function Fn is by data field functional FSWith inverting domain functionalCollectively form：Its Linear Operators LbAct on the production of contrast source
Raw scattered wave field；djRepresent the scattered field received, ηSWithIt is expressed as follows respectively：
7. a kind of finite difference contrast source full waveform inversion method based on wavelet iterative estimate according to claim 6, its
It is characterised by, the step-length of nth iteration in the step S4
Wherein the conjugated gradient direction of nth iteration is pj,n,It is gradient direction, adopts and update contrast source with the following method:
8. a kind of finite difference contrast source full waveform inversion method based on wavelet iterative estimate according to claim 7,
Characterized in that, the step S6 medium waves are in the spread speed calculating method of medium：
R ∈ D, v (r) therein and vb(r) it is respectively spread speed of the ripple in the spread speed and background model of medium.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201710231308.0A CN106950596B (en) | 2017-04-11 | 2017-04-11 | A kind of finite difference comparison source full waveform inversion method based on wavelet iterative estimate |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201710231308.0A CN106950596B (en) | 2017-04-11 | 2017-04-11 | A kind of finite difference comparison source full waveform inversion method based on wavelet iterative estimate |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN106950596A true CN106950596A (en) | 2017-07-14 |
| CN106950596B CN106950596B (en) | 2019-02-26 |
Family
ID=59474578
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201710231308.0A Expired - Fee Related CN106950596B (en) | 2017-04-11 | 2017-04-11 | A kind of finite difference comparison source full waveform inversion method based on wavelet iterative estimate |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN106950596B (en) |
Cited By (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN108508482A (en) * | 2018-05-25 | 2018-09-07 | 中国海洋石油集团有限公司 | A kind of subterranean fracture seismic scattering response characteristic analogy method |
| CN109239771A (en) * | 2018-08-10 | 2019-01-18 | 杭州电子科技大学 | A kind of elastic wave imaging method based on non-homogeneous background media |
| CN109541681A (en) * | 2017-09-22 | 2019-03-29 | 中国海洋大学 | A kind of waveform inversion method of streamer seismic data and a small amount of OBS data aggregate |
| CN112444870A (en) * | 2019-08-30 | 2021-03-05 | 中国石油化工股份有限公司 | Multi-scale finite difference contrast source inversion method based on backscattering theory |
Citations (14)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN103091711A (en) * | 2013-01-24 | 2013-05-08 | 中国石油天然气集团公司 | Method and device for full-wave-shape inversion |
| CN103135132A (en) * | 2013-01-15 | 2013-06-05 | 中国科学院地质与地球物理研究所 | Hybrid-domain full wave form inversion method of central processing unit (CPU)/graphics processing unit (GPU) synergetic parallel computing |
| CN103713315A (en) * | 2012-09-28 | 2014-04-09 | 中国石油化工股份有限公司 | Seismic anisotropy parameter full waveform inversion method and device |
| US20150012221A1 (en) * | 2013-07-08 | 2015-01-08 | Reeshidev Bansal | Full-Wavefield Inversion of Primaries and Multiples in Marine Environment |
| CN104570082A (en) * | 2013-10-29 | 2015-04-29 | 中国石油化工股份有限公司 | Extraction method for full waveform inversion gradient operator based on green function characterization |
| CN105005076A (en) * | 2015-06-02 | 2015-10-28 | 中国海洋石油总公司 | Seismic wave full waveform inversion method based on least square gradient update speed model |
| CN105319581A (en) * | 2014-07-31 | 2016-02-10 | 中国石油化工股份有限公司 | Efficient time domain full waveform inversion method |
| CN105445798A (en) * | 2014-08-21 | 2016-03-30 | 中国石油化工股份有限公司 | Full waveform inversionmethod and system based on gradient processing |
| CN105467444A (en) * | 2015-12-10 | 2016-04-06 | 中国石油天然气集团公司 | An elastic wave full-waveform inversion method and apparatus |
| CN105549080A (en) * | 2016-01-20 | 2016-05-04 | 中国石油大学(华东) | Undulating ground surface waveform inversion method based on auxiliary coordinate system |
| CN105676277A (en) * | 2015-12-30 | 2016-06-15 | 中国石油大学(华东) | Full waveform joint inversion method for improving high steep structure velocity inversion efficiency |
| US20160216389A1 (en) * | 2015-01-23 | 2016-07-28 | Advanced Geophysical Technology Inc. | Beat tone full waveform inversion |
| CN106501852A (en) * | 2016-10-21 | 2017-03-15 | 中国科学院地质与地球物理研究所 | A kind of multiple dimensioned full waveform inversion method of three-dimensional acoustic wave equation arbitrarily-shaped domain and device |
| CN106526674A (en) * | 2016-11-14 | 2017-03-22 | 中国石油化工股份有限公司 | Three-dimensional full waveform inversion energy weighted gradient preprocessing method |
-
2017
- 2017-04-11 CN CN201710231308.0A patent/CN106950596B/en not_active Expired - Fee Related
Patent Citations (14)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN103713315A (en) * | 2012-09-28 | 2014-04-09 | 中国石油化工股份有限公司 | Seismic anisotropy parameter full waveform inversion method and device |
| CN103135132A (en) * | 2013-01-15 | 2013-06-05 | 中国科学院地质与地球物理研究所 | Hybrid-domain full wave form inversion method of central processing unit (CPU)/graphics processing unit (GPU) synergetic parallel computing |
| CN103091711A (en) * | 2013-01-24 | 2013-05-08 | 中国石油天然气集团公司 | Method and device for full-wave-shape inversion |
| US20150012221A1 (en) * | 2013-07-08 | 2015-01-08 | Reeshidev Bansal | Full-Wavefield Inversion of Primaries and Multiples in Marine Environment |
| CN104570082A (en) * | 2013-10-29 | 2015-04-29 | 中国石油化工股份有限公司 | Extraction method for full waveform inversion gradient operator based on green function characterization |
| CN105319581A (en) * | 2014-07-31 | 2016-02-10 | 中国石油化工股份有限公司 | Efficient time domain full waveform inversion method |
| CN105445798A (en) * | 2014-08-21 | 2016-03-30 | 中国石油化工股份有限公司 | Full waveform inversionmethod and system based on gradient processing |
| US20160216389A1 (en) * | 2015-01-23 | 2016-07-28 | Advanced Geophysical Technology Inc. | Beat tone full waveform inversion |
| CN105005076A (en) * | 2015-06-02 | 2015-10-28 | 中国海洋石油总公司 | Seismic wave full waveform inversion method based on least square gradient update speed model |
| CN105467444A (en) * | 2015-12-10 | 2016-04-06 | 中国石油天然气集团公司 | An elastic wave full-waveform inversion method and apparatus |
| CN105676277A (en) * | 2015-12-30 | 2016-06-15 | 中国石油大学(华东) | Full waveform joint inversion method for improving high steep structure velocity inversion efficiency |
| CN105549080A (en) * | 2016-01-20 | 2016-05-04 | 中国石油大学(华东) | Undulating ground surface waveform inversion method based on auxiliary coordinate system |
| CN106501852A (en) * | 2016-10-21 | 2017-03-15 | 中国科学院地质与地球物理研究所 | A kind of multiple dimensioned full waveform inversion method of three-dimensional acoustic wave equation arbitrarily-shaped domain and device |
| CN106526674A (en) * | 2016-11-14 | 2017-03-22 | 中国石油化工股份有限公司 | Three-dimensional full waveform inversion energy weighted gradient preprocessing method |
Non-Patent Citations (4)
| Title |
|---|
| A ABUBAKAR 等: ""A finite-difference contrast source inversion method"", 《INVERSION PROBLEMS》 * |
| 李庆洋 等: ""不依赖子波、基于包络的WFI初始模型建立方法研究"", 《地球物理学报》 * |
| 段晓亮 等: ""基于逆散射理论的地震波速度正则化反演"", 《物理学报》 * |
| 邵婕 等: ""地震波散射理论及应用研究进展"", 《地球物理学进展》 * |
Cited By (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN109541681A (en) * | 2017-09-22 | 2019-03-29 | 中国海洋大学 | A kind of waveform inversion method of streamer seismic data and a small amount of OBS data aggregate |
| CN108508482A (en) * | 2018-05-25 | 2018-09-07 | 中国海洋石油集团有限公司 | A kind of subterranean fracture seismic scattering response characteristic analogy method |
| CN109239771A (en) * | 2018-08-10 | 2019-01-18 | 杭州电子科技大学 | A kind of elastic wave imaging method based on non-homogeneous background media |
| CN109239771B (en) * | 2018-08-10 | 2020-01-31 | 杭州电子科技大学 | elastic wave imaging method based on non-uniform background medium |
| CN112444870A (en) * | 2019-08-30 | 2021-03-05 | 中国石油化工股份有限公司 | Multi-scale finite difference contrast source inversion method based on backscattering theory |
Also Published As
| Publication number | Publication date |
|---|---|
| CN106950596B (en) | 2019-02-26 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN108345031B (en) | Full waveform inversion method for elastic medium active source and passive source mixed mining seismic data | |
| CN105388518B (en) | A kind of centroid frequency and earthquake inversion of quality factor method in the united well of Frequency spectrum ratio | |
| CN103238158B (en) | Utilize the marine streamer data source inverting simultaneously that mutually related objects function is carried out | |
| CN106950596B (en) | A kind of finite difference comparison source full waveform inversion method based on wavelet iterative estimate | |
| CN110031896B (en) | Seismic random inversion method and device based on multi-point geostatistics prior information | |
| CN106526674A (en) | Three-dimensional full waveform inversion energy weighted gradient preprocessing method | |
| CN110058302A (en) | A kind of full waveform inversion method based on pre-conditional conjugate gradient accelerating algorithm | |
| US20150247940A1 (en) | History matching of time-lapse crosswell data using ensemble kalman filtering | |
| RU2008150484A (en) | ELECTROMAGNETIC SURVEY | |
| Alielahi et al. | Applying a time-domain boundary element method for study of seismic ground response in the vicinity of embedded cylindrical cavity | |
| KR20090075843A (en) | Iterative inversion of data from simultaneous geophysical sources | |
| CN103105623B (en) | Data waveform processing method in seismic exploration | |
| CN105093278B (en) | Full waveform inversion gradient operator extracting method based on the main energy-optimised algorithm of excitation | |
| CN108645994A (en) | A kind of geology stochastic inversion methods and device based on Multiple-Point Geostatistics | |
| CN108387933A (en) | A kind of method, apparatus and system of definitely interval quality factors | |
| CN107894618B (en) | A kind of full waveform inversion gradient preprocess method based on model smoothing algorithm | |
| CN106569262B (en) | Background velocity model reconstruction method under low frequency seismic data missing | |
| CN107462924A (en) | A kind of absolute wave impedance inversion method independent of well-log information | |
| CN110187384A (en) | Bayes's time-lapse seismic difference inversion method and device | |
| CN113156493A (en) | Time-frequency domain full-waveform inversion method and device using normalized seismic source | |
| Yin et al. | Improving horizontal resolution of high-frequency surface-wave methods using travel-time tomography | |
| CN105911584A (en) | Implicit staggered-grid finite difference elastic wave numerical simulation method and device | |
| CN115600373A (en) | Viscous anisotropic medium qP wave simulation method, system, equipment and application | |
| KR20160012922A (en) | Seismic imaging apparatus and method using iterative direct waveform inversion | |
| Yablokov et al. | Uncertainty quantification of multimodal surface wave inversion using artificial neural networks |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant | ||
| CF01 | Termination of patent right due to non-payment of annual fee | ||
| CF01 | Termination of patent right due to non-payment of annual fee |
Granted publication date: 20190226 Termination date: 20200411 |