CN105676277A - Full waveform joint inversion method for improving high steep structure velocity inversion efficiency - Google Patents

Abstract

The invention discloses a full waveform joint inversion method for improving high steep structure velocity inversion efficiency, and belongs to the petroleum geophysical exploration field. The method comprises the steps: establishing an observation system by inputting the initial velocity field, the shot record and the focus wavelet; utilizing traditional reverse time migration to obtain the imaging result which is taken as a reflection coefficient model, and applying linear forward modeling and linear wave field back propagation; calculating the gradient directions of two prism waveforms and adding the gradient directions together to calculate the inverse gradient direction of the prism waveform, using the linear search method or a parabolic fitting method to calculate the updated step, and using the calculated gradient direction and step to update the velocity; using the conventional full waveform inversion method to update the velocity again, and determining whether or not to satisfy the error condition; if not, using the updated velocity as the input velocity, and updating the velocity again; and if so, outputting the inversion result. The full waveform joint inversion method for improving high steep structure velocity inversion efficiency reduces the dependency about whether the high steep structure information is missing in the initial model, and improves the inversion accuracy and efficiency for the high steep structure.

Description

A kind of Full wave shape joint inversion method improving high-dip structure velocity inversion efficiency

Technical field

The invention belongs to field of petroleum geophysical exploration, be specifically related to a kind of Full wave shape joint inversion method improving high-dip structure velocity inversion efficiency.

Background technology

Conventional primary reflection exploration, is limited by observation system, to the steep dip subsurface structure illumination deficiencies such as tomography, salt dome flank, high-dip structure accurately image difficulty. The prism ripple information that high-dip structure produces, it is possible to be effectively increased the lighting level of high-dip structure, improve high-dip structure imaging effect. Therefore make full use of prism ripple information the imaging capability improving high-dip structure is very important.

Due to full waveform inversion high accuracy and high-resolution feature so that it is become a kind of powerful of velocity modeling, be increasingly becoming the focus of research. Full waveform inversion is the process of a nonlinear data fitting, carrys out undated parameter model by reducing the difference between observation data and prediction data, and this process repeats down in an iterative manner, until data difference is sufficiently small. Full waveform inversion amount of calculation is very big, one of computational efficiency problem research emphasis becoming full waveform inversion.

Summary of the invention

For the above-mentioned technical problem existed in prior art, the present invention proposes a kind of Full wave shape joint inversion method improving high-dip structure velocity inversion efficiency, the dependency whether high-dip structure information in initial model is lacked can be reduced, improve the inversion accuracy to high-dip structure and efficiency.

To achieve these goals, the present invention adopts the following technical scheme that

A kind of Full wave shape joint inversion method improving high-dip structure velocity inversion efficiency, comprises the steps:

Step 1: the input speed of initial velocity field, field inspection big gun record and source wavelet, and set up observation system;

Step 2: use the imaging formula that following formula represents to obtain imaging results:

$I=\underset{x_s}{\Σ}\underset{t}{\Σ}uw/u^2$

Wherein, I represents imaging results, and u and w represents main story wave field and anti-pass wave field, w=L respectively^*(R^*(p_obs)), L^*For anti-pass operator, R^*Represent and will be defined in data residual error spatial spread on cymoscope to the whole model space, p_obsFor field inspection big gun record;

Step 3: using the imaging results in step 2 as reflection coefficient;

Step 4: apply linear forward simulation and linear wave field anti-pass, calculates the gradient direction of two kinds of prism ripple prism1 and prism2 respectively;

$g{\left(v\right)}_{prism1}=\frac{2}{v^3}\underset{x_s}{\Σ}\underset{t}{\Σ}\frac{{\∂}^{2}B\left(u\right)}{\∂t^2}{w}^{\′}$

$g{\left(v\right)}_{prism2}=\frac{2}{v^3}\underset{x_s}{\Σ}\underset{t}{\Σ}\frac{{\∂}^{2}u}{\∂t^2}{B}^{*}\left(w^{\′}\right)$

Wherein, g (v)_prism1With g (v)_prism2Representing the gradient direction of prism1 and the gradient direction of prism2 respectively, v is the speed of velocity field, and w' represents the anti-pass wave field of field inspection big gun record and the earthquake record residual difference of numerical simulation, x_sFor hypocentral location, t is the time;

Step 5: two gradient directions of gained in step 4 are added the gradient direction that can try to achieve prism waveform inversion;

$g{\left(v\right)}_{prism}=\frac{2}{v^3}\underset{x_s}{\Σ}\underset{t}{\Σ}\frac{{\∂}^{2}B\left(u\right)}{\∂t^2}{w}^{\′}+\frac{2}{v^3}\underset{x_s}{\Σ}\underset{t}{\Σ}\frac{{\∂}^{2}u}{\∂t^2}{B}^{*}\left(w^{\′}\right)$

Wherein, g (v)_prismRepresent the gradient direction of prism waveform inversion;

Step 6: ask for renewal step-length α with linear search method or parabolic approximating method;

Step 7: utilize the gradient direction asked for and update step-length renewal speed;

v_k=v_k-1+α_kg(v_k-1)_prism

Wherein, subscript k represents the number of times of iteration;

