CN106054244B - The LPF of window multiple dimensioned full waveform inversion method when blocking - Google Patents
The LPF of window multiple dimensioned full waveform inversion method when blocking Download PDFInfo
- Publication number
- CN106054244B CN106054244B CN201610429648.XA CN201610429648A CN106054244B CN 106054244 B CN106054244 B CN 106054244B CN 201610429648 A CN201610429648 A CN 201610429648A CN 106054244 B CN106054244 B CN 106054244B
- Authority
- CN
- China
- Prior art keywords
- window
- full waveform
- inversion
- lpf
- waveform inversion
- 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.)
- Expired - Fee Related
Links
- 238000000034 method Methods 0.000 title claims abstract description 64
- 230000000903 blocking effect Effects 0.000 title claims abstract description 44
- 230000006870 function Effects 0.000 claims description 31
- 230000008569 process Effects 0.000 claims description 14
- 238000012545 processing Methods 0.000 claims description 14
- 238000009795 derivation Methods 0.000 claims description 6
- 238000012986 modification Methods 0.000 claims description 6
- 230000004048 modification Effects 0.000 claims description 6
- 238000001914 filtration Methods 0.000 claims description 5
- 238000001615 p wave Methods 0.000 claims description 5
- 238000004422 calculation algorithm Methods 0.000 claims description 3
- 238000005070 sampling Methods 0.000 claims description 3
- 238000004088 simulation Methods 0.000 claims 4
- 238000012360 testing method Methods 0.000 abstract description 14
- 230000000694 effects Effects 0.000 description 5
- 238000005516 engineering process Methods 0.000 description 4
- 230000008901 benefit Effects 0.000 description 3
- 238000009499 grossing Methods 0.000 description 3
- 230000009286 beneficial effect Effects 0.000 description 2
- 238000004364 calculation method Methods 0.000 description 2
- 238000010276 construction Methods 0.000 description 2
- 238000011161 development Methods 0.000 description 2
- 235000013399 edible fruits Nutrition 0.000 description 2
- 239000002360 explosive Substances 0.000 description 2
- 230000006872 improvement Effects 0.000 description 2
- HPTJABJPZMULFH-UHFFFAOYSA-N 12-[(Cyclohexylcarbamoyl)amino]dodecanoic acid Chemical compound OC(=O)CCCCCCCCCCCNC(=O)NC1CCCCC1 HPTJABJPZMULFH-UHFFFAOYSA-N 0.000 description 1
- 241001269238 Data Species 0.000 description 1
- 230000001133 acceleration Effects 0.000 description 1
- 238000002939 conjugate gradient method Methods 0.000 description 1
- 238000012937 correction Methods 0.000 description 1
- 230000007423 decrease Effects 0.000 description 1
- 230000003247 decreasing effect Effects 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 238000010586 diagram Methods 0.000 description 1
- 238000011156 evaluation Methods 0.000 description 1
- 230000005284 excitation Effects 0.000 description 1
- 238000003384 imaging method Methods 0.000 description 1
- 230000010354 integration Effects 0.000 description 1
- 238000010606 normalization Methods 0.000 description 1
- 239000003208 petroleum Substances 0.000 description 1
- 238000003672 processing method Methods 0.000 description 1
- 238000000926 separation method Methods 0.000 description 1
- 230000006641 stabilisation Effects 0.000 description 1
- 238000011105 stabilization Methods 0.000 description 1
- 238000013517 stratification Methods 0.000 description 1
- 230000007704 transition Effects 0.000 description 1
- 238000012795 verification Methods 0.000 description 1
- 239000002699 waste material Substances 0.000 description 1
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
- G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
Landscapes
- Engineering & Computer Science (AREA)
- Remote Sensing (AREA)
- Physics & Mathematics (AREA)
- Life Sciences & Earth Sciences (AREA)
- Acoustics & Sound (AREA)
- Environmental & Geological Engineering (AREA)
- Geology (AREA)
- General Life Sciences & Earth Sciences (AREA)
- General Physics & Mathematics (AREA)
- Geophysics (AREA)
- Geophysics And Detection Of Objects (AREA)
Abstract
The present invention relates to a kind of multiple dimensioned full waveform inversion method of LPF of window during block.This method carries out full waveform inversion first in low-frequency range, with window refutation strategy when blocking, and obtained inversion result is used for Mid Frequency inverting as initial model.Then using the inversion result of Mid Frequency as initial model, for Conventional Time domain full waveform inversion and final inversion result is obtained.Window refutation strategy has saved a large amount of forward modeling times when blocking.LPF Multi-scale inversion strategy, geological data can be divided into different frequency range and carry out inverting.Window is applied in combination with LPF when will block, and can obtain preferable inversion result in the case of missing low frequency and strong Gauusian noise jammer.Model test results show, the combined strategy alleviates cycle slip phenomenon while computational efficiency is improved, and when blocking window the multiple dimensioned full waveform inversion method of LPF in inversion accuracy, computational efficiency, low frequency rely on, anti-noise ability etc. is better than Conventional Time domain full waveform inversion.
Description
Technical field
The present invention relates to a kind of data processing method of seismic prospecting, and window and LPF are multiple dimensioned anti-when especially blocking
Drill strategy.Whole calculating process carries out full waveform inversion using super random epicentre coding strategy, it is therefore an objective to accelerates whole Full wave shape
The process of inverting.
Background technology：
With the development of petroleum industry and deepening continuously for exploration and development, lithology exploration is gradually moved towards from the construction exploration phase
Stage, its difficulties in exploration are gradually increased.In order to comply with this requirement, seismic data processing technology constantly improve, in the meantime entirely
Waveform inversion develops rapidly, and as the study hotspot of current geophysics circle.In the 1980s, Tarantola is carried
Time-domain full waveform inversion reason method is gone out, and object function is passed through into residual error anti-pass wave field and main story ripple to model parameter derivation
Field computing cross-correlation is obtained, and Jacobi matrixes are asked for so as to avoid.In the 1990s, Pratt promotes full waveform inversion
Frequency domain is arrived, proposition only needs several discrete frequencies to can be obtained by high-precision inversion result, and low frequency is to high frequency
Refutation strategy can solve the problems, such as to be absorbed in local minimum.Shin et al. develops on the basis of frequency domain full waveform inversion
Laplace-Fourier domains full waveform inversion, this method can obtain a high accuracy in the case where lacking low frequency
Initial model.Wang et al. (2011) proposes the acceleration full waveform inversion based on CUDA programmings on GPU, and this method effectively adds
The fast computational efficiency of full waveform inversion.
Bunks (1995) proposes the multiple dimensioned full waveform inversion based on LPF, and this method is pointed out first to observing data
LPF is carried out, then the source wavelet of estimation is handled to same filter, can solve Full wave shape to a certain extent
The cycle slip problem of inverting.But carry out full waveform inversion using the multi-scale method and still there are problems that, such as：Stratum is
One viscoelastic medium, and Wave equation forward modeling is a complicated process, the meeting during underground propagation when seismic wave
There are frequency shift effect or other complexity changes.Forward modeling is carried out using the source wavelet after LPF, it may appear that forward modeling is remembered
Record and the unmatched problem of observational record, good solution is not provided on this problem Bunks et al..Chen and
Wu et al. handles observation data and forward modeling data using time attentuating filter simultaneously, i.e., be multiplied by geological data one with
The exponential function that time increases and decayed, so first drills shallow-layer information, ignores infrastructure.As inverting be gradually reduced declining
The influence of subtraction function so that inverting is by shallow-layer to deep layer Layer by layer inversion.This method is greatly reduced the inverted parameters of each step,
Solves the problems, such as the cycle slip of full waveform inversion to a certain extent.But this method during forward modeling when wasting substantial amounts of
Between, because the forward record that forward modeling obtains decays to zero substantially multiplied by with an attenuation function, the deep information.This method is without profit
The deep information is used, but wasting the plenty of time carrys out forward modeling deep wave field, largely reduces the calculating effect of full waveform inversion
Rate.
The content of the invention：
The purpose of the present invention is to be directed to data mismatch problem after above-mentioned existing LPF, there is provided window when one kind is blocked
The multiple dimensioned full waveform inversion method of LPF.
The present invention solution be：Observational record and forward modeling are handled using Butterworth most smoothing low-pass filters simultaneously
Record, so effectively prevent the inversion result that factor data mismatches and produces mistake, and this method is more suitable for actual seismic data
Handling process, accurate inversion result can be obtained.For time-domain full waveform inversion computational efficiency it is low the problem of, the present invention carries
Go out window full waveform inversion strategy when blocking, i.e., when inverting is just started, it is anti-only to intercept the progress of real data forward part time
Drill, so window data carry out Data Matching when also only needing corresponding to forward modeling, wave field when so being blocked by forward modeling in window come
Instead of forward modeling All Time wave field, the time required for forward modeling highly shortened.With the progress of inverting, gradually increase inverting
Time window length, window refutation strategy can reduce the number of each step inverted parameters when this is blocked, while add the stabilization of program
The reliability of property and inversion result.
The purpose of the present invention is achieved through the following technical solutions：
The core of the multiple dimensioned full waveform inversion method of the LPF of window is complete on MATLAB 2013a platforms when blocking
Into.LPF is carried out to observation data and forward modeling data using Butterworth filter, it is most flat by setting Butterworth
The different cut frequency of filter slide, to realize time-domain Multi-scale inversion strategy.Low pass filter can effectively filter out earthquake
Radio-frequency component in data, retain the low-frequency component in geological data, because low-frequency component contains subsurface structure macroscopic information, have
Beneficial to the foundation of initial model, while avoid the generation of cycle slip phenomenon.Window refutation strategy when blocking, that is, select different length when
Between window be used for full waveform inversion to extract observation data, while according to time span progress forward modeling corresponding to window when this, so
All time datas of each step all forward modelings can be avoided, have largely saved the calculating time.
Window refutation strategy and LPF Multi-scale inversion strategy all have compacting cycle slip and improve inverting in itself when blocking
The effect of precision.In order to obtain more stable inverting energy flow, more accurate inversion result and more efficient inversion method, this hair
Bright is that further the two is used in combination with (as shown in Figure 1).
The multiple dimensioned full waveform inversion method of the LPF of window, comprises the following steps when blocking：
MATLAB2013a softwares are installed, and the seismic data processing software Crews instruments that MATLAB language is write are installed a,
Bag；
B, the geological data of collection is pre-processed, and observational record and observation system input computer is subjected to the time
Domain full waveform inversion；
C, source wavelet is estimated from observational record, the source wavelet obtained using estimation carries out forward modeling, and (present invention utilizes
Ricker wavelet replaces the source wavelet in real data)；
D, linear increment initial model is established, initial value of the rate pattern as inverting, passes through the renewal of computation model
Direction, realize from initial model gradually to the close process of true model；
E, object function is constructed according to the principle of least square：Wherein v represents underground P-wave
Degree, dobsFor the observational record of collection, dcalIn the forward record that rate pattern is obtained by forward modeling.Seeking the gradient of object function
During to object function both ends on p wave interval velocity v derivations, obtain pressure gradient expression formula：
Wherein pbRepresent by with focus fsThe wave field obtained to underground propagation, pfForward modeling wave field is represented, s represents focus
Number, r represent the number of wave detector.
F, the time span of window when blocking is chosen, intercepts observational record forward part seismic signal, and records window length when blocking
Degree, time span of the time span of window as forward modeling when this is blocked；
G, according to the dynamic law of ACOUSTIC WAVE EQUATION, the source wavelet that is obtained using the observation system of input, estimation and just
Beginning rate pattern carries out seismic wave field forward modeling.The geological data when step only needs the forward modeling to block in window, it is not necessary to will be all
Time geological data whole forward modeling comes out.Store main story wavefield data and forward record.
The LPF multiple dimensioned full waveform inversion method relevant parameter of window when setting is blocked as requested, including model are big
Small nz × nx, grid is away from dx, dz, forward modeling time window length It, time sampling interval dt, Ricker wavelet dominant frequency f0, Butterworth filter
The cut frequency f of ripple devicecut, the maximum iteration itermax of each super random epicentre, update gradient minimum value gtol, target
Function precision tol, the maximum vmax and minimum value vmin of inversion result；
Observational record and forward record when h, to blocking in window while carry out LPF, it is ensured that with identical low pass filtered
Ripple device processing data.After filtering, obtain handling observational record and handle forward record；
I, processing observational record and processing forward record are made the difference, obtains residual error data；
J, using residual error data as with focus, acted on anti-pass operator with focus, obtain the residual of the model space
Poor anti-pass wave field；
K, main story wave field does cross-correlation with residual error anti-pass wave field, obtains model modification gradient；
L, with super-memory gradient method optimized algorithm computation model more new direction, and found by Wolfe convergence criterions it is suitable
Step-length；
M, judge it is enough to meet end condition, the multiple dimensioned full waveform inversion of LPF of window when output is blocked if meeting
Inversion method result.If being unsatisfactory for end condition, the initial model using inversion result as next circulation, Step d is returned to；
N, the initial velocity model using the output result of previous step as Conventional Time domain full waveform inversion, then carry out normal
Advise time-domain full waveform inversion.Final inversion result is exported after iteration 200 times.
Beneficial effect：Pass through the improvement to original theory method and the integration of relative earthquake full waveform inversion technology, this hair
The bright technology such as window and LPF when will successfully block has been applied in full waveform inversion.
The multiple dimensioned full waveform inversion method of LPF of window, this method when blocking proposed can not only be alleviated because of low frequency
Lack the cycle slip problem brought, and the time for having saved forward modeling of high degree.Significantly carried using window refutation strategy when blocking
The high computational efficiency and inversion accuracy of full waveform inversion.Observation data and forward modeling data are carried out using Butterworth filter
LPF, radio-frequency component is filtered out, retain low-frequency component, Butterworth most smoothing filter is different to be blocked frequently by setting
Rate, to realize the time-domain Multi-scale inversion strategy of frequency from low to high, the cycle slip phenomenon of full waveform inversion is effectively suppressed.
The refutation strategy improves the precision of full waveform inversion to a certain extent, and is more suitable for actual seismic data processing, while this hair
The LPF refutation strategy of bright proposition without being filtered to seismic wavelet, effectively avoid factor data mismatch and caused by it is anti-
Drill mistake.Focus all uses super random epicentre coding strategy in whole calculating process, it is therefore an objective to is accelerating the same of inverting speed
When suppress crosstalk noise inside super big gun.Compared with the full waveform inversion of Conventional Time domain, the mistake of inverting subsurface velocities of the present invention
Window and LPF refutation strategy when blocking have been used in journey, has solved problems with：
Window and LPF strategy combination, are applied in full waveform inversion when the 1, will block first, can quickly obtain
One high-precision initial velocity model, the initial velocity model have recovered the macroscopic information of Marmousi models substantially, profit
Conventional full wave shape inverting is carried out with the initial velocity model, cycle slip problem can be alleviated.
2nd, time histories sample will be blocked to be incorporated into full waveform inversion so that each forward modeling need not calculate the ripple at all moment
, wave field when blocking in window is only calculated, highly shortened the calculating time required for forward modeling.
Window LPF Multi-scale inversion strategy can effectively reduce the unknown number number of inverting each time when the 3rd, blocking, according to
The gradient and model modification direction that cross-correlation is tried to achieve are more accurate.Due to during model modification iteration, it is necessary to pass through
Wolfe linear search finds suitable iteration step length, in the case of gradient direction or descent direction inaccuracy, Wolfe lines
Property search can waste the substantial amounts of time, and window LPF refutation strategy proposed by the present invention when blocking accurate can obtain mould
The gradient direction of type renewal, the substantial amounts of calculating time is saved so during linear search.
4th, the inverting of conventional full wave shape needs to store main story wave field and residual error anti-pass wave field, and storing the two wave fields needs largely
Calculator memory.Window refutation strategy proposed by the present invention when blocking, it is only necessary to which the main story wave field and residual error when storing in window are anti-
Wave field is passed, has saved the memory headroom of required computer to a certain extent, especially in real data processing.
5th, super random epicentre coding strategy can effectively suppress the crosstalk noise inside super big gun, and this method is ensureing inverting essence
Reduce the number of forward modeling while spending to a certain extent, be effectively shortened the calculating time of full waveform inversion.
Model test results show, in the case of missing low-frequency data and very strong Gauusian noise jammer, utilize
The full waveform inversion of window LPF strategy can also obtain preferable inversion result when blocking.The refutation strategy can be realized
By the Layer by layer inversion method of shallow-layer to deep layer, and can is significantly improved, alleviates simultaneously low-frequency range in terms of computational efficiency
The cycle slip phenomenon caused by lacking low-frequency data.It can be seen that window LPF refutation strategy proposed by the present invention when blocking can
A high-precision initial velocity model is provided for conventional full wave shape inverting.
Brief description of the drawings
The LPF of window multiple dimensioned full waveform inversion method flow chart when Fig. 1 is blocked.
Fig. 2 Ricker wavelets
(left side) Ricker wavelet oscillogram, (right side) Ricker wavelet spectrogram.
Fig. 3 Butterworths most smoothing filter figure.
Fig. 4 rate patterns
(left side) linear increment initial velocity model figure, (right side) true velocity model figure.
Fig. 5 filter shapes figure and spectrogram
(a) 15Hz LPFs oscillogram and true seismic wave figure,
(b) 25Hz LPFs oscillogram and true seismic wave figure,
(c) 15Hz LPFs, 25Hz LPFs and true earthquake record spectrogram,
(d) 10Hz high-pass filterings oscillogram,
(e) 10Hz high-pass filterings spectrogram.
The super random epicentre observational records of Fig. 6
(left side) is without data of making an uproar, (right side) noisy data.
Window schematic diagram when Fig. 7 is blocked.
Fig. 8 observational record LPFs
(left side) 15Hz LPF observational records,
(right side) 25Hz LPF observational records.
Fig. 9 Conventional Times domain full waveform inversion result figure.
The multiple dimensioned full waveform inversion result figure of Figure 10 LPFs
(a) 15Hz LPFs full waveform inversion result figure,
(b) 25Hz LPFs full waveform inversion result figure,
(c) the final inversion result figure of LPF.
Window full waveform inversion result figure when Figure 11 is blocked.
Window LPF multiple dimensioned full waveform inversion result figure when Figure 12 is blocked.
Figure 13 Conventional Times domain full waveform inversion result figure (missing low-frequency test).
The multiple dimensioned full waveform inversion result figure of Figure 14 LPFs (missing low-frequency test).
(a) 15Hz LPFs full waveform inversion result figure,
(b) 25Hz LPFs full waveform inversion result figure,
(c) the final inversion result figure of LPF.
Window full waveform inversion result figure (missing low-frequency test) when Figure 15 is blocked.
The multiple dimensioned full waveform inversion result figure of window LPF (missing low-frequency test) when Figure 16 is blocked.
(a) window+15Hz LPFs full waveform inversion result figure when blocking,
(b) 25Hz LPFs full waveform inversion result figure,
(c) window LPF final inversion result figure when blocking.
Figure 17 Conventional Times domain full waveform inversion result figure (is tested) containing Gaussian noise.
Window full waveform inversion result figure (being tested containing Gaussian noise) when Figure 18 is blocked.
The multiple dimensioned full waveform inversion result figure of Figure 19 LPFs (is tested) containing Gaussian noise.
The multiple dimensioned full waveform inversion result figure of window LPF (being tested containing Gaussian noise) when Figure 20 is blocked.
Embodiment
Below in conjunction with the accompanying drawings with example to the further detailed description of the present invention
The multiple dimensioned full waveform inversion method of the LPF of window, comprises the following steps when blocking：
A, program is that completion is write under MATLAB2013a software frames, it is necessary first to MATLAB2013a softwares are installed,
And the seismic data processing software Crews kits that MATLAB language is write are installed；
B, the geological data of collection is needed to be pre-processed first, and observational record and observation system input computer are carried out
Time-domain full waveform inversion.Wherein pretreatment includes：
B1, multiple attenuation, eliminate reverberation and suppress the seismic wave that ghosting etc. is unable to forward modeling.
B2, underfrequency protection denoising, the compensation of missing seismic channel, relative amplitude preserved processing etc..
B3, define observation system.
C, source wavelet is estimated from observational record, the source wavelet is used for forward modeling, and (present invention is replaced using Ricker wavelet
Source wavelet (Fig. 2) in real data)；
D, linear increment initial model (Fig. 4 is left) is established, initial value of the rate pattern as inverting, passes through computation model
More new direction, to realize from initial model gradually to the close process of true model；
E, object function is constructed according to the principle of least square：Wherein v represents underground P-wave
Degree, dobsFor the observational record of collection, dcalIn the forward record that rate pattern is obtained by forward modeling.Seeking the gradient of object function
During to object function both ends on p wave interval velocity v derivations, obtain pressure gradient expression formula：
Wherein pbRepresent by with focus fsThe wave field obtained to underground propagation, pfForward modeling wave field is represented, s represents focus
Number, r represent the number of wave detector.
F, window (Fig. 7) when blocking is chosen, intercepts observational record forward part seismic signal, and records and blocks time window length, should
Time span of the time span of window as forward modeling when blocking；
G, according to the dynamic law of ACOUSTIC WAVE EQUATION, the source wavelet that is obtained using the observation system of input, estimation and just
Beginning rate pattern carries out seismic wave field forward modeling.The geological data when step only needs the forward modeling to block in window, it is not necessary to will be all
Time geological data whole forward modeling comes out.Store main story wavefield data and forward record.
The LPF multiple dimensioned full waveform inversion method relevant parameter of window when setting is blocked as requested, including model are big
Small nz × nx, grid is away from dx, dz, forward modeling time window length It, time sampling interval dt, Ricker wavelet dominant frequency f0, Butterworth filter
Cut frequency (the f of ripple devicecut), the maximum iteration itermax of each super random epicentre, update gradient minimum value gtol, mesh
Scalar functions precision tol, the maximum vmax and minimum value vmin of inversion result；
Super random epicentre coding：Time-domain full waveform inversion takes a substantial amount of time during single-shot forward modeling, someone
It is proposed using while excite multiple focus (super big gun) to replace single-shot explosive source, can so shorten required for forward modeling when
Between, but carry out full waveform inversion using explosive source simultaneously and bring crosstalk noise to inversion result again.In order to improve all-wave
The computational efficiency of shape inverting suppresses this software of crosstalk noise and utilizes super random epicentre coding strategy simultaneously.
The Ricker wavelet that the super gun excitation of random phase encoding earth's surface goes out, it is desirable to phase is positive and negative to be randomly assigned, while each
Single-shot in super big gun excites delay different, and delay maximum is no more than 400ms.If random phase encoding information is：
Wherein s' ∪ s={ 1,2 ... Ns, the number of element is fixed in s, and it represents contained single-shot number in a super big gun
Number.Single-shot quantity can adjust mixed proportion according to the actual requirements in super big gun.Element size is randomly assigned in s.i∈N+Represent
Random positive integer.Super random epicentre coding can be expressed as：
Spf=DyA Γ TLFf (4)
Wherein Dy represents that (meaning of dynamic coding is exactly in supermemory gradient method often iteration 10 times to dynamic focus coding function
Random phase encoding information is re-replaced afterwards, finally achievable earth's surface all standing.This method can effectively increase earth's surface big gun
Coverage density, while weaken influence of the crosstalk noise to full waveform inversion, finally realize the huge leap forward in computational efficiency.),
A represents random amplitude coding function, and Γ represents random phase encoding function, and T represents random delay coding function, and L represents random
Position encoded function, F represent random Ricker wavelet dominant frequency coding function, f source functions.
Observational record and forward record when h, to blocking in window while carry out LPF, it is ensured that with identical low pass filtered
Ripple device (Fig. 3) processing data.After filtering, obtain handling observational record and handle forward record；
Butterworth LPF：Carry out the system function of approximate wave filter with Butterworth function.Butterworth filters
Device is the wave filter that there is most smoothness properties to define according to amplitude versus frequency characte in passband.
The low pass squared magnitude function of Butterworth filter represents
The principal character of Butterworth filter is as follows：
H1, to all N,
H2, to all N,I.e.
H3、|Ha(jΩ)|2It is Ω monotonically decreasing function；
H4、|Ha(jΩ)|2With order N increase closer to ideal low-pass filter；
The amplitude versus frequency characte of wave filter becomes to become better and better with filter order N increase, in cut-off frequency ΩcPlace
There is the value in more frequency band areas in the case that functional value is always 1/2, in passband close to 1；The more rapid convergence in stopband
In zero.
I, processing observational record and processing forward record are made the difference, obtains residual error data；
J, using residual error data as with focus, acted on anti-pass operator with focus, obtain the residual of the model space
Poor anti-pass wave field；
K, main story wave field does cross-correlation with residual error anti-pass wave field and obtains model modification gradient；
The object function of full waveform inversion is commonly defined as forward modeling data with observing the L2 functionals of data residual error, and form is such as
Under：
Wherein DcalIt is forward record, DcalRepresent the observational record in actual field.Object function is obtained to speed v derivations
Arrive：
The present invention is for apparent elaboration theoretical method, it is assumed that density is constant, then only needs to carry out inverting to speed.
Chang Midu ACOUSTIC WAVE EQUATIONs are：
Wherein u is acoustic wavefield, and f is source function.V represents the acoustic velocity of underground medium.
ACOUSTIC WAVE EQUATION is simplified to the form for being write as operator：
Lu=s (9)
Wherein forward modeling operator：And to velocity derivatives：
Equation (9) both sides obtain (Pratt et al., 1998) to speed v derivations simultaneously：
Obtain wave field is to subsurface velocities v derivative：
Equation (11) is substituted into equation (7) and obtained：
Wherein L-1Represent adjoint operator or be anti-pass operator, L-1[Dcal-Dobs] represent residual error anti-pass wave field.fs=Dcal-
DobsReferred to as with focus.
L, with super-memory gradient method optimized algorithm computation model more new direction, and found by Wolfe convergence criterions it is suitable
Step-length.
Supermemory gradient method：School is carried out to the information of current gradient with current gradient and the before information of multi-step gradient
Just, to accelerate the convergence rate of object function.Its iteration form is：
vk+1=vk+αkdk (13)
Wherein vkRepresent iteration kth portion rate pattern, αkRepresent step-length, dkRepresent descent direction.
Supermemory gradient method descent direction can be represented by following form：
As k ＞ mem, βk-i+1∈[0,sk i], (i=1 ... mem)
Here
Wherein, 0 ＜ ρ ＜ 1, mem ＞ 0 are a positive integers, and we term it memory degree by mem.Different memory degree target letters
Number convergence rate and model inversion precision are different, and Memory Gradient quantity, which is crossed, can at most cause gradient direction correction excessive, target letter
Number decline is slack-off, and Memory Gradient quantity deteriorates to conjugate gradient method at least.Obtained by experimental verification, in seismic full-field shape inverting
In best memory degree substantially in mem=4,5.Supermemory gradient method has certain advantage, but when gradient remembers excessive, calculating
Slow, memory requirements increase.So advantage that the suitable memory degree competence exertion supermemory gradient method of selection is maximum.
M, judge it is enough to meet end condition, the multiple dimensioned full waveform inversion of LPF of window when output is blocked if meeting
Inversion method result.If being unsatisfactory for end condition, the initial model using inversion result as next circulation, Step d is returned to.
N, the initial velocity model using the output result of previous step as Conventional Time domain full waveform inversion, then carry out normal
Advise time-domain full waveform inversion.Final inversion result is exported after iteration 200 times.
Embodiment 1
Required, MATLAB2013a is installed under the flagship edition systems of Windows 7, and install according to exploration
The seismic data processing software Crews kits that MATLAB language is write.
The present invention is tested using Marmousi p wave interval velocity models, and Marmousi rate patterns have complicated structure
Make, trickle stratification and deep oil storage Tibetan construction, be well suited for carrying out theoretical method test.Original Marmousi rate patterns compared with
Greatly, it is contemplated that hardware memory and calculating time factor, the present invention carries out vacuating processing to archetype, after vacuating
The LPF of window multiple dimensioned full waveform inversion method test when Marmousi rate patterns are blocked.True model (Fig. 4
It is right) and linear increment initial model (Fig. 4 is left).
Model parameter is as follows：
The LPF of window multiple dimensioned full waveform inversion method test parameter (time-domain) when table 1 blocks
| Model size | Grid spacing | Lateral separation | Longitudinal depth | Velocity interval | Dominant frequency of seismic wavelet |
| 69*384 | 12.5m | 2.4km | 0.8625km | 1.5-4km/s | 22Hz |
Window refutation strategy parameter is as follows when blocking：
Observational record is handled using window when blocking, time window length It=0.2 is blocked in setting:0.1:2.1s, that is, select
Window 0-0.2s when starting is blocked is selected, the time span increase 0.1s of window when then being blocked for 10 times per iteration, and it is constant to originate for 0 moment
(Fig. 7), the total length of window is 2.1s when blocking.The process totally 20 windows when blocking, full waveform inversion iteration 200 times altogether.
Time span when forward modeling only needs forward modeling to corresponding to, so can largely save the calculating time required for forward modeling.Cut
Window refutation strategy is successively carried out inverting, is greatly reduced the parameter of inverting each time, effectively avoided by shallow-layer to deep layer when disconnected
The generation of cycle slip phenomenon, while enhance the stability of program.
LPF Multi-scale inversion policing parameter is as follows：
In order to overcome the problems such as being absorbed in local minimum and cycle slip during full waveform inversion, and also to improve imaging matter
Amount, full waveform inversion is carried out with LPF multi-scale method, sets two cut frequency of 15Hz and 25Hz.Two cut frequency
Whole full waveform inversion process is divided into three frequency ranges：
Ⅰ、0-15Hz：Low-frequency range.The macroscopic information of subsurface model can be obtained by carrying out inverting in low-frequency range, in low-frequency range
Initial model of the inversion result as next frequency band inverting, so gradually rate pattern accurately can be finally inversed by.
But low-frequency information is often lacked in real data, land geological data lacks below 10Hz data, and LPF simply exists
Competence exertion acts in the case of low-frequency data being present, the effect of the multiple dimensioned full waveform inversion of LPF when lacking low-frequency data
Fruit is not fine (Figure 13).So the present invention in low-frequency range combine window refutation strategy when blocking eliminate because lack low frequency and
Caused cycle slip phenomenon.
Ⅱ、0-25Hz：Mid Frequency.The frequency range primarily serves the effect of an intermediate transition, and the inversion result of low-frequency range is made
For the initial velocity model of Mid Frequency.
Ⅲ、0-end：Full frequency band.The initial velocity model of the inversion result Whole frequency band of midband, the step can obtain
The inversion result of time-domain Full wave shape.
Based on above parameter, full waveform inversion is carried out under the test environment shown in table 2.
The computing power test environment of table 2
The full waveform inversion method of table 3 contrasts
As can be drawn from Table 3, Conventional Time domain full waveform inversion and the LPF of window during proposed by the present invention block are more
Yardstick full waveform inversion has obvious difference (model error, in terms of calculating the time).The meter of Conventional Time domain full waveform inversion
Evaluation time is 101.3 minutes, and the multiple dimensioned full waveform inversion of LPF of window only needs 43.6 minutes cans to complete when blocking
Whole calculating process.Window is low when from the angle combination inversion result figure (Fig. 9,10,11,12) of model error it can be seen that blocking
The multiple dimensioned full waveform inversion inversion result precision of pass filter is higher, and the rate results that inverting obtains are more accurate, and Conventional Time
Even there is mistake in the result of domain full waveform inversion, is fallen far short with true model.
Based on above parameter, various factors test is carried out under the test environment shown in table 2：
The LPF of window multiple dimensioned full waveform inversion method test result when being blocked under the different affecting factors of table 4
Note：Model inversion normalization error calculation formula is in table 3 and table 4：Wherein vtrue generations
Table true velocity model, vk represent inversion result.
From table 4 and Figure 13, as can be seen that missing low frequency has a serious impact to full waveform inversion in 14,15.Work as earthquake
When low frequency is lacked in data, Conventional Time domain full waveform inversion result (Figure 13), the multiple dimensioned full waveform inversion of LPF
As a result (Figure 14) and manifest error all occurs in window full waveform inversion result (Figure 15) when blocking, model inversion error is substantially inclined
It is high.Simultaneously as the factor of missing low frequency causes model modification gradient calculation result incorrect so that carrying out, Wolfe is linear
The substantial amounts of time is wasted during step-size in search, this explains why identical iterations, when lacking low frequency
The calculating time of full frequency band full waveform inversion will be considerably longer than by waiting the calculating time of full waveform inversion.In order to preferably solve because
The inverting Problem-Error that missing low-frequency component is brought, window refutation strategy and LPF multi-scale strategy knot when the present invention will block
Close and use, as can be seen from Figure 16, combination strategy causes the result of full waveform inversion to obtain obvious improvement.Model error
It also illustrate that, the inversion result that the strategy obtains is more nearly with true model.
From table 4 and Figure 17,18,19,20 as can be seen that when Gaussian noise is contained in geological data, inverting knot
Fruit becomes unintelligible.Conventional Time domain full waveform inversion (Figure 17) and when blocking window full waveform inversion (Figure 18) all occur it is bright
Aobvious mistake, although window refutation strategy accelerates inverting process to a certain extent when illustrating to block, its anti-noise ability compared with
It is weak.The multiple dimensioned full waveform inversion of LPF (Figure 19) has obvious advantage compared with the above two, it was demonstrated that LPF is multiple dimensioned
Strategy has stronger anti-noise ability.Window refutation strategy and LPF Multi-scale inversion strategy are used in combination with when will block
(Figure 20), the defects of window anti-noise ability is weak when not only accelerating the process of inverting, and compensate for blocking.
Fig. 1 is the flow chart of whole refutation process, from flow chart it can be seen that anti-using window when blocking in low-frequency range first
Drill strategy and carry out full waveform inversion, window inversion result is used for next step intermediate frequency as initial model when then low-frequency range being blocked again
Section inverting, finally by the initial model of the inversion result of Mid Frequency Conventional Time domain full waveform inversion the most.Window is low when this is blocked
The multiple dimensioned full waveform inversion strategy of pass filter, in inversion accuracy, computational efficiency, low frequency dependence, anti-noise ability etc. better than normal
Advise time-domain full waveform inversion.
Claims (1)
1. a kind of multiple dimensioned full waveform inversion method of LPF of window when blocking, it is characterised in that comprise the following steps：
MATLAB2013a softwares are installed, and the seismic data processing software Crews kits that MATLAB language is write are installed a,；
B, the geological data of collection is pre-processed, and observational record and observation system is inputted into computer, carry out time-domain
Full waveform inversion；
C, source wavelet is estimated from observational record, the source wavelet that the estimation obtains is used for forward simulation；
D, linear increment initial velocity model is established, initial value of the rate pattern as inverting, passes through the renewal of computation model
Direction, realize the process gradually approached from initial model to true model；
E, object function is built according to principle of least square method：
Wherein v represents p wave interval velocity, dobsFor the observational record of collection, dcalIn the forward record that rate pattern is obtained by forward modeling,
In the gradient procedure of object function is sought, the p wave interval velocity v derivations to object function both ends, pressure gradient expression formula is obtained：
<mrow>
<mfrac>
<mrow>
<mo>&part;</mo>
<mi>E</mi>
<mrow>
<mo>(</mo>
<mi>v</mi>
<mo>)</mo>
</mrow>
</mrow>
<mrow>
<mo>&part;</mo>
<mi>v</mi>
</mrow>
</mfrac>
<mo>=</mo>
<mo>-</mo>
<mfrac>
<mn>2</mn>
<msup>
<mi>v</mi>
<mn>3</mn>
</msup>
</mfrac>
<munder>
<mo>&Sigma;</mo>
<mi>s</mi>
</munder>
<munder>
<mo>&Sigma;</mo>
<mi>r</mi>
</munder>
<mfrac>
<mrow>
<msup>
<mo>&part;</mo>
<mn>2</mn>
</msup>
<msub>
<mi>p</mi>
<mi>f</mi>
</msub>
</mrow>
<mrow>
<mo>&part;</mo>
<msup>
<mi>t</mi>
<mn>2</mn>
</msup>
</mrow>
</mfrac>
<mo>*</mo>
<msub>
<mi>p</mi>
<mi>b</mi>
</msub>
</mrow>
Wherein pbRepresent by with focus fsThe wave field obtained to underground propagation, pfForward simulation wave field is represented, s represents focus
Number, r represent the number of wave detector；
F, the time span of window when blocking is chosen, intercepts observational record forward part seismic signal, and records and blocks time window length, should
Time span of the time span of window as forward simulation when blocking；
G, according to the dynamic law of ACOUSTIC WAVE EQUATION, the source wavelet obtained using the observation system of input, estimation and initial speed
Spend model and carry out seismic wave field forward modeling, store main story wavefield data and forward record；
The LPF multiple dimensioned full waveform inversion method relevant parameter of window when setting is blocked as requested, including model size nz
× nx, grid is away from dx, dz, forward simulation time window length It, time sampling interval dt, Ricker wavelet dominant frequency f0, Butterworth filter
The cut frequency f of ripple devicecut, the maximum iteration itermax of each super random epicentre, update gradient minimum value gtol, target
Function precision tol, the maximum vmax and minimum value vmin of inversion result；
Observational record and forward record when h, to blocking in window while carry out LPF, it is ensured that with identical low pass filter
Processing data, after filtering, obtain handling observational record and handle forward record；
I, processing observational record and processing forward record are made the difference, obtains residual error data；
J, using residual error data as with focus, acted on anti-pass operator with focus, the residual error for obtaining the model space is anti-
Pass wave field；
K, main story wave field does cross-correlation with anti-pass residual error wave field, obtains model modification gradient；
L, with super-memory gradient method optimized algorithm computation model more new direction, and pass through Wolfe convergence criterions and find suitable step-length；
M, judge whether to meet end condition, the multiple dimensioned full waveform inversion method of LPF of window when output is blocked if meeting
Inversion result；If being unsatisfactory for end condition, the initial model using inversion result as next circulation, Step d is returned to；
N, the initial velocity model using m steps output result as Conventional Time domain full waveform inversion, Conventional Time is then carried out
Domain full waveform inversion, final inversion result is exported after iteration 200 times.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201610429648.XA CN106054244B (en) | 2016-06-16 | 2016-06-16 | The LPF of window multiple dimensioned full waveform inversion method when blocking |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201610429648.XA CN106054244B (en) | 2016-06-16 | 2016-06-16 | The LPF of window multiple dimensioned full waveform inversion method when blocking |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN106054244A CN106054244A (en) | 2016-10-26 |
| CN106054244B true CN106054244B (en) | 2017-11-28 |
Family
ID=57169034
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201610429648.XA Expired - Fee Related CN106054244B (en) | 2016-06-16 | 2016-06-16 | The LPF of window multiple dimensioned full waveform inversion method when blocking |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN106054244B (en) |
Families Citing this family (10)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN107450102A (en) * | 2017-07-28 | 2017-12-08 | 西安交通大学 | Multiple dimensioned full waveform inversion method based on the controllable envelope generating operator of resolution ratio |
| CN109541681B (en) * | 2017-09-22 | 2020-07-17 | 中国海洋大学 | Wave inversion method for combining streamer seismic data and small amount of OBS data |
| CN107765308B (en) * | 2017-10-12 | 2018-06-26 | 吉林大学 | Reconstruct low-frequency data frequency domain full waveform inversion method based on convolution thought Yu accurate focus |
| CN108549100B (en) * | 2018-01-11 | 2019-08-02 | 吉林大学 | The multiple dimensioned full waveform inversion method of time-domain for opening up frequency based on non-linear high order |
| CN109143336B (en) * | 2018-08-03 | 2019-09-13 | 中国石油大学(北京) | A method of overcome the period in full waveform inversion to jump |
| CN112578431B (en) * | 2019-09-27 | 2024-04-09 | 中国石油化工股份有限公司 | Method and system for storing full waveform inversion wave field optimization in finite state |
| CN111722287B (en) * | 2020-06-19 | 2021-10-08 | 南京大学 | Seismic phase characteristic identification waveform inversion method based on progressive data assimilation method |
| US11467299B2 (en) | 2020-12-16 | 2022-10-11 | Saudi Arabian Oil Company | Full waveform inversion velocity guided first arrival picking |
| CN113552625B (en) * | 2021-06-21 | 2022-05-13 | 中国地质科学院地球物理地球化学勘查研究所 | Multi-scale full waveform inversion method for conventional land-domain seismic data |
| CN116540297B (en) * | 2023-05-06 | 2024-03-08 | 中国科学院地质与地球物理研究所 | Full waveform inversion method, system and equipment for elastic wave seismic data |
Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN103499835A (en) * | 2013-10-13 | 2014-01-08 | 中国石油集团西北地质研究所 | Method for inverting near-surface velocity model by utilizing preliminary waveforms |
| CN105093278A (en) * | 2014-05-16 | 2015-11-25 | 中国石油化工股份有限公司 | Extraction method for full waveform inversion gradient operator based on excitation main energy optimization algorism |
| CN105549080A (en) * | 2016-01-20 | 2016-05-04 | 中国石油大学(华东) | Undulating ground surface waveform inversion method based on auxiliary coordinate system |
| CN105467444B (en) * | 2015-12-10 | 2017-11-21 | 中国石油天然气集团公司 | A kind of elastic wave full waveform inversion method and device |
Family Cites Families (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US8437998B2 (en) * | 2010-09-27 | 2013-05-07 | Exxonmobil Upstream Research Company | Hybrid method for full waveform inversion using simultaneous and sequential source method |
-
2016
- 2016-06-16 CN CN201610429648.XA patent/CN106054244B/en not_active Expired - Fee Related
Patent Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN103499835A (en) * | 2013-10-13 | 2014-01-08 | 中国石油集团西北地质研究所 | Method for inverting near-surface velocity model by utilizing preliminary waveforms |
| CN105093278A (en) * | 2014-05-16 | 2015-11-25 | 中国石油化工股份有限公司 | Extraction method for full waveform inversion gradient operator based on excitation main energy optimization algorism |
| CN105467444B (en) * | 2015-12-10 | 2017-11-21 | 中国石油天然气集团公司 | A kind of elastic wave full waveform inversion method and device |
| CN105549080A (en) * | 2016-01-20 | 2016-05-04 | 中国石油大学(华东) | Undulating ground surface waveform inversion method based on auxiliary coordinate system |
Non-Patent Citations (1)
| Title |
|---|
| 基于特征能量梯度算子的时间域全波形反演;张文祥 等;《地球物理学进展》;20151231;第30卷(第1期);第363-371页 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN106054244A (en) | 2016-10-26 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN106054244B (en) | The LPF of window multiple dimensioned full waveform inversion method when blocking | |
| CN105467444B (en) | A kind of elastic wave full waveform inversion method and device | |
| CN103713315B (en) | A kind of seismic anisotropy parameter full waveform inversion method and device | |
| CN107422379A (en) | Multiple dimensioned seismic full-field shape inversion method based on local auto-adaptive convexification method | |
| CN106526674A (en) | Three-dimensional full waveform inversion energy weighted gradient preprocessing method | |
| CN107505654A (en) | Full waveform inversion method based on earthquake record integration | |
| CN108549100B (en) | The multiple dimensioned full waveform inversion method of time-domain for opening up frequency based on non-linear high order | |
| CN109407151B (en) | Time-domain full waveform inversion method based on wave field local correlation time shift | |
| CN104570082B (en) | Extraction method for full waveform inversion gradient operator based on green function characterization | |
| CN110058302A (en) | A kind of full waveform inversion method based on pre-conditional conjugate gradient accelerating algorithm | |
| CN106908835A (en) | Band limit Green's function filters multiple dimensioned full waveform inversion method | |
| CN107765302A (en) | Inversion method when time-domain single-frequency waveform independent of source wavelet is walked | |
| CN105388518A (en) | Centroid frequency and spectral ratio integrated borehole seismic quality factor inversion method | |
| CN105549079A (en) | Method and device for establishing full-waveform inversion model for geophysics parameters | |
| CN105891888A (en) | Multi-domain frequency-division parallel multi-scale full-waveform inversion method | |
| CN107765308B (en) | Reconstruct low-frequency data frequency domain full waveform inversion method based on convolution thought Yu accurate focus | |
| CN104965223B (en) | Viscoelastic acoustic wave full-waveform inversion method and apparatus | |
| CN111025387B (en) | Pre-stack earthquake multi-parameter inversion method for shale reservoir | |
| CN107894618B (en) | A kind of full waveform inversion gradient preprocess method based on model smoothing algorithm | |
| CN111045077B (en) | Full waveform inversion method of land seismic data | |
| CN105093278A (en) | Extraction method for full waveform inversion gradient operator based on excitation main energy optimization algorism | |
| CN106033124A (en) | Multi-seismic resource sticky sound least square reverse time migration method based on stochastic optimization | |
| CN104808243A (en) | Prestack seismic Bayesian inversion method and prestack seismic Bayesian inversion device | |
| CN112327358A (en) | Acoustic seismic data forward modeling method in viscous medium | |
| CN104391324A (en) | Seismic trace set dynamic correction stretching correction pre-processing technology before AVO inversion depending on frequency |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| C06 | Publication | ||
| PB01 | Publication | ||
| C10 | Entry into 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: 20171128 Termination date: 20180616 |