CN106054244A - Low-pass filter multi-scale full waveform inversion method of cut-off time window - Google Patents

Abstract

The invention relates to a low-pass filter multi-scale full waveform inversion method of a cut-off time window. The method comprises the steps of firstly, in a low band, performing full waveform inversion using a cut-off time window inversion strategy to obtain an inversion result serving as an initial model for intermediate band inversion; and then performing full waveform inversion of a conventional time domain using the intermediate band inversion result as an initial model to obtain a final inversion result. The cut-off time window inversion strategy saves a lot forward time. The low-pass filter multi-scale inversion strategy can divide seismic data to different bands for inversion. A good inversion result can be obtained by using the cut-off time window and the low-pass filter together under the condition of low frequency deficiency and strong Gaussian noise interference. Model test results show that the combined strategy relieves a cycle-slip phenomenon while improving the calculation efficiency, and the low-pass filter multi-scale full waveform inversion method of the cut-off time window is superior to the full waveform inversion of the conventional time domain on the aspects of inversion precision, calculation efficiency, low-frequency dependency, anti-noise capability and the like.

Description

The low-pass filtering of window multiple dimensioned full waveform inversion method when blocking

Technical field

The present invention relates to the data processing method of a kind of seismic prospecting, when especially blocking, window and low-pass filtering are multiple dimensioned instead Drill strategy.Whole calculating process utilizes super random epicentre coding strategy to carry out full waveform inversion, it is therefore an objective to accelerate whole Full wave shape The process of inverting.

Background technology:

Along with the development of petroleum industry and deepening continuously of exploration and development, gradually move towards lithology exploration from the structure exploration phase In the stage, its difficulties in exploration is gradually increased.In order to comply with this requirement, seismic data processing technology is the most perfect, the most entirely Waveform inversion develops rapidly, and becomes the study hotspot of current geophysics circle.The eighties in 20th century, Tarantola carries Go out time domain full waveform inversion reason method, and model parameter derivation has been passed through residual error anti-pass wave field and main story ripple by object function Field computing cross-correlation obtains, thus avoids and ask for Jacobi matrix.The nineties in 20th century, full waveform inversion is promoted by Pratt Arrive frequency domain, proposed to have only to several discrete frequency and can be obtained by high-precision inversion result, and low frequency is to high frequency The refutation strategy problem that can solve to be absorbed in local minimum.Shin et al. development on the basis of frequency domain full waveform inversion Laplace-Fourier territory full waveform inversion, the method can obtain a high accuracy in the case of disappearance low frequency Initial model.Wang et al. (2011) proposes acceleration full waveform inversion based on CUDA programming on GPU, and the method adds effectively The fast computational efficiency of full waveform inversion.

Bunks (1995) proposes multiple dimensioned full waveform inversion based on low-pass filtering, and the method is pointed out first to observation data Carry out low-pass filtering, then the source wavelet estimated is processed in order to same filter, Full wave shape can be to some extent solved The cycle slip problem of inverting.But utilize this multi-scale method to carry out full waveform inversion and still there are some problems, such as: stratum is One viscoelastic medium, and Wave equation forward modeling is a complicated process, when seismic wave meeting during underground propagation Frequency shift effect or other complicated changes occur.Source wavelet after utilizing low-pass filtering is just drilled, it may appear that just drilling note Record and the unmatched problem of observational record, do not provide good solution about this problem Bunks et al..Chen and Wu et al. utilizes time attentuating filter to process observation data and just drilling data simultaneously, be i.e. multiplied by geological data one along with The exponential function that time increases and decays, the most first drills shallow-layer information, ignores infrastructure.Along with inverting carries out being gradually reduced declining The impact of subtraction function so that inverting by shallow-layer to deep layer Layer by layer inversion.The method is greatly reduced the inverted parameters of each step, Solve to a certain extent full waveform inversion cycle slip problem.When but the method wastes substantial amounts of during just drilling Between, because just drilling the forward record obtained to be multiplied by an attenuation function again, the deep information decays to zero substantially.The method does not has profit Using the deep information, but wasting the plenty of time is just drilling deep wave field, largely reduces the calculating effect of full waveform inversion Rate.

Summary of the invention:

It is an object of the invention to for data mismatch problem after above-mentioned existing low-pass filtering, it is provided that window when one is blocked The multiple dimensioned full waveform inversion method of low-pass filtering.