Step 8: utilize conventional full wave shape inversion method renewal speed;

v_k+1=v_k+α_kg(v_k)

Wherein g (v_k) for the gradient direction of kth time iteration conventional full wave shape inverting;

Step 9: judge v_k+1With v_kDifference whether meet error condition;

If: judged result is v_k+1With v_kDifference be unsatisfactory for error condition, then perform step 2;

Or judged result is v_k+1With v_kDifference meet error condition, then perform step 10;

Step 10: the rate results of output inverting.

Preferably, in step 4, specifically include

Step 4.1: objective function:

$E\left(v\right)=\frac{1}{2}\underset{x_r}{\Σ}\underset{x_s}{\Σ}\underset{t}{\Σ}{\left(Ru\right(t,x_r,x_s)-p{\left(t,x_r,x_s\right)}_{obs})}^2---\left(1\right);$

Wherein, Ru (t, x_r,x_s) and p (t, x_r,x_s)_obsThe respectively earthquake record of numerical simulation and field inspection big gun record, x_sAnd x_rRepresenting shot point and geophone station respectively, t is the time;

Step 4.2: by object function variation, obtains variation expression formula:

$\begin{array}{c}\mathrm{\δ}E\left(v\right)=\frac{1}{2}\mathrm{\δ}\[(Ru-p_{obs})\·(Ru-p_{obs})\]\\ =\[\mathrm{\δ}(Ru-p_{obs})\]\·(Ru-p_{obs})\\ =\left(R\mathrm{\δ}u\right)\·(Ru-p_{obs})\end{array}---\left(2\right);$

Step 4.3: definition two dimension Acoustic Wave-equation:

$\frac{1}{v^2}\frac{{\∂}^{2}u}{\∂t^2}-\▿\·\▿u=s(t,x_s)---\left(3\right);$

Equation (3) is deformed:

$\frac{1}{{(v+\mathrm{\δ}v)}^2}\frac{{\∂}^2(u+\mathrm{\δ}u)}{\∂t^2}-\▿\·\▿(u+\mathrm{\δ}u)=s(t,x_s)---\left(4\right);$

Wherein, s (t, x_s) represent focus item,Represent Laplace operator;

Equation (4) is carried out Taylor expansion:

$(\frac{1}{v^2}-\frac{2\mathrm{\δ}v}{v^3})\frac{{\∂}^2(u+\mathrm{\δ}u)}{\∂t^2}-\▿\·\▿(u+\mathrm{\δ}u)=s(t,x_s)---\left(5\right);$

Launch to omit higher order term to formula (5):

$\frac{1}{v^2}\frac{{\∂}^{2}u}{\∂t^2}-\▿\·\▿u+\frac{1}{v^2}\frac{{\∂}^2\mathrm{\δ}u}{\∂t^2}-\▿\·\▿\mathrm{\δ}u-\frac{2\mathrm{\δ}v}{v^3}\frac{{\∂}^{2}u}{\∂t^2}=s(t,x_s)---\left(6\right);$

Equation (4) and equation (6) are subtracted each other, obtain:

$\frac{1}{v^2}\frac{{\∂}^2\mathrm{\δ}u}{\∂t^2}-\▿\·\▿\mathrm{\δ}u=\frac{2\mathrm{\δ}v}{v^3}\frac{{\∂}^{2}u}{\∂t^2}---\left(7\right);$

Can obtain further:

$\mathrm{\δ}u=L\left(\frac{2\mathrm{\δ}v}{v^3}\frac{{\∂}^{2}u}{\∂t^2}\right)---\left(8\right);$

Wherein, L represents positive communication process;

Step 4.4: equation (8) is substituted into variation expression formula (2), can obtain:

$\begin{array}{c}\mathrm{\δ}E\left(v\right)=\left(R\right(L\left(\frac{2\mathrm{\δ}v}{v^3}\frac{{\∂}^{2}u}{\∂t^2}\right)\left)\right)\·(Ru-p_{obs})\\ =\left(L\right(\frac{2\mathrm{\δ}v}{v^3}\frac{{\∂}^{2}u}{\∂t^2}\left)\right)\·R^{*}(Ru-p_{obs})\\ =\left(\frac{2\mathrm{\δ}v}{v^3}\frac{{\∂}^{2}u}{\∂t^2}\right)\·L^{*}\left(R^{*}\right(Ru-p_{obs}\left)\right)\end{array}---\left(9\right);$

Wherein, R^*(Ru-p_obs) represent data residual error spatial spread to whole space, L^*(R^*(Ru-p_obs)) represent that the residue wave field inverse time propagates;

Step 4.5: make the data residual error be:

Δ p=Ru-p_obs(10);

Step 4.6: by main story wave field u and data residual delta p=Ru-p_obsIt is decomposed into following two parts:

u(t,x_r,x_s)=u₁(t,x_r,x_s)+u₂(t,x_r,x_s) (11);

Δp(t,x_r,x_s)=Δ p₁(t,x_r,x_s)+Δp₂(t,x_r,x_s) (12);

