CN105319589A - Full-automatic three-dimensional tomography inversion method using a local event slope - Google Patents

Abstract

The invention provides a full-automatic three-dimensional tomography inversion method using a local event slope and belongs to the seismic imaging and inversion field in oil and gas exploration and development. The method includes the following steps that: (1) an event slope is picked up in an automated manner; (2) an event slope coordinate transformation identical equation is utilized to perform quality monitoring, data points (S, R, PS, PR, t) pick are picked up in a filtering manner; and (3), a regularized three-dimensional tomography equation is established and solved.

Description

A kind of fully automatic stereo chromatography conversion method utilizing local lineups slope

Technical field

The invention belongs to the seismic imaging in oil-gas exploration and development and inverting field, be specifically related to a kind of fully automatic stereo chromatography conversion method utilizing local lineups slope, the enforcement of this technology does not need to carry out manual intervention, and has degree of precision.

Background technology

Estimating that the speed (especially its lower wave number part, also claims macroscopic velocity) of underground medium distributes is one of key problem of seismic prospecting, all depends on a rate pattern accurately to underground structural offset imaging, predicting oil/gas reservoir parameter.Seismic Tomography has higher inversion accuracy, and be the important method estimating subsurface velocity model, wherein most widely used is reflection tomographic.Reflection tomography has the advantages such as efficiency is high, easy realization, but the lineups that its needs pick up along the entire profile continuous distribution on pre-stack seismic road collection, this operation (particularly under three-dimensional situation) very consuming time, and be almost impossible mission for low signal-to-noise ratio earthquake data before superposition.The macroscopic velocity estimation technique that an other class is used widely in the industry is migration velocity analysis (MVA).The common imaging gather (CIG) that the principle of tomographic inversion generates with migration imaging by such technology combines, and reduces the requirement of velocity estimation technique to earthquake data before superposition signal to noise ratio (S/N ratio), has higher practical value.But such technology needs repeatedly to implement pre-stack depth migration and man-machine interaction pickup common imaging gather, this makes such technology very poor efficiency under three-dimensional situation, reduces its practical value to a certain extent.

Summary of the invention

The object of the invention is to solve the difficult problem existed in above-mentioned prior art, a kind of fully automatic stereo chromatography conversion method utilizing local lineups slope is provided, improve precision and the automaticity of velocity estimation technique.Utilize local lineups slope to carry out this operation of continuous pickup that Tomography Velocity inverting avoids lineups, comparatively reflection tomography is adapted to low signal-to-noise ratio earthquake data before superposition more in theory; Meanwhile, except when reflection seismic is walked, lineups slope is used for retrain subsurface velocity model, alleviates the ray multipath phenomenon of reflection tomographic, there is higher precision and inverting stability.In addition, the invention belongs to the category of tomography, only need an enforcement data pickup operation, therefore compared with migration velocity analysis, there is higher counting yield.In addition, by introducing " lineups slope coordinate transform identical relation " this quality monitoring condition, present invention achieves a kind of tomographic inversion technology of full-automation.

The present invention is achieved by the following technical solutions:

Utilize a fully automatic stereo chromatography conversion method for local lineups slope, comprising:

(1) robotization pickup lineups slope;

(2) utilize lineups slope coordinate transform identical relation to carry out quality monitoring, filter and pick up out data point (S, R, P _s, P _r, t) _pick;

The three-dimensional chromatography equation of regularization foundation with solve.

Described step (1) is achieved in that

Pre stack data is divided common-shot-gather of hanking, common detector gather and common offset road collection, then utilize structure tensor analytical approach to scan the optimal partial lineups slope obtaining these three each sampled points of road collection data, be denoted as P respectively _r, P _s, P _m.

Described step (2) is achieved in that

Described lineups slope coordinate transform identical relation is as follows:

P _S+P _R＝P _M

If the error of lineups slope coordinate transform identical relation meets error range given in advance, be then effectively pickup, preserve the data point (S, R, the P that pick up and obtain _s, P _r, t) _pick;

Described error absolute error ε ₁=P _s+ P _r-P _mor relative error ε ₂=(P _s+ P _r-P _m)/P _mrepresent.

The three-dimensional chromatography equation of regularization in described step (3) is as shown in (7) formula:

$\left[\begin{array}{c}\mathrm{DF}\\ {\mathrm{\ϵ}}_{d}I\\ {\mathrm{\ϵ}}_{C_1}{D}_1\\ {\mathrm{\ϵ}}_{C_2}{D}_2\end{array}\right]\mathrm{\Δm}=\left[\begin{array}{c}\mathrm{D\Δd}\\ 0\\ 0\\ 0\end{array}\right]---\left(7\right)$

(7) in formula, D is data precondition matrix, is used for the data space component of different dimensions of the three-dimensional tomographic inversion of weighting;

