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

Full-automatic stereo chromatographic inversion method using local event slope

Technical Field

The invention belongs to the field of seismic imaging and inversion in oil-gas exploration and development, and particularly relates to a full-automatic three-dimensional chromatographic inversion method utilizing local in-phase axis slope.

Background

Estimating the velocity distribution of the subsurface medium (especially its low wavenumber portions, also known as macroscopic velocities) is one of the core problems of seismic exploration, and migration imaging of subsurface structures, and prediction of hydrocarbon reservoir parameters, all rely on an accurate velocity model. The seismic tomography technology has higher inversion accuracy and is an important method for estimating the subsurface velocity model, wherein the reflection tomography is most widely applied. Reflection tomography has the advantages of high efficiency, ease of implementation, etc., however it requires picking the in-phase axes on the prestack seismic gathers, which are distributed continuously along the entire profile, an operation that is time consuming (especially in three dimensions) and a task that is almost impossible to accomplish for low signal-to-noise ratio prestack seismic data. Another class of macroscopic velocity estimation techniques that is widely used in the industry is the offset velocity analysis (MVA). The technology combines the principle of chromatographic inversion with common imaging point gather (CIG) generated by migration imaging, reduces the requirement of the velocity estimation technology on the signal-to-noise ratio of the pre-stack seismic data, and has high practical value. However, the technology needs to implement prestack depth migration for multiple times and pick up a common imaging point gather through human-computer interaction, so that the technology is very inefficient under the three-dimensional condition, and the practical value of the technology is reduced to a certain extent.

Disclosure of Invention

The invention aims to solve the problems in the prior art, provides a full-automatic stereo tomographic inversion method using local in-phase axis slope, and improves the precision and automation degree of the speed estimation technology. The operation of continuous pickup of the same-phase axis is avoided by utilizing the local same-phase axis slope to carry out chromatographic velocity inversion, and the method is theoretically more suitable for low signal-to-noise ratio prestack seismic data than reflection chromatography; meanwhile, besides the travel time of the reflection earthquake, the same-phase axis slope is used for constraining the underground velocity model, the ray multipath phenomenon of reflection tomography is relieved, and the method has higher precision and inversion stability. In addition, the invention belongs to the field of tomography, only needs to carry out one data pickup operation, and therefore, has higher calculation efficiency compared with the offset velocity analysis. In addition, the invention realizes a full-automatic chromatographic inversion technology by introducing the quality monitoring condition of 'same-phase axis slope coordinate transformation identity'.

The invention is realized by the following technical scheme:

a full-automatic stereo tomographic inversion method using local event slope, comprising:

the method comprises the steps of automatically picking up the slope of a same phase shaft;

⑵ quality monitoring is performed by using the same phase axis slope coordinate transformation identity, filtering and picking data points (S, R, P)_S,P_R,t)_pick；

Establishing and solving a regularization stereo chromatography equation.

The step (1) is realized by the following steps:

the prestack data is divided into a common shot point gather, a common survey point gather and a common offset gather, then the optimal local homophase axis slope of each sampling point of the three gather data is obtained by scanning by using a structure tensor analysis method, and the prestack data is divided intoRespectively recorded as P_R，P_S，P_M。

The step (2) is realized by the following steps:

the same phase axis slope coordinate transformation identity is as follows:

P_S+P_R＝P_M

if the error of the homophase axis slope coordinate transformation identity satisfies the predetermined error range, the data points (S, R, P) are effectively picked up and stored_S,P_R,t)_pick；

Absolute error for said error₁＝P_S+P_R-P_MOr relative error₂＝(P_S+P_R-P_M)/P_MAnd (4) showing.

The regularized stereo chromatography equation in the step (3) is shown as the formula (7):

$\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 the formula, D is a data precondition matrix used for weighting data space components of different dimensions of the stereo tomography inversion;

f is in m_nThe Frechet derivative matrix of the value is N_data×N_modelIn which N is_dataIs the dimensionality of the data components used for the stereo tomographic inversion;

i is a unit matrix with the same dimension as the model space

All positive real numbers between 0 and 1 are set by a user;

D₁、D₂respectively carrying out first-order difference operators on the updating quantity of the B spline basis function coefficient in the transverse direction and the longitudinal direction;

