阅读说明：本技术 基于深度加权的地震斜率层析成像方法 (Depth weighting-based seismic slope tomography method ) 是由 张倩锋 张建中 杨华臣 于 2021-06-11 设计创作，主要内容包括：本发明涉及一种基于深度加权的地震斜率层析成像方法,属于地震资料处理的速度建模技术领域。其解决了由于复杂构造区域深部反射同相轴微弱,有效信号较少而导致深部速度建模效果差的问题,在保证浅层速度精度的同时,提高了深部地层速度的反演精度。本发明包括如下步骤：对地震数据进行预处理；拾取反射波的走时和斜率,形成观测数据；建立初始速度模型和初始射线段参数,形成模型空间；射线段参数和速度模型联合反演；对反演结果进行质量控制,如果反演结果不收敛,则继续进行反演迭代,如果反演结果收敛,则输出速度模型和射线段参数。(The invention relates to a depth weighting-based seismic slope tomography method, and belongs to the technical field of velocity modeling of seismic data processing. The method solves the problem of poor deep velocity modeling effect caused by weak deep reflection phase axis and less effective signals in a complex construction region, and improves the inversion accuracy of deep stratum velocity while ensuring the accuracy of shallow velocity. The invention comprises the following steps: preprocessing the seismic data; picking up travel time and slope of the reflected wave to form observation data; establishing an initial velocity model and initial ray segment parameters to form a model space; joint inversion of ray segment parameters and a velocity model; and (4) performing quality control on the inversion result, if the inversion result is not converged, continuing inversion iteration, and if the inversion result is converged, outputting a velocity model and ray segment parameters.)

1. A seismic slope tomography method based on depth weighting is characterized by comprising the following steps:

the method comprises the following steps: preprocessing seismic data;

step two: picking up the travel time and the slope of the reflected wave to form observation data:

s1: picking up reflected wave travel time, a wave detection point end and a shot point end slope in a common shot domain gather and a common wave detection point domain gather;

s2: obtaining the position and elevation information of shot points and wave detection points from a field observation system file, and finally obtaining observation data d as follows:

， （1）

wherein the content of the first and second substances,x _s 、z _s respectively the position and elevation of the shot point,p _sx is the slope of the horizontal direction of the shot point end,x _r 、z _r respectively the position and elevation of the demodulator probe,p _rx for the horizontal direction slope at the end of the demodulator probe,t _sr in order to reflect the wave travel time,Nsubscript to pick up total number of observations obtainednIs shown asn(ii) observed data;

step three: establishing an initial ray segment parameter and an initial velocity model, and obtaining a model space m as follows:

， （2）

wherein m is _ray As initial ray segment parameter, m _v An initial velocity model;

s1: the initial ray segment parameters are expressed as:

， （3）

wherein the content of the first and second substances,x _c 、z _c in order to be able to reflect the position of the point,φ _s andφ _r respectively the angles of the ray emitted from the reflection point to the shot point and the wave detection point,t _s andt _r respectively, the single-pass travel time of the ray from the reflection point to the shot point and the wave detection point,Nis the number of ray pairs, subscriptnIs shown asnA pair of rays;

s2: the initial velocity model is expressed as:

， （4）

wherein M is the total number of speed nodes describing the whole discrete speed model;v _m =(v ₀ +kz)km/swherein, in the step (A),v ₀ as the velocity of the earth's surface,kin order to have a velocity vertical gradient,zis depth;

step four: joint inversion of ray segment parameters and a velocity model:

s1: the depth weighted slope tomographic objective function is established as follows:

， （5）

wherein d is _cal To calculate the data, d _obs To observe data, m _ray As a parameter of the ray segment, m _v Is a velocity model, λ is the damping coefficient, L is the Laplace operator of velocity perturbation, m _{v_prior} Being a prior velocity model, C_dIs the inverse of the data covariance matrix, C_mFor the inverse of the model covariance matrix, superscriptTRepresents a transpose of a matrix;