F is at m _nthe Fr é chet Jacobian matrix of place's value is a N _data× N _modelmatrix, wherein N _datait is the dimension of the data component used by three-dimensional tomographic inversion;

I is a unit matrix identical with model space dimension

be the arithmetic number between 0 to 1, be set by the user;

D ₁, D ₂to the first difference operator of B-spline basis function coefficient update amount on horizontal and vertical respectively;

To each data point (s, r, p _sx, p _rx, t) _i, Δ d is the data residual error of inverting, Δ d (m)=d-f (m), wherein:

d＝(S,R,P _S,P _R,t) _i,i＝1,2,...,n _data(3)

m＝[(x ₀,z ₀,θ _s,θ _r) _i,i＝1,2,...,n _data；v(x,z)](4)

F represents the non-linear positive operator of three-dimensional chromatography, and namely reflection spot is respectively to the initial value ray tracing in shot point, geophone station direction.

Described step (3) is completed by iteration, specific as follows:

If the model component of kth time iteration is m _k, then the process of kth+1 iteration is:

(31) model component m, is utilized _kforward modelling data component f (m _k) and set up the three-dimensional chromatography equation of described regularization;

(32) the model modification amount Δ m that the three-dimensional chromatography equation of regularization obtains kth+1 iteration, is solved _k+1;

(33) the model component m that model modification obtains kth+1 iteration, is carried out _k+1=m _k+ Δ m _k+1;

(34) if data residual error meets user-defined threshold values, then proceed to step (35), if do not meet, then return step (31) and enter next iteration;

(35) exit iteration, export final model component.

Data residual error in described step (34) is two norms of the residual error of the data component of forward modelling and the data component of pickup.

Compared with prior art, the invention has the beneficial effects as follows:

The first, whole tomographic inversion iterative process is full-automatic, does not need to carry out manual intervention, which increases the efficiency of whole velocity estimation process, decrease labor workload, shorten treatment cycle;

The second, be used for retraining rate pattern owing to introducing lineups slope, the rate pattern that inverting obtains has higher precision, meets the requirement of follow-up migration imaging.

Accompanying drawing explanation

Fig. 1 true velocity model;

The part big gun road collection that Fig. 2 wave equation finite difference is just being drilled;

The full-automatic pickup result of Fig. 3 offset distance-2km is fitted on common offset road collection;

The initial velocity model that Fig. 4 is given;

The reflection point position distribution that Fig. 5 is initial;

The rate pattern of Fig. 6 the 6th iteration;

The reflection point position distribution of Fig. 7 the 6th iteration;

The rate pattern of Fig. 8 the 17th iteration;

The reflection point position distribution of Fig. 9 the 17th iteration.

Figure 10, using the rate pattern of the 17th iteration as input, implements the migrated section that Gaussian beam pre-stack depth migration obtains;

Figure 11, using the rate pattern of the 17th iteration as input, implements the common imaging gather that Gaussian beam pre-stack depth migration obtains;

Figure 12 Figure 11 is the partial enlarged drawing of about 150 at horizontal ordinate.

The step block diagram of Figure 13 this method.

Figure 14 is by s ₀to s ₁transmission process.

Embodiment

Below in conjunction with accompanying drawing, the present invention is described in further detail:

The invention belongs to the chromatography imaging technique under linearized inversion theoretical frame, its main technology contents and flow process are: first on Gong Bao road collection, common detector gather, common offset road collection, utilize structure tensor analytical technology to calculate lineups slope, fully automatically filter out according to " lineups slope coordinate transform identical relation " this quality monitoring condition on these three road collection the data point met the demands, the data as inverting input; Then given initial model parameter, carries out ray tracing and is just drilling and calculating computed tomography equation, solve the model modification amount that chromatography equation obtains current iteration; Finally carry out model modification and enter next iteration until meet convergent requirement given in advance.This process full automation, does not carry out extra manual intervention.

As shown in figure 13, step of the present invention is as follows:

(1) the robotization pickup of lineups slope.

Pre stack data is divided common-shot-gather of hanking, common detector gather and common offset road collection, utilize structure tensor analytical technology to scan the optimal partial lineups slope obtaining these three each sampled points of road collection data, be denoted as P respectively _r, P _s, P _m.

For the estimation of seismic event slope local, the present invention adopts the computing method of local structure tensor (please refer to document: LucasJ.vanVlietandPietW.Verbeek.1995, EstimatorsforOrientationandAnisotropyinDigitizedImages).If I is two-dimension earthquake image, in two-dimensional image I, the structure tensor of representation space directional information is defined by image gradient value, structure tensor represents the change direction in region and the variable quantity size along change direction, and seismic strata texture and tomography texture are determined by local each point azimuth information variation relation.Introduce the fuzzy local detail of Gaussian function, make structure tensor highlight the complicacy of signal in region.To two dimensional image, structure tensor is the matrix of a 2*2:

$G=\left(\begin{array}{cc}<{g_x}^2>& <g_x{g}_y>\\ <g_x{g}_y>& <{g_y}^2>\end{array}\right)---\left(1\right)$

Wherein g _xwith g _yrepresent seismic image in the horizontal direction with the gradient of vertical direction, <> represents the smooth filtering of dimensional Gaussian.For positive semidefinite matrix G, eigen vector can by solving | and G-λ I|=0 obtains:

$G=\left[\begin{array}{cc}{v}_1& v_2\end{array}\right]\left[\begin{array}{cc}{\mathrm{\&lambda;}}_1& 0\\ 0& {\mathrm{\&lambda;}}_2\end{array}\right]\left[\begin{array}{cc}{v_1}^T& {v_2}^T\end{array}\right]---\left(2\right)$

λ ₁: eigenvalue of maximum, tensor energy is at first characteristic tensor direction v ₁energy.

λ ₂: minimal eigenvalue, tensor energy is at second characteristic tensor direction v ₂energy.

(λ ₁-λ ₂)/λ ₁: the linearity, the consistance of reflection local direction.

Proper vector describes the directivity of image local linear structure, for each point of image, and proper vector v ₁be orthogonal to the main structure direction of image, proper vector v ₂be parallel to the main structure direction of image.

In the present invention to certain in earthquake data before superposition a bit, the vector v solved by structure tensor technology ₂for v ₂=(e ₁, e ₂), then e ₂/ e ₁be the local lineups slope of this point.P _r, P _s, P _mcorrespond respectively to the lineups slope that earthquake data before superposition Gong Bao road collection, common detector gather and common offset road are concentrated.

(2) picks up data point is filtered in " lineups slope coordinate transform identical relation " quality monitoring.

To each reflection axle point in the same way of position, reflection horizon in common offset section, examine position according to information, big gun when walking, utilize above-mentioned structure tensor analytical technology to calculate the slope P concentrated in big gun road _r, the slope P that geophone station road is concentrated _s, with the slope P in common offset section _mcompare, if P _s+ P _r=P _merror (this error absolute error ε of this identical relation ₁=P _s+ P _r-P _mor relative error ε ₂=(P _s+ P _r-P _m)/P _mdescribe) meet error range (this scope (i.e. absolute error ε given in advance ₁with relative error ε ₂scope) given by User Defined, these two parameter-dependents are in the quality of geological data, higher then these two error upper limits of quality of geological data are given less, otherwise it is given relatively large, user can carry out quality monitoring by the data point of picking up out being shown, if picks up data point is crossed the error upper limit to be amplified at least and is again picked up), be then effectively pickup, preserve pickup result.Input introductory die component in Figure 13 is given by user, and its concrete element is defined by formula (4).

Regularization chromatography equation foundation with solve.

Three-dimensional tomographic inversion belongs to the speed estimation method of data fitting class, namely finds one as the model vector m of (4) formula definition, makes the data d that forward modelling obtains _modthe error minimize of the picks up data d that=f (m) and formula (4) define.F represents the non-linear positive operator of three-dimensional chromatography, and namely reflection spot is respectively to the initial value ray tracing in shot point, geophone station direction.The concrete criterion that whether correct three-dimensional tomographic inversion weigh model is: consider that right reflection examined by a big gun, if the coordinate of reflection spot, incident angle and reflection angle and rate pattern correct time, by reflection spot respectively with incident angle and reflection angle to shot point direction and geophone station direction emergent ray, when arriving inspection surface, two transmitted ray just arrive respectively and reach correct shot point and geophone station with correct exit direction, correct two way travel time.

The three-dimensional chromatographic data component that the present invention adopts and model component are respectively:

d＝(S,R,P _S,P _R,t) _i,i＝1,2,...,n _data(3)

m＝[(x ₀,z ₀,θ _s,θ _r) _i,i＝1,2,...,n _data；v(x,z)](4)

In formula (3), data are respectively the n picking up and arrive _dataindividual data point (i.e. the result of step (2)), each data point is specifically respectively shot point horizontal ordinate, geophone station horizontal ordinate, the horizontal component of shot point place slowness vector, the horizontal component of geophone station place slowness vector and two way travel time.

Equally, each data point corresponds to a reflection/diffraction process, corresponds to a group model parameter (x of (4) formula in a model ₀, z ₀, θ _s, θ _r) _i, i.e. reflection spot horizontal ordinate, ordinate, incident angle, reflection angle.

Weigh data fitting errors with two norms, what three-dimensional tomographic inversion was summed up as following functional asks extreme-value problem:

$S\left(m\right)=\frac{1}{2}{\left|\right|d-f\left(m\right)\left|\right|}_2^2=\frac{1}{2}{\mathrm{\&Delta;d}}^T\left(m\right)C_D^{-1}\mathrm{\&Delta;d}\left(m\right)---\left(5\right)$