for each data point (s, r, p)_sx,p_rx,t)_iΔ d is the inverted data residual, Δ d (m) d-f (m), where:

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 nonlinear positive operator of the stereo chromatography, namely the initial ray tracing from the reflection point to the shot point and the wave detection point respectively.

The step (3) is accomplished by iteration, specifically as follows:

let the model component of the kth iteration be m_kThen, the process of the (k + 1) th iteration is:

(31) using the model component m_kForward calculated data component f (m)_k) Establishing the regularized stereo chromatography equation;

(32) solving the regularization three-dimensional chromatography equation to obtain the model updating quantity delta m of the (k + 1) th iteration_k+1；

(33) Carrying out model updating to obtain a model component m of the (k + 1) th iteration_k+1＝m_k+Δm_k+1；

(34) If the data residual meets the user-defined threshold, the step (35) is carried out, and if the data residual does not meet the user-defined threshold, the step (31) is returned to enter the next iteration;

(35) and exiting iteration and outputting a final model component.

The data residual in said step (34) is a two-norm of the residuals of the forward calculated data component and the picked-up data component.

Compared with the prior art, the invention has the beneficial effects that:

firstly, the whole chromatographic inversion iteration process is full-automatic, and manual intervention is not needed, so that the efficiency of the whole speed estimation process is improved, the manual workload is reduced, and the processing period is shortened;

secondly, because the in-phase axis slope is introduced to constrain the velocity model, the velocity model obtained by inversion has higher precision and meets the requirement of subsequent offset imaging.

Drawings

FIG. 1 true velocity model;

FIG. 2 is a partial shot gather of finite difference forward modeling of the wave equation;

FIG. 3 shows a full auto pick-up result with an offset of-2 km, which is attached to a common offset gather;

FIG. 4 shows an initial velocity model;

FIG. 5 initial reflection point location distribution;

FIG. 6 speed model for iteration 6;

FIG. 7 shows a distribution of reflection point positions for iteration 6;

FIG. 8 velocity model for iteration 17;

fig. 9 reflection point position distribution of iteration 17.

FIG. 10 is a migration profile obtained by performing a Gaussian beam prestack depth migration using the velocity model of the 17 th iteration as input;

FIG. 11 shows a common imaging point gather obtained by performing Gaussian beam prestack depth migration using the 17 th iteration velocity model as input;

fig. 12 is a partial enlarged view of fig. 11 on the abscissa of about 150.

FIG. 13 is a block diagram of the steps of the method.

FIG. 14 is a graph formed by s₀To s₁The transmission process of (1).

Detailed Description

The invention is described in further detail below with reference to the accompanying drawings:

the invention belongs to a tomography technology under a linearized inversion theory framework, and the main technical content and the flow are as follows: firstly, calculating the slope of the same phase axis on a common shot gather, a common geophone point gather and a common offset gather by using a structure tensor analysis technology, and filtering out data points meeting the requirements in a full-automatic manner according to the quality monitoring condition of 'coordinate transformation identity of the same phase axis slope' on the three gathers to be used as inverted data input; then, giving initial model parameters, performing ray tracing forward modeling and calculating a chromatography equation, and solving the chromatography equation to obtain the model updating quantity of the iteration; and finally, updating the model and entering next iteration until the preset convergence requirement is met. The process is fully automated without additional human intervention.

As shown in fig. 13, the steps of the present invention are as follows:

the method comprises the following steps of automatically picking up the slope of a same phase axis.

The prestack data is divided into a common shot point gather, a common survey point gather and a common offset gather, the optimal local homophase axis slope of each sampling point of the three gather data is obtained by scanning by utilizing a structure tensor analysis technology and is respectively recorded as P_R，P_S，P_M。

For the estimation of the local slope of the seismic event, the invention adopts the calculation method of the local structure tensor (please refer to Lucas J. van Vlietan PietW. Verbeek.1995, EstimatorsforOrientationand dAnisotrypinDigizedImages). And setting I as a two-dimensional seismic image, wherein a structure tensor representing spatial direction information in the two-dimensional image I is defined by an image gradient value, the structure tensor represents the change direction of a region and the variation along the change direction, and seismic stratum textures and fault textures are determined by the variation relation of azimuth information of local points. Introducing a Gaussian function blurs local details so that the structure tensor highlights the complexity of the signal in the region. For two-dimensional images, the structure tensor is a 2 x 2 matrix:

$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 is_xAnd g_yRepresenting the gradient of the seismic image in both the horizontal and vertical directions,<·>representing a two-dimensional gaussian smooth filter. For a semi-positive definite matrix G, eigenvalues and eigenvectors may be obtained by solving | G- λ I | ═ 0:

$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)$

λ₁: maximum eigenvalue, tensor energy in the first eigentensor direction v₁The energy of (a).

λ₂: minimum eigenvalue, tensor energy in the second eigentensor direction v₂The energy of (a).

(λ₁-λ₂)/λ₁: linearity, reflecting the consistency of local directions.

The feature vector describes the directionality of the local linear structure of the image, the feature vector v for each point of the image₁Normal to the main structural direction of the image, eigenvector v₂Parallel to the main structural direction of the image.

In the invention, a certain point in the seismic data before the superposition is subjected to structural tensionVector v solved by quantity technique₂Is v is₂＝(e₁,e₂) Then e₂/e₁I.e. the local in-phase axis slope at that point. P_R，P_S，P_MRespectively corresponding to the same-phase axis slopes in the pre-stack seismic data common shot gather, common survey point gather and common offset gather.

And the quality monitoring and filtering method comprises the following steps of 'homophasic axis slope coordinate transformation identity' and picking data points.

Calculating the slope P in the shot trace set for each reflection coaxial point of the reflection horizon in the common offset profile by utilizing the structure tensor analysis technology according to travel time information and shot point positions_RSlope P of trace concentration of the wave-detecting point_SSlope P in common offset profile_MMaking a comparison if P_S+P_R＝P_MError of this identity (absolute error for the error)₁＝P_S+P_R-P_MOr relative error₂＝(P_S+P_R-P_M)/P_MDescription) satisfies a predetermined error range (the range (i.e., absolute error)₁And relative error₂The range of the error is given by a user in a self-defined mode, the two parameters depend on the quality of the seismic data, the upper error limits are given to be smaller when the quality of the seismic data is higher, otherwise, the given value is relatively larger, the user can monitor the quality by displaying the picked data points, and if the number of the picked data points is too small, the upper error limits are amplified and picked again), the picking result is effectively picked, and the picking result is stored. The input initial modulus component in fig. 13 is given by the user, and its specific elements are defined by equation (4).

And establishing and solving a regularized analytic equation.

The stereo chromatographic inversion belongs to a speed estimation method of data fitting class, namely a model vector m defined by formula (4) is found, so that forward calculation is carried out to obtain data d_modThe error of the picked-up data d defined by the equation (4) is minimized. f represents the nonlinear positive operator of the stereo chromatography, i.e. reflection point to shot point, detection respectivelyInitial ray tracing of the point direction. The specific criterion for judging whether the stereo chromatography inversion measurement model is correct is as follows: considering the reflection of a shot-checking pair, if the coordinate, the incident angle, the reflection angle and the velocity model of the reflection point are correct, the reflection point emits rays to the shot point direction and the wave detection point direction respectively according to the incident angle and the reflection angle, and when the two transmission rays reach the observation surface, the two transmission rays reach the correct shot point and wave detection point respectively just in the correct emitting direction and the correct two-way travel.

The three-dimensional chromatography data component and the model component adopted by the invention are respectively as follows:

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)

the data in the formula (3) are n picked up respectively_dataAnd (3) data points (namely the result of the step (2)), wherein each data point is specifically a shot point abscissa, a demodulator probe abscissa, a horizontal component of a slowness vector at a shot point, a horizontal component of a slowness vector at a demodulator probe and a two-way travel time.

Similarly, each data point corresponds to a reflection/diffraction process, in the model corresponding to a set of model parameters (x) of equation (4)₀,z₀,θ_s,θ_r)_iI.e. the abscissa, ordinate, angle of incidence, angle of reflection of the reflection point.

Measuring data fitting errors by using two norms, and solving the three-dimensional chromatography inversion as an extreme value solving problem of the following functional:

$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 matrixAlso known as a data covariance matrix, serves to weight the different data components.

Due to inherent ill-conditioned nature and multiple solution of the inverse problem, the invention introduces additional model prior constraint for ensuring the reliability of the inversion result, and the error functional of the transformation formula (5) 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)$

whereinIs a small positive real number (The three parameters are distributed between 0 and 1, the three parameters also need to be given by a user through test customization, and when the initial model (the form is defined by formula (4)) of inversion is better, the given parameters are smaller when the data quality is higher, and conversely, the given parameters are larger, the three parameters are used for weighting the weight between the data fitting item and the regularization item; m is_rFor a given prior model, the form of the prior model is an inversion model component defined by a formula (4), the prior model needs to be given by user definition when the inversion is implemented, the prior model is given by a plurality of methods, and the prior model m is taken by the invention_r＝0。The introduction of this term makes the amount of model update per iteration small for the damping term.Andconstraining the smoothness of the velocity model in the transverse and longitudinal directions, respectively, D₁、D₂Respectively, the updating amount of the B spline base function coefficient is in the transverse directionFor a sequence of 7 numbers, the first order difference operator is defined as the matrix of 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),$

similarly, for a sequence of N points, the first order difference operator is a (N-1) XN matrix in the form of the above equation.

According to the local optimization theory of inversion, the optimization problem described in formula (6) can be equivalent to iterative solution of the following linear equation system, and the aim of convergence of inversion solution is achieved through multiple iterations, wherein the equation system (7) is called a chromatography equation of stereo chromatography 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) Where Δ d is the inverted data residual, Δ d (m) d-f (m); i is a unit matrix with the same dimension as the model spaceWherein N is_modelIs the dimension of the model component used in the stereo tomographic inversion, F is the Frechet derivative matrix, which is an N_data×N_modelIn which N is_dataIs the dimensionality of the data components used for stereo tomographic inversion, D is a data preconditions matrix used to weight the different dimensional data space components of the stereo tomographic inversion, for each data point (s, r, p)_sx,p_rx,t)_iD is defined by the form:

wherein the invention provides $({\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 a transmission ray tracing operator (mature in the field of geophysical prospecting for oil and gas), and m is a model component defined by formula (4); f is in m_nA Frechet derivative matrix of values, whereinF_ijThe elements in the ith row and the jth column of the Frechet derivative matrix are expressed as follows:

the calculation of the Frechet derivative of the stereo tomography is based on the ray perturbation theory (Fara & Madariaga,1987) and the properties of the ray propagation matrix (Gilbert & Bacus,1996), and the invention deduces the Frechet derivative of the stereo tomography inversion under a ray center coordinate system in the coordinate system by considering the compact form of the ray perturbation theory under the coordinate system. In the positive problem of the stereo tomography, the reflection process of S (shot point) → X (reflection point) → R (detection point) is divided into two transmission processes of X → S and X → R, and therefore, the transmission process from the subsurface to the surface shown in fig. 14 may be considered when calculating the friechet derivative of the stereo tomography.

First, we discuss the Frechet derivative formula derivation for stereo chromatography. Obviously, in Frechet derivativesIn the specification, the following are:

$\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_RThe method comprises the following steps:

$\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 ofAndthe two methods are completely consistent and respectively correspond to the transmission process from the reflection point to the shot point and the direction from the detection point.

As shown in FIG. 14, consider a channel s₀To s₁In a transmission process ofAnd (3) fully differentiating two sides:

$\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 calculate from A (3)

The calculation of (b) can be obtained by ray perturbation theory. The ray perturbation theory essentially establishes a linear relationship between the perturbation amount of the observed ray field and the perturbation amount of the initial ray field and the velocity, and the linear relationship can be converted into (S, P)_S) The amount of disturbance of (a) and (x)₀,z₀,θ_sV) linear relationship of the disturbance variable.

The ray perturbation theory is simplest in form under a ray center coordinate system, so the ray perturbation theory is firstly applied under the ray center coordinate system, and then (delta S, delta P) is deduced according to the geometric relation between rectangular coordinates and ray center coordinates_S) And (Δ x)₀,Δz₀,Δθ_sΔ v) of the linear relation to obtain