s2: the nonlinear problem of the formula (5) is linearized, and an inversion equation set based on depth-weighted slope tomography is derived and represented as follows:

， （6）

wherein, C_dIs the inverse of the data covariance matrix, C_mIs the inverse of the model covariance matrix, w _{_vd}Is a depth weighting coefficient matrix, G _{_rayn}Is a ray segment parameter kernel function matrix, G _{_vn}For the velocity kernel function matrix, Δ n _ray For the unknowns, Δ n, of the ray parameters to be solved for _v For the unknown quantity of velocity to be solved, Δ d is the difference between the calculated data and the observed data,μa weight factor between the control ray segment parameter kernel function and the speed kernel function;

s3: the kernel function matrix in the formula (6) is obtained by solving a paraxial ray first-order approximation formula;

s4: solving the formula (6) by adopting a least square QR algorithm, and obtaining the model update quantity as follows:

， （7）

in the formula,. DELTA.m _ray For ray segment parameter update quantity, Δ m _v In order to update the volume of the velocity model, _d_ _v a correction factor matrix that is a weighting factor for the velocity nodes,μa weight factor between the control ray segment parameter kernel function and the speed kernel function;

s5: updating the model space;

step five: and (4) performing quality control on the inversion result, if the inversion result is not converged, continuing inversion iteration, and if the inversion result is converged, outputting a velocity model and ray segment parameters.

2. The depth-weighted seismic slope tomography method as claimed in claim 1, wherein in the fourth step, the velocity kernel is weighted, and the depth-weighted inversion equation set is established as follows:

， （8）

wherein, C_dIs the inverse of the data covariance matrix, C_mIs the inverse of the model covariance matrix, w _{_vd}Is a depth weighting coefficient matrix, G _{_rayn}Is a ray segment parameter kernel function matrix, G _{_vn}For the velocity kernel function matrix, Δ n _ray For the unknowns, Δ n, of the ray parameters to be solved for _v For the unknown quantity of velocity to be solved, Δ d is the difference between the calculated data and the observed data,μa weight factor between the control ray segment parameter kernel function and the speed kernel function; the element value calculation formula of the depth weighting matrix is as follows:

， （9）

wherein the content of the first and second substances,w _im _ _v is shown asiData tomThe depth weighting coefficients of the kernel functions of the individual velocity nodes,k _m is shown asiA reflected ray hasTo closemThe depth of each speed node is greater than the depth of each speed node,N _i is shown withiThe total number of velocity nodes associated with each reflected ray,αis a constant; obtaining a correction factor matrix of velocity node weighting coefficients from depth weighting coefficients and data from different depths _d_ _v This is a diagonal matrix whose elemental values are expressed as:

， （10）

in the formula (I), the compound is shown in the specification, _m is as followsmThe velocity correction factor for each velocity node,w _im is as followsiData tomThe depth weighting factor for each velocity node,N _m to pass throughmTotal number of ray pairs for each velocity node.

Technical Field

The invention relates to a depth weighting-based seismic slope tomography method, and belongs to the technical field of velocity modeling of seismic data processing.

Background

Seismic slope tomography was proposed by Fred ric Billette in 1998, which combines the travel time of the reflected waves and the slope information of the local coherence events to build a model of the subsurface velocity and estimate the location of the reflection points. Compared with the reflected wave time-lapse tomography, the slope tomography introduces slope information to constrain a ray propagation path, which is beneficial to reducing inversion multi-solution; meanwhile, the method does not need to establish the one-to-one correspondence relationship between the continuous reflection homophase axis and the underground interface, only needs to pick up the travel time and the slope of the local coherent homophase axis, and is more suitable for low signal-to-noise ratio data processing in a complex structure area. The slope tomography solves the problem that geological horizons need to be divided in advance, so that the method has higher practical value and speed modeling effect than time-lapse tomography, and is successfully applied to actual data processing at present.