Wherein, Δ d (m)=d-f (m).Diagonal matrix also referred to as data covariance matrix, play weighting different pieces of information component.

The pathosis intrinsic due to indirect problem and multi-solution, the present invention introduces the prior-constrained reliability for ensureing inversion result of extra model, and the error functional of transformation (5) formula is as follows:

$\begin{array}{c}S\left(m\right)={\left|\right|d-f\left(m\right)\left|\right|}_2^2+{\mathrm{\&epsiv;}}_d^2{\left|\right|m-m_r\left|\right|}_2^2+\\ {\mathrm{\&epsiv;}}_{C_1}^2{\left|\right|D_1(m-m_r)\left|\right|}_2^2+{\mathrm{\&epsiv;}}_{C_2}^2{\left|\right|D_2(m-m_r)\left|\right|}_2^2\end{array}---\left(6\right)$

Wherein for less arithmetic number ( these three parameter distribution are between 0 to 1, these three parameters also need user to pass through to test self-defined providing, when the initial model (form is as formula (4) definition) of inverting is better, higher then these parameters of data quality are given less, otherwise given larger), for the weight between weighted data matching item and regularization term; m _rfor the prior model that certain is given, the form of prior model is the inverse model component defined by formula (4), needs User Defined to provide prior model when implementing inverting, and the method for given prior model is a lot, and the present invention gets prior model m _r=0. for damping term, this introducing makes the model modification amount of each iteration less. with retrain the smooth degree of rate pattern on horizontal and vertical respectively, D ₁, D ₂to the first difference operator of B-spline basis function coefficient update amount on horizontal and vertical respectively.One is had to the ordered series of numbers of 7 numbers, its first order difference operator definitions is the matrix of following 6 × 7,

$\left(\begin{array}{ccccccc}1& -1& 0& 0& 0& 0& 0\\ 0& 1& -1& 0& 0& 0& 0\\ 0& 0& 1& -1& 0& 0& 0\\ 0& 0& 0& 1& -1& 0& 0\\ 0& 0& 0& 0& 1& -1& 0\\ 0& 0& 0& 0& 0& 1& -1\end{array}\right),$

In like manner, for the ordered series of numbers having N number of point, first difference operator is that shape is as above formula (N-1) × N matrix).

Theoretical according to the suboptimization of inverting, (6) optimization problem that formula is portrayed can be equivalent to the iterative of following system of linear equations, reached the object of inversion solution convergence by successive ignition, wherein system of equations (7) is called the chromatography equation of three-dimensional tomographic inversion.

$\left[\begin{array}{c}\mathrm{DF}\\ {\mathrm{\&epsiv;}}_{d}I\\ {\mathrm{\&epsiv;}}_{C_1}{D}_1\\ {\mathrm{\&epsiv;}}_{C_2}{D}_2\end{array}\right]\mathrm{\&Delta;m}=\left[\begin{array}{c}\mathrm{D\&Delta;d}\\ 0\\ 0\\ 0\end{array}\right]---\left(7\right)$

(7), in formula, Δ d is the data residual error of inverting, Δ d (m)=d-f (m); I is a unit matrix identical with model space dimension wherein N _modelbe the dimension of the model component used by three-dimensional tomographic inversion, F is Fr é chet Jacobian matrix, is a N _data× N _modelmatrix, wherein N _databe the dimension of the data component used by three-dimensional tomographic inversion, D is data precondition matrix, is used for the data space component of different dimensions of the three-dimensional tomographic inversion of weighting, to each data point (s, r, p _sx, p _rx, t) _i, D is by following formal definition:

wherein the present invention gets $({\mathrm{\&sigma;}}_s,{\mathrm{\&sigma;}}_r,{\mathrm{\&sigma;}}_{p_{\mathrm{sr}}},{\mathrm{\&sigma;}}_{p_{\mathrm{rx}}},{\mathrm{\&sigma;}}_t)=(\mathrm{1.0,1.0},{10}^6,{10}^3).$

Wherein f represents that transmitted ray follows the trail of operator (comparatively ripe in oil gas field of geophysical exploration), and m is the model component that formula (4) defines; F is at m _nthe Fr é chet Jacobian matrix of place's value, wherein f _ijrepresent the element of the i-th row jth row of Fr é chet Jacobian matrix, specific as follows:

The calculating of the Fr é chet derivative of three-dimensional chromatography is based on ray-tracing diffracted (Fara & Madariaga, 1987) with character (the Gilbert & Bacus of ray propogator matrix, 1996), consider the succinct form of ray-tracing diffracted under ray center coordinate system, the present invention derives and obtains the Fr é chet derivative of three-dimensional tomographic inversion under this coordinate system.S (shot point) → X (reflection spot) → R (geophone station) this reflection process is split as these two transmission process of X → S and X → R by the direct problem of three-dimensional chromatography, therefore considers the transmission process of the ground shown in Figure 14 down to earth's surface when calculating the Fr é chet derivative of three-dimensional chromatography.