The perturbation coordinate q and the corresponding slowness vector component p have the following components according to the ray perturbation theory:

$\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)$ called slave s₀Point propagation to s₁Ray of a pointThe propagation matrix is a 2 × 2 matrix, and is denoted 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₀))^Tand (Q)₂(s₁,s₀),P₂(s₁,s₀))^TThe initial value of each paraxial ray tracing system is (1,0)^TAnd (0,1)^TPhysically corresponding to the normalized plane wave seismic source and the normalized point seismic source, respectively. The paraxial ray tracing system under the ray center coordinate system is as follows:

$\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)$

the ray propagation matrix has the following properties:

∏(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 from A (4), A (5), A (7) and 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}$

considering the geometrical relationship between the radial center coordinate system and the rectangular coordinate system, there are:

$\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^ξThe horizontal coordinate in the rectangular coordinate system and the component of the slowness vector in the coordinate direction are shown as a.1, and theta and α are respectively included angles between the outgoing transmitted ray and the arrival on the ground surface and the vertical axis, A (10) and A (11) are substituted into A (9) to 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)$

the derivative can be determined from A (12) and A (13)Which corresponds to the one to be soughtAnd

the present invention utilizes the least square QR method (LSQR) (please refer to the literature: PaigeCC, SaundersMA. LSQR: Analgorithmfor sparse and paraSequestresses. ACMTransactionon sparse software,1982,8(1):43-71) to solve the matrix equation set (7) (where Δ m is an unknown quantity and others are known quantities), which is an iterative method that can efficiently solve large-scale sparse matrices in the least square sense.

The two steps of steps (1) and (2) provide a data component d required by the stereo tomography inversion, and the right-end term Δ d (m) ═ d-f (m) of the equation system (7) is calculated by using d in each iteration of the inversion, wherein f (m) is obtained by using the forward modeling of the model in the previous iteration.

An embodiment of the method of the invention is as follows:

fig. 1 is a real underground medium velocity model, the velocity model is 17 km long in the transverse direction and 3.72 km deep in the longitudinal direction, fig. 2 is (partial) data of a surface observation shot gather obtained by solving an acoustic wave equation through finite difference numerical simulation using the velocity model of fig. 1, fig. 3 is a graph in which all travel time and slope of the picked data points can be well attached to seismic data (a common offset gather with an offset of 2 km) by using partial inversion data points automatically extracted in steps (1) and (2), which shows that the automatic picking process of the invention has higher precision, and the picked data points are subsequently used as input in step (3). FIG. 4 is a velocity model for the initial constant gradient given in this example, with the shallowest and deepest velocities of the velocity model being 2800m/s and 6000m/s, respectively, and the velocity gradient being 0.860215s^-1FIG. 5 is a superimposed graph of the initial velocity model and the reflection point position, in which the black short line segment is obtained from the initialized reflection point position and the local formation dip converted from the initialized incidence and reflection angles, and the conversion formula isWherein theta is_sIs the angle of incidence, θ_rIs the angle of reflection, θ_dipIs the local formation dip angle, and the initial model is used as the output of the step (3)And (6) adding.

After the retrieval of the inverted data component is finished, initializing the model component, and then starting the iterative update of the inversion, wherein fig. 6 is a velocity model output after 6 iterations, and fig. 7 is the distribution of the positions of the reflection points output after 6 iterations, so that the shallow layer of the velocity model is effectively recovered, and the positions of the reflection points are gradually converged to the position of a real reflection interface. Fig. 8 and 9 are velocity model and reflection point position distribution diagrams after 17 iterations, respectively, and it can be seen that as the iteration number increases, the inversion gradually approaches to the direction of the real solution. The velocity model obtained in the 17 th iteration is used as an input to perform gaussian beam prestack depth migration, and the obtained migration overlay cross section is shown in fig. 10, which shows that the migration image is relatively focused. Fig. 11 is an angle domain common imaging point gather output by gaussian beam prestack depth migration, and fig. 12 is a partial enlarged view of fig. 11, which shows that the common imaging point gather is leveled, so that the velocity model obtained by inversion has higher precision, meets the requirement of migration imaging, and can provide a high-quality construction profile for subsequent seismic interpretation.