For actual observation of seismic data, on one hand, the energy of seismic waves from deep strata is weaker due to the spherical diffusion and medium absorption and attenuation effects of the seismic waves; on the other hand, the signal-to-noise ratio of deep seismic data from a region with complex and various geological structures and large speed change is often low, and for example, in a hidden type buried hill under the control of multi-stage structure motion in the middle-ancient world, the in-screen homomorphic axis of the buried hill is unclear and is quite disordered. Therefore, the effective deep reflection travel time and slope data which can be picked up are often less, so that the deep stratum ray coverage is insufficient, and the inversion speed effect is poor. For seismic data with less deep effective reflection data, an accurate deep velocity field cannot be established by the conventional slope tomography method, and particularly, when the difference between the initial velocity and the real velocity is large, the accuracy of the established depth velocity field is very low, and the requirement of velocity modeling cannot be met.

Disclosure of Invention

Aiming at the defects in the prior art, the invention provides a seismic slope tomography method based on depth weighting.

The invention is realized by adopting the following technical scheme: the invention relates to a depth weighting-based seismic slope tomography method, which comprises the following steps:

the method comprises the following steps: preprocessing seismic data;

step two: picking up the travel time and the slope of the reflected wave to form observation data:

s1: picking up reflected wave travel time, a wave detection point end and a shot point end slope in a common shot domain gather and a common wave detection point domain gather;

s2: obtaining the position and elevation information of shot points and wave detection points from a field observation system file, and finally obtaining observation data d as follows:

wherein x is_s、z_sAre respectively asLocation and elevation of shot, p_sxIs the horizontal slope, x, of the shot point end_r、z_rRespectively the position and elevation of the probe point, p_rxTo detect the horizontal slope at the point end, t_srWhen the time is the reflected wave travel time, N is the total number of the observed data obtained by picking up, and a subscript N represents the nth observed data;

step three: establishing an initial ray segment parameter and an initial velocity model, and obtaining a model space m as follows:

m＝(m_ray，m_v)， (12)

wherein m is_rayAs initial ray segment parameter, m_vAn initial velocity model;

s1: the initial ray segment parameters are expressed as:

wherein x is_c、z_cIn order to be able to reflect the position of the point,andrespectively the angle of ray emitted from the reflection point to the shot point and the detection point, t_sAnd t_rRespectively, the single-pass travel time of the ray from the reflection point to the shot point and the wave detection point, wherein N is the number of ray pairs, and a subscript N represents the nth ray pair;

s2: the initial velocity model is expressed as:

wherein M is the total number of speed nodes describing the whole discrete speed model; v. of_m＝(v₀+ kz) km/s, where v₀The earth surface velocity, k the velocity vertical gradient and z the depth;

step four: joint inversion of ray segment parameters and a velocity model:

s1: the depth weighted slope tomographic objective function is established as follows:

wherein d is_calTo calculate the data, do_bsTo observe data, m_rayAs a parameter of the ray segment, m_vIs a velocity model, λ is the damping coefficient, L is the Laplace operator of velocity perturbation, m_{v_prior}Being a prior velocity model, C_dIs the inverse of the data covariance matrix, C_mThe inverse of the model covariance matrix is used, and the superscript T represents the transposition of the matrix;

s2: the nonlinear problem of the formula (5) is linearized, and an inversion equation set based on depth-weighted slope tomography is derived and represented as follows:

wherein, C_dIs the inverse of the data covariance matrix, C_mIs the inverse of the model covariance matrix, w_d__vIs a depth weighting coefficient matrix, G_n__rayIs a ray segment parameter kernel function matrix, G_n__vFor the velocity kernel function matrix, Δ n_rayFor the unknowns, Δ n, of the ray parameters to be solved for_vThe unknown quantity to be solved for the speed is delta d, the difference between the calculated data and the observed data is delta d, and mu is a weight factor between a control ray section parameter kernel function and a speed kernel function;

s3: the kernel function matrix in the formula (6) is obtained by solving a paraxial ray first-order approximation formula;

s4: solving the formula (6) by adopting a least square QR algorithm, and obtaining the model update quantity as follows:

in the formula,. DELTA.m_rayFor ray segment parameter update quantity, Δ m_vAs a velocity modelThe amount of the update is such that,a correction factor matrix of the velocity node weighting coefficient is used, and mu is a weighting factor between the control ray section parameter kernel function and the velocity kernel function;

s5: updating the model space;

step five: and (4) performing quality control on the inversion result, if the inversion result is not converged, continuing inversion iteration, and if the inversion result is converged, outputting a velocity model and ray segment parameters.

The invention has the beneficial effects that: by adopting the seismic slope tomography method based on depth weighting, the shallow velocity precision is ensured, and meanwhile, the deep velocity modeling effect is improved, so that the integrally established velocity model quality is improved, and an effective method is provided for seismic data velocity modeling with low signal-to-noise ratio and fewer deep effective signals.

Drawings

FIG. 1 is a flow chart of the present invention;

FIG. 2 is a diagram of a true velocity model of the present invention;

FIG. 3 is a true velocity model smoothing diagram of the present invention;

FIG. 4 is a diagram of a velocity model built for non-depth weighted slope tomography;

FIG. 5 is a diagram of a velocity model established by the present invention;

fig. 6 is a graph of velocity curves at 7km versus non-depth-weighted slope tomography and initial velocity curves established by the present invention.

Detailed Description

In order to make the purpose and technical solution of the present invention more apparent, the present invention is further described in detail below with reference to the accompanying drawings.

As shown in the flow chart of FIG. 1, the depth weighting-based seismic slope tomography method of the invention comprises the following steps:

the first embodiment is as follows:

the theoretical model test of the present invention is explained and illustrated below with reference to specific embodiments.

In order to further explain the realization idea and the realization process of the method and prove the effectiveness of the method, a theoretical model is used for testing and is compared with the result of the seismic slope tomography method without depth weighting.

S1: and taking the theoretical speed model as a real speed model. The real speed model has a transverse width of 14km and a depth of 4 km. The discretization was performed using a square cell grid, with a grid size of 50m, see fig. 2, and for comparison, the theoretical velocity model is shown smoothly in fig. 3.

S2: an observation system: the scattering points are uniformly distributed on each speed interface from each scattering pointDirections of 10 °, 20 °, 30 ° and 35 ° are emitted toward the earth surface, and the rays are traced to the earth surface, so that 6876 sets of observation data are obtained, wherein the deep data are small.

S3: and establishing an initial speed model v which is 2 km/s.

S4: and establishing initial ray segment parameters, and optimizing inversion ray segment parameters on the initial velocity model.

S5: joint inversion of the velocity model and ray segment parameters; taking the initial velocity model and the optimized ray segment parameters as initial conditions, performing ray tracing forward, and continuously iterating and solving an inversion equation set to obtain an updated velocity model and ray segment parameters; s6: and judging whether the inversion result meets the precision requirement, if not, performing mesh unit subdivision, if so, outputting a final speed model and ray parameters, and finally, inverting and establishing the speed as shown in figure 5.

Fig. 4 is a velocity model built by depth-weighted slope tomography, and fig. 6 is a comparison graph of a velocity curve built by the invention with an inversion velocity curve and an initial velocity curve which are not depth-weighted at the position where x is 6 km. Comparing fig. 4, 5 and 6, it can be seen that the depth-weighted seismic slope tomography method improves the deep velocity, fits better with the real velocity, and improves the modeling effect of the deep velocity.

The seismic slope tomography method based on depth weighting increases the influence of deep data on deep stratum velocity modeling, and improves the deep velocity modeling effect while ensuring the precision of shallow velocity. The method has the advantages of simple calculation, easy realization, less time consumption, strong adaptability and high reliability of the inversion result.

The invention can be widely applied to the seismic slope tomography occasions.

The above description is only for the purpose of illustrating the preferred embodiments of the present invention and is not to be construed as limiting the invention, but rather the intention is to cover all modifications, equivalents, and alternatives falling within the spirit and scope of the invention as defined by the appended claims.