The Fr é chet Derivative Formula that first three-dimensional chromatography is discussed is derived.Obviously, at Fr é chet derivative in, have:

$\frac{\&PartialD;(S,P_S)}{{\&PartialD;\mathrm{\&theta;}}_r}=0,\frac{\&PartialD;(R,P_R)}{{\&PartialD;\mathrm{\&theta;}}_s}=0---A\left(1\right)$

Consider t=t _s+ t _r, have:

$\frac{\&PartialD;t}{\&PartialD;(x_0,z_0,{\mathrm{\&theta;}}_s,{\mathrm{\&theta;}}_r,v)}=\frac{{\&PartialD;t}_S}{\&PartialD;(x_0,z_0,{\mathrm{\&theta;}}_s,{\mathrm{\&theta;}}_r,v)}+\frac{{\&PartialD;t}_R}{\&PartialD;(x_0,z_0,{\mathrm{\&theta;}}_s,{\mathrm{\&theta;}}_r,v)}---A\left(2\right)$

Wherein, $\frac{{\&PartialD;t}_S}{{\&PartialD;\mathrm{\&theta;}}_r}=0,\frac{{\&PartialD;t}_R}{{\&PartialD;\mathrm{\&theta;}}_s}=0.$

Analysis is known with method of asking be on all four, both correspond to respectively by the transmission process of reflection spot to shot point and geophone station direction.

As shown in figure 14, consider that one by s ₀to s ₁transmission process, have both sides total differential obtains:

$\mathrm{\&Delta;}{t}_s=\frac{\sin \mathrm{\&theta;}}{v\left(s_0\right)}\mathrm{\&Delta;x}-\frac{\cos \mathrm{\&theta;}}{v\left(s_0\right)}\mathrm{\&Delta;z}-\underset{s_0}{\overset{s_1}{\&Integral;}}\frac{\mathrm{\&Delta;v}}{v^2}\mathrm{ds}---A\left(3\right)$

Can be calculated by A (3)

calculating can be in the hope of by ray-tracing diffracted.Ray-tracing diffracted is the linear relationship establishing the observation disturbance quantity of ray field and the disturbance quantity of initial ray field and speed in essence, can converse (S, P by this linear relationship _s) disturbance quantity and (x ₀, z ₀, θ _s, the linear relationship of disturbance quantity v).

Ray-tracing diffracted form under ray center coordinate system is the simplest, therefore first under ray center coordinate system, applies ray-tracing diffracted, then derives (Δ S, Δ P according to the geometric relationship of rectangular coordinate and ray center coordinate _s) and (Δ x ₀, Δ z ₀, Δ θ _s, Δ v) linear relation, try to achieve

Disturbance coordinate q and corresponding slowness vector component p, has according to ray-tracing diffracted:

$\left(\begin{array}{c}{\mathrm{\&Delta;q}}_1\\ {\mathrm{\&Delta;p}}_1\end{array}\right)=\mathrm{\&Pi;}(s_1,s_0)\left(\begin{array}{c}{\mathrm{\&Delta;q}}_0\\ {\mathrm{\&Delta;p}}_0\end{array}\right)+\underset{s_0}{\overset{s_1}{\&Integral;}}\mathrm{\&Pi;}(s_1,s)\mathrm{\&Delta;w}(s,\mathrm{\&Delta;v}\left(s\right))\mathrm{ds}---A\left(4\right)$

Wherein: $\mathrm{\&Delta;w}={(\frac{\&PartialD;\mathrm{\&Delta;H}}{\&PartialD;p},-\frac{\&PartialD;\mathrm{\&Delta;H}}{\&PartialD;q})}^T={(0,-\frac{1}{v^2}\frac{\&PartialD;\mathrm{\&Delta;v}}{\&PartialD;q}+\frac{1}{v^3}\frac{\&PartialD;v}{\&PartialD;q}\mathrm{\&Delta;v})}^T\mathrm{\&Pi;}(s_1,s_0)$ Be called from s ₀point propagates into s ₁the ray propogator matrix of point, is the matrix of one 2 × 2, is designated as:

$\mathrm{\&Pi;}(s_1,s_0)=\left(\begin{array}{cc}{Q}_1(s_1,s_0)& Q_2(s_1,s_0)\\ P_1(s_1,s_0)& P_2(s_1,s_0)\end{array}\right)---A\left(5\right)$

(Q ₁(s ₁, s ₀), P ₁(s ₁, s ₀)) ^twith (Q ₂(s ₁, s ₀), P ₂(s ₁, s ₀)) ^tbe respectively Paraxial ray-gracing system initial value for (1,0) ^twith (0,1) ^tsolution, physically correspond respectively to normalized Plane wave source and single source.Paraxial ray-gracing system under ray center coordinate system is:

$\frac{d}{\mathrm{ds}}\left(\begin{array}{c}\mathrm{\&Delta;q}\\ \mathrm{\&Delta;p}\end{array}\right)=\stackrel{\&OverBar;}{S}\left(\begin{array}{c}\mathrm{\&Delta;q}\\ \mathrm{\&Delta;p}\end{array}\right),\stackrel{\&OverBar;}{S}=\left(\begin{array}{cc}0& v\left(s\right)\\ -v^{-2}\frac{{\&PartialD;}^{2}v}{{\&PartialD;q}^2}& 0\end{array}\right)---A\left(6\right)$

Ray propogator matrix has following character:

∏(s ₁,s)＝∏(s ₁,s ₀)∏(s ₀,s)＝∏(s ₁,s ₀)∏ ^-1(s,s ₀)A(7)

${\mathrm{\&Pi;}}^{-1}(s,s_0)=\left(\begin{array}{cc}{P}_2(s,s_0)& {-Q}_2(s,s_0)\\ {-P}_1(s,s_0)& Q_1(s,s_0)\end{array}\right)---A\left(8\right)$

Obtained by A (4), A (5), A (7), A (8):

$\left(\begin{array}{c}\mathrm{\&Delta;}{q}_1\\ {\mathrm{\&Delta;p}}_1\end{array}\right)=\mathrm{\&Pi;}(s_1,s_0)[\left(\begin{array}{c}{\mathrm{\&Delta;q}}_0\\ {\mathrm{\&Delta;p}}_0\end{array}\right)+\left(\begin{array}{c}\mathrm{\&Delta;}\stackrel{\&OverBar;}{q}\left(\mathrm{\&Delta;v}\right)\\ \mathrm{\&Delta;}\stackrel{\&OverBar;}{p}\left(\mathrm{\&Delta;v}\right)\end{array}\right)]---A\left(9\right)$

Wherein,

$\mathrm{\&Delta;}\stackrel{\&OverBar;}{q}\left(\mathrm{\&Delta;v}\right)=-\underset{s_0}{\overset{s_1}{\&Integral;}}{Q}_2(s,s_0)(-\frac{1}{v^2}\frac{\&PartialD;\mathrm{\&Delta;v}}{\&PartialD;q}+\frac{1}{v^3}\frac{\&PartialD;v}{\&PartialD;q}\mathrm{\&Delta;v})\mathrm{ds}$

$\mathrm{\&Delta;}\stackrel{\&OverBar;}{p}\left(\mathrm{\&Delta;v}\right)=\underset{s_0}{\overset{s_1}{\&Integral;}}{Q}_1(s,s_0)(-\frac{1}{v^2}\frac{\&PartialD;\mathrm{\&Delta;v}}{\&PartialD;q}+\frac{1}{v^3}\frac{\&PartialD;v}{\&PartialD;q}\mathrm{\&Delta;v})\mathrm{ds}$

Consider the geometric relationship of ray center coordinate system and rectangular coordinate system, have:

$\mathrm{\&Delta;\&xi;}=\frac{1}{\cos \mathrm{\&alpha;}}\mathrm{\&Delta;}{q}_1,\mathrm{\&Delta;}{p}^{\mathrm{\&xi;}}=\frac{\cos \mathrm{\&alpha;}}{v\left(s_1\right)}\mathrm{\&Delta;\&alpha;}=\cos \mathrm{\&alpha;\&Delta;}{p}_1---A\left(10\right)$

Δq ₀＝cosθΔx ₀-sinθΔz ₀， ${\mathrm{\&Delta;p}}_0=\frac{\sin \left(\mathrm{\&Delta;\&theta;}\right)}{v\left(s_0\right)}\&ap;\frac{1}{v\left(s_0\right)}\mathrm{\&Delta;\&theta;}---A\left(11\right)$

Wherein ξ, p ^ξbe respectively the component of slowness vector on horizontal ordinate under rectangular coordinate system and this coordinate direction, as figure A.1 shown in, with when arriving earth's surface and the angle of the longitudinal axis when θ, α are respectively transmitted ray outgoing.A (10), A (11) are substituted into A (9) obtain:

$\begin{array}{c}\mathrm{\&Delta;\&xi;}=\frac{\cos \mathrm{\&theta;}}{\cos \mathrm{\&alpha;}}{Q}_1(s_1,s_0)\mathrm{\&Delta;}{x}_0\\ -\frac{\sin \mathrm{\&theta;}}{\cos \mathrm{\&alpha;}}{Q}_1(s_1,s_0)\mathrm{\&Delta;}{z}_0\\ +\frac{1}{v\left(s_0\right)\cos \mathrm{\&alpha;}}{Q}_2(s_1,s_0)\mathrm{\&Delta;\&theta;}\\ +\frac{1}{\cos \mathrm{\&alpha;}}[Q_1(s_1,s_0)\mathrm{\&Delta;}\stackrel{\&OverBar;}{q}\left(\mathrm{\&Delta;v}\right)+Q_2(s_1,s_0)\mathrm{\&Delta;}\stackrel{\&OverBar;}{p}\left(\mathrm{\&Delta;v}\right)]\end{array}---A\left(12\right)$