The solution of the present invention is: utilize Butterworth smoothing low-pass filters process observational record simultaneously and just drilling Record, so effectively prevent factor data and does not mates and produce wrong inversion result, and the method is more suitable for actual seismic data Handling process, can obtain accurate inversion result.For the problem that time domain full waveform inversion computational efficiency is low, the present invention carries Going out window full waveform inversion strategy when blocking, i.e. just starting inverting when, only the intercepting real data forward part time is carried out instead Drill, so the time window data also having only to just to drill correspondence carry out Data Matching, so by just drill block time window in wave field come Replace just drilling All Time wave field, highly shortened and just drill the required time.Along with the carrying out of inverting, it is gradually increased inverting Time window length, when this blocks, window refutation strategy can reduce the number of each step inverted parameters, adds stablizing of program simultaneously Property and the reliability of inversion result.

It is an object of the invention to be achieved through the following technical solutions:

When blocking, the core of the multiple dimensioned full waveform inversion method of the low-pass filtering of window is complete on MATLAB 2013a platform Become.Butterworth filter is utilized to observation data and just to drill data and carry out low-pass filtering, the most flat by setting Butterworth The cut frequency that filter slide is different, realizes time domain Multi-scale inversion strategy.Low pass filter can effectively filter out earthquake Radio-frequency component in data, retains the low-frequency component in geological data, owing to low-frequency component contains subsurface structure macroscopic information, has It is beneficial to the foundation of initial model, avoids the generation of cycle slip phenomenon simultaneously.Window refutation strategy when blocking, i.e. select different length time Between window extract observation data for full waveform inversion, just drill according to the time span that window time this is corresponding, so simultaneously Each step can be avoided the most just to drill all time datas, saved the calculating time largely.

When blocking, window refutation strategy and low-pass filtering Multi-scale inversion strategy itself all have compacting cycle slip and improve inverting The effect of precision.In order to obtain more stable inverting energy flow process, more accurate inversion result and more efficient inversion method, this Bright is further the two to be used in combination with (as shown in Figure 1).

The low-pass filtering of window multiple dimensioned full waveform inversion method when blocking, comprises the following steps:

A, installation MATLAB2013a software, and the seismic data processing software Crews instrument that MATLAB language is write is installed Bag；

B, the geological data gathered is carried out pretreatment, and observational record and observation system are inputted computer carry out the time Territory full waveform inversion；

C, from observational record, estimate source wavelet, utilize and estimate that the source wavelet that obtains just drills that (present invention utilizes Ricker wavelet replaces the source wavelet in real data)；

D, setting up linear increment initial model, this rate pattern is as the initial value of inverting, by the renewal of computation model Direction, it is achieved by initial model gradually to the process that true model is close；

E, according to the principle of least square construct object function:Wherein v represents underground P-wave Degree, d_obsFor the observational record gathered, d_calAt rate pattern by just drilling the forward record obtained.Seeking the gradient of object function During to object function two ends about p wave interval velocity v derivation, obtain gradient expression formula:

$\frac{\∂E\left(v\right)}{\∂v}=-\frac{2}{v^3}\underset{s}{\Σ}\underset{r}{\Σ}\frac{{\∂}^2{p}_f}{\∂t^2}*p_b---\left(1\right)$

Wherein p_bRepresent by with focus f_sThe wave field obtained to underground propagation, p_fRepresenting and just drilling wave field, s represents focus Number, r represents the number of cymoscope.

F, choose the time span of window when blocking, intercept observational record forward part seismic signal, and it is long to record window when blocking Degree, when this blocks, the time span of window is as the time span just drilled；

G, dynamic law according to ACOUSTIC WAVE EQUATION, utilize the observation system of input, estimate the source wavelet that obtains and just Beginning rate pattern carries out seismic wave field just drills.This step has only to just drill the geological data in window when blocking, it is not necessary to will be all Time geological data is the most just performed.Storage main story wavefield data and forward record.

Set the low-pass filtering multiple dimensioned full waveform inversion method relevant parameter of window when blocking as requested, big including model Little nz × nx, grid, away from dx, dz, is just drilling time window length It, time sampling interval dt, Ricker wavelet dominant frequency f₀, Butterworth is filtered Cut frequency f of ripple device_cut, maximum iteration time itermax of each super random epicentre, update gradient minima gtol, target Function precision tol, the maximum vmax of inversion result and minima vmin；

H, the observational record in window when blocking and forward record are carried out low-pass filtering simultaneously, it is ensured that by identical low pass filtered Ripple device processes data.After filtering, obtain processing observational record and processing forward record；

I, by process observational record and process forward record do difference, obtain residual error data；

J, using residual error data as with focus, act on on focus with anti-pass operator, obtain the residual of the model space Difference anti-pass wave field；

K, main story wave field and residual error anti-pass wave field do cross-correlation, obtain model modification gradient；