Wherein, u₁(t,x_r,x_s) and u₂(t,x_r,x_s) represent first order reflection wave field and second order prism wave field respectively; Δ p₁(t,x_r,x_s) and Δ p₂(t,x_r,x_s) represent the data residual error of data residual sum second order prism wave field of first order reflection wave field respectively;

Step 4.7: equation (10), (11), (12) are substituted into equation (9):

Step 4.8: ask for the object function gradient to rate pattern:

$\begin{array}{c}g\left(v\right)=\frac{\mathrm{\δ}E\left(v\right)}{\mathrm{\δ}v}=\frac{2}{v^3}\underset{x_s}{\mathrm{\Σ}}\underset{t}{\mathrm{\Σ}}\frac{{\∂}^2{u}_1}{\∂t^2}{L}^{*}\left(R^{*}\right({\mathrm{\Δp}}_1)\\ +\frac{2}{v^3}\underset{x_s}{\mathrm{\Σ}}\underset{t}{\mathrm{\Σ}}\frac{{\∂}^2{u}_1}{\∂t^2}{L}^{*}\left(R^{*}\right({\mathrm{\Δp}}_2)+\frac{2}{v^3}\underset{x_s}{\mathrm{\Σ}}\underset{t}{\mathrm{\Σ}}\frac{{\∂}^2{u}_2}{\∂t^2}{L}^{*}\left(R^{*}\right({\mathrm{\Δp}}_1)\end{array}---\left(14\right);$

Make g (v)_reflect, g (v)_prism1With g (v)_prism2Represent the gradient of primary reflection respectively, the gradient of prism1 and the gradient of prism2, then

$g{\left(v\right)}_{reflect}=\frac{2}{v^3}\underset{x_s}{\Σ}\underset{t}{\Σ}\frac{{\∂}^2{u}_1}{\∂t^2}{L}^{*}\left(R^{*}\right({\mathrm{\Δp}}_1)---\left(15\right);$

$g{\left(v\right)}_{prism1}=\frac{2}{v^3}\underset{x_s}{\Σ}\underset{t}{\Σ}\frac{{\∂}^2{u}_2}{\∂t^2}{L}^{*}\left(R^{*}\right({\mathrm{\Δp}}_1)---\left(16\right);$

$g{\left(v\right)}_{prism2}=\frac{2}{v^3}\underset{x_s}{\Σ}\underset{t}{\Σ}\frac{{\∂}^2{u}_1}{\∂t^2}{L}^{*}\left(R^{*}\right({\mathrm{\Δp}}_2)---\left(17\right);$

By first order reflection wave field u₁(t,x_r,x_s) represented by main story wave field u, second order prism wave field u₂(t,x_r,x_s) represented by linear forward simulation B (u), the L of the data residual error of first order reflection wave field^*(R^*(Δp₁) represented by anti-pass wave field w', the L of the data residual error of second order prism wave field^*(R^*(Δp₂) by the linear wavelength anti-pass B of anti-pass wave field w'^*(w') represent, do corresponding variable replacement, then

$g{\left(v\right)}_{prism1}=\frac{2}{v^3}\underset{x_s}{\Σ}\underset{t}{\Σ}\frac{{\∂}^{2}B\left(u\right)}{\∂t^2}{w}^{\′}---\left(18\right);$

$g{\left(v\right)}_{prism2}=\frac{2}{v^3}\underset{x_s}{\Σ}\underset{t}{\Σ}\frac{{\∂}^{2}u}{\∂t^2}{B}^{*}\left(w^{\′}\right)---\left(19\right).$

The Advantageous Effects that the present invention brings:

The present invention proposes a kind of Full wave shape joint inversion method improving high-dip structure velocity inversion efficiency, compared with prior art, a kind of Full wave shape joint inversion method improving high-dip structure velocity inversion efficiency, prism ripple information can be made full use of, the salt side wing carries out speed accurately update, therefore, in describing the vertical boundary of salt dome flank, the present invention plays very big effect;The present invention is low to the more conventional full waveform inversion method of degree of dependence of initial model, can when initial model high steep information is less even lack relatively accurately inverting there is the model of high-dip structure, improve the inversion accuracy to high-dip structure and efficiency, provide migration velocity field accurately for high accuracy formation method.

Accompanying drawing explanation

Fig. 1 is the FB(flow block) of a kind of Full wave shape joint inversion method improving high-dip structure velocity inversion efficiency of the present invention.

Fig. 2 is the true velocity model that the present invention uses.

Fig. 3 is the initial velocity model that the present invention uses.

Fig. 4 is the inversion speed of the 10th iteration using the present invention to obtain.

Fig. 5 is the inversion speed of the 10th iteration using conventional full wave shape inversion method to obtain.

Fig. 6 is the final inversion speed using the present invention to obtain.

Fig. 7 is the inversion speed of the identical iterations using conventional full wave shape inversion method to obtain.

Fig. 8 is inversion speed curve.

Detailed description of the invention

Below in conjunction with accompanying drawing and detailed description of the invention, the present invention is described in further detail:

The flow chart (as shown in Figure 1) of a kind of Full wave shape joint inversion method utilizing prism ripple information to improve high-dip structure velocity inversion efficiency, comprises the steps:

Input initial velocity field, big gun record and source wavelet, and set up observation system;

Use tradition reverse-time migration to obtain imaging results, use reverse-time migration imaging results as reflectivity model, apply linear forward simulation and linear wave field anti-pass;

Calculate the gradient direction of two kinds of prism ripple prism1 and prism2, and this two parts gradient is added the gradient direction that can try to achieve prism waveform inversion, ask for renewal step-length with linear search method or parabolic approximating method, utilize the gradient direction and step-length renewal speed asked for;

Utilize conventional full wave shape inversion method renewal speed again, it may be judged whether meet error condition, if the error condition of being unsatisfactory for, utilize the speed updated again to carry out speed renewal as input speed, if meeting error condition, the rate results of output inverting.

Concretely comprise the following steps:

Step 1: the input speed of initial velocity field, field inspection big gun record and source wavelet, and set up observation system;

Step 2: use the imaging formula that following formula represents to obtain imaging results:

$I=\underset{x_s}{\Σ}\underset{t}{\Σ}uw/u^2$

Wherein I represents imaging results, and u and w represents main story wave field and anti-pass wave field, w=L respectively^*(R^*(p_obs)), L^*For anti-pass operator, R^*Represent and will be defined in data residual error spatial spread on cymoscope to the whole model space, p_obsFor field inspection big gun record;

Step 3: using reverse-time migration imaging results as reflectivity model, apply linear forward simulation and linear wave field anti-pass, the two process is represented by B and B* respectively;

Step 4: calculate the gradient direction of two kinds of prism ripple prism1 and prism2;

Characteristic according to prism ripple can be classified as two types, the first is called prism1, be characterized as being seismic wave first occurs reflection at lenticular body flank, reflection occur again subsequently and propagate to ground and be detected device and receive in deposition interface, is therefore referred to as IF prism wave reflection; The second is called prism2, is that seismic wave first reflects in deposition interface subsequently in the reflection of lenticular body flank and is propagated back to ground receiver formation;

Utilize full waveform inversion to obtain formation velocity, be generally the minimum solving object function by full waveform inversion problem definition, adopt the object function of the L2 norm form shown in following formula:

$E\left(v\right)=\frac{1}{2}\underset{x_r}{\Σ}\underset{x_s}{\Σ}\underset{t}{\Σ}{\left(Ru\right(t,x_r,x_s)-p{\left(t,x_r,x_s\right)}_{obs})}^2---\left(1\right)$

Wherein, Ru (t, x_r,x_s) and p (t, x_r,x_s)_obsThe respectively earthquake record of numerical simulation and field inspection numerical value, x_sAnd x_rRepresenting shot point and geophone station respectively, t is the time;

The variation of object function can be expressed as:

$\begin{array}{c}\mathrm{\δ}E\left(v\right)=\frac{1}{2}\mathrm{\δ}\[(Ru-p_{obs})\·(Ru-p_{obs})\]\\ =\[\mathrm{\δ}(Ru-p_{obs})\]\·(Ru-p_{obs})\\ =\left(R\mathrm{\δ}u\right)\·(Ru-p_{obs})\end{array}---\left(2\right)$

A given velocity disturbance δ v, it can cause the disturbance δ u of a main story wave field, this relation v+ δ v → u+ δ u to meet two dimension Acoustic Wave-equation (3), brings equation into and can obtain (4);

$\frac{1}{v^2}\frac{{\∂}^{2}u}{\∂t^2}-\▿\·\▿u=s(t,x_s)---\left(3\right)$

$\frac{1}{{(v+\mathrm{\δ}v)}^2}\frac{{\∂}^2(u+\mathrm{\δ}u)}{\∂t^2}-\▿\·\▿(u+\mathrm{\δ}u)=s(t,x_s)---\left(4\right)$

Wherein, s (t, x_s) represent focus item,Represent Laplace operator, equation (4) carried out Taylor and launches approximate:

$(\frac{1}{v^2}-\frac{2\mathrm{\δ}v}{v^3})\frac{{\∂}^2(u+\mathrm{\δ}u)}{\∂t^2}-\▿\·\▿(u+\mathrm{\δ}u)=s(t,x_s)---\left(5\right)$

Launch to omit higher order term to formula (5):

$\frac{1}{v^2}\frac{{\∂}^{2}u}{\∂t^2}-\▿\·\▿u+\frac{1}{v^2}\frac{{\∂}^2\mathrm{\δ}u}{\∂t^2}-\▿\·\▿\mathrm{\δ}u-\frac{2\mathrm{\δ}v}{v^3}\frac{{\∂}^{2}u}{\∂t^2}=s(t,x_s)---\left(6\right)$

Equation (4) and equation (6) are subtracted each other and can obtain:

$\frac{1}{v^2}\frac{{\∂}^2\mathrm{\δ}u}{\∂t^2}-\▿\·\▿\mathrm{\δ}u=\frac{2\mathrm{\δ}v}{v^3}\frac{{\∂}^{2}u}{\∂t^2}---\left(7\right)$

Can obtain further:

$\mathrm{\δ}u=L\left(\frac{2\mathrm{\δ}v}{v^3}\frac{{\∂}^{2}u}{\∂t^2}\right)---\left(8\right)$

Wherein, L represents positive communication process, and the variation expression formula (2) of object function can be rewritten as:

$\begin{array}{c}\mathrm{\δ}E\left(v\right)=\left(R\right(L\left(\frac{2\mathrm{\δ}v}{v^3}\frac{{\∂}^{2}u}{\∂t^2}\right)\left)\right)\·(Ru-p_{obs})\\ =\left(L\right(\frac{2\mathrm{\δ}v}{v^3}\frac{{\∂}^{2}u}{\∂t^2}\left)\right)\·R^{*}(Ru-p_{obs})\\ =\left(\frac{2\mathrm{\δ}v}{v^3}\frac{{\∂}^{2}u}{\∂t^2}\right)\·L^{*}\left(R^{*}\right(Ru-p_{obs}\left)\right)\end{array}---\left(9\right)$

Wherein, R^*(Ru-p_obs) represent and will be defined in data residual error spatial spread on cymoscope to whole space, L^*(R^*(Ru-p_obs)) represent that the residue wave field inverse time propagates.

In order to simplify representation, the data residual error is made to be:

Δ p=Ru-p_obs(10)

In order to prism ripple is easily separated, main story wave field u and data residual delta p=Ru-p_obsCan be analyzed to following two parts:

u(t,x_r,x_s)=u₁(t,x_r,x_s)+u₂(t,x_r,x_s)(11)

Δp(t,x_r,x_s)=Δ p₁(t,x_r,x_s)+Δp₂(t,x_r,x_s)(12)

Wherein, u₁(t,x_r,x_s) and u₂(t,x_r,x_s) represent first order reflection wave field and secondary prism wave field respectively; Δ p₁(t,x_r,x_s) and Δ p₂(t,x_r,x_s) represent the data residual error of data residual sum second order prism wave field of first order reflection wave field respectively.

Equation (9) can be further rewritten as:

The gradient of rate pattern is by cost functional:

$\begin{array}{c}g\left(v\right)=\frac{\mathrm{\δ}E\left(v\right)}{\mathrm{\δ}v}=\frac{2}{v^3}\underset{x_s}{\mathrm{\Σ}}\underset{t}{\mathrm{\Σ}}\frac{{\∂}^2{u}_1}{\∂t^2}{L}^{*}\left(R^{*}\right({\mathrm{\Δp}}_1)\\ +\frac{2}{v^3}\underset{x_s}{\mathrm{\Σ}}\underset{t}{\mathrm{\Σ}}\frac{{\∂}^2{u}_1}{\∂t^2}{L}^{*}\left(R^{*}\right({\mathrm{\Δp}}_2)+\frac{2}{v^3}\underset{x_s}{\mathrm{\Σ}}\underset{t}{\mathrm{\Σ}}\frac{{\∂}^2{u}_2}{\∂t^2}{L}^{*}\left(R^{*}\right({\mathrm{\Δp}}_1)\end{array}---\left(14\right)$

Make g (v)_reflect, g (v)_prism1With g (v)_prism2Represent the gradient of primary reflection respectively, the gradient of prism1 and the gradient of prism2, then

$g{\left(v\right)}_{reflect}=\frac{2}{v^3}\underset{x_s}{\Σ}\underset{t}{\Σ}\frac{{\∂}^2{u}_1}{\∂t^2}{L}^{*}\left(R^{*}\right({\mathrm{\Δp}}_1)---\left(15\right)$

$g{\left(v\right)}_{prism1}=\frac{2}{v^3}\underset{x_s}{\Σ}\underset{t}{\Σ}\frac{{\∂}^2{u}_2}{\∂t^2}{L}^{*}\left(R^{*}\right({\mathrm{\Δp}}_1)---\left(16\right)$

$g{\left(v\right)}_{prism2}=\frac{2}{v^3}\underset{x_s}{\Σ}\underset{t}{\Σ}\frac{{\∂}^2{u}_1}{\∂t^2}{L}^{*}\left(R^{*}\right({\mathrm{\Δp}}_2)---\left(17\right);$

First order reflection wave field u₁(t,x_r,x_s) can be represented by main story wave field u, second order prism wave field u₂(t,x_r,x_s) can linear forward simulation B (u) try to achieve, the L of the data residual error of first order reflection wave field^*(R^*(Δp₁) can by the anti-pass wave field w' of field inspection big gun record Yu the earthquake record residual difference of numerical simulation, the L of the data residual error of second order prism wave field^*(R^*(Δp₂) can by the linear wavelength anti-pass B of w'^*(w') try to achieve, do corresponding variable replacement, then

$g{\left(v\right)}_{prism1}=\frac{2}{v^3}\underset{x_s}{\Σ}\underset{t}{\Σ}\frac{{\∂}^{2}B\left(u\right)}{\∂t^2}{w}^{\′}---\left(18\right);$

$g{\left(v\right)}_{prism2}=\frac{2}{v^3}\underset{x_s}{\Σ}\underset{t}{\Σ}\frac{{\∂}^{2}u}{\∂t^2}{B}^{*}\left(w^{\′}\right)---\left(19\right).$

Can be directly separating and utilize prism shape information by asking for of two kinds of prism waveform gradient directions, without the big gun record isolated by the big gun record of input containing only prism ripple;

Step 5: two gradient directions of gained in step 4 are added the gradient direction that can try to achieve prism waveform inversion;

$g{\left(v\right)}_{prism}=\frac{2}{v^3}\underset{x_s}{\Σ}\underset{t}{\Σ}\frac{{\∂}^{2}B\left(u\right)}{\∂t^2}{w}^{\′}+\frac{2}{v^3}\underset{x_s}{\Σ}\underset{t}{\Σ}\frac{{\∂}^{2}u}{\∂t^2}{B}^{*}\left(w^{\′}\right)---\left(20\right)$

Wherein g (v)_prismRepresent the gradient direction of prism waveform inversion;

Step 6: ask for renewal step-length α with linear search method or parabolic approximating method;

Step 7: utilize the gradient direction asked for and update step-length renewal speed;

v_k=v_k-1+α_kg(v_k-1)_prism(21)

Wherein subscript k represents the number of times of iteration;

Step 8: utilize conventional full wave shape inversion method renewal speed;

v_k+1=v_k+α_kg(v_k)(22)

Wherein g (v_k) for the gradient direction of kth time iteration conventional full wave shape inverting;

Step 9: judge v_k+1With v_kDifference whether meet error condition;

If: judged result is v_k+1With v_kDifference be unsatisfactory for error condition, then perform step 2;

Or judged result is v_k+1With v_kDifference meet error condition, then perform step 10;

Step 10: the rate results of output inverting.

The present invention utilize prism ripple information improve high-dip structure velocity inversion efficiency Full wave shape joint inversion method the more conventional full waveform inversion method of degree of dependence of initial model is low, can when initial model high steep information is less even lack relatively accurately inverting there is the model of high-dip structure, improve the inversion accuracy to high-dip structure and efficiency.

Application experiment

The present invention utilizes prism ripple information to improve the Full wave shape joint inversion method of high-dip structure velocity inversion efficiency, is applied to international standard Marmousi model data, achieves and desirably calculate effect. Fig. 2 is Marmousi true velocity model;Fig. 3 is Marmousi initial velocity model; Use tradition reverse-time migration to obtain imaging results, use reverse-time migration imaging results as reflectivity model, apply linear forward simulation and linear wave field anti-pass; Calculate the gradient direction of prism waveform inversion, ask for renewal step-length with linear search method or parabolic approximating method, utilize the gradient direction and step-length renewal speed asked for; Utilize conventional full wave shape inverting (FWI) method renewal speed again, judge whether to meet error condition, if the error condition of being unsatisfactory for, utilize the speed updated again to carry out speed renewal as input speed, obtain, by 10 iteration, the inversion result (as shown in Figure 4) utilizing the Full wave shape joint inversion method renewal speed that prism ripple information improves high-dip structure velocity inversion efficiency to obtain the 10th iteration; Iterate until meeting error condition, the rate results (as shown in Figure 6) of output inverting. Compared to the inversion result (as shown in Figure 5) of tradition the 10th iteration of full waveform inversion, improve, by prism ripple information, the inversion result (as shown in Figure 4) of the 10th iteration that the Full wave shape joint inversion method of high-dip structure velocity inversion efficiency obtains the inversion speed of the high-dip structure such as major fault is considerably more rapid; Inversion result (as shown in Figure 7) compared to the tradition identical iterations of full waveform inversion, the overall inversion result of the inversion result (as shown in Figure 6) that obtains of Full wave shape joint inversion method of high-dip structure velocity inversion efficiency is improved also more preferably by prism ripple information, giving more high frequency detail accurately, relatively initial model (as shown in Figure 3) has had bigger improvement. Fig. 8 is inversion speed curve, and the inversion speed curve of gained of the present invention is closer to true velocity curve.

A kind of Full wave shape joint inversion method utilizing prism ripple information to improve high-dip structure velocity inversion efficiency is the invention provides for this, develop and a kind of can reduce the dependency whether high-dip structure information in initial model is lacked, improve the velocity modeling method of the inversion accuracy to high-dip structure and efficiency, provide migration velocity field accurately for high accuracy formation method.

Certainly, described above is not limitation of the present invention, and the present invention is also not limited to the example above, the change made in the essential scope of the present invention of those skilled in the art, remodeling, interpolation or replacement, also should belong to protection scope of the present invention.

Claims (2)

1. the Full wave shape joint inversion method improving high-dip structure velocity inversion efficiency, it is characterised in that: carry out in accordance with the following steps:

Step 1: the input speed of initial velocity field, field inspection big gun record and source wavelet, and set up observation system;

Step 2: use the imaging formula that following formula represents to obtain imaging results:

$I=\underset{x_s}{\Σ}\underset{t}{\Σ}uw/u^2$

Wherein I represents imaging results, and u and w represents main story wave field and anti-pass wave field, w=L respectively^*(R^*(p_obs)), L^*For anti-pass operator, R^*Represent data residual error spatial spread to whole space, p_obsFor field inspection big gun record;

Step 3: using the imaging results in step 2 as reflection coefficient;

Step 4: apply linear forward simulation and linear wave field anti-pass, calculates the gradient direction of two kinds of prism ripple prism1 and prism2 respectively;

$g{\left(v\right)}_{prism1}=\frac{2}{v^3}\underset{x_s}{\Σ}\underset{t}{\Σ}\frac{{\∂}^{2}B\left(u\right)}{\∂t^2}{w}^{\′}$

$g{\left(v\right)}_{prism2}=\frac{2}{v^3}\underset{x_s}{\Σ}\underset{t}{\Σ}\frac{{\∂}^{2}u}{\∂t^2}{B}^{*}\left(w^{\′}\right)$

Wherein g (v)_prism1With g (v)_prism2Representing the gradient direction of prism1 and the gradient direction of prism2 respectively, v is the speed of velocity field, and w' represents the anti-pass wave field of field inspection big gun record and the earthquake record residual difference of numerical simulation, x_sFor hypocentral location, t is the time;

Step 5: two gradient directions of gained in step 4 are added the gradient direction that can try to achieve prism waveform inversion;

$g{\left(v\right)}_{prism}=\frac{2}{v^3}\underset{x_s}{\Σ}\underset{t}{\Σ}\frac{{\∂}^{2}B\left(u\right)}{\∂t^2}{w}^{\′}+\frac{2}{v^3}\underset{x_s}{\Σ}\underset{t}{\Σ}\frac{{\∂}^{2}u}{\∂t^2}{B}^{*}\left(w^{\′}\right)$

Wherein g (v)_prismRepresent the gradient direction of prism waveform inversion;

Step 6: ask for renewal step-length α with linear search method or parabolic approximating method;

Step 7: utilize the gradient direction asked for and update step-length renewal speed;

v_k=v_k-1+α_kg(v_k-1)_prism

Wherein, subscript k represents the number of times of iteration;

Step 8: utilize conventional full wave shape inversion method renewal speed;

v_k+1=v_k+α_kg(v_k)

Wherein, g (v_k) for the gradient direction of kth time iteration conventional full wave shape inverting;

Step 9: judge v_k+1With v_kDifference whether meet error condition;

If: judged result is v_k+1With v_kDifference be unsatisfactory for error condition, then perform step 2;

Or judged result is v_k+1With v_kDifference meet error condition, then perform step 10;

Step 10: the rate results of output inverting.

2. the Full wave shape joint inversion method of raising high-dip structure velocity inversion efficiency according to claim 1, it is characterised in that: in step 4, specifically include

Step 4.1: objective function:

$E\left(v\right)=\frac{1}{2}\underset{x_r}{\Σ}\underset{x_s}{\Σ}\underset{t}{\Σ}{\left(Ru\right(t,x_r,x_s)-p{\left(t,x_r,x_s\right)}_{obs})}^2---\left(1\right);$

Wherein, Ru (t, x_r,x_s) and p (t, x_r,x_s)_obsThe respectively earthquake record of numerical simulation and field inspection big gun record, x_sAnd x_rRepresenting shot point and geophone station respectively, t is the time;

Step 4.2: by object function variation, obtains variation expression formula:

$\begin{array}{c}\mathrm{\δ}E\left(v\right)=\frac{1}{2}\mathrm{\δ}\[(Ru-p_{obs})\·(Ru-p_{obs})\]\\ =\[\mathrm{\δ}(Ru-p_{obs})\]\·(Ru-p_{obs})\\ =\left(R\mathrm{\δ}u\right)\·(Ru-p_{obs})\end{array}---\left(2\right);$

Step 4.3: definition two dimension Acoustic Wave-equation:

$\frac{1}{v^2}\frac{{\∂}^{2}u}{\∂t^2}-\▿\·\▿u=s(t,x_s)---\left(3\right);$

Equation (3) is deformed:

$\frac{1}{{(v+\mathrm{\δ}v)}^2}\frac{{\∂}^2(u+\mathrm{\δ}u)}{\∂t^2}-\▿\·\▿(u+\mathrm{\δ}u)=s(t,x_s)---\left(4\right);$

Wherein, s (t, x_s) represent focus item, represent Laplace operator;

Equation (4) is carried out Taylor expansion:

$(\frac{1}{v^2}-\frac{2\mathrm{\δ}v}{v^3})\frac{{\∂}^2(u+\mathrm{\δ}u)}{\∂t^2}-\▿\·\▿(u+\mathrm{\δ}u)=s(t,x_s)---\left(5\right);$

Launch to omit higher order term to formula (5):

$\frac{1}{v^2}\frac{{\∂}^{2}u}{\∂t^2}-\▿\·\▿u+\frac{1}{v^2}\frac{{\∂}^2\mathrm{\δ}u}{\∂t^2}-\▿\·\▿\mathrm{\δ}u-\frac{2\mathrm{\δ}v}{v^3}\frac{{\∂}^{2}u}{\∂t^2}=s(t,x_s)---\left(6\right);$

Equation (4) and equation (6) are subtracted each other, obtain:

$\frac{1}{v^2}\frac{{\∂}^2\mathrm{\δ}u}{\∂t^2}-\▿\·\▿\mathrm{\δ}u=\frac{2\mathrm{\δ}v}{v^3}\frac{{\∂}^{2}u}{\∂t^2}---\left(7\right);$

Can obtain further:

$\mathrm{\δ}u=L\left(\frac{2\mathrm{\δ}v}{v^3}\frac{{\∂}^{2}u}{\∂t^2}\right)---\left(8\right);$

Wherein, L represents positive communication process;

Step 4.4: equation (8) is substituted into variation expression formula (2), can obtain:

$\begin{array}{c}\mathrm{\δ}E\left(v\right)=\left(R\right(L\left(\frac{2\mathrm{\δ}v}{v^3}\frac{{\∂}^{2}u}{\∂t^2}\right)\left)\right)\·(Ru-p_{obs})\\ =\left(L\right(\frac{2\mathrm{\δ}v}{v^3}\frac{{\∂}^{2}u}{\∂t^2}\left)\right)\·R^{*}(Ru-p_{obs})\\ =\left(\frac{2\mathrm{\δ}v}{v^3}\frac{{\∂}^{2}u}{\∂t^2}\right)\·L^{*}\left(R^{*}\right(Ru-p_{obs}\left)\right)\end{array}---\left(9\right);$

Wherein, R^*(Ru-p_obs) represent data residual error spatial spread to whole space, L^*(R^*(Ru-p_obs)) represent that the residue wave field inverse time propagates;

Step 4.5: make the data residual error be:

Δ p=Ru-p_obs(10);

Step 4.6: by main story wave field u and data residual delta p=Ru-p_obsIt is decomposed into following two parts:

u(t,x_r,x_s)=u₁(t,x_r,x_s)+u₂(t,x_r,x_s) (11);

Δp(t,x_r,x_s)=Δ p₁(t,x_r,x_s)+Δp₂(t,x_r,x_s) (12);

Wherein, u₁(t,x_r,x_s) and u₂(t,x_r,x_s) represent first order reflection wave field and second order prism wave field respectively; Δ p₁(t,x_r,x_s) and Δ p₂(t,x_r,x_s) represent the data residual error of data residual sum second order prism wave field of first order reflection wave field respectively;

Step 4.7: equation (10), (11), (12) are substituted into equation (9):

Step 4.8: ask for the object function gradient to rate pattern:

$\begin{array}{c}g\left(v\right)=\frac{\mathrm{\δ}E\left(v\right)}{\mathrm{\δ}v}=\frac{2}{v^3}\underset{x_s}{\mathrm{\Σ}}\underset{t}{\mathrm{\Σ}}\frac{{\∂}^2{u}_1}{\∂t^2}{L}^{*}\left(R^{*}\right({\mathrm{\Δp}}_1)\\ +\frac{2}{v^3}\underset{x_s}{\mathrm{\Σ}}\underset{t}{\mathrm{\Σ}}\frac{{\∂}^2{u}_1}{\∂t^2}{L}^{*}\left(R^{*}\right({\mathrm{\Δp}}_2)+\frac{2}{v^3}\underset{x_s}{\mathrm{\Σ}}\underset{t}{\mathrm{\Σ}}\frac{{\∂}^2{u}_2}{\∂t^2}{L}^{*}\left(R^{*}\right({\mathrm{\Δp}}_1)\end{array}---\left(14\right);$

Make g (v)_reflect, g (v)_prism1With g (v)_prism2Represent the gradient of primary reflection respectively, the gradient of prism1 and the gradient of prism2, then

$g{\left(v\right)}_{reflect}=\frac{2}{v^3}\underset{x_s}{\Σ}\underset{t}{\Σ}\frac{{\∂}^2{u}_1}{\∂t^2}{L}^{*}\left(R^{*}\right({\mathrm{\Δp}}_1)---\left(15\right);$

$g{\left(v\right)}_{prism1}=\frac{2}{v^3}\underset{x_s}{\Σ}\underset{t}{\Σ}\frac{{\∂}^2{u}_2}{\∂t^2}{L}^{*}\left(R^{*}\right({\mathrm{\Δp}}_1)---\left(16\right);$

$g{\left(v\right)}_{prism2}=\frac{2}{v^3}\underset{x_s}{\Σ}\underset{t}{\Σ}\frac{{\∂}^2{u}_1}{\∂t^2}{L}^{*}\left(R^{*}\right({\mathrm{\Δp}}_2)---\left(17\right);$

By first order reflection wave field u₁(t,x_r,x_s) represented by main story wave field u, second order prism wave field u₂(t,x_r,x_s) represented by linear forward simulation B (u), the L of the data residual error of first order reflection wave field^*(R^*(Δp₁) represented by anti-pass wave field w', the L of the data residual error of second order prism wave field^*(R^*(Δp₂) by the linear wavelength anti-pass B of anti-pass wave field w'^*(w') represent, do corresponding variable replacement, then

$g{\left(v\right)}_{prism1}=\frac{2}{v^3}\underset{x_s}{\Σ}\underset{t}{\Σ}\frac{{\∂}^{2}B\left(u\right)}{\∂t^2}{w}^{\′}---\left(18\right);$

$g{\left(v\right)}_{prism2}=\frac{2}{v^3}\underset{x_s}{\Σ}\underset{t}{\Σ}\frac{{\∂}^{2}u}{\∂t^2}{B}^{*}\left(w^{\′}\right)---\left(19\right).$