CN105572728A - Reverse illumination speed inversion method based on least-square objective functional - Google Patents

Abstract

A reverse illumination speed inversion method based on a least-square objective functional is disclosed. The method comprises the following steps of calculating a reverse illumination matrix; according to top boundary and bottom boundary positions of a current layer position, using a bending degree of a co-phase shaft on an angle domain common imaging point gather to calculate a residual depth; under least-square sense, through searching a speed model, using a minimum difference of actual observation travel time and forward modeling travel time to construct the least-square objective functional. The method has advantages that the method can be compatible with current processing software and implementation is convenient; and calculating efficiency and calculating precision are high. In the method, the angle domain common imaging point gather without illusion is used to calculate the residual depth so that the calculating precision is high; ray tracing from a bottom interface of a reflecting layer to a top interface is used to calculate the reverse illumination matrix and reflection ray tracking with tedious calculating is not needed so that the calculating efficiency is high. By using the method, a requirement to initial speed model precision is low, an inversed non-linear degree is low and adaptability is high.

Description

Based on the antidromic illumination velocity inversion method of least square cost functional

Technical field

The invention belongs to oil-gas exploration seismic data process field, is that a kind of earthquake data before superposition that is applicable to is for realizing the method for Depth Domain velocity modeling.

Prior art

Along with the develop rapidly of seismic exploration technique and computer level, pre-stack depth migration has become main force's technology of complex hydrocarbon exploratory area high precision imaging.But, pre-stack depth migration is very strong to the susceptibility of velocity accuracy, if Depth Domain rate pattern is inaccurate, offset method accurate, advanced so more all can not carry out accurate imaging to underground structure, therefore, a good Depth Domain velocity modeling method seems most important.

At present, existing Depth Domain velocity modeling method mainly contains three kinds: one utilizes time domain speed to carry out time and depth transfer, but all there is larger error due to conversion formula, the general method as acquisition initial depth territory rate pattern; The second is full waveform inversion, and its modeling accuracy is high, but it is high to initial velocity accuracy requirement, in refutation process, nonlinear degree is strong, and needs to estimate source wavelet etc., and calculated amount is huge, the rarely seen example that should make good use of, is difficult to popularization and carries out large-scale commercial production; The third is migration velocity analysis, and its advantage is moderate accuracy, and counting yield can accept, and current most of business softwares all adopt the method to carry out Depth Domain velocity modeling.But migration velocity analysis uses offset domain common image gathers, produce in the process of this road collection and there is multipath imaging problem, larger to follow-up velocity modeling process influence.

Summary of the invention

The object of the invention is for the problems referred to above, propose a kind of antidromic illumination velocity inversion method based on least square cost functional, more and more higher requirement seismic processing technique proposed with the increase meeting oil-gas exploration and development difficulty.Pre-stack depth migration is as main force's technology of high precision imaging, and it key of large-scale promotion application can be the accurate foundation of Depth Domain rate pattern.

The present invention is under Least squares inversion framework, set up the cost functional of Depth Domain rate pattern inverting, adopt current reflective layer bottom boundary to interface, top antidromic illumination matrix by the residue depth assignment that produces in communication process to the grid of wave path process to obtain residual velocity, iterative inversion is repeatedly to obtain high accuracy depth domain rate pattern.The angle domain common image gathers that the present invention applies without illusion asks for the residue degree of depth, and reflection horizon bottom boundary asks for antidromic illumination matrix to the ray tracing at interface, top, simple to operate, counting yield and computational accuracy high, be beneficial to large-scale promotion application in actual production.

Technical scheme of the present invention comprises:

(1) antidromic illumination matrix is calculated, interface, top and the bottom boundary of current layer position is determined in pre-stack depth migration result, the starting point of ray tracing fixes on the bottom boundary of current layer, given ray exit direction, given fixed step size carries out ray tracing, and the path addition in each grid is obtained the total length of ray in each grid, total length in each grid is lined up smoothly according to grid, obtain the antidromic illumination vector of a ray, the antidromic illumination vector of each ray is lined up according to ray order, obtains antidromic illumination matrix;

(2) according to interface, top and the bottom boundary position of current layer position, utilize the degree of crook of lineups on angle domain common image gathers to calculate the residue degree of depth;

(3) under least square meaning, by finding a rate pattern, such that the difference of actual observation whilst on tour and forward simulation whilst on tour is minimum builds least square cost functional.

Such scheme comprises further:

The raypath computing formula of step (1) described calculating antidromic illumination matrix is:

$p_{i,j}={\∫}_{l_1}^{l_2}\left|\right|n_0[1+\frac{l}{v}(\mathrm{\λ}\·n_0)-\frac{\mathrm{\λl}}{v}-\frac{n_0}{2v^2}{l}^2[{\mathrm{\λ}}^2-{(\mathrm{\λ}\·n_0)}^2\left]\right||\mathrm{dl}---\left(1\right)$

Wherein, p _i,jbe the path in a grid, dl is fixed step size, l ₁, l ₂be respectively in grid and start and terminate path; n ₀be ray exit direction vector, v is grid local velocity, and λ is velocity gradient

The spread pattern of antidromic illumination matrix P is:

$P=\left(\begin{array}{cccccc}{p}_{\mathrm{1,1}}& p_{\mathrm{1,2}}& ...& p_{1,j}& ...& p_{1,N}\\ p_{\mathrm{2,1}}& p_{\mathrm{2,2}}& ...& p_{2,j}& ...& p_{2,N}\\ .& .& .& .& .& .\\ .& .& .& .& .& .\\ .& .& .& .& .& .\\ p_{i,1}& p_{i,2}& ...& p_{i,j}& ...& p_{i,N}\\ .& .& .& .& .& .\\ .& .& .& .& .& .\\ .& .& .& .& .& .\\ p_{M,1}& p_{M,2}& ...& p_{M,j}& ...& p_{M,N}\end{array}\right)---\left(2\right)$

Wherein, i ∈ [1, M] represents ray sequence number, and M is ray sum; J ∈ [1, N] represents grid numbering, and N is grid sum.

Step (2) described calculating residue degree of depth method: the reference mark of first selected velocity inverting, then utilizes the maximum energy value of common imaging gather road, energy cross-correlation method pick-up angles territory collection lineups, obtain the residue depth delta z of different angles:

$\left\{\begin{array}{c}\mathrm{energy}\left({\mathrm{\τ}}_{\max }\right)=\underset{\mathrm{\τ}=1}{\overset{{\mathrm{kd}}_{\mathrm{down}}-{\mathrm{kd}}_{\mathrm{up}}+1}{\max }}[\underset{z=d_{\mathrm{up}}}{\overset{d_{\mathrm{down}}}{\mathrm{\Σ}}}f\left(z\right)*f(z+\mathrm{\τ})]\\ \mathrm{\Δz}=z_{{\mathrm{\τ}}_{\max }}-z_0\end{array}\right.---\left(3\right)$

Wherein, energy (τ _max) be cross-correlation maximum energy value, τ is depth-sampling point amount of movement, and span is [1, kd _down-kd _up+ 1]; F (z) for the degree of depth be the offset data energy value at z place, the span of z is top, the end degree of depth [d _up, d _down]; Max [] is max function, z ₀for the lineups degree of depth that zero degree in angle domain common image gathers is corresponding.

Step (3) builds least square cost functional and is expressed from the next:

$E\left(s\right)=\frac{1}{2}{\left|\right|\mathrm{P\Δs}-2\mathrm{s\Δ}z\sin {\mathrm{\θ}}_1\cos {\mathrm{\θ}}_2\left|\right|}^2---\left(4\right)$

Wherein, E (s) is least square cost functional, and s represents grid slowness, is the inverse of Grid Velocity; Δ s represents slowness disturbance, for carrying out velocity inversion renewal; P represents antidromic illumination matrix, and Δ z represents the residue degree of depth; θ ₁for ray emergence angle, θ ₂for bottom boundary inclination angle, current layer position; || || the L2 mould of representation vector

Its Grid Velocity inverting more new formula is:

$\frac{1}{v_n}=\frac{1}{v_{n-1}}+\mathrm{\Δs}---\left(5\right)$

Wherein, v _n-1and v _nbe respectively the velocity vector that n-th (n>=1) inverting upgrades front and back.

Antidromic illumination velocity inversion method based on least square cost functional of the present invention, its technical advantage is mainly manifested in following three aspects:

The first, method and current process software can be good at compatibility and implement conveniently.The Depth Domain rate pattern that the time domain rate pattern that the method utilizes current conventional process software to provide is become by time and depth transfer inputs as initial velocity model, the pre-stack depth migration module when pre-processing software is utilized to obtain angle domain common image gathers, compatible good, it is convenient to implement;

The second, counting yield and computational accuracy high.The method application asks for the residue degree of depth without the angle domain common image gathers of illusion, and computational accuracy is high; Application reflection horizon bottom boundary asks for antidromic illumination matrix to the ray tracing at interface, top, and do not need to calculate loaded down with trivial details reflected ray and follow the trail of, counting yield is high;

Three, the requirement of method to initial velocity model precision is low, and the nonlinear degree of inverting is weak, strong adaptability.

Accompanying drawing explanation

Fig. 1 is based on the antidromic illumination velocity inversion method operational flowchart of least square cost functional

The complicated rate pattern of Fig. 2

The input initial velocity model of Fig. 3 complex model

Fig. 4 antidromic illumination inversion speed model

The pre-stack depth migration result that Fig. 5 inversion speed model is corresponding

The input initial velocity model of Fig. 6 real data

The rate pattern that the inverting of Fig. 7 real data antidromic illumination obtains

The pre-stack depth migration result that Fig. 8 Real data inversion rate pattern is corresponding

Embodiment

Below by specific embodiment, technical scheme of the present invention is described further.

The detailed technology operation steps of this invention as shown in Figure 1.This invention major technique key point is following three: (1), calculating antidromic illumination matrix; (2) the residue degree of depth, is calculated; (3), least square cost functional is built.

(1) antidromic illumination matrix is calculated

Interface, top and the bottom boundary of current layer position is determined in pre-stack depth migration result, the starting point of ray tracing fixes on the bottom boundary of current layer, given ray exit direction, given fixed step size carries out ray tracing, and the path addition in each grid is obtained the total length of ray in each grid, total length in each grid is lined up smoothly according to grid, obtain the antidromic illumination vector of a ray, the antidromic illumination vector of each ray is lined up according to ray order, obtains antidromic illumination matrix.Raypath computing formula is:

$p_{i,j}={\∫}_{l_1}^{l_2}\left|\right|n_0[1+\frac{l}{v}(\mathrm{\λ}\·n_0)-\frac{\mathrm{\λl}}{v}-\frac{n_0}{2v^2}{l}^2[{\mathrm{\λ}}^2-{(\mathrm{\λ}\·n_0)}^2\left]\right||\mathrm{dl}---\left(1\right)$

Wherein, p _i,jbe the path in a grid, dl is fixed step size, l ₁, l ₂be respectively in grid and start and terminate path; n ₀be ray exit direction vector, v is grid local velocity, and λ is velocity gradient.

The spread pattern of antidromic illumination matrix P is:

$P=\left(\begin{array}{cccccc}{p}_{\mathrm{1,1}}& p_{\mathrm{1,2}}& ...& p_{1,j}& ...& p_{1,N}\\ p_{\mathrm{2,1}}& p_{\mathrm{2,2}}& ...& p_{2,j}& ...& p_{2,N}\\ .& .& .& .& .& .\\ .& .& .& .& .& .\\ .& .& .& .& .& .\\ p_{i,1}& p_{i,2}& ...& p_{i,j}& ...& p_{i,N}\\ .& .& .& .& .& .\\ .& .& .& .& .& .\\ .& .& .& .& .& .\\ p_{M,1}& p_{M,2}& ...& p_{M,j}& ...& p_{M,N}\end{array}\right)---\left(2\right)$

Wherein, i ∈ [1, M] represents ray sequence number, and M is ray sum; J ∈ [1, N] represents grid numbering, and N is grid sum.

(2) the residue degree of depth is calculated

According to interface, top and the bottom boundary position of current layer position, utilize the degree of crook of lineups on angle domain common image gathers to calculate the residue degree of depth.First the reference mark (according to accuracy requirement, the amount doesn't matter) of selected velocity inverting, then utilizes the maximum energy value of common imaging gather road, energy cross-correlation method pick-up angles territory collection lineups, obtains the residue depth delta z of different angles:

$\left\{\begin{array}{c}\mathrm{energy}\left({\mathrm{\τ}}_{\max }\right)=\underset{\mathrm{\τ}=1}{\overset{{\mathrm{kd}}_{\mathrm{down}}-{\mathrm{kd}}_{\mathrm{up}}+1}{\max }}[\underset{z=d_{\mathrm{up}}}{\overset{d_{\mathrm{down}}}{\mathrm{\Σ}}}f\left(z\right)*f(z+\mathrm{\τ})]\\ \mathrm{\Δz}=z_{{\mathrm{\τ}}_{\max }}-z_0\end{array}\right.---\left(3\right)$

Wherein, energy (τ _max) be cross-correlation maximum energy value, τ is depth-sampling point amount of movement, and span is [1, kd _down-kd _up+ 1]; F (z) for the degree of depth be the offset data energy value at z place, the span of z is top, the end degree of depth [d _up, d _down]; Max [] is max function, z ₀for the lineups degree of depth that zero degree in angle domain common image gathers is corresponding.

(3) least square cost functional is built

Under least square meaning, can searching rate pattern be passed through, make the difference of actual observation whilst on tour and forward simulation whilst on tour minimum come establishing target functional:

$E\left(s\right)=\frac{1}{2}{\left|\right|\mathrm{P\Δs}-2\mathrm{s\Δ}z\sin {\mathrm{\θ}}_1\cos {\mathrm{\θ}}_2\left|\right|}^2---\left(4\right)$

Wherein, E (s) is least square cost functional, and s represents grid slowness, is the inverse of Grid Velocity; Δ s represents slowness disturbance, for carrying out velocity inversion renewal; P represents antidromic illumination matrix, and Δ z represents the residue degree of depth; θ ₁for ray emergence angle, θ ₂for bottom boundary inclination angle, current layer position; || || the L2 mould of representation vector.

Its Grid Velocity inverting more new formula is:

$\frac{1}{v_n}=\frac{1}{v_{n-1}}+\mathrm{\Δs}---\left(5\right)$

Wherein, v _n-1and v _nbe respectively the velocity vector that n-th (n>=1) inverting upgrades front and back.

Below for complex model and certain exploratory area real data processing procedure, the effect of this invention is described.

Fig. 2-Fig. 5 is the test result of this invention to complex model.Fig. 2 is real complicated rate pattern, and discrete grid block is horizontal, longitudinal sampling number is 801*638, and corresponding sampling interval is respectively 10m and 20m, and this model exists the complex structures such as high steep dip, tomography and thin layer.

This invention to the implementation step of complex model is:

First, the initial rate pattern of input shown in Fig. 3, antidromic illumination matrix is calculated: the bottom boundary starting point of ray tracing being fixed on ground floor according to summary of the invention (1), given ray exit direction and fixed step size carry out ray tracing, path in each grid is added and obtains the total length of ray in each grid, total length in each grid is lined up smoothly according to grid, obtain the antidromic illumination vector of a ray, the antidromic illumination vector of each ray is lined up according to ray order, obtains antidromic illumination matrix;

The residue degree of depth is calculated: according to 80 reference mark of data sampling number selected velocity inverting according to summary of the invention (2), then utilize the maximum energy value of common imaging gather road, energy cross-correlation method pick-up angles territory collection lineups, obtain the residue degree of depth of different angles;

Then the least square cost functional as shown in formula (4) is built, formula (5) is utilized to carry out velocity inversion renewal, ground floor iteration ends starts the inverting of the second layer, until all layer positions iteration ends exports final inversion speed model (as shown in Figure 4).Can find out: rate pattern (Fig. 4) same to true model (Fig. 1) goodness of fit that inverting obtains is very high, the display of pre-stack depth migration imaging result in Figure 5, visible inverse model precision is higher, can obtain high-quality imaging section for carrying out skew.

Fig. 6-Fig. 8 is the test result of this invention to real data.Real data discrete grid block is horizontal, longitudinal sampling number is 901*1001, and corresponding sampling interval is respectively 15m and 8m.Implementation step illustrates: the initial velocity model of input shown in Fig. 6, first antidromic illumination matrix is calculated: the bottom boundary starting point of ray tracing being fixed on ground floor, given ray exit direction and fixed step size carry out ray tracing, path in each grid is added and obtains the total length of ray in each grid, total length in each grid is lined up smoothly according to grid, obtain the antidromic illumination vector of a ray, the antidromic illumination vector of each ray is lined up according to ray order, obtains antidromic illumination matrix; Secondly the residue degree of depth is calculated: according to 90 reference mark of data sampling number selected velocity inverting, then utilize the maximum energy value of common imaging gather road, energy cross-correlation method pick-up angles territory collection lineups, obtain the residue degree of depth of different angles; Then the least square cost functional as shown in formula (4) is built, formula (5) is utilized to carry out velocity inversion renewal, ground floor iteration ends starts the inverting of the second layer, until all layer positions iteration ends exports final inversion speed model, display in the figure 7.Can see: from shallow-layer to deep layer, speed details is obtained for well portrays, and the position, pre-stack depth migration result middle level shown in Fig. 8 is clearly demarcated, tomography is clear, and the velocity inversion precision demonstrating this invention is higher.

Claims (4)

1., based on an antidromic illumination velocity inversion method for least square cost functional, it is characterized in that comprising:

(1) antidromic illumination matrix is calculated, interface, top and the bottom boundary of current layer position is determined in pre-stack depth migration result, the starting point of ray tracing fixes on the bottom boundary of current layer, given ray exit direction, given fixed step size carries out ray tracing, and the path addition in each grid is obtained the total length of ray in each grid, total length in each grid is lined up smoothly according to grid, obtain the antidromic illumination vector of a ray, the antidromic illumination vector of each ray is lined up according to ray order, obtains antidromic illumination matrix;

(2) according to interface, top and the bottom boundary position of current layer position, utilize the degree of crook of lineups on angle domain common image gathers to calculate the residue degree of depth;

(3) under least square meaning, by finding a rate pattern, such that the difference of actual observation whilst on tour and forward simulation whilst on tour is minimum builds least square cost functional.

2. the antidromic illumination velocity inversion method based on least square cost functional according to claim 1, is characterized in that, the raypath computing formula of step (1) described calculating antidromic illumination matrix is:

$p_{i,j}={\∫}_{l_1}^{l_2}\left|\right|n_0[1+\frac{1}{v}(\mathrm{\λ}\·n_0)]-\frac{\mathrm{\λ}1}{v}-\frac{n_0}{2v^2}{l}^2[{\mathrm{\λ}}^2-{(\mathrm{\λ}\·n_0)}^2]\left|\right|d1---\left(1\right)$

Wherein, p _i,jbe the path in a grid, dl is fixed step size, l ₁, l ₂be respectively in grid and start and terminate path; n ₀be ray exit direction vector, v is grid local velocity, and λ is velocity gradient

The spread pattern of antidromic illumination matrix P is:

$P=\left(\begin{array}{cccccc}{p}_{\mathrm{1,1}}& p_{\mathrm{1,2}}& ...& p_{1,j}& ...& p_{1,N}\\ p_{\mathrm{2,1}}& p_{\mathrm{2,2}}& ...& p_{2,j}& ...& p_{2,N}\\ .& .& .& .& .& .\\ .& .& .& .& .& .\\ .& .& .& .& .& .\\ p_{i,1}& p_{i,2}& ...& p_{i,j}& ...& p_{i,N}\\ .& .& .& .& .& .\\ .& .& .& .& .& .\\ .& .& .& .& .& .\\ p_{M,1}& p_{M,2}& ...& p_{M,j}& ...& p_{M,N}\end{array}\right)---\left(2\right)$

Wherein, i ∈ [1, M] represents ray sequence number, and M is ray sum; J ∈ [1, N] represents grid numbering, and N is grid sum.

3. the antidromic illumination velocity inversion method based on least square cost functional according to claim 2, it is characterized in that, step (2) described calculating residue degree of depth method: the reference mark of first selected velocity inverting, then utilize the maximum energy value of common imaging gather road, energy cross-correlation method pick-up angles territory collection lineups, obtain the residue depth delta z of different angles:

$\left\{\begin{array}{c}\mathrm{energy}\left({\mathrm{\τ}}_{\max }\right)=\underset{\mathrm{\τ}=1}{\overset{{\mathrm{kd}}_{\mathrm{down}}-{\mathrm{kd}}_{\mathrm{up}}+1}{\max }}[\underset{z=d_{\mathrm{up}}}{\overset{d_{\mathrm{down}}}{\mathrm{\Σ}}}f\left(z\right)*f(z+\mathrm{\τ})]\\ \mathrm{\Δz}=z_{{\mathrm{\τ}}_{\max }}-z_0\end{array}\right.---\left(3\right)$

Wherein, energy (τ _max) be cross-correlation maximum energy value, τ is depth-sampling point amount of movement, and span is [1, kd _down-kd _up+ 1]; F (z) for the degree of depth be the offset data energy value at z place, the span of z is top, the end degree of depth [d _up, d _down]; Max [] is max function, z ₀for the lineups degree of depth that zero degree in angle domain common image gathers is corresponding.

4. the antidromic illumination velocity inversion method based on least square cost functional according to claim 1,2 or 3, is characterized in that, step (3) builds least square cost functional and is expressed from the next:

$E\left(s\right)=\frac{1}{2}{\left|\right|\mathrm{P\Δs}-2\mathrm{s\Δ}z\sin {\mathrm{\θ}}_1\cos {\mathrm{\θ}}_2\left|\right|}^2---\left(4\right)$

Wherein, E (s) is least square cost functional, and s represents grid slowness, is the inverse of Grid Velocity; Δ s represents slowness disturbance, for carrying out velocity inversion renewal; P represents antidromic illumination matrix, and Δ z represents the residue degree of depth; θ ₁for ray emergence angle, θ ₂for bottom boundary inclination angle, current layer position; || || the L2 mould of representation vector

Its Grid Velocity inverting more new formula is:

$\frac{1}{v_n}=\frac{1}{v_{n-1}}+\mathrm{\Δs}---\left(5\right)$

Wherein, v _n-1and v _nbe respectively the velocity vector that n-th (n>=1) inverting upgrades front and back.