L, use super-memory gradient method optimized algorithm computation model more new direction, and found suitably by Wolfe convergence criterion Step-length；

M, judgement are enough to meet end condition, if meeting, and the multiple dimensioned full waveform inversion of low-pass filtering of window when output is blocked Inversion method result.If being unsatisfactory for end condition, the initial model that inversion result is circulated as the next one, return Step d；

N, using the output result of previous step as the initial velocity model of Conventional Time territory full waveform inversion, often then carry out Rule time domain full waveform inversion.Iteration exports final inversion result 200 times afterwards.

Beneficial effect: by the improvement of original theory method and the integration of relative earthquake full waveform inversion technology, this Bright when successfully will block the technology such as window and low-pass filtering be applied in full waveform inversion.

The low-pass filtering of window multiple dimensioned full waveform inversion method when blocking proposed, the method can not only be alleviated because of low frequency The disappearance cycle slip problem brought, and the saving of the high degree time just drilled.When utilization is blocked, window refutation strategy significantly carries The computational efficiency of high full waveform inversion and inversion accuracy.Butterworth filter is utilized to observation data and just to drill data and carry out Low-pass filtering, filters radio-frequency component, retains low-frequency component, by setting different the blocking frequently of Butterworth smoothing filter Rate, realizes frequency time domain Multi-scale inversion strategy from low to high, has effectively suppressed the cycle slip phenomenon of full waveform inversion. This refutation strategy improves the precision of full waveform inversion to a certain extent, and is more suitable for actual seismic data and processes, simultaneously this The low-pass filtering refutation strategy of bright proposition without seismic wavelet is filtered, be prevented effectively from factor data do not mate and cause anti- Drill mistake.During whole calculating, focus all uses super random epicentre coding strategy, it is therefore an objective to accelerating the same of inverting speed Time suppress the crosstalk noise within super big gun.Compared with the full waveform inversion of Conventional Time territory, the mistake of inverting subsurface velocities of the present invention Window and low-pass filtering refutation strategy when Cheng Zhongyong blocks, solve problems with:

When 1, first will block, window and low-pass filtering strategy combination, be applied in full waveform inversion, can quickly obtain One high-precision initial velocity model, this initial velocity model has recovered the macroscopic information of Marmousi model substantially, profit Carry out conventional full wave shape inverting by this initial velocity model, cycle slip problem can be alleviated.

2, time histories sample will be blocked to be incorporated in full waveform inversion so that just drilling the ripple that need not calculate all moment every time , the wave field only calculated when blocking in window, highly shortened and just drill the required calculating time.

When 3, blocking, window low-pass filtering Multi-scale inversion strategy can effectively reduce the unknown number number of inverting each time, according to Gradient and model modification direction that cross-correlation is tried to achieve are more accurate.Owing to, during model modification iteration, needing to pass through Wolfe linear search find suitable iteration step length, when gradient direction or descent direction inaccurate in the case of, Wolfe line Property search can waste the substantial amounts of time, and window low-pass filtering refutation strategy when blocking that the present invention proposes accurate can obtain mould The gradient direction that type updates, so saves the substantial amounts of calculating time during linear search.

4, the inverting of conventional full wave shape needs to store main story wave field and residual error anti-pass wave field, and storage the two wave field needs a large amount of Calculator memory.Window refutation strategy when blocking that the present invention proposes, it is only necessary to during storage, the main story wave field in window and residual error are anti- Biography wave field, has saved the memory headroom of required computer, to a certain extent especially in real data processes.

5, super random epicentre coding strategy can effectively suppress the crosstalk noise within super big gun, and the method is ensureing inverting essence Decrease the number of times just drilled while degree to a certain extent, be effectively shortened the calculating time of full waveform inversion.

Model test results shows, even if in the case of disappearance low-frequency data and the strongest Gauusian noise jammer, utilizes When blocking, the full waveform inversion of window low-pass filtering strategy can also obtain preferable inversion result.This refutation strategy can realize Low-frequency range, by the Layer by layer inversion method of shallow-layer to deep layer, can be significantly improved in terms of computational efficiency again, alleviate simultaneously The cycle slip phenomenon caused because of disappearance low-frequency data.Can be seen that the present invention proposes when blocking, window low-pass filtering refutation strategy can A high-precision initial velocity model is provided for conventional full wave shape inverting.

Accompanying drawing explanation

The low-pass filtering of window multiple dimensioned full waveform inversion method flow chart when Fig. 1 blocks.

Fig. 2 Ricker wavelet