$\begin{array}{c}\mathrm{\&Delta;}{p}^{\mathrm{\&xi;}}=\cos {\mathrm{\&alpha;}\cos \mathrm{\&theta;P}}_1(s_1,s_0)\mathrm{\&Delta;}{x}_0\\ -{\cos \mathrm{\&alpha;}\sin \mathrm{\&theta;P}}_1(s_1,s_0)\mathrm{\&Delta;}{z}_0\\ +\frac{\cos \mathrm{\&alpha;}}{v\left(s_0\right)}{P}_2(s_1,s_0)\mathrm{\&Delta;\&theta;}\\ +\cos \mathrm{\&alpha;}{[P}_1(s_1,s_0)\mathrm{\&Delta;}\stackrel{\&OverBar;}{q}\left(\mathrm{\&Delta;v}\right)+P_2(s_1,s_0)\mathrm{\&Delta;}\stackrel{\&OverBar;}{p}\left(\mathrm{\&Delta;v}\right)]\end{array}---A\left(13\right)$

Derivative can be tried to achieve according to A (12), A (13) it is corresponding to be asked with

The present invention utilizes least square QR method (LSQR) (to please refer to document: PaigeCC, SaundersMA.LSQR:Analgorithmforsparselinearequationsandsp arseleastsquares.ACMTransactiononMathematicalSoftware, 1982,8 (1): 43-71) (in system of equations (7), Δ m is unknown quantity to solution matrix system of equations (7), other are all known quantities), the method is a kind of method of iteration, can under least square meaning Solving Large Scale Sparse matrix efficiently.

The data component d that step (1), (2) this two step provide three-dimensional tomographic inversion to need, the often wheel iteration d of inverting carrys out right-hand vector Δ d (the m)=d-f (m) of group of equations (7), and wherein f (m) utilizes the forward modeling of last round of iteration to simulate to obtain.

An embodiment of the inventive method is as follows:

Fig. 1 is real underground medium rate pattern, this rate pattern laterally long 17 kilometers, longitudinally dark 3.72 kilometers, Fig. 2 utilizes the rate pattern of Fig. 1 to carry out finite difference numerical simulation to solve observation big gun road, earth's surface collection (part) data that ACOUSTIC WAVE EQUATION obtains, Fig. 3 utilizes step (1), (2) the part invert data point automatically extracted, when the data point that pickup obtains is walked, slope better can be fitted with geological data (offset distance is the common offset road collection of 2km), illustrate that automatic Picking flow process of the present invention has higher precision, the follow-up input as step (3) of picks up data point.Fig. 4 is the rate pattern of the initial normal gradient that this example is given, and the most shallow, the bottommost layer speed of rate pattern are respectively 2800m/s and 6000m/s, and velocity gradient is 0.860215s ^-1, Fig. 5 is that in initial velocity model and reflection point position superimposed figure, figure, black short line segment obtains according to initialized reflection point position and the Local Layer inclination angle that conversed by initialized incidence, reflection angle, and reduction formula is wherein θ _sincident angle, θ _rreflection angle, θ _dipbe Local Layer inclination angle, initial model is as the input of step (3).

The data component pickup of inverting is complete, the good model component of initialization, the iteration that just can start inverting upgrades, Fig. 6 is the rate pattern exported after 6 iteration, Fig. 7 is the reflection point position distribution exported after 6 iteration, can find out that the shallow-layer of rate pattern is effectively restored, the position of reflection spot is also progressively to true reflecting interface everywhere convergent.Fig. 8, Fig. 9 are rate pattern after 17 iteration and reflection point position distribution plan respectively, can find out the increase along with iterations, and inverting is progressively to the true direction approximation separated.Using the rate pattern that the 17th time iteration obtains as input, implement Gaussian beam pre-stack depth migration, the migration stack section obtained as shown in Figure 10, can find out that migrated image comparatively focuses on.Figure 11 is the angle domain common image gathers that Gaussian beam pre-stack depth migration exports, Figure 12 is the partial enlarged drawing of Figure 11, can find out that common imaging gather is evened up, illustrate that the rate pattern that inverting obtains has higher precision, meet the requirement of migration imaging, a high-quality tectonic profile can be provided for follow-up seismic interpretation.

