CN108680957A - Local cross-correlation time-frequency domain Phase-retrieval method based on weighting - Google Patents
Local cross-correlation time-frequency domain Phase-retrieval method based on weighting Download PDFInfo
- Publication number
- CN108680957A CN108680957A CN201810490671.9A CN201810490671A CN108680957A CN 108680957 A CN108680957 A CN 108680957A CN 201810490671 A CN201810490671 A CN 201810490671A CN 108680957 A CN108680957 A CN 108680957A
- Authority
- CN
- China
- Prior art keywords
- frequency domain
- time
- weighting
- data
- phase
- 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 85
- 238000012360 testing method Methods 0.000 claims abstract description 15
- 230000009466 transformation Effects 0.000 claims description 16
- 238000004088 simulation Methods 0.000 claims description 15
- 238000004422 calculation algorithm Methods 0.000 claims description 12
- 238000005457 optimization Methods 0.000 claims description 10
- 239000011159 matrix material Substances 0.000 claims description 6
- 238000000844 transformation Methods 0.000 claims description 6
- 238000001914 filtration Methods 0.000 claims description 4
- 238000012986 modification Methods 0.000 claims description 4
- 230000004048 modification Effects 0.000 claims description 4
- 230000008901 benefit Effects 0.000 abstract description 5
- 230000006870 function Effects 0.000 description 44
- 230000008569 process Effects 0.000 description 8
- 238000010276 construction Methods 0.000 description 4
- 230000008859 change Effects 0.000 description 3
- 238000010606 normalization Methods 0.000 description 3
- 238000001228 spectrum Methods 0.000 description 3
- 238000004458 analytical method Methods 0.000 description 2
- 244000145845 chattering Species 0.000 description 2
- 238000000354 decomposition reaction Methods 0.000 description 2
- 230000005624 perturbation theories Effects 0.000 description 2
- 238000010587 phase diagram Methods 0.000 description 2
- 238000011160 research Methods 0.000 description 2
- 241001269238 Data Species 0.000 description 1
- 230000003044 adaptive effect Effects 0.000 description 1
- 238000013459 approach Methods 0.000 description 1
- 238000004364 calculation method Methods 0.000 description 1
- 238000012937 correction Methods 0.000 description 1
- 238000005314 correlation function Methods 0.000 description 1
- 238000013461 design Methods 0.000 description 1
- 230000000368 destabilizing effect Effects 0.000 description 1
- 238000010586 diagram Methods 0.000 description 1
- 230000000694 effects Effects 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 238000007306 functionalization reaction Methods 0.000 description 1
- 230000006872 improvement Effects 0.000 description 1
- 230000010354 integration Effects 0.000 description 1
- 238000005259 measurement Methods 0.000 description 1
- 230000003071 parasitic effect Effects 0.000 description 1
- 238000011056 performance test Methods 0.000 description 1
- 238000000614 phase inversion technique Methods 0.000 description 1
- 230000007480 spreading Effects 0.000 description 1
- 238000010998 test method Methods 0.000 description 1
- 238000003325 tomography Methods 0.000 description 1
- 238000005303 weighing 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/30—Analysis
- G01V1/307—Analysis for determining seismic attributes, e.g. amplitude, instantaneous phase or frequency, reflection strength or polarity
-
- 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/36—Effecting static or dynamic corrections on records, e.g. correcting spread; Correlating seismic signals; Eliminating effects of unwanted energy
-
- 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/36—Effecting static or dynamic corrections on records, e.g. correcting spread; Correlating seismic signals; Eliminating effects of unwanted energy
- G01V1/364—Seismic filtering
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V2210/00—Details of seismic processing or analysis
- G01V2210/30—Noise handling
- G01V2210/32—Noise reduction
- G01V2210/324—Filtering
-
- 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/63—Seismic attributes, e.g. amplitude, polarity, instant phase
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 local cross-correlation time-frequency domain Phase-retrieval method based on weighting, is that the time-frequency domain phase information of seismic data is introduced into mutually related objects function, is referred to as the local cross-correlation time-frequency domain Phase-retrieval method based on weighting.Dependence of the full waveform inversion to initial velocity model is alleviated by the introducing of phase information.Weight factor is added in time-frequency domain object function simultaneously, the stability of noise resisting ability and inverting can be enhanced well.In the local cross-correlation time-frequency domain Phase-retrieval method of the weighting of low-frequency range, it can obtain a good initial velocity model, then in the local cross-correlation time-frequency domain Phase-retrieval method for using it for high frequency band weighting, high-resolution inversion result may finally be obtained.It lacks at the same time in the test of low frequency component and strong Gaussian Background noise and shows that the local cross-correlation time-frequency domain Phase-retrieval method based on weighting has many advantages, such as stronger anti-noise ability and do not depend on initial velocity model.
Description
Technical field
The present invention relates to a kind of geophysical exploration methods, are built especially with seismic data time-frequency domain phase information
The object function of rate pattern inverting, and energy is carried out to phase object function using the amplitude information of time-frequency domain analogue data
The local cross-correlation time-frequency domain Phase-retrieval for the problems such as destabilizing factor and Waveform Matching misplace in refutation process is eliminated in correction
Method.
Background technology：
Full waveform inversion is a least squares problem, by the residual error for constantly reducing observation data and analogue data
It is updated to initial velocity model is previously set, the final purpose for realizing high-precision inverting.With many existing inverting sides
Method is compared, and in theory, full waveform inversion method is to obtain one of the algorithm of subsurface velocity model most accurate method.However, i.e.
Make full waveform inversion that there is high-resolution huge advantage, but is still restricted by many unfavorable factors, such as missing low frequency component,
Strong nonlinearity object function and inaccurate initial model.For these reasons, the object function of full waveform inversion is usually restrained
To local minimum point, and there is serious cycle skipping problem in real data test.In order to overcome these obstacles, many to grind
Study carefully personnel to be attempted from the following aspects：Waveform Matching, inverting and Phase-retrieval when walking.
In order to expand analogue data and observe data between waveform can matching range, have conducted in-depth research,
And propose many new methods.Bunks et al. (1995.Multiscale seismic waveform
inversion.Geophysics,60(5):1457-1473) by the method for filtering by the high fdrequency component and low frequency of seismic data
Component separates, and realizing the first fast Construction of A Model of inverting macroscopic view, inverting details rate pattern constructs again.However, low in seismic data
Frequency component is very faint, and is always subjected to pollution (the Crosshole seismic waveform tomography-of noise
I.Strategy for real data application.Geophysical Journal of the Royal
Astronomical Society,2006,167(1):1224-1236).Importantly, ultralow frequency focus cost is high, and
And it is unavailable usually in seismic prospecting.Therefore, when not having low frequency component in seismic data, multi-scale strategy cannot be well
The macrostructure of resume speed model.Earthquake enveloping method (2014, Seismic envelope inversion and
modulation signal model:Geophysics 79WA13-WA24.) by using non-linear Hilbert transform energy
It is enough that low frequency component is extracted from original earthquake data.Hu et al. (Adaptive multi-step Full Waveform
Inversion based on Waveform Mode Decomposition.Journal of Applied Geophysics,
2017,139:A kind of method of waveform mode decomposition 195-210.) is proposed to demodulate ultralow frequency component, and is deriving companion
With in focal process do not use chain rule.It in this way can be so that ultralow frequency component be perfectly preserved at in source.Also very
Mostly based on the innovative approach of object function come spreading wave form can matched range, such as accumulated energies norm (Accumulated
energy norm for full waveform inversion of marine data.Journal of Applied
Geophysics, 2017,147.), optimal transport target function (Optimal Transport for Seismic Full
Waveform Inversion [J] .2016,14 (8)), second-order time integrates wave field object function (Full waveform
inversion of the second-order time integral wavefield.Chinese Journal of
Geophysics, 2016,59 (6)) and normalization integral method (The normalized integration method:an
alternative to full waveform inversion.2012.)。
Seimic travel time information and subsurface velocities value (Computerized geophysical in a linear relationship
tomography.Proceedings of the IEEE,1979,67(7):1065-1073.).Therefore, Zoumaling tunnel is not
There are problems that cycle skipping, such logical method is often to assume come inversion speed model based on ray tracing theory and high frequency.So
And shown according to the achievement in research of forefathers, information can only reflect background velocity model when walking.When the construction of rate pattern is complicated
When, the model resolution of Zoumaling tunnel method is limited.Luo et al., (Wave-equation traveltime
inversion.Geophysics,1991,56(5):645-653.) and Zhou et al. (Acoustic wave-equation
traveltime and waveform inversion of crosshole seismic data.Geophysics,1995,
60(3):Inversion method 765-773.) is proposed when wave equation is walked come the letter between inverting when connecting full waveform inversion and walking
Ditch.However, by this method for complicated rate pattern when, such as Marmousi rate patterns, when wave equation is walked inverting obtain
Rate pattern cannot still meet the inverting of conventional full wave shape to initial model requirement.Later, Luo et al. (Full-traveltime
inversion[J].Geophysics,2016,81(5):R261-R274. it is multiple to solve) to have also been proposed inversion method when walking entirely
The inversion problem of parasitic mode type, and do not need the seismic data of accurate initial velocity model or ultralow frequency.Recklessly et al.
(Single-frequency waveform-based multiscale wave-equation traveltime
inversion[C]//Seg Technical Program Expanded.2017:1692-1696) by time domain single-frequency waveform with
Wave equation Travel Time Inversion is combined, and carries out inverting step by step from low frequency component to high fdrequency component, may finally be obtained preferably
Result.
Phase property in seismic data and information when walking are closely related, have compared with amplitude information and Construction of A Model
Better linear relationship (Theoretical background for continental-and global-scale full-
waveform inversion in the time–frequency domain.Geophysical Journal of the
Royal Astronomical Society,2008,175(2):665-685.).Therefore can be believed by detaching phase and amplitude
It ceases to design more linear inversion methods.Many object functions about pure phase information have heretofore been proposed.Such as：
Fichtner et al. proposes a kind of time-frequency domain Phase-retrieval method；Et al. (Misfit functions for
full waveform inversion based on instantaneous phase and envelope
measurements.Geophysical Journal International,2011,185(2):845-870.) propose base
In the instantaneous phase inversion method of Hilbert transform；Choi et al. (Frequency-domain waveform inversion
using the phase derivative.Geophysical Journal International,2013,195(3):
Phase information, inversion result 1904-1916.) has been used to show that phase information is ground in frequency domain or Laplce's Fourier domain
Shake the preferable selection of inverting.Normalized crosscorrelation object function (Application of multi-source simultaneously
waveform inversion to marine streamer data using the global correlation
norm.Geophysical Prospecting,2012,60(4):748-758.) similar with Phase-retrieval method.It, which is calculated, sees
Similitude between measured data and analogue data, compared with conventional full waveform inversion object function, mutually related objects function tool
Have with the better linear relationship of subsurface velocity model, reduce dependence of the existing inversion method to initial velocity model.
Invention content：
The purpose of the present invention is in view of the above shortcomings of the prior art, provide during a kind of alleviation full waveform inversion because scarce
It loses low frequency component and has the local cross-correlation time-frequency domain Phase-retrieval based on weighting of stronger dependence to initial velocity model
Method.
Idea of the invention is that：The phase information of time-frequency domain seismic data is applied in mutually related objects function, to delay
Dependence of the full waveform inversion method to initial velocity model is solved, and then alleviates the cycle skipping occurred during full waveform inversion and asks
Topic, finally obtains high-precision full waveform inversion result.Entire calculating process uses minute gun paralleling tactic, it is therefore an objective to accelerate whole speed
Spend the process of model inversion.
The purpose of the present invention is what is be achieved through the following technical solutions：
Local cross-correlation time-frequency domain Phase-retrieval method based on weighting, includes the following steps：
A, seismic data is imported into computer, it is pre-processed using seismic data processing software, and prepares
The local cross-correlation time-frequency domain Phase-retrieval being weighted；The present invention utilizes the ground that forward simulation generates on true velocity model
Shake data replace the seismic data of practical field acquisition to carry out numerical value test；
B, linear increment initial velocity model is established, initial value of the initial velocity model as inverting；
C, observation system, source wavelet and linear increment initial velocity model are input to the local cross-correlation based on weighting
In time-frequency domain Phase-retrieval algorithm, and Wave equation forward modeling is carried out, obtain analogue data, and stores main story wave field；
D, low-pass filtering is carried out at the same time to observation data and analogue data, filters out seismic data high frequency components.
E, Gabor transformation is carried out by road to observation data and analogue data and obtains the observation data and simulation number of time-frequency domain
According to；
F, the object function of the local cross-correlation time-frequency domain Phase-retrieval method based on weighting is built：Wherein indicate WuIt indicates weight factor, and setsIndicate observation data,Indicate analogue data time-frequency domain phase,Indicate Gaussian window
The central instant of function,It indicates local angular frequency, and meets
G, the partial derivative of object function：And it utilizes
The partial derivative of object function is transmitted to time-frequency domain analogue data in the partial derivative of speed parameter by chain rule；
H, the pressure gradient expression formula of time-frequency domain is transformed into time and space domain using Gabor inverse transformations：The gradient of so much big gun data can be expressed as：
I, the local cross-correlation time-frequency domain Phase-retrieval method based on weighting is extracted according to the pressure gradient expression formula of object function
Adjoint focus：And continued through finite difference numerical simulation
Method travels to rate pattern space, obtains anti-pass wave field；
J, the update gradient of rate pattern is obtained using zero-lag cross-correlation to main story wave field and anti-pass wave field；
K, L-BFGS optimization algorithms computation model more new direction is utilized：vk+1=vk-αkHkgk, wherein vkFor rate pattern
Parameter, αkFor the step-length that Wolfe linear search obtains, HkFor approximate Hessian inverses of a matrix matrix, gkFor model modification gradient.
And suitable step-length is found by Wolfe convergence criterions；
L, judge whether to meet end condition, the local cross-correlation time-frequency domain Phase-retrieval based on weighting is exported if meeting
Methods and results.If being unsatisfactory for end condition, current inversion result is continued as the initial velocity model recycled next time；
M, the inversion result of low-frequency range is further subjected to the part based on weighting mutually using the seismic data of high band
Time-frequency domain Phase-retrieval is closed, and exports final inversion result.
Advantageous effect：The phase information of seismic data have with the preferable linear relationship of subsurface velocity model, utilize phase
Information builds object function, largely reduces dependence of the full waveform inversion to initial velocity model.When therefore utilizing
Frequency domain phase mutually related objects function, even if in the case where initial velocity model is poor or seismic data lacks low frequency, most
It can also obtain high-precision inversion result eventually.
But the phase change of seismic data is unrelated with amplitude variations, phase information cannot quantify strong reflection stratum well
With the difference on weak reflection stratum.If not being distinguish with weak reflection to strong reflection, this will likely make strong reflection interface
Corresponding gradient gradient numerical value having the same corresponding with weak reflecting interface, rate pattern renewal amount is asymmetric, causes entire
Inverting flow process may be updated to the direction of mistake.Therefore the present invention introduces weight factor in the object function of Phase-retrieval,
To alleviate above-mentioned problem.
The introducing of weight factor also solves interference of the noise to earthquake data phase information well simultaneously.Because local
It shakes in data there are when faint noise jamming, these faint noises have on the phase diagram of seismic wave and effective information
Equal weight, this may be such that the phase of analogue data constantly to the phase matched of faint noise, cause entire inverting to
The direction of mistake updates, or even cannot get any information of subsurface model.Due to not containing noise contribution in analogue data, by mould
Quasi- data are added to as weight factor in time-frequency domain cross-correlation phase object function can suppress noise to Phase-retrieval well
Influence.
Description of the drawings
Local cross-correlation time-frequency domain Phase-retrieval method flow diagrams of the Fig. 1 based on weighting.
The adjoint hypocentral plot that Fig. 2 distinct methods are calculated
(a) the adjoint hypocentral plot of conventional full wave shape inverting,
(b) the adjoint hypocentral plot of normalized crosscorrelation full waveform inversion,
(c) the adjoint hypocentral plot of the local cross-correlation time-frequency domain Phase-retrieval without weight factor,
(d) the adjoint hypocentral plot of the local cross-correlation time-frequency domain Phase-retrieval weighted.
Fig. 3 rate pattern figures
(a) true velocity model figure,
(b) initial velocity model figure.
Fig. 4 conventional full wave shape inversion results figure (missing low-frequency test)
(a) the inversion result figure of 6-15Hz waveform components,
(b) using (a) as initial velocity model, the inversion result figure of 6-40Hz waveform components.
Fig. 5 normalized crosscorrelation full waveform inversions result figure (missing low-frequency test)
(a) the inversion result figure of 6-15Hz waveform components,
(b) using (a) as initial velocity model, the inversion result figure of 6-40Hz waveform components.
The local cross-correlation time-frequency domain Phase-retrieval result figure (missing low-frequency test) of Fig. 6 weightings
(a) the inversion result figure of 6-15Hz waveform components,
(b) using (a) as initial velocity model, the inversion result figure of 6-40Hz waveform components.
The single track sectional view (being derived from Fig. 4,5,6 final inversion results) of Fig. 7 distinct methods inversion results
(a) conventional full wave shape inversion result figure,
(b) normalized crosscorrelation full waveform inversion result figure,
(c) the local cross-correlation time-frequency domain Phase-retrieval result figure weighted.
The model error curves figure that Fig. 8 distinct methods inversion result changes with iterations.
The adjoint hypocentral plot (strong Gauusian noise jammer test) that Fig. 9 distinct methods are calculated
(a) the adjoint hypocentral plot of conventional full wave shape inverting,
(b) the adjoint hypocentral plot of normalized crosscorrelation full waveform inversion,
(c) the adjoint hypocentral plot of the local cross-correlation time-frequency domain Phase-retrieval without weight factor,
(d) the adjoint hypocentral plot of the local cross-correlation time-frequency domain Phase-retrieval weighted.
Figure 10 distinct methods inversion results figure (strong Gauusian noise jammer test)
(a) conventional full wave shape inversion result figure,
(b) normalized crosscorrelation full waveform inversion result figure,
(c) the local cross-correlation time-frequency domain Phase-retrieval result figure weighted.
Specific implementation mode
Below in conjunction with the accompanying drawings with example to further detailed description of the invention.
(1) obtains time-frequency domain seismic data by Gabor transformation：Its basic thought is by seismic signal plus a sliding
Gauss function.When using Gauss function come when intercepting unstable signal, the signal at this in window can regard stable state as
, then the signal after adding window is subjected to Fourier transformation, and constantly move window function along the time axis, it carries out in Fu successively
Leaf transformation, the process are the time frequency analysis to time-domain seismic signal.With the progress of analysis, the frequency spectrum that signal can be obtained exists
Changing rule on time shaft is called the time-frequency spectrum of the signal, by observing the feature of time-frequency spectrum, more can intuitively reflect letter
The rule that number frequency changes over time.
(2) builds cross-correlation phase object function in time-frequency domain：Seismic data Phase-retrieval actually constantly utilizes
The phase information of analogue data goes the phase information of matching observation data, the present invention to believe using the time-frequency domain phase of seismic data
Breath, and a time-frequency domain object function is built in conjunction with mutually related objects function.Cross-correlation phase object function can be better
The phase portrait degree for weighing two data proves rate pattern closer to true velocity mould when phase portrait degree increases
Type proves that inversion result updates iteration to the direction of mistake when phase portrait degree reduces.By judging that object function increases
Or reduce and can control Local Optimization Algorithm to being correctly oriented update iteration so that initial velocity model is gradually to true
Rate pattern is close, the final result for realizing high-precision inverting.
(3) weight constraints are added to time-frequency domain phase object function in：The phase change of seismic data is unrelated with amplitude variations,
This point makes the phase information in seismic data at strong reflection interface and weak reflecting interface value having the same, this will likely make
Obtain the corresponding gradient in strong reflection interface gradient numerical value having the same corresponding with weak reflecting interface so that rate pattern renewal amount
Asymmetry causes entire inverting flow process that may be updated to the direction of mistake.Simultaneously when there are faint noises in seismic data
When interference, these faint noises have the weight equal with effective information on the phase diagram of seismic data, this may
The phase of analogue data can be made constantly to the phase matched of faint noise, cause entire inverting to be carried out to the direction of mistake, very
To any effective information that cannot get subsurface velocity model at all.Therefore it is needed wherein while designed phase object function
In addition a weight term, to avoid the generation of the above problem.So that object function is more stablized during inverting.
(4) utilizes the gradient of chain rule derived object function：In order to seek partial derivative of the object function to speed parameter,
The partial derivative for obtaining analogue data to speed parameter using chain rule in time-frequency domain is needed, could be combined in this way with state
Method obtains the update gradient of subsurface velocity model.
(5) combinations Gabor inverse transformations obtain time-domain with focus：Time-frequency domain analogue data has been obtained using chain rule
About the partial derivative of speed parameter, in order to be further combined with time finite element method forward simulation, need anti-using Gabor
The seismic data for becoming time-frequency domain of changing commanders is transformed into time-domain.The adjoint shake similar with conventional full wave shape inverting can be obtained in this way
Source, therefore this method has implementing procedure similar with conventional full wave shape inverting.
(6) is utilized obtains anti-pass wave field with focus anti-pass, and carries out zero-lag cross-correlation with main story wave field and obtain target
The update gradient of function：The method that most original seeks gradient is solved using perturbation theory, by each to subsurface model
After point increases small disturbance, then carries out forward simulation and obtain corresponding analogue data.The perturbation theory does not utilize chain method
First derivative then is approximate, but its calculation amount is huge.By this paper by taking test model (69 × 192) as an example, each iteration needs just
Simulation 26496 times, forward simulation about 15 seconds are drilled, then a model modification needs 110.4 hours.Certainly this does not also include
Optimization algorithm step-length chooses improper the case where needing to recalculate.After Adjoint State Method, iteration once need 4 times just
Drill the time of simulation about 1 minute.So seeking the gradient of object function by Adjoint State Method, full waveform inversion is substantially reduced
Gradient calculates the required time.
(7) obtains the more new direction of rate pattern using L-BFGS optimization algorithms：L-BFGS optimization algorithms effectively overcome
Explicitly Hessian matrixes and its inverse difficulty are sought, with super-linear convergence, computational efficiency is high, and committed memory is small,
The advantages that precision is high is relatively suitble to solve large-scale nonlinear optimization problem.
Local cross-correlation time-frequency domain Phase-retrieval method based on weighting, includes the following steps：
A, MATLAB2013a copyrighted softwares are installed, and the seismic data processing software write based on MATLAB language is installed
Crews kits.Seismic data is imported into computer simultaneously, it is pre-processed using seismic data processing software,
And it is ready for the local cross-correlation time-frequency domain Phase-retrieval of weighting；The present invention utilizes the forward simulation on true velocity model to produce
Raw seismic data replaces the seismic data of practical field acquisition to carry out numerical value test；
B, linear increment initial velocity model is established, initial value of the initial velocity model as inverting, such as Fig. 3 (b) institutes
Show.Linear increment initial velocity model is substantially free of any tectonic information of true velocity model (Fig. 3 a), in this way can be fine
Verify the stability and robustness of the method for the present invention in ground；
C, observation system, source wavelet and linear initial velocity model are input to the local cross-correlation time-frequency based on weighting
In the Phase-retrieval algorithm of domain, and Wave equation forward modeling is carried out, obtain analogue data, and stores main story wave field.Wherein Chang Mi
The two-dimentional second order ACOUSTIC WAVE EQUATION of degree can be expressed as：
Wherein x, z representation space reference axis, v indicate that the parameter of rate pattern, p indicate the pressure field of Acoustic Wave Propagation.Pass through
Finite difference method solves above-mentioned Acoustic Wave-equation；
D, low-pass filtering is carried out at the same time to observation data and analogue data, filters out seismic data high frequency components：Due to complete
Waveform inversion method is a strong nonlinearity optimization problem, and can not usually search global minimum, leads to serious cycle
Chattering.So it is unpractical that full waveform inversion method, which is directly applied to seismic data inversion,.The present invention uses Bart
Seismic data is divided into different dimensions in frequency and carries out full waveform inversion respectively by Butterworth filter.When ground of the utilization compared with low-frequency range
Shake data carry out inverting when, local minimum is less, can utilize low-frequency range large scale seismic data obtain one compared with
Good rate pattern, then carry out the inverting of high band.Based on this theory, the earthquake of low-frequency band (6-15Hz components) is used first
Data carry out one macro-velocity model of inverting, then as the initial velocity model of high frequency band (6-40Hz components) into one
Step obtains final inversion result；
E, Gabor transformation is carried out by road to observation data and analogue data, obtains the observation data and simulation number of time-frequency domain
According to.Wherein the Fourier transformation of seismic observation data and analogue data can be expressed as：
Seismic observation data and analogue data Gabor transformation can be expressed as：
WhereinIndicate the central instant of Gauss function；Indicate that local angular frequency, d (t) indicate that time-domain observes number
According to,Indicate that frequency domain observes data,Indicating that time-frequency domain observes data, u (t) indicates time-domain analogue data,Indicate frequency domain analogue data,Indicate time-frequency domain analogue data.F [] indicates that Fourier transformation, G [] indicate
Gabor transformation；
F, the object function of the local cross-correlation time-frequency domain Phase-retrieval method based on weighting is built：
Wherein indicate WuIt indicates weight factor, and setsWithSimulation number is indicated respectively
, in the phase of time-frequency domain, and meet according to observation data
G, can and it utilize chain rule by the partial derivative of object function the partial derivative of speed parameter about object function
It is transmitted to partial derivative of the analogue data to speed parameter：
Wherein haveAndBring above formula into：
H, the pressure gradient expression formula of time-frequency domain is transformed into time-domain using Gabor inverse transformations, wherein seismic observation data
It can be expressed as with analogue data inverse Fourier transform：
Seismic observation data and analogue data Gabor inverse transformations can be expressed as：
And there are the Gabor transformations of conjugate complex number：
It substitutes the above to obtain：
Final gradient formula is obtained according to Gabor reconstructed formulas simultaneously：
I, it seeks with focus (as shown in Figure 2)：
And the method for being continued through finite difference numerical simulation travels to rate pattern space, obtains anti-pass wave field；
J, the update gradient of rate pattern is obtained using the method for zero-lag cross-correlation to main story wave field and anti-pass wave field, it is right
Answer the gradient of more big gun seismic datas that can be expressed as：
Wherein PfIndicate main story wave field,Indicate anti-pass wave field；
K, L-BFGS optimization algorithms computation model more new direction is utilized：vk+1=vk-αkHkgk, wherein vkFor rate pattern
Parameter, αkFor the step-length that Wolfe linear search obtains, HkFor approximate Hessian inverses of a matrix matrix, gkFor model modification gradient.
And suitable step-length is found by Wolfe convergence criterions；
L, judge whether to meet end condition (as shown in Figure 8), when exporting the local cross-correlation based on weighting if meeting
Frequency domain Phase-retrieval methods and results.If being unsatisfactory for end condition, current inversion result is continued first as what is recycled next time
Beginning rate pattern；
M, the inversion result of low-frequency range is further subjected to the part based on weighting mutually using the seismic data of high band
Time-frequency domain Phase-retrieval is closed, and exports final inversion result (as shown in Figure 6).
Embodiment 1
It is that full waveform inversion method is most difficult to the construction restored according to the numerical experimentation of front, on the left of Marmousi models.
Therefore, we intercept Marmousi models left part to test method proposed by the present invention, and with other two kinds of sides
Method is compared (conventional full wave shape inversion method and normalization normalized crosscorrelation full waveform inversion method).Initial velocity model
As shown in Figure 3b, it shares 36 focus to be uniformly distributed on the surface, each focus corresponds to 192 wave detectors.Seismic data is remembered
The record time is 2s, time interval 1ms.While in order to prove the local cross-correlation time-frequency domain Phase-retrieval method based on weighting
Validity utilizes butterworth high pass filter to cut off Ricker wavelet 6Hz components below in testing.In numerical value test,
We are updated initial model to calculate descent direction using L-BFGS optimization algorithms.
Model parameter is as follows：
Table 1：Local cross-correlation time-frequency domain Phase-retrieval method test parameter based on weighting
| Model size | Grid spacing | Lateral distance | Longitudinal depth | Velocity interval | Dominant frequency of seismic wavelet |
| 69*192 | 12.5m | 2.4km | 0.8625km | 1.5-4km/s | 15Hz |
Based on the above parameter, under the test environment shown in table 2, the local cross-correlation time-frequency domain phase based on weighting is carried out
Inversion method is tested.
2 computer performance test environment of table
By Fig. 4 (a) and Fig. 5 (a) comparisons as can be seen that normalized crosscorrelation full waveform inversion result is better than conventional all-wave
The result of shape inverting.There is one group of exceptional value on the top of Marmousi models in Fig. 4 (a), this is by conventional full wave shape inverting
Method will appear apparent cycle skipping problem when lacking low frequency component.Due to normalized crosscorrelation full waveform inversion side
Method has used mutually related objects function and normalization seismic data to carry out inversion speed model, mainly measures observation data and simulation
Similitude between data.Therefore, it has and the stronger linear relationship of subsurface velocity model than conventional full waveform inversion.From
As can be seen that the local cross-correlation time-frequency domain Phase-retrieval result of weighting has prodigious improvement in Fig. 6 a and Fig. 7, especially exist
Horizontal tectonics on the left of Marmousi.Since phase cross-correlation function has and the better linear relationship of subsurface velocity model, pole
Dependence of the full waveform inversion to initial velocity model is reduced in big degree.In addition, introducing weight in time-frequency domain object function
The factor not only increases the stability of refutation process, and has big advantage to alleviating cycle skipping problem.Numerical experimentation
The result shows that the local cross-correlation time-frequency domain Phase-retrieval method based on weighting can be effectively relieved to initial velocity model according to
Rely.
When not having low frequency component in earthquake record data, there is serious cycle skipping in conventional Full wave shape result
(Fig. 4), inverse model is gradually deviated from true model.As shown in figure 8, the fitting difference of conventional full wave shape inversion result is always continuous
Increase.However, as iterations increase, the models fitting of the local cross-correlation time-frequency domain Phase-retrieval method based on weighting is poor
It is being gradually reduced.Demonstrate,proved the local cross-correlation time-frequency domain Phase-retrieval method based on weighting have do not depend on initial velocity model
Advantage, and with the principal frequency component that seismic data is continuously improved, high-resolution inversion result may finally be obtained.
Conventional full wave shape inversion result is can be seen that from Figure 10 a to differ greatly with true velocity model, this is mainly due to
Lack low frequency component, and by the severe jamming of strong Gaussian Background noise.Figure 10 b are compared discovery, normalizing with Figure 10 a
It is more a little better than conventional full wave shape inversion result to change cross-correlation full waveform inversion result, is primarily due to mutually related objects function tool
There is the ability for resisting Gaussian noise.However, since signal-to-noise ratio is too low, most of significant wave information is all submerged in strong background noise
In, cause final inversion result still to deviate true velocity model.In conjunction with Fig. 9 d and Figure 10 c it is found that strong Gaussian Background noise
The local cross-correlation time-frequency domain Phase-retrieval method of weighting is almost no impact.It demonstrates, the part based on weighting is mutually
Related time-frequency domain Phase-retrieval method can decay with the strong Gaussian Background noise in source, with the week for overcoming low signal-to-noise ratio to bring
Phase chattering.
Claims (1)
1. a kind of local cross-correlation time-frequency domain Phase-retrieval method based on weighting, includes the following steps：
A, seismic data is imported into computer, it is pre-processed using seismic data processing software, and is ready for
The local cross-correlation time-frequency domain Phase-retrieval of weighting；The present invention utilizes the earthquake number that forward simulation generates on true velocity model
Numerical value test is carried out according to the seismic data of practical field acquisition is replaced；
B, linear increment initial velocity model is established, initial value of the initial velocity model as inverting；
C, observation system, source wavelet and linear increment initial velocity model are input to the local cross-correlation time-frequency based on weighting
In the Phase-retrieval algorithm of domain, and Wave equation forward modeling is carried out, obtain analogue data, and stores main story wave field；
D, low-pass filtering is carried out at the same time to observation data and analogue data, filters out seismic data high frequency components；
E, Gabor transformation is carried out by road to observation data and analogue data and obtains the observation data and analogue data of time-frequency domain；
F, the object function of the local cross-correlation time-frequency domain Phase-retrieval method based on weighting is built：Wherein indicate WuIt indicates weight factor, and sets
Indicate observation data,Indicate analogue data time-frequency domain phase,Indicate the central instant of Gauss function,Table
Show local angular frequency, and meets
G, the partial derivative of object function：And utilize chain method
The partial derivative of object function is then transmitted to time-frequency domain analogue data in the partial derivative of speed parameter；
H, the pressure gradient expression formula of time-frequency domain is transformed into time and space domain using Gabor inverse transformations：The gradient of so much big gun data can be expressed as：
I, the companion of the local cross-correlation time-frequency domain Phase-retrieval method based on weighting is extracted according to the pressure gradient expression formula of object function
With focus：And the method for being continued through finite difference numerical simulation
Rate pattern space is traveled to, anti-pass wave field is obtained；
J, the update gradient of rate pattern is obtained using zero-lag cross-correlation to main story wave field and anti-pass wave field；
K, L-BFGS optimization algorithms computation model more new direction is utilized：vk+1=vk-αkHkgk, wherein vkFor the parameter of rate pattern,
αkFor the step-length that Wolfe linear search obtains, HkFor approximate Hessian inverses of a matrix matrix, gkFor model modification gradient, and lead to
It crosses Wolfe convergence criterions and finds suitable step-length；
L, judge whether to meet end condition, the local cross-correlation time-frequency domain Phase-retrieval method based on weighting is exported if meeting
If as a result, being unsatisfactory for end condition, current inversion result is continued as the initial velocity model recycled next time；
When the inversion result of low-frequency range further m, being carried out the local cross-correlation based on weighting using the seismic data of high band
Frequency domain Phase-retrieval, and export final inversion result.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201810490671.9A CN108680957B (en) | 2018-05-21 | 2018-05-21 | Local cross-correlation time-frequency domain Phase-retrieval method based on weighting |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201810490671.9A CN108680957B (en) | 2018-05-21 | 2018-05-21 | Local cross-correlation time-frequency domain Phase-retrieval method based on weighting |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN108680957A true CN108680957A (en) | 2018-10-19 |
| CN108680957B CN108680957B (en) | 2019-08-13 |
Family
ID=63807347
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201810490671.9A Expired - Fee Related CN108680957B (en) | 2018-05-21 | 2018-05-21 | Local cross-correlation time-frequency domain Phase-retrieval method based on weighting |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN108680957B (en) |
Cited By (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN111999764A (en) * | 2020-05-20 | 2020-11-27 | 中国矿业大学 | Method for constructing least square reverse time migration under salt based on time-frequency domain objective function |
Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US9470811B2 (en) * | 2014-11-12 | 2016-10-18 | Chevron U.S.A. Inc. | Creating a high resolution velocity model using seismic tomography and impedance inversion |
| CN106908835A (en) * | 2017-03-01 | 2017-06-30 | 吉林大学 | Band limit Green's function filters multiple dimensioned full waveform inversion method |
| CN107422379A (en) * | 2017-07-27 | 2017-12-01 | 中国海洋石油总公司 | Multiple dimensioned seismic full-field shape inversion method based on local auto-adaptive convexification method |
| CN107765302A (en) * | 2017-10-20 | 2018-03-06 | 吉林大学 | Inversion method when time-domain single-frequency waveform independent of source wavelet is walked |
| CN107843925A (en) * | 2017-09-29 | 2018-03-27 | 中国石油化工股份有限公司 | A kind of reflection wave inversion method based on orrection phase place |
-
2018
- 2018-05-21 CN CN201810490671.9A patent/CN108680957B/en not_active Expired - Fee Related
Patent Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US9470811B2 (en) * | 2014-11-12 | 2016-10-18 | Chevron U.S.A. Inc. | Creating a high resolution velocity model using seismic tomography and impedance inversion |
| CN106908835A (en) * | 2017-03-01 | 2017-06-30 | 吉林大学 | Band limit Green's function filters multiple dimensioned full waveform inversion method |
| CN107422379A (en) * | 2017-07-27 | 2017-12-01 | 中国海洋石油总公司 | Multiple dimensioned seismic full-field shape inversion method based on local auto-adaptive convexification method |
| CN107843925A (en) * | 2017-09-29 | 2018-03-27 | 中国石油化工股份有限公司 | A kind of reflection wave inversion method based on orrection phase place |
| CN107765302A (en) * | 2017-10-20 | 2018-03-06 | 吉林大学 | Inversion method when time-domain single-frequency waveform independent of source wavelet is walked |
Non-Patent Citations (2)
| Title |
|---|
| 杨贺龙 等: "波场相位相关时移全波形反演", 《石油物探》 * |
| 梁煌 等: "互相关与最小二乘加权目标函数全波形反演", 《世界地质》 * |
Cited By (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN111999764A (en) * | 2020-05-20 | 2020-11-27 | 中国矿业大学 | Method for constructing least square reverse time migration under salt based on time-frequency domain objective function |
| CN111999764B (en) * | 2020-05-20 | 2021-04-13 | 中国矿业大学 | Method for constructing least square reverse time migration under salt based on time-frequency domain objective function |
Also Published As
| Publication number | Publication date |
|---|---|
| CN108680957B (en) | 2019-08-13 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Chen et al. | Elastic least-squares reverse time migration via linearized elastic full-waveform inversion with pseudo-Hessian preconditioning | |
| CN103293552B (en) | A kind of inversion method of Prestack seismic data and system | |
| CN106405651B (en) | Full waveform inversion initial velocity model construction method based on logging matching | |
| CN109669212B (en) | Seismic data processing method, stratum quality factor estimation method and device | |
| CN109407151B (en) | Time-domain full waveform inversion method based on wave field local correlation time shift | |
| CA2726462A1 (en) | Estimation of soil properties using waveforms of seismic surface waves | |
| CN101201409B (en) | Method for revising earthquake data phase | |
| CN109459789B (en) | Time-domain full waveform inversion method based on amplitude decaying and linear interpolation | |
| CN104820242B (en) | A kind of road collection amplitude towards prestack inversion divides compensation method | |
| Yin et al. | Improving horizontal resolution of high-frequency surface-wave methods using travel-time tomography | |
| Zhou et al. | Passive surface-wave waveform inversion for source-velocity joint imaging | |
| Pan et al. | Individual and joint inversions of shallow-seismic Rayleigh and love waves: full-waveform inversion versus random-objective waveform inversion | |
| CN111999764B (en) | Method for constructing least square reverse time migration under salt based on time-frequency domain objective function | |
| CN108680957B (en) | Local cross-correlation time-frequency domain Phase-retrieval method based on weighting | |
| Yu et al. | Application of weighted early-arrival waveform inversion to shallow land data | |
| Biondi et al. | Elastic-parameter estimation by combining full-waveform inversion by model extension and target-oriented elastic inversion | |
| Cai et al. | Data weighted full-waveform inversion with adaptive moment estimation for near-surface seismic refraction data | |
| Sava et al. | Interferometric imaging condition for wave-equation migration | |
| Yu et al. | An application of multiscale early arrival waveform inversion to shallow seismic data | |
| Chai et al. | Q-compensated acoustic impedance inversion of attenuated seismic data: Numerical and field-data experiments | |
| CN111208568B (en) | Time domain multi-scale full waveform inversion method and system | |
| Zhang et al. | Frequency extension and robust full-waveform inversion based on n th power operation | |
| Zhang et al. | Autoencoded elastic wave-equation traveltime inversion: Toward reliable near-surface tomogram | |
| Kurzmann et al. | Real data applications of seismic full waveform inversion | |
| Huang et al. | Data-assimilated time-lapse visco-acoustic full-waveform inversion: Theory and application for injected CO2 plume monitoring |
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 |
Granted publication date: 20190813 Termination date: 20200521 |
|
| CF01 | Termination of patent right due to non-payment of annual fee |