(left) Ricker wavelet oscillogram, (right) Ricker wavelet spectrogram.

Fig. 3 Butterworth smoothing filter figure.

Fig. 4 rate pattern

(left) linear increment initial velocity model figure, (right) true velocity model figure.

Fig. 5 filter shape figure and spectrogram

(a) 15Hz low-pass filtering oscillogram and true seismic wave figure,

(b) 25Hz low-pass filtering oscillogram and true seismic wave figure,

(c) 15Hz low-pass filtering, 25Hz low-pass filtering and true earthquake record spectrogram,

(d) 10Hz high-pass filtering oscillogram,

(e) 10Hz high-pass filtering spectrogram.

Fig. 6 surpasses random epicentre observational record

(left) without making an uproar data, (right) noisy data.

Window schematic diagram when Fig. 7 blocks.

Fig. 8 observational record low-pass filtering

(left) 15Hz low-pass filtering observational record,

(right) 25Hz low-pass filtering observational record.

Fig. 9 Conventional Time territory full waveform inversion result figure.

Figure 10 low-pass filtering multiple dimensioned full waveform inversion result figure

(a) 15Hz low-pass filtering full waveform inversion result figure,

(b) 25Hz low-pass filtering full waveform inversion result figure,

(c) low-pass filtering final inversion result figure.

Window full waveform inversion result figure when Figure 11 blocks.

Window low-pass filtering multiple dimensioned full waveform inversion result figure when Figure 12 blocks.

Figure 13 Conventional Time territory full waveform inversion result figure (disappearance low-frequency test).

Figure 14 low-pass filtering multiple dimensioned full waveform inversion result figure (disappearance low-frequency test).

(a) 15Hz low-pass filtering full waveform inversion result figure,

(b) 25Hz low-pass filtering full waveform inversion result figure,

(c) low-pass filtering final inversion result figure.

Window full waveform inversion result figure (disappearance low-frequency test) when Figure 15 blocks.

Window low-pass filtering multiple dimensioned full waveform inversion result figure (disappearance low-frequency test) when Figure 16 blocks.

Window+15Hz low-pass filtering full waveform inversion result figure when () blocks a,

(b) 25Hz low-pass filtering full waveform inversion result figure,

Window low-pass filtering final inversion result figure when () blocks c.

Figure 17 Conventional Time territory full waveform inversion result figure (is tested containing Gaussian noise).

Window full waveform inversion result figure (testing containing Gaussian noise) when Figure 18 blocks.

Figure 19 low-pass filtering multiple dimensioned full waveform inversion result figure (is tested containing Gaussian noise).

Window low-pass filtering multiple dimensioned full waveform inversion result figure (testing containing Gaussian noise) when Figure 20 blocks.

Detailed description of the invention

Detailed description further to the present invention with example below in conjunction with the accompanying drawings

The low-pass filtering of window multiple dimensioned full waveform inversion method when blocking, comprises the following steps:

A, program are to have write under MATLAB2013a software frame, it is necessary first to install MATLAB2013a software, And the seismic data processing software Crews tool kit that MATLAB language is write is installed；

First b, the geological data of collection need to carry out pretreatment, and observational record and observation system input computer are carried out Time domain full waveform inversion.Wherein pretreatment includes:

The seismic wave that B1, multiple attenuation, elimination reverberation and compacting ghosting etc. can not just drilled.

B2, underfrequency protection denoising, disappearance seismic channel compensation, relative amplitude preserved processing etc..

B3, definition observation system.

C, estimating source wavelet from observational record, (present invention utilizes Ricker wavelet to replace to this source wavelet for just drilling Source wavelet (Fig. 2) in real data)；

D, setting up linear increment initial model (Fig. 4 left), this rate pattern, as the initial value of inverting, passes through computation model More new direction, realize by initial model gradually to the process that true model is close；

E, according to the principle of least square construct object function:Wherein v represents underground P-wave Degree, d_obsFor the observational record gathered, d_calAt rate pattern by just drilling the forward record obtained.Seeking the gradient of object function During to object function two ends about p wave interval velocity v derivation, obtain gradient expression formula:

$\frac{\∂E\left(v\right)}{\∂v}=-\frac{2}{v^3}\underset{s}{\Σ}\underset{r}{\Σ}\frac{{\∂}^2{p}_f}{\∂t^2}*p_b---\left(2\right)$

Wherein p_bRepresent by with focus f_sThe wave field obtained to underground propagation, p_fRepresenting and just drilling wave field, s represents focus Number, r represents the number of cymoscope.

F, choose window when blocking (Fig. 7), intercept observational record forward part seismic signal, and record and block time window length, should When blocking, the time span of window is as the time span just drilled；