The present invention calculates seismic event slope accurately by adopting structure tensor analysis, and combines " seismic event slope coordinate transform identical relation " this Quality Control condition, achieves a kind of full automatic three-dimensional tomographic inversion technology.Except when reflection seismic waves is walked, lineups slope is incorporated in reflection seismic tomography by this technology, enhance precision and the stability of tomography, the whole flow process of this technology does not need to carry out manual intervention, there is higher precision and automaticity, can provide an accurate subsurface velocity model efficiently for seismic exploration on a large scale, be follow-up migration imaging and reservoir prediction service.

Technique scheme is one embodiment of the present invention, for those skilled in the art, on the basis that the invention discloses application process and principle, be easy to make various types of improvement or distortion, and the method be not limited only to described by the above-mentioned embodiment of the present invention, therefore previously described mode is just preferred, and does not have restrictive meaning.

Claims (6)

1. utilize a fully automatic stereo chromatography conversion method for local lineups slope, it is characterized in that: described method comprises:

(1) robotization pickup lineups slope;

(2) utilize lineups slope coordinate transform identical relation to carry out quality monitoring, filter and pick up out data point (S, R, P _s, P _r, t) _pick;

The three-dimensional chromatography equation of regularization foundation with solve.

2. the fully automatic stereo chromatography conversion method utilizing local lineups slope according to claim 1, is characterized in that: described step (1) is achieved in that

Pre stack data is divided common-shot-gather of hanking, common detector gather and common offset road collection, then utilize structure tensor analytical approach to scan the optimal partial lineups slope obtaining these three each sampled points of road collection data, be denoted as P respectively _r, P _s, P _m.

3. the fully automatic stereo chromatography conversion method utilizing local lineups slope according to claim 2, is characterized in that: described step (2) is achieved in that

Described lineups slope coordinate transform identical relation is as follows:

P _S+P _R＝P _M

If the error of lineups slope coordinate transform identical relation meets error range given in advance, be then effectively pickup, preserve the data point (S, R, the P that pick up and obtain _s, P _r, t) _pick;

Described error absolute error ε ₁=P _s+ P _r-P _mor relative error ε ₂=(P _s+ P _r-P _m)/P _mrepresent.

4. the fully automatic stereo chromatography conversion method utilizing local lineups slope according to claim 3, is characterized in that: the three-dimensional chromatography equation of the regularization in described step (3) is as shown in (7) formula:

$\left[\begin{array}{c}\mathrm{DF}\\ {\mathrm{\&epsiv;}}_{d}I\\ {\mathrm{\&epsiv;}}_{C_1}{D}_1\\ {\mathrm{\&epsiv;}}_{C_2}{D}_2\end{array}\right]\mathrm{\&Delta;m}=\left[\begin{array}{c}\mathrm{D\&Delta;d}\\ 0\\ 0\\ 0\end{array}\right]---\left(7\right)$

(7) in formula, D is data precondition matrix, is used for the data space component of different dimensions of the three-dimensional tomographic inversion of weighting;

F is at m _nthe Fr é chet Jacobian matrix of place's value is a N _data× N _modelmatrix, wherein N _datathe dimension of the data component used by three-dimensional tomographic inversion, N _modelit is the dimension of the model component used by three-dimensional tomographic inversion;

I is a unit matrix identical with model space dimension

be the arithmetic number between 0 to 1, be set by the user;

D ₁, D ₂to the first difference operator of B-spline basis function coefficient update amount on horizontal and vertical respectively;

To each data point (s, r, p _sx, p _rx, t) _i, Δ d is the data residual error of inverting, Δ d (m)=d-f (m), wherein:

d＝(S,R,P _S,P _R,t) _i,i＝1,2,...,n _data(3)

m＝[(x ₀,z ₀,θ _s,θ _r) _i,i＝1,2,...,n _data；v(x,z)](4)

F represents the non-linear positive operator of three-dimensional chromatography, and namely reflection spot is respectively to the initial value ray tracing in shot point, geophone station direction.

5. the fully automatic stereo chromatography conversion method utilizing local lineups slope according to claim 4, is characterized in that: described step (3) is completed by iteration, specific as follows:

If the model component of kth time iteration is m _k, then the process of kth+1 iteration is:

(31) model component m, is utilized _kforward modelling data component f (m _k) and set up the three-dimensional chromatography equation of described regularization;

(32) the model modification amount Δ m that the three-dimensional chromatography equation of regularization obtains kth+1 iteration, is solved _k+1;

(33) the model component m that model modification obtains kth+1 iteration, is carried out _k+1=m _k+ Δ m _k+1;

(34) if data residual error meets user-defined threshold values, then proceed to step (35), if do not meet, then return step (31) and enter next iteration;

(35) exit iteration, export final model component.

6. the fully automatic stereo chromatography conversion method utilizing local lineups slope according to claim 5, is characterized in that: the data residual error in described step (34) is two norms of the residual error of the data component of forward modelling and the data component of pickup.