The invention adopts the structure tensor analysis to calculate the earthquake homophase slope with high precision, and combines the quality control condition of 'earthquake homophase slope coordinate transformation identity', thus realizing a full-automatic stereo chromatographic inversion technique. Besides the travel time of the reflected seismic wave, the technology introduces the slope of the same phase axis into the reflection seismic tomography, enhances the precision and stability of the tomography, does not need manual intervention in the whole process of the technology, has higher precision and automation degree, can efficiently provide an accurate underground velocity model for large-scale seismic exploration, and provides service for subsequent migration imaging and reservoir prediction.

The above-described embodiment is only one embodiment of the present invention, and it will be apparent to those skilled in the art that various modifications and variations can be easily made based on the application and principle of the present invention disclosed in the present application, and the present invention is not limited to the method described in the above-described embodiment of the present invention, so that the above-described embodiment is only preferred, and not restrictive.

Claims (6)

1. A full-automatic stereo chromatographic inversion method using local event slope is characterized in that: the method comprises the following steps:

the method comprises the steps of automatically picking up the slope of a same phase shaft;

⑵ quality monitoring is performed by using the same phase axis slope coordinate transformation identity, filtering and picking data points (S, R, P)_S,P_R,t)_pick；

Establishing and solving a regularization stereo chromatography equation.

2. The full-automatic tomographic inversion method using local event slope according to claim 1, wherein: the step (1) is realized by the following steps:

the prestack data is divided into a common shot point gather, a common survey point gather and a common offset gather, then the optimal local homophase axis slope of each sampling point of the three gather data is obtained by scanning by using a structure tensor analysis method and is respectively recorded as P_R，P_S，P_M。

3. The full-automatic tomographic inversion method using local event slope according to claim 2, wherein: the step (2) is realized by the following steps:

the same phase axis slope coordinate transformation identity is as follows:

P_S+P_R＝P_M

if the error of the homophase axis slope coordinate transformation identity satisfies the predetermined error range, the data points (S, R, P) are effectively picked up and stored_S,P_R,t)_pick；

Absolute error for said error₁＝P_S+P_R-P_MOr relative error₂＝(P_S+P_R-P_M)/P_MAnd (4) showing.

4. The full-automatic tomographic inversion method using local event slope according to claim 3, wherein: the regularized stereo chromatography equation in the step (3) is shown as the formula (7):

$\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 the formula, D is a data precondition matrix used for weighting data space components of different dimensions of the stereo tomography inversion;

f is in m_nThe Frechet derivative matrix of the value is N_data×N_modelIn which N is_dataIs the dimensionality, N, of the data component used for the stereo tomographic inversion_modelIs the dimensionality of the model components used for the stereo tomographic inversion;

i is a unit matrix with the same dimension as the model space

All positive real numbers between 0 and 1 are set by a user;

D₁、D₂respectively carrying out first-order difference operators on the updating quantity of the B spline basis function coefficient in the transverse direction and the longitudinal direction;

for each data point (s, r, p)_sx,p_rx,t)_iΔ d is the inverted data residual, Δ d (m) d-f (m), where:

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 nonlinear positive operator of the stereo chromatography, namely the initial ray tracing from the reflection point to the shot point and the wave detection point respectively.

5. The full-automatic tomographic inversion method using local event slope according to claim 4, wherein: the step (3) is accomplished by iteration, specifically as follows:

let the model component of the kth iteration be m_kThen, the process of the (k + 1) th iteration is:

(31) using the model component m_kForward calculated data component f (m)_k) Establishing the regularized stereo chromatography equation;

(32) solving the regularization three-dimensional chromatography equation to obtain the model updating quantity delta m of the (k + 1) th iteration_k+1；

(33) Carrying out model updating to obtain a model component m of the (k + 1) th iteration_k+1＝m_k+Δm_k+1；

(34) If the data residual meets the user-defined threshold, the step (35) is carried out, and if the data residual does not meet the user-defined threshold, the step (31) is returned to enter the next iteration;

(35) and exiting iteration and outputting a final model component.

6. The full-automatic tomographic inversion method using local event slope according to claim 5, wherein: the data residual in said step (34) is a two-norm of the residuals of the forward calculated data component and the picked-up data component.