G, dynamic law according to ACOUSTIC WAVE EQUATION, utilize the observation system of input, estimate the source wavelet that obtains and just Beginning rate pattern carries out seismic wave field just drills.This step has only to just drill the geological data in window when blocking, it is not necessary to will be all Time geological data is the most just performed.Storage main story wavefield data and forward record.

Set the low-pass filtering multiple dimensioned full waveform inversion method relevant parameter of window when blocking as requested, big including model Little nz × nx, grid, away from dx, dz, is just drilling time window length It, time sampling interval dt, Ricker wavelet dominant frequency f₀, Butterworth is filtered Cut frequency (the f of ripple device_cut), maximum iteration time itermax of each super random epicentre, update gradient minima gtol, mesh Scalar functions precision tol, the maximum vmax of inversion result and minima vmin；

Super random epicentre encodes: time domain full waveform inversion takes a substantial amount of time during single big gun is just drilled, and has people Propose to utilize and excite multiple focus (super big gun) to replace single gun excitation focus simultaneously, so can shorten when just drilling required Between, but utilize explosive source simultaneously to carry out full waveform inversion and bring crosstalk noise to inversion result again.In order to improve all-wave The computational efficiency of shape inverting is suppressed this software of crosstalk noise simultaneously and is utilized super random epicentre coding strategy.

The Ricker wavelet that the super gun excitation in random phase encoding earth's surface goes out, it is desirable to the positive and negative random assortment of phase place, the most each Single gun excitation time delay in super big gun is different, and time delay maximum is less than 400ms.If random phase encoding information is:

$\mathrm{\Γ}\left(n\right)=\left\{\begin{array}{c}0;n=s^{\′}\\ {(-1)}^i;n=s,i\∈N^{+}\end{array}\right.---\left(3\right)$

Wherein s' ∪ s={1,2 ... N_s, in s, the number of element is fixed, and it represents contained single big gun number in a super big gun Number.In super big gun, single big gun quantity can adjust mixed proportion according to the actual requirements.Element size random assortment in s.i∈N⁺Represent Random positive integer.Super random epicentre coding can be expressed as:

Spf=Dy A Γ T L F f (4)

Wherein Dy represents that (meaning of dynamic coding is exactly the every iteration of supermemory gradient method 10 times to dynamic focus coding function The most again change random phase encoding information, finally can realize earth's surface all standing.The method can increase earth's surface big gun effectively Coverage density, weaken the crosstalk noise impact on full waveform inversion simultaneously, 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 represents random Ricker wavelet dominant frequency coding function, f source function.

H, the observational record in window when blocking and forward record are carried out low-pass filtering simultaneously, it is ensured that by identical low pass filtered Ripple device (Fig. 3) processes data.After filtering, obtain processing observational record and processing forward record；

Butterworth LPF: approximate the system function of wave filter with Butterworth function.Butterworth filters Device is the wave filter in passband according to amplitude-frequency characteristic with smoothness properties definition.

The low pass squared magnitude function of Butterworth filter represents

${\left|H_a\left(j\mathrm{\Ω}\right)\right|}^2=\frac{1}{1+{(\mathrm{\Ω}/{\mathrm{\Ω}}_c)}^{2N}},N=1,2,\mathrm{......}---\left(5\right)$

The principal character of Butterworth filter is as follows:

H1, to all of N,

H2, to all of N,I.e.

H3、|H_a(jΩ)|²It it is the monotonically decreasing function of Ω；

H4、|H_a(jΩ)|²Along with order N increase and closer to ideal low-pass filter；

The amplitude-frequency characteristic of wave filter becomes to become better and better along with the increase of filter order N, at cut-off frequency Ω_cPlace In the case of functional value is always 1/2, there is the value in more frequency band district in passband close to 1；More rapid convergence in stopband In zero.

I, by process observational record and process forward record do difference, obtain residual error data；

J, using residual error data as with focus, act on on focus with anti-pass operator, obtain the residual of the model space Difference anti-pass wave field；

K, main story wave field and residual error anti-pass wave field do cross-correlation and obtain model modification gradient；

The object function of full waveform inversion is commonly defined as the L2 functional just drilling data with observation data residual error, and form is such as Under:

$E\left(v\right)=\frac{1}{2}\left|\right|D_{cal}-D_{obs}|{|}^2---\left(6\right)$

Wherein D_calIt is forward record, D_calRepresent the observational record in actual field.Speed v derivation is obtained by object function Arrive:

$\frac{\∂E\left(v\right)}{\∂v}=\frac{\∂D_{cal}}{\∂v}\[D_{cal}-D_{obs}\]---\left(7\right)$

The present invention is for apparent elaboration theoretical method, it is assumed that density is constant, then have only to speed is carried out inverting. Chang Midu ACOUSTIC WAVE EQUATION is:

$\frac{{\∂}^{2}u}{\∂x^2}+\frac{{\∂}^{2}u}{\∂z^2}-\frac{1}{v^2}\frac{{\∂}^{2}u}{\∂t^2}=f---\left(8\right)$

Wherein u is acoustic wavefield, and f is source function.V represents the acoustic velocity of underground medium.

ACOUSTIC WAVE EQUATION is simplified the form of operator write as:

Lu=s (9)

The most just calculating son:And to velocity derivatives:

Speed v derivation is obtained (Pratt et al., 1998) by equation (9) both sides simultaneously:

$\frac{\∂L}{\∂v}u+L\frac{\∂u}{\∂v}=0---\left(10\right)$

Obtaining wave field to the derivative of underground speed v is:

$\frac{\∂u}{\∂v}=L^{-1}\frac{\∂L}{\∂v}u---\left(11\right)$

Equation (11) is substituted in equation (7) and obtains:

$\frac{\∂E\left(v\right)}{\∂v}=-\frac{2}{v^3}*\frac{{\∂}^{2}u}{\∂t^2}*L^{-1}\[D_{cal}-D_{obs}\]---\left(12\right)$

Wherein L^-1Expression adjoint operator or referred to as anti-pass operator, L^-1[D_cal-D_obs] represent residual error anti-pass wave field.f_s=D_cal- D_obsIt is referred to as with focus.

L, use super-memory gradient method optimized algorithm computation model more new direction, and found suitably by Wolfe convergence criterion Step-length.

Supermemory gradient method: the information of current gradient is carried out school by the information of current gradient and multi-step gradient before Just, to accelerate the convergence rate of object function.Its iteration form is:

v_k+1=v_k+α_kd_k (13)

Wherein v_kRepresent iteration kth portion rate pattern, α_kRepresent step-length, d_kRepresent descent direction.

Supermemory gradient method descent direction can be represented by following form:

$d_k=\left\{\begin{array}{cc}-g_k,& k<mem\\ -{\mathrm{\&beta;}}_k{g}_k-\underset{i=2}{\overset{mem}{\&Sigma;}}{\mathrm{\&beta;}}_{k-i+1}{g}_{k-i+1},& k\&GreaterEqual;mem\end{array}\right.---\left(14\right)$

As k ＞ mem, β_k-i+1∈[0,s_k ⁱ], (i=1 ... mem)

Here

Wherein, 0 ＜ ρ ＜ 1, mem ＞ 0 are positive integers, and mem is we term it memory degree.Different memory degree target letters Number convergence rate and model inversion precision are different, and Memory Gradient quantity is crossed and can be caused gradient direction correction at most excessively, target letter Number declines slack-off, and Memory Gradient quantity deteriorates to conjugate gradient method at least.Obtain through experimental verification, in seismic full-field shape inverting In best memory degree substantially at mem=4,5.Supermemory gradient method has certain advantage, but when gradient memory is too much, calculates Slowing, memory requirements increases.So selecting the advantage that suitably memory degree competence exertion supermemory gradient method is maximum.

M, judgement are enough to meet end condition, if meeting, and the multiple dimensioned full waveform inversion of low-pass filtering of window when output is blocked Inversion method result.If being unsatisfactory for end condition, the initial model that inversion result is circulated as the next one, return Step d.

N, using the output result of previous step as the initial velocity model of Conventional Time territory full waveform inversion, often then carry out Rule time domain full waveform inversion.Iteration exports final inversion result 200 times afterwards.

Embodiment 1

According to exploration requirement, MATLAB2013a is installed under Windows 7 Ultimate system, and installs The seismic data processing software Crews tool kit that MATLAB language is write.

The present invention utilizes Marmousi p wave interval velocity model to test, and Marmousi rate pattern has the structure of complexity Make, trickle layer manages and structure is hidden in deep oil storage, is well suited for carrying out theoretical method test.Original Marmousi rate pattern is relatively Greatly, it is contemplated that hardware memory and calculate time factor, archetype is taken out dilute process by the present invention, to smoke dilute after The low-pass filtering of window multiple dimensioned full waveform inversion method test when Marmousi rate pattern blocks.True model (Fig. 4 Right) and linear increment initial model (Fig. 4 is left).

Model parameter is as follows:

The low-pass filtering of window multiple dimensioned full waveform inversion method test parameter (time domain) when table 1 blocks

Model size	Mesh spacing	Lateral separation	Longitudinal degree of depth	Velocity interval	Dominant frequency of seismic wavelet
						69*384	12.5m	2.4km	0.8625km	1.5-4km/s	22Hz

When blocking, window refutation strategy parameter is as follows:

Observational record utilizes window when blocking process, sets and block time window length It=0.2:0.1:2.1s, i.e. select Window 0-0.2s when selecting initial blocking, when then every iteration is blocked for 10 times, the time span of window increases 0.1s, and initial 0 moment is constant (Fig. 7), until when blocking the total length of window be 2.1s.This process totally 20 windows when blocking, full waveform inversion iteration 200 times altogether. Just drilling and having only to just drill the time span to during correspondence, so can largely save and just drill the required calculating time.Cut Time disconnected, window refutation strategy is by shallow-layer to deep layer, successively carries out inverting, greatly reduces the parameter of inverting each time, be prevented effectively from The generation of cycle slip phenomenon, enhances the stability of program simultaneously.

Low-pass filtering Multi-scale inversion policing parameter is as follows:

In order to be absorbed in the problem such as local minimum and cycle slip during overcoming full waveform inversion, and also to improve into picture element Amount, carries out full waveform inversion with low-pass filtering 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:

I, 0-15Hz: low-frequency range.Carry out inverting in low-frequency range and can obtain the macroscopic information of subsurface model, in low-frequency range Inversion result, so can gradually rate pattern accurately inverting out as the initial model of next frequency band inverting. But often lacking low-frequency information in real data, land geological data disappearance 10Hz data below, low-pass filtering simply exists There is competence exertion effect in the case of low-frequency data, the effect of the multiple dimensioned full waveform inversion of low-pass filtering when lacking low-frequency data Fruit is not very well (Figure 13).So the present invention combines window refutation strategy when blocking in low-frequency range and eliminates because of disappearance low frequency The cycle slip phenomenon caused.

II, 0-25Hz: Mid Frequency.This frequency range primarily serves the effect of an intermediate transition, and the inversion result of low-frequency range is made Initial velocity model for Mid Frequency.

III, 0-end: full frequency band.The initial velocity model of the inversion result Whole frequency band of midband, this step i.e. can obtain The inversion result of time domain Full wave shape.

Based on above parameter, under the test environment shown in table 2, carry out full waveform inversion.

Table 2 computing power test environment

Table 3 full waveform inversion method contrasts

As can be drawn from Table 3, the low-pass filtering of window when blocking that Conventional Time territory full waveform inversion and the present invention propose is many Yardstick full waveform inversion has obvious difference (model error calculates time aspect).The meter of Conventional Time territory full waveform inversion Evaluation time is 101.3 minutes, and when blocking, the multiple dimensioned full waveform inversion of the low-pass filtering of window only needs just can complete for 43.6 minutes Whole calculating process.Combine inversion result figure (Fig. 9,10,11,12) from the angle of model error and can be seen that the low of window when blocking Pass filter multiple dimensioned full waveform inversion inversion result precision is higher, and the rate results that inverting obtains is more accurate, and Conventional Time The result of territory full waveform inversion even occurs in that mistake, falls far short with true model.

Based on above parameter, under the test environment shown in table 2, carry out various factors test:

The low-pass filtering multiple dimensioned full waveform inversion method test result of window when blocking under table 4 different affecting factors

Note: in table 3 and table 4, model inversion normalization error calculation formula is:Wherein vtrue generation Table true velocity model, vk represents inversion result.

From table 4 and Figure 13, it can be seen that full waveform inversion is had a serious impact by disappearance low frequency in 14,15.Work as earthquake The when that data lacking low frequency, Conventional Time territory full waveform inversion result (Figure 13), the multiple dimensioned full waveform inversion of low-pass filtering Result (Figure 14) and when blocking window full waveform inversion result (Figure 15) all occur in that manifest error, model inversion error is the most inclined High.Simultaneously as the factor of disappearance low frequency causes model modification gradient calculation result incorrect so that carrying out, Wolfe is linear Waste the substantial amounts of time during step-size in search, this explains the most identical iterations, when lacking low frequency Wait the calculating time of calculating time full frequency band the to be considerably longer than full waveform inversion of full waveform inversion.In order to preferably solve because of The inverting Problem-Error that disappearance low-frequency component brings, window refutation strategy and low-pass filtering multi-scale strategy knot when the present invention will block Closing and use, as can be seen from Figure 16, this combination strategy makes the result of full waveform inversion obtain obvious improvement.Model error Also illustrate that, the inversion result that this strategy obtains is more nearly with true model.

From table 4 and Figure 17, it can be seen that when in geological data containing Gaussian noise when, inverting is tied in 18,19,20 Fruit becomes unintelligible.Conventional Time territory full waveform inversion (Figure 17) and when blocking window full waveform inversion (Figure 18) all occur in that bright Aobvious mistake, although window refutation strategy accelerates inverting process to a certain extent when illustrating to block, but its anti-noise ability is relatively Weak.The multiple dimensioned full waveform inversion of low-pass filtering (Figure 19) has obvious advantage compared with the above two, it was demonstrated that low-pass filtering is multiple dimensioned Strategy has stronger anti-noise ability.When will block, window refutation strategy and low-pass filtering Multi-scale inversion strategy are used in combination with (Figure 20), the process of inverting is not only accelerated, and the defect that when compensate for blocking, window anti-noise ability is weak.

Fig. 1 is the flow chart of whole refutation process, can be seen that first window is anti-when low-frequency range utilizes and blocks from flow chart Drill strategy and carry out full waveform inversion, when the most again low-frequency range being blocked window inversion result as initial model for next step intermediate frequency Section inverting, finally by the initial model of the inversion result of Mid Frequency Conventional Time territory full waveform inversion the most.When this blocks, window is low Pass filter multiple dimensioned full waveform inversion strategy, is better than often at aspects such as inversion accuracy, computational efficiency, low frequency dependence, anti-noise abilities Rule time domain full waveform inversion.

Claims (1)

1. the multiple dimensioned full waveform inversion method of low-pass filtering of window when blocking, it is characterised in that comprise the following steps:

A, installation MATLAB2013a software, and the seismic data processing software Crews tool kit that MATLAB language is write is installed；

B, the geological data gathered is carried out pretreatment, and observational record and observation system are inputted computer, carry out time domain Full waveform inversion；

C, from observational record estimate source wavelet, the source wavelet that this estimation obtains is for forward simulation；

D, setting up linear increment initial model, this rate pattern is as the initial value of inverting, by the more new direction of computation model, Realize the process gradually approached by initial model to true model；

E, according to principle of least square method build object function:

Wherein v represents p wave interval velocity, d_obsFor the observational record gathered, d_calAt rate pattern by just drilling the forward record obtained, In the gradient procedure seeking object function, the p wave interval velocity v derivation to object function two ends, obtain gradient expression formula:

$\frac{\&part;E\left(v\right)}{\&part;v}=-\frac{2}{v^3}\underset{s}{\&Sigma;}\underset{r}{\&Sigma;}\frac{{\&part;}^2{p}_f}{\&part;t^2}*p_b$

Wherein p_bRepresent by with focus f_sThe wave field obtained to underground propagation, p_fRepresenting forward simulation wave field, s represents focus Number, r represents the number of cymoscope；

F, choose the time span of window when blocking, intercept observational record forward part seismic signal, and record and block time window length, should When blocking, the time span of window is as the time span of forward simulation；

G, dynamic law according to ACOUSTIC WAVE EQUATION, utilize the observation system of input, estimate the source wavelet that obtains and initial speed Degree model carries out seismic wave field just drills, storage main story wavefield data and forward record；

Set the low-pass filtering multiple dimensioned full waveform inversion method relevant parameter of window when blocking 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 f₀, Butterworth is filtered Cut frequency f of ripple device_cut, the maximum iteration time iter max of each super random epicentre, update gradient minima gtol, mesh Scalar functions precision tol, the maximum v max of inversion result and minima v min；

H, the observational record in window when blocking and forward record are carried out low-pass filtering simultaneously, it is ensured that with identical low pass filter Process data, after filtering, obtain processing observational record and processing forward record；

I, by process observational record and process forward record do difference, obtain residual error data；

J, using residual error data as with focus, act on on focus with anti-pass operator, the residual error obtaining the model space is anti- Pass wave field；

K, main story wave field and anti-pass residual error wave field do cross-correlation, obtain model modification gradient；

L, use super-memory gradient method optimized algorithm computation model more new direction, and find suitable step-length by Wolfe convergence criterion；

M, judge whether to meet end condition, if meeting, the multiple dimensioned full waveform inversion method of low-pass filtering of window when output is blocked Inversion result；If being unsatisfactory for end condition, the initial model that inversion result is circulated as the next one, return Step d；

N, using m step output result as the initial velocity model of Conventional Time territory full waveform inversion, then carry out Conventional Time Territory full waveform inversion, iteration exports final inversion result 200 times afterwards.