CN102116870A - Elastic wave gaussian beam pre-stack depth migration technology - Google Patents

Abstract

The invention provides an elastic wave gaussian beam pre-stack depth migration technology, which comprises the following steps of: 1) according to a selected gaussian beam primary width, dividing a single slot record into a series of stacked areas by a gaussian window, wherein the center of each area is a beam center; 2) decomposing a multi-component seismic record into different wave-formed local plane wave components by utilizing local slant stack; 3) respectively test-firing the elastic dynamic gaussian beams from hypocenter and beam centers and calculating the travel time and amplitude information of each gaussian beam in an effective range; 4) according to a PP wave and PS wave-deduced imaging formula, calculating an corresponding image value of each beam center by utilizing the travel time, the amplitude information of the gaussian beam, the polar information and surface local inclined information; and 5) repeating the calculations of steps 3) and 4), thereby acquiring an imaging result of the single slot, wherein the final migration result is the stack of all imaging results of the single slot. The method can be used for repeatedly imaging the multiple waves, is higher in precision, has calculation efficiency and flexibility and is a multi-component seismic imaging method having higher application value.

Description

Elastic wave Gaussian beam pre-stack depth migration technology

Affiliated technical field

The invention belongs to the multi-component seismic data process field, is a kind of Depth Domain formation method that is applicable to the multiwave multicomponent earthquake data.

Background technology

Along with the raising of physical prospecting technical equipment, many ripples (many components) seismic prospecting is subjected to the attention of industry member day by day, and has progressed into suitability for industrialized production.The multi-component seismic technology has demonstrated special advantages and huge application potential in the fine description of the petroleuon-gas prediction of nonuniformity reservoir, hydrocarbon-bearing pool and dynamic monitoring.Yet the existing formation method that is applicable to multi-component earthquake data is perfect not enough at present, often has the contradiction between imaging precision and the counting yield.

In traditional multi-component earthquake data processing procedure, often at first geological data is divided into vertical (P) ripple and horizontal stroke (S) ripple by the wave field separation method, utilize traditional sound wave offset imaging method then, imaging is not carried out in P ripple and S wavelength-division after separating.Yet above-mentioned disposal route does not make full use of the vectorial property of elasticity wave field; In addition, often the energy of different wave modes can not be separated fully in the process of wave field separation, non-remaining type wave energy often causes the much noise in the imaging results to be crosstalked, and then has a strong impact on imaging effect.

Based on the elastic wave seismic imaging method of vector wave field extrapolation, be the effective way that addresses the above problem.Because conventional one way wave imaging method can't correctly be described the elastic wave propagation process of coupling, existing elastic wave formation method all carries out in time domain, it roughly can be divided into two classes: a class is the elastic wave reverse-time migration, it is based on the elastic wave wave equation, multi-component earthquake data is carried out the inverse time extrapolation as boundary condition, and image-forming condition or simple crosscorrelation image-forming condition were asked for the imaging value when utilization excited then.Sun and McMechan (1986) at first are applied to reverse-time migration in the elastic wave skew, after this, Chang and McMechan (1994), Botelho and Stoffa (1991), Yan and Sava (2008) and Du Qizhen and Qin Tong people such as (2009) promote it and be applied in the migration imagings of three-dimensional and anisotropic medium.Another kind of is elastic wave Kirchhoff skew, its when calculating the walking of this type ripple and transformed wave, amplitude and polarity information, track superposes to the wave field energy and asks for the imaging value during subsequently along the walking of selected wave mode.Kuo and Dai (1984) have at first proposed this method based on the far-field approximation of displacement Green function in the homogeneous isotropic medium.After this, Sena and Toksoz (1993), Hokstad (2000) and Druzhinin people such as (2003) are developed it.Above-mentioned two class methods, have separately advantage and deficiency, though its Elastic Wave reverse-time migration has high imaging precision, its calculated amount is very huge, and its imaging results insoluble reversal of poles problem during often with a large amount of low-frequency disturbance and transformed wave imaging.Elastic wave Kirchhoff skew, though have flexible and efficient characteristics, it is difficult to multiple reflection to cause its imaging precision to underground complex geological structure not high to carrying out imaging.

Summary of the invention

The object of the present invention is to provide a kind of elastic wave Gaussian beam pre-stack depth migration technology, it is a kind of elastic wave seismic imaging method of taking into account precision and efficient, with the multi-component earthquake data is input, utilize the amplitude of elastic oscillation mechanics Gaussian beam, when walking and polarity information carry out the continuation of vector wave field, and utilize the vector projection of wave field to extract the imaging value of ripple in length and breadth.

Technical scheme that the present invention adopts is: elastic wave Gaussian beam pre-stack depth migration technology, and performing step is:

(1) according to the Gaussian beam original width of selecting, by Gaussian window single shot record is divided into a series of overlapping areas, each regional center is a beam center;

(2) utilize the local dip stack to be decomposed into the part plan wave component of different wave modes the multi-component seismic record;

(3) respectively from focus and beam center trial fire elastokinetics Gaussian beam, when calculating walking in every Gaussian beam effective range, amplitude information;

(4) according to the imaging formula of derivation PP ripple and PS ripple, when utilizing the walking of Gaussian beam, the local obliquity information on amplitude, polarity information and the face of land calculates the imaging value of certain corresponding beam center;

(5) repeat the computing in (3), (4) step, obtain the imaging results of single big gun, final migration result is the stack of all single big gun imaging results.

Wherein the computing method in the step (3) are:

On the described inspection surface of S, suppose by focus x _sExcite acceptance point x _rThe two component elastic wave seismologic records that receive are u _i(x _rω), the reverse elastic wave displacement field u of continuation then _m(x; x _rω) can represent by the Kirchhoff-Helmholtz integration:

$u_m(x;x_r;\mathrm{\ω})={\∫}_S{\mathrm{dx}}_r[t_i(x_r;\mathrm{\ω})G_{\mathrm{im}}^{*}(x;x_r;\mathrm{\ω})-u_i(x_r;\mathrm{\ω}){\mathrm{\Σ}}_{\mathrm{im}}^{*}(x;x_r;\mathrm{\ω})]---\left(1\right)$

Wherein * represents complex conjugate; G _Lm(x; x _rω) be displacement Green tensor, it represents x _rThe component of the m direction of the x place displacement that the l of place direction specific body force is caused; t _i(x _r) be x _rPlace's stress; ∑ _Im(x; x _r) be stress Green tensor;

If hypothesis S is the free face of land, then according to free stress boundary condition t (x; ω)=0, thus formula (1) becomes:

$u_m(x;x_r;\mathrm{\ω})=-{\∫}_S{\mathrm{dx}}_r{u}_i(x_r;\mathrm{\ω}){\mathrm{\Σ}}_{\mathrm{im}}^{*}(x;x_r;\mathrm{\ω})---\left(2\right)$

And get G _Lm(x; x _rHigh-frequency approximation ω) is separated, substitution formula (2),

$u_m(x;x_r;\mathrm{\ω})=u_m^p(x;x_r;\mathrm{\ω})+u_m^s(x;x_r;\mathrm{\ω})$

$=-\mathrm{i\ω}\underset{v}{\mathrm{\Σ}}{\∫}_S{\mathrm{dx}}_r\mathrm{\ρ}\left(x_r\right)U_m^{v*}(x;x_r;\mathrm{\ω})[u_1(x_r;\mathrm{\ω})W_1^v\left(x_r\right)+u_2(x_r;\mathrm{\ω})W_2^v\left(x_r\right)]$

(3)

Formula (3) is the elastic wave wave field extrapolation formula of decoupling zero, wherein,

Be respectively displacement field u _m(x; x _rP ripple ω) and S wave component;

Weights Has following form

$\left\{\begin{array}{c}{W}_1^p\left(x_r\right)=2v_s^2\left(x_r\right)p_2^p\left(x_r\right)e_1^p\left(x_r\right)\\ W_2^p\left(x_r\right)=2v_s^2\left(x_r\right)p_2^p\left(x_r\right)e_2^p\left(x_r\right)+\left(\frac{v_p^2\left(x_r\right)-2v_s^2\left(x_r\right)}{v_p\left(x_r\right)}\right)\\ W_1^s\left(x_r\right)=v_s^2\left(x_r\right)p_2^s\left(x_r\right)e_1^s\left(x_r\right)+v_s^2\left(x_r\right)p_1^s\left(x_r\right)e_2^s\left(x_r\right)\\ W_2^s\left(x_r\right)=-2v_s^2\left(x_r\right)p_1^s\left(x_r\right)e_1^s\left(x_r\right)\end{array}\right.---\left(4\right)$

Utilize the superposition integral of elastokinetics Gaussian beam to come displacement vector in the calculating formula (3)

Obtain the elasticity wave field extrapolation formula represented based on Gaussian beam;

Be the plane before the primary wave of Gaussian beam, the multi-component seismic record is decomposed into the part plan ripple of different wave modes difference inceptive directions, utilize corresponding Gaussian beam to carry out wave field extrapolation then, its specific implementation process is as follows:

Seismologic record is added a series of Gaussian windows, it is divided into a series of overlapping regional areas, the Gaussian window function has following character

$\frac{\mathrm{\ΔL}}{\sqrt{2\mathrm{\π}}{w}_0}\sqrt{\left|\frac{\mathrm{\ω}}{{\mathrm{\ω}}_r}\right|}\underset{L}{\mathrm{\Σ}}\exp [-|\frac{\mathrm{\ω}}{{\mathrm{\ω}}_r}\left|\frac{(x_r-L)}{2w_0^2}\right]\≈1---\left(5\right)$

Wherein, x=L is the center line of Gaussian window, also is the horizontal coordinate of beam center L; Δ L is its horizontal interval, and formula (5) substitution formula (3) is got

$u_m(x;x_r;\mathrm{\ω})=-\frac{\mathrm{i\ω\ΔL}}{\sqrt{2\mathrm{\π}}{w}_0}\sqrt{\left|\frac{\mathrm{\ω}}{{\mathrm{\ω}}_r}\right|}\underset{v}{\mathrm{\Σ}}\underset{L}{\mathrm{\Σ}}{\∫}_S{\mathrm{dx}}_r\mathrm{\ρ}\left(x_r\right)U_m^{v*}(x;x_r;\mathrm{\ω})[u_1(x_r;\mathrm{\ω})W_1^v\left(x_r\right)$ (6)

$+u_2(x_r;\mathrm{\ω})W_2^v\left(x_r\right)\left]\exp \right[-\left|\frac{\mathrm{\ω}}{{\mathrm{\ω}}_r}\right|\frac{(x_r-L)}{2w_0^2}]$

Introduce the phase shift correction factor, will

Gaussian beam by the outgoing of L place is represented

$U_m^v(x;x_r;\mathrm{\ω})\≈{\mathrm{\Ψ}}^v\∫{\hat{u}}_m^v(x;L;\mathrm{\ω})\exp [-\mathrm{i\ω}{p}_1^v\left(L\right)(x_r-L)]\frac{{\mathrm{dp}}_1^v\left(L\right)}{p_2^v\left(L\right)}---\left(7\right)$

With following formula substitution formula (6), make ρ (x _r) ≈ ρ (L),

And the exchange integral order, just can obtain the P ripple displacement of reverse continuation And S ripple displacement

$u_m^p(x;x_r;\mathrm{\ω})=-\frac{\mathrm{\ΔL\ω}}{4\mathrm{\π}}\underset{L}{\mathrm{\Σ}}\∫\frac{d{p}_1^p\left(L\right)}{p_2^p\left(L\right)}\sqrt{\mathrm{\ρ}\left(L\right)}{\hat{u}}_m^{p*}(x;L;\mathrm{\ω})---\left(8\right)$

$\×[W_1^p\left(L\right)D_1^p(L;p_1^p;\mathrm{\ω})+W_2^p\left(L\right)D_2^p(L;p_1^p;\mathrm{\ω})]$

$u_m^s(x;x_r;\mathrm{\ω})=-\frac{\mathrm{\ΔL\ω}}{4\mathrm{\π}}\underset{L}{\mathrm{\Σ}}\∫\frac{d{p}_1^s\left(L\right)}{p_2^s\left(L\right)}\sqrt{\mathrm{\ρ}\left(L\right)}{\hat{u}}_m^{s*}(x;L;\mathrm{\ω})---\left(9\right)$

$\×[W_1^s\left(L\right)D_1^s(L;p_1^s;\mathrm{\ω})+W_2^s\left(L\right)D_2^s(L;p_1^s;\mathrm{\ω})]$

Wherein, Be windowing local dip stack to different wave mode multi-component seismic records

$D_n^v(L;P_1^v;\mathrm{\ω})=\sqrt{\frac{\left|\mathrm{\ω}\right|}{2\mathrm{\π}}}{\∫}_S{\mathrm{dx}}_r{u}_n(x_r;\mathrm{\ω})\exp [\mathrm{i\ω}{p}_1^v\left(L\right)(x_r-L)-|\frac{\mathrm{\ω}}{{\mathrm{\ω}}_r}\left|\frac{{(x_r-L)}^2}{2w_0^2}\right]---\left(10\right)$

Weights Be at this moment

$W_1^p\left(L\right)=2{\mathrm{\γ}}^2\left(L\right)p_2^p\left(L\right)e_1^p\left(L\right),$ $W_2^p\left(L\right)=2{\mathrm{\γ}}^2\left(L\right)p_2^p\left(L\right)e_2^p\left(L\right)+\left(\frac{1-2{\mathrm{\γ}}^2\left(L\right)}{v_p\left(L\right)}\right)$ $W_1^s\left(L\right)=p_2^s\left(L\right)e_1^s\left(L\right)+p_1^s\left(L\right)e_2^s\left(L\right),$ $W_2^s\left(L\right)=-2p_1^s\left(L\right)e_1^s\left(L\right),$ $\mathrm{\γ}\left(L\right)=\frac{v_s\left(L\right)}{v_p\left(L\right)}$

(11)

Wherein in the step (4), according to the Claerbout image-forming principle, at first with focus displacement wave field

Represent by the elastokinetics Gaussian beam

$U_m^p(x;x_s;\mathrm{\ω})\≈\frac{i}{4\mathrm{\π}{v}_p^2\left(x_s\right)}\sqrt{\frac{{\mathrm{\ω}}_r{w}_0^2}{\mathrm{\ρ}\left(x_s\right)}}\∫\frac{d{p}_1^p\left(x_s\right)}{p_2^p\left(x_s\right)}{\hat{u}}_m^p(x;x_s;\mathrm{\ω})---\left(12\right)$

Next according to P ripple and S wave propagation characteristics, convolution (8), the reception wave field of the reverse continuation that formula (9) is represented is defined as follows imaging formula

$I^{\mathrm{pp}}\left(x\right)={\∫U}_2^{p*}(x;x_s;\mathrm{\ω})u_2^p(x;x_r;\mathrm{\ω})\mathrm{d\ω}$

$=\frac{\mathrm{\ΔL}\sqrt{{\mathrm{\ω}}_r{w}_0^2}}{16{\mathrm{\π}}^2}\underset{L}{\mathrm{\Σ}}\∫\mathrm{d\ω}\frac{\mathrm{i\ω}}{v_p^2\left(x_s\right)}\sqrt{\frac{\mathrm{\ρ}\left(L\right)}{\mathrm{\ρ}\left(x_s\right)}}\∫\∫\frac{d{p}_1^p\left(x_s\right)d{p}_1^p\left(L\right)}{p_2^p\left(x_s\right)p_2^p\left(L\right)}---\left(13\right)$

$\×{\hat{u}}_2^{p*}(x;x_s;\mathrm{\ω}){\hat{u}}_1^{p*}(x;L;\mathrm{\ω})[W_1^p\left(L\right)D_1^p(L;p_1^p;\mathrm{\ω})+W_2^p\left(L\right)D_2^p(L;p_1^p;\mathrm{\ω})]$

$I^{\mathrm{ps}}\left(x\right)=\frac{\mathrm{\ΔL}\sqrt{{\mathrm{\ω}}_r{w}_0^2}}{16{\mathrm{\π}}^2}\underset{L}{\mathrm{\Σ}}\∫\mathrm{d\ω}\frac{\mathrm{i\ω}}{v_p^2\left(x_s\right)}\sqrt{\frac{\mathrm{\ρ}\left(L\right)}{\mathrm{\ρ}\left(x_s\right)}}\∫\∫\frac{d{p}_1^p\left(x_s\right)d{p}_1^s\left(L\right)}{p_2^p\left(x_s\right)p_2^s\left(L\right)}\mathrm{sgn}\left(\mathrm{\α}\right)---\left(14\right)$

$\×{\hat{u}}_2^{p*}(x;x_s;\mathrm{\ω}){\hat{u}}_1^{s*}(x;L;\mathrm{\ω})[W_1^s\left(L\right)D_1^s(L;p_1^s;\mathrm{\ω})+W_2^s\left(L\right)D_2^s(L;p_1^s;\mathrm{\ω})]$

Wherein, I ^Pp(x) be the single big gun imaging of P-P value, I ^Ps(x) be the single big gun imaging of P-S value, all big guns records are offset and superpose, just obtain final elastic wave imaging results.

The invention has the beneficial effects as follows: this method not only can have higher precision, and take into account counting yield and dirigibility to multiple reflection to carrying out imaging, is a kind of multi-component seismic formation method with higher using value.

Description of drawings

Fig. 1 is the flat bed model skew tentative calculation design sketch of embodiment.(a) horizontal component single shot record among the figure; (b) vertical component single shot record; (c) vertical component record P-P ripple list big gun imaging results; (d) the single big gun imaging results of the P-P of two components record; (e) horizontal component record P-S ripple list big gun imaging results; (f) the P-S imaging results that writes down through two components behind the polarity correction; (g) P-P ripple imaging stack result; (h) P-S ripple imaging stack result.

Fig. 2 is the polarization figure of P-S transformed wave at the reflecting interface place.

Fig. 3 is the complex model skew The trial result figure of embodiment.(a) p wave interval velocity model wherein; (b) S wave velocity model; (c) X component single shot record; (d) Z component single shot record; (e) P-P ripple imaging stack result; (f) P-S ripple imaging stack result.

Embodiment

Elastic wave Gaussian beam pre-stack depth migration technology, performing step is:

(1) according to the Gaussian beam original width of selecting, by Gaussian window single shot record is divided into a series of overlapping areas, each regional center is a beam center;

(2) utilize the local dip stack to be decomposed into the part plan wave component of different wave modes the multi-component seismic record;

(3) respectively from focus and beam center trial fire elastokinetics Gaussian beam, when calculating walking in every Gaussian beam effective range, amplitude information;

(4) according to the imaging formula of derivation PP ripple and PS ripple, when utilizing the walking of Gaussian beam, the local obliquity information on amplitude, polarity information and the face of land calculates the imaging value of certain corresponding beam center;

(5) repeat the computing in (3), (4) step, obtain the imaging results of single big gun, final migration result is the stack of all single big gun imaging results.

Wherein the computing method in the step (3) are:

On the described inspection surface of S, suppose by focus x _sExcite acceptance point x _rThe two component elastic wave seismologic records that receive are u _i(x _rω), the reverse elastic wave displacement field u of continuation then _m(x; x _rω) can represent by the Kirchhoff-Helmholtz integration:

$u_m(x;x_r;\mathrm{\ω})={\∫}_S{\mathrm{dx}}_r[t_i(x_r;\mathrm{\ω})G_{\mathrm{im}}^{*}(x;x_r;\mathrm{\ω})-u_i(x_r;\mathrm{\ω}){\mathrm{\Σ}}_{\mathrm{im}}^{*}(x;x_r;\mathrm{\ω})]---\left(1\right)$

Wherein * represents complex conjugate; G _Lm(x; x _rω) be displacement Green tensor, it represents x _rThe component of the m direction of the x place displacement that the l of place direction specific body force is caused; t _i(x _r) be x _rPlace's stress; ∑ _Im(x; x _r) be stress Green tensor;

If hypothesis S is the free face of land, then according to free stress boundary condition t (x; ω)=0, thus formula (1) becomes:

$u_m(x;x_r;\mathrm{\ω})=-{\∫}_S{\mathrm{dx}}_r{u}_i(x_r;\mathrm{\ω}){\mathrm{\Σ}}_{\mathrm{im}}^{*}(x;x_r;\mathrm{\ω})---\left(2\right)$

And get G _Lm(x; x _rHigh-frequency approximation ω) is separated, substitution formula (2),

$u_m(x;x_r;\mathrm{\ω})=u_m^p(x;x_r;\mathrm{\ω})+u_m^s(x;x_r;\mathrm{\ω})$

$=-\mathrm{i\ω}\underset{v}{\mathrm{\Σ}}{\∫}_S{\mathrm{dx}}_r\mathrm{\ρ}\left(x_r\right)U_m^{v*}(x;x_r;\mathrm{\ω})[u_1(x_r;\mathrm{\ω})W_1^v\left(x_r\right)+u_2(x_r;\mathrm{\ω})W_2^v\left(x_r\right)]$

(3)

Formula (3) is the elastic wave wave field extrapolation formula of decoupling zero, wherein,

Be respectively displacement field u _m(x; x _rP ripple ω) and S wave component;

Weights

Has following form

$\left\{\begin{array}{c}{W}_1^p\left(x_r\right)=2v_s^2\left(x_r\right)p_2^p\left(x_r\right)e_1^p\left(x_r\right)\\ W_2^p\left(x_r\right)=2v_s^2\left(x_r\right)p_2^p\left(x_r\right)e_2^p\left(x_r\right)+\left(\frac{v_p^2\left(x_r\right)-2v_s^2\left(x_r\right)}{v_p\left(x_r\right)}\right)\\ W_1^s\left(x_r\right)=v_s^2\left(x_r\right)p_2^s\left(x_r\right)e_1^s\left(x_r\right)+v_s^2\left(x_r\right)p_1^s\left(x_r\right)e_2^s\left(x_r\right)\\ W_2^s\left(x_r\right)=-2v_s^2\left(x_r\right)p_1^s\left(x_r\right)e_1^s\left(x_r\right)\end{array}\right.---\left(4\right)$

Utilize the superposition integral of elastokinetics Gaussian beam to come displacement vector in the calculating formula (3)

Obtain the elasticity wave field extrapolation formula represented based on Gaussian beam;

Be the plane before the primary wave of Gaussian beam, the multi-component seismic record is decomposed into the part plan ripple of different wave modes difference inceptive directions, utilize corresponding Gaussian beam to carry out wave field extrapolation then, its specific implementation process is as follows:

Seismologic record is added a series of Gaussian windows, it is divided into a series of overlapping regional areas, the Gaussian window function has following character

$\frac{\mathrm{\ΔL}}{\sqrt{2\mathrm{\π}}{w}_0}\sqrt{\left|\frac{\mathrm{\ω}}{{\mathrm{\ω}}_r}\right|}\underset{L}{\mathrm{\Σ}}\exp [-|\frac{\mathrm{\ω}}{{\mathrm{\ω}}_r}\left|\frac{(x_r-L)}{2w_0^2}\right]\≈1---\left(5\right)$

Wherein, x=L is the center line of Gaussian window, also is the horizontal coordinate of beam center L; Δ L is its horizontal interval, and formula (5) substitution formula (3) is got

$u_m(x;x_r;\mathrm{\ω})=-\frac{\mathrm{i\ω\ΔL}}{\sqrt{2\mathrm{\π}}{w}_0}\sqrt{\left|\frac{\mathrm{\ω}}{{\mathrm{\ω}}_r}\right|}\underset{v}{\mathrm{\Σ}}\underset{L}{\mathrm{\Σ}}{\∫}_S{\mathrm{dx}}_r\mathrm{\ρ}\left(x_r\right)U_m^{v*}(x;x_r;\mathrm{\ω})[u_1(x_r;\mathrm{\ω})W_1^v\left(x_r\right)$ (6)

$+u_2(x_r;\mathrm{\ω})W_2^v\left(x_r\right)\left]\exp \right[-\left|\frac{\mathrm{\ω}}{{\mathrm{\ω}}_r}\right|\frac{(x_r-L)}{2w_0^2}]$

Introduce the phase shift correction factor, will Gaussian beam by the outgoing of L place is represented

$U_m^v(x;x_r;\mathrm{\ω})\≈{\mathrm{\Ψ}}^v\∫{\hat{u}}_m^v(x;L;\mathrm{\ω})\exp [-\mathrm{i\ω}{p}_1^v\left(L\right)(x_r-L)]\frac{{\mathrm{dp}}_1^v\left(L\right)}{p_2^v\left(L\right)}---\left(7\right)$

With following formula substitution formula (6), make ρ (x _r) ≈ ρ (L),

And the exchange integral order, just can obtain the P ripple displacement of reverse continuation

And S ripple displacement

$u_m^p(x;x_r;\mathrm{\ω})=-\frac{\mathrm{\ΔL\ω}}{4\mathrm{\π}}\underset{L}{\mathrm{\Σ}}\∫\frac{d{p}_1^p\left(L\right)}{p_2^p\left(L\right)}\sqrt{\mathrm{\ρ}\left(L\right)}{\hat{u}}_m^{p*}(x;L;\mathrm{\ω})---\left(8\right)$

$\×[W_1^p\left(L\right)D_1^p(L;p_1^p;\mathrm{\ω})+W_2^p\left(L\right)D_2^p(L;p_1^p;\mathrm{\ω})]$

$u_m^s(x;x_r;\mathrm{\ω})=-\frac{\mathrm{\ΔL\ω}}{4\mathrm{\π}}\underset{L}{\mathrm{\Σ}}\∫\frac{d{p}_1^s\left(L\right)}{p_2^s\left(L\right)}\sqrt{\mathrm{\ρ}\left(L\right)}{\hat{u}}_m^{s*}(x;L;\mathrm{\ω})---\left(9\right)$

$\×[W_1^s\left(L\right)D_1^s(L;p_1^s;\mathrm{\ω})+W_2^s\left(L\right)D_2^s(L;p_1^s;\mathrm{\ω})]$

Wherein,

Be windowing local dip stack to different wave mode multi-component seismic records

$D_n^v(L;P_1^v;\mathrm{\ω})=\sqrt{\frac{\left|\mathrm{\ω}\right|}{2\mathrm{\π}}}{\∫}_S{\mathrm{dx}}_r{u}_n(x_r;\mathrm{\ω})\exp [\mathrm{i\ω}{p}_1^v\left(L\right)(x_r-L)-|\frac{\mathrm{\ω}}{{\mathrm{\ω}}_r}\left|\frac{{(x_r-L)}^2}{2w_0^2}\right]---\left(10\right)$

Weights

Be at this moment

$W_1^p\left(L\right)=2{\mathrm{\γ}}^2\left(L\right)p_2^p\left(L\right)e_1^p\left(L\right),$ $W_2^p\left(L\right)=2{\mathrm{\γ}}^2\left(L\right)p_2^p\left(L\right)e_2^p\left(L\right)+\left(\frac{1-2{\mathrm{\γ}}^2\left(L\right)}{v_p\left(L\right)}\right)$ $W_1^s\left(L\right)=p_2^s\left(L\right)e_1^s\left(L\right)+p_1^s\left(L\right)e_2^s\left(L\right),$ $W_2^s\left(L\right)=-2p_1^s\left(L\right)e_1^s\left(L\right),$ $\mathrm{\γ}\left(L\right)=\frac{v_s\left(L\right)}{v_p\left(L\right)}$

(11)

Wherein in the step (4), according to the Claerbout image-forming principle, at first with focus displacement wave field

Represent by the elastokinetics Gaussian beam

$U_m^p(x;x_s;\mathrm{\ω})\≈\frac{i}{4\mathrm{\π}{v}_p^2\left(x_s\right)}\sqrt{\frac{{\mathrm{\ω}}_r{w}_0^2}{\mathrm{\ρ}\left(x_s\right)}}\∫\frac{d{p}_1^p\left(x_s\right)}{p_2^p\left(x_s\right)}{\hat{u}}_m^p(x;x_s;\mathrm{\ω})---\left(12\right)$

Next according to P ripple and S wave propagation characteristics, convolution (8), the reception wave field of the reverse continuation that formula (9) is represented is defined as follows imaging formula

$I^{\mathrm{pp}}\left(x\right)={\∫U}_2^{p*}(x;x_s;\mathrm{\ω})u_2^p(x;x_r;\mathrm{\ω})\mathrm{d\ω}$

$=\frac{\mathrm{\ΔL}\sqrt{{\mathrm{\ω}}_r{w}_0^2}}{16{\mathrm{\π}}^2}\underset{L}{\mathrm{\Σ}}\∫\mathrm{d\ω}\frac{\mathrm{i\ω}}{v_p^2\left(x_s\right)}\sqrt{\frac{\mathrm{\ρ}\left(L\right)}{\mathrm{\ρ}\left(x_s\right)}}\∫\∫\frac{d{p}_1^p\left(x_s\right)d{p}_1^p\left(L\right)}{p_2^p\left(x_s\right)p_2^p\left(L\right)}---\left(13\right)$

$\×{\hat{u}}_2^{p*}(x;x_s;\mathrm{\ω}){\hat{u}}_1^{p*}(x;L;\mathrm{\ω})[W_1^p\left(L\right)D_1^p(L;p_1^p;\mathrm{\ω})+W_2^p\left(L\right)D_2^p(L;p_1^p;\mathrm{\ω})]$

$I^{\mathrm{ps}}\left(x\right)=\frac{\mathrm{\ΔL}\sqrt{{\mathrm{\ω}}_r{w}_0^2}}{16{\mathrm{\π}}^2}\underset{L}{\mathrm{\Σ}}\∫\mathrm{d\ω}\frac{\mathrm{i\ω}}{v_p^2\left(x_s\right)}\sqrt{\frac{\mathrm{\ρ}\left(L\right)}{\mathrm{\ρ}\left(x_s\right)}}\∫\∫\frac{d{p}_1^p\left(x_s\right)d{p}_1^s\left(L\right)}{p_2^p\left(x_s\right)p_2^s\left(L\right)}\mathrm{sgn}\left(\mathrm{\α}\right)---\left(14\right)$

$\×{\hat{u}}_2^{p*}(x;x_s;\mathrm{\ω}){\hat{u}}_1^{s*}(x;L;\mathrm{\ω})[W_1^s\left(L\right)D_1^s(L;p_1^s;\mathrm{\ω})+W_2^s\left(L\right)D_2^s(L;p_1^s;\mathrm{\ω})]$

Wherein, I ^Pp(x) be the single big gun imaging of P-P value, I ^Ps(x) be the single big gun imaging of P-S value, all big guns records are offset and superpose, just obtain final elastic wave imaging results.

At first test correctness of the present invention by three layers of layer-cake model.Model meshes is 301 * 401, all is 10m to sampling interval in length and breadth, and two flat reflectors lay respectively at 1400m and 2700m place.Along with the increase of the degree of depth, p wave interval velocity is respectively 2500m/s, 2900m/s, and 3200m/s, the S wave velocity is respectively 1500m/s, 1800m/s, 2600m/s, density is set at constant.Simulation elastic wave record dominant frequency is 30Hz, has 51 big guns, every big gun 151 roads, and track pitch is 20m, and it is 1500 that time-sampling is counted, and sampling interval is 2ms.Fig. 1 a, Fig. 1 b is respectively the X of the 26th big gun single shot record, Z component.Fig. 1 c is only to the single big gun imaging results of the P-P of Z component record, can obviously see the crosstalk noise (arrow indication place) that the S wave component in the Z component record is caused, Fig. 1 d is that the single big gun imaging results of the P-P that obtains after the weighting is write down to two components in employing formula (15), (23), and this moment, crosstalk noise obtained effective compacting.Fig. 1 e is the single big gun imaging results of the P-S of X component record, Fig. 1 f is employing formula (15), (23) to the single big gun imaging results of P-S after the weighting of two components record, can see, than Fig. 1 e, crosstalk noise among Fig. 1 f has obtained effective compacting (arrow indication place), and the reversal of poles phenomenon in the imaging results has obtained effective correction (oval sign place).Fig. 1 g, 1h are respectively the P-P after all big gun skew stacks, and the P-S imaging results can both all reflect the true form at interface clearly, and because the S velocity of wave propagation is slower, make the P-S imaging results have higher resolution relatively.

When transformed wave is carried out imaging, different and the reversal of poles phenomenon in the imaging section that causes in transformed wave polarization direction often appears, have a strong impact on the effect of stacking image.At the problems referred to above, the present invention proposes a kind of bearing calibration.Fig. 2 has shown P-S conversion wave propagation process, can see because incident P ripple has the incident angle (α of distinct symbols ₁＞0, α ₂＜0, with counterclockwise for just), make O ₁The P-S transformed wave that point produces has the displacement component of positive horizontal direction (just to be), O to the right ₂The P-S transformed wave that produces has the displacement component of negative horizontal direction, finally makes R ₁, R ₂The X component seismic record that point receives has opposite polarity.The polarity that is to say conversion S ripple in the X component record is positive and negative relevant with P ripple incident angle, therefore can be according to positive and negative directly imaging results correction of P ripple incident angle α, and the created symbol function obtains the P-S transformed wave imaging formula shown in the formula (14) for this reason.

Utilize a complex structure model to test the imaging effect of the present invention to the complicated geological model.Rate pattern such as Fig. 3 a, shown in the 3b, its sizing grid is 1000 * 550, all is 10m to sampling interval in length and breadth, 5m is set at constant with density equally.It is 7 layers that model is divided into, and wherein comprises tomography, complex structure such as hollow.Simulation elastic wave record dominant frequency is 25Hz, has 200 big guns, every big gun 361 roads, and track pitch is 10m, and it is 1200 that time-sampling is counted, and sampling interval is 2.1ms.Fig. 3 c, 3d are respectively the X of the 1st big gun single shot record, and Z component can be seen because the wave field information of underground structure complexity in causing writing down is also very complicated.Fig. 3 e, 3f are respectively final P-P, the P-S imaging results, and the structure situation of rate pattern of can having seen reflection that both are all good, the tomography in the model, bottom hollow and that rise and fall have all obtained imaging preferably, have proved the validity of this paper method.

Claims (3)

1. elastic wave Gaussian beam pre-stack depth migration technology is characterized in that performing step is:

(1) according to the Gaussian beam original width of selecting, by Gaussian window single shot record is divided into a series of overlapping areas, each regional center is a beam center;

(2) utilize the local dip stack to be decomposed into the part plan wave component of different wave modes the multi-component seismic record;

(3) respectively from focus and beam center trial fire elastokinetics Gaussian beam, when calculating walking in every Gaussian beam effective range, amplitude information;

(4) according to the imaging formula of derivation PP ripple and PS ripple, when utilizing the walking of Gaussian beam, the local obliquity information on amplitude, polarity information and the face of land calculates the imaging value of certain corresponding beam center;

(5) repeat the computing in (3), (4) step, obtain the imaging results of single big gun, final migration result is the stack of all single big gun imaging results.

2. elastic wave Gaussian beam pre-stack depth migration technology according to claim 1 is characterized in that the computing method in the step (3) are:

On the described inspection surface of S, suppose by focus x _sExcite acceptance point x _rThe two component elastic wave seismologic records that receive are u _i(x _rω), the reverse elastic wave displacement field u of continuation then _m(x; x _rω) can represent by the Kirchhoff-Helmholtz integration:

$u_m(x;x_r;\mathrm{\ω})={\∫}_S{\mathrm{dx}}_r[t_i(x_r;\mathrm{\ω})G_{\mathrm{im}}^{*}(x;x_r;\mathrm{\ω})-u_i(x_r;\mathrm{\ω}){\mathrm{\Σ}}_{\mathrm{im}}^{*}(x;x_r;\mathrm{\ω})]---\left(1\right)$

Wherein * represents complex conjugate; G _Lm(x; x _rω) be displacement Green tensor, it represents x _rThe component of the m direction of the x place displacement that the l of place direction specific body force is caused; t _i(x _r) be x _rPlace's stress; ∑ _Im(x; x _r) be stress Green tensor;

If hypothesis S is the free face of land, then according to free stress boundary condition t (x; ω)=0, thus formula (1) becomes:

$u_m(x;x_r;\mathrm{\ω})=-{\∫}_S{\mathrm{dx}}_r{u}_i(x_r;\mathrm{\ω}){\mathrm{\Σ}}_{\mathrm{im}}^{*}(x;x_r;\mathrm{\ω})---\left(2\right)$

And get G _Lm(x; x _rHigh-frequency approximation ω) is separated, substitution formula (2),

$u_m(x;x_r;\mathrm{\ω})=u_m^p(x;x_r;\mathrm{\ω})+u_m^s(x;x_r;\mathrm{\ω})$

$=-\mathrm{i\ω}\underset{v}{\mathrm{\Σ}}{\∫}_S{\mathrm{dx}}_r\mathrm{\ρ}\left(x_r\right)U_m^{v*}(x;x_r;\mathrm{\ω})[u_1(x_r;\mathrm{\ω})W_1^v\left(x_r\right)+u_2(x_r;\mathrm{\ω})W_2^v\left(x_r\right)]$

(3)

Formula (3) is the elastic wave wave field extrapolation formula of decoupling zero, wherein,

Be respectively displacement field u _m(x; x _rP ripple ω) and S wave component;

Weights

Has following form

$\left\{\begin{array}{c}{W}_1^p\left(x_r\right)=2v_s^2\left(x_r\right)p_2^p\left(x_r\right)e_1^p\left(x_r\right)\\ W_2^p\left(x_r\right)=2v_s^2\left(x_r\right)p_2^p\left(x_r\right)e_2^p\left(x_r\right)+\left(\frac{v_p^2\left(x_r\right)-2v_s^2\left(x_r\right)}{v_p\left(x_r\right)}\right)\\ W_1^s\left(x_r\right)=v_s^2\left(x_r\right)p_2^s\left(x_r\right)e_1^s\left(x_r\right)+v_s^2\left(x_r\right)p_1^s\left(x_r\right)e_2^s\left(x_r\right)\\ W_2^s\left(x_r\right)=-2v_s^2\left(x_r\right)p_1^s\left(x_r\right)e_1^s\left(x_r\right)\end{array}\right.---\left(4\right)$

Utilize the superposition integral of elastokinetics Gaussian beam to come displacement vector in the calculating formula (3)

Obtain the elasticity wave field extrapolation formula represented based on Gaussian beam;

Be the plane before the primary wave of Gaussian beam, the multi-component seismic record is decomposed into the part plan ripple of different wave modes difference inceptive directions, utilize corresponding Gaussian beam to carry out wave field extrapolation then, its specific implementation process is as follows:

Seismologic record is added a series of Gaussian windows, it is divided into a series of overlapping regional areas, the Gaussian window function has following character

$\frac{\mathrm{\ΔL}}{\sqrt{2\mathrm{\π}}{w}_0}\sqrt{\left|\frac{\mathrm{\ω}}{{\mathrm{\ω}}_r}\right|}\underset{L}{\mathrm{\Σ}}\exp [-|\frac{\mathrm{\ω}}{{\mathrm{\ω}}_r}\left|\frac{(x_r-L)}{2w_0^2}\right]\≈1---\left(5\right)$

Wherein, x=L is the center line of Gaussian window, also is the horizontal coordinate of beam center L; Δ L is its horizontal interval, and formula (5) substitution formula (3) is got

$u_m(x;x_r;\mathrm{\ω})=-\frac{\mathrm{i\ω\ΔL}}{\sqrt{2\mathrm{\π}}{w}_0}\sqrt{\left|\frac{\mathrm{\ω}}{{\mathrm{\ω}}_r}\right|}\underset{v}{\mathrm{\Σ}}\underset{L}{\mathrm{\Σ}}{\∫}_S{\mathrm{dx}}_r\mathrm{\ρ}\left(x_r\right)U_m^{v*}(x;x_r;\mathrm{\ω})[u_1(x_r;\mathrm{\ω})W_1^v\left(x_r\right)$ (6)

$+u_2(x_r;\mathrm{\ω})W_2^v\left(x_r\right)\left]\exp \right[-\left|\frac{\mathrm{\ω}}{{\mathrm{\ω}}_r}\right|\frac{(x_r-L)}{2w_0^2}]$

Introduce the phase shift correction factor, will

Gaussian beam by the outgoing of L place is represented

$U_m^v(x;x_r;\mathrm{\ω})\≈{\mathrm{\Ψ}}^v\∫{\hat{u}}_m^v(x;L;\mathrm{\ω})\exp [-\mathrm{i\ω}{p}_1^v\left(L\right)(x_r-L)]\frac{{\mathrm{dp}}_1^v\left(L\right)}{p_2^v\left(L\right)}---\left(7\right)$

With following formula substitution formula (6), make ρ (x _r) ≈ ρ (L),

And the exchange integral order, just can obtain the P ripple displacement of reverse continuation

And S ripple displacement

$u_m^p(x;x_r;\mathrm{\ω})=-\frac{\mathrm{\ΔL\ω}}{4\mathrm{\π}}\underset{L}{\mathrm{\Σ}}\∫\frac{d{p}_1^p\left(L\right)}{p_2^p\left(L\right)}\sqrt{\mathrm{\ρ}\left(L\right)}{\hat{u}}_m^{p*}(x;L;\mathrm{\ω})---\left(8\right)$

$\×[W_1^p\left(L\right)D_1^p(L;p_1^p;\mathrm{\ω})+W_2^p\left(L\right)D_2^p(L;p_1^p;\mathrm{\ω})]$

$u_m^s(x;x_r;\mathrm{\ω})=-\frac{\mathrm{\ΔL\ω}}{4\mathrm{\π}}\underset{L}{\mathrm{\Σ}}\∫\frac{d{p}_1^s\left(L\right)}{p_2^s\left(L\right)}\sqrt{\mathrm{\ρ}\left(L\right)}{\hat{u}}_m^{s*}(x;L;\mathrm{\ω})---\left(9\right)$

$\×[W_1^s\left(L\right)D_1^s(L;p_1^s;\mathrm{\ω})+W_2^s\left(L\right)D_2^s(L;p_1^s;\mathrm{\ω})]$

Wherein,

Be windowing local dip stack to different wave mode multi-component seismic records

$D_n^v(L;P_1^v;\mathrm{\ω})=\sqrt{\frac{\left|\mathrm{\ω}\right|}{2\mathrm{\π}}}{\∫}_S{\mathrm{dx}}_r{u}_n(x_r;\mathrm{\ω})\exp [\mathrm{i\ω}{p}_1^v\left(L\right)(x_r-L)-|\frac{\mathrm{\ω}}{{\mathrm{\ω}}_r}\left|\frac{{(x_r-L)}^2}{2w_0^2}\right]---\left(10\right)$

Weights

Be at this moment

$W_1^p\left(L\right)=2{\mathrm{\γ}}^2\left(L\right)p_2^p\left(L\right)e_1^p\left(L\right),$ $W_2^p\left(L\right)=2{\mathrm{\γ}}^2\left(L\right)p_2^p\left(L\right)e_2^p\left(L\right)+\left(\frac{1-2{\mathrm{\γ}}^2\left(L\right)}{v_p\left(L\right)}\right)$ $W_1^s\left(L\right)=p_2^s\left(L\right)e_1^s\left(L\right)+p_1^s\left(L\right)e_2^s\left(L\right),$ $W_2^s\left(L\right)=-2p_1^s\left(L\right)e_1^s\left(L\right),$ $\mathrm{\γ}\left(L\right)=\frac{v_s\left(L\right)}{v_p\left(L\right)}$

(11)。

3. elastic wave Gaussian beam pre-stack depth migration technology according to claim 1 is characterized in that in the step (4), according to the Claerbout image-forming principle, at first with focus displacement wave field

Represent by the elastokinetics Gaussian beam

$U_m^p(x;x_s;\mathrm{\ω})\≈\frac{i}{4\mathrm{\π}{v}_p^2\left(x_s\right)}\sqrt{\frac{{\mathrm{\ω}}_r{w}_0^2}{\mathrm{\ρ}\left(x_s\right)}}\∫\frac{d{p}_1^p\left(x_s\right)}{p_2^p\left(x_s\right)}{\hat{u}}_m^p(x;x_s;\mathrm{\ω})---\left(12\right)$

Next according to P ripple and S wave propagation characteristics, convolution (8), the reception wave field of the reverse continuation that formula (9) is represented is defined as follows imaging formula

$I^{\mathrm{pp}}\left(x\right)={\∫U}_2^{p*}(x;x_s;\mathrm{\ω})u_2^p(x;x_r;\mathrm{\ω})\mathrm{d\ω}$

$=\frac{\mathrm{\ΔL}\sqrt{{\mathrm{\ω}}_r{w}_0^2}}{16{\mathrm{\π}}^2}\underset{L}{\mathrm{\Σ}}\∫\mathrm{d\ω}\frac{\mathrm{i\ω}}{v_p^2\left(x_s\right)}\sqrt{\frac{\mathrm{\ρ}\left(L\right)}{\mathrm{\ρ}\left(x_s\right)}}\∫\∫\frac{d{p}_1^p\left(x_s\right)d{p}_1^p\left(L\right)}{p_2^p\left(x_s\right)p_2^p\left(L\right)}---\left(13\right)$

$\×{\hat{u}}_2^{p*}(x;x_s;\mathrm{\ω}){\hat{u}}_1^{p*}(x;L;\mathrm{\ω})[W_1^p\left(L\right)D_1^p(L;p_1^p;\mathrm{\ω})+W_2^p\left(L\right)D_2^p(L;p_1^p;\mathrm{\ω})]$

$I^{\mathrm{ps}}\left(x\right)=\frac{\mathrm{\ΔL}\sqrt{{\mathrm{\ω}}_r{w}_0^2}}{16{\mathrm{\π}}^2}\underset{L}{\mathrm{\Σ}}\∫\mathrm{d\ω}\frac{\mathrm{i\ω}}{v_p^2\left(x_s\right)}\sqrt{\frac{\mathrm{\ρ}\left(L\right)}{\mathrm{\ρ}\left(x_s\right)}}\∫\∫\frac{d{p}_1^p\left(x_s\right)d{p}_1^s\left(L\right)}{p_2^p\left(x_s\right)p_2^s\left(L\right)}\mathrm{sgn}\left(\mathrm{\α}\right)---\left(14\right)$

$\×{\hat{u}}_2^{p*}(x;x_s;\mathrm{\ω}){\hat{u}}_1^{s*}(x;L;\mathrm{\ω})[W_1^s\left(L\right)D_1^s(L;p_1^s;\mathrm{\ω})+W_2^s\left(L\right)D_2^s(L;p_1^s;\mathrm{\ω})]$

Wherein, I ^Pp(x) be the single big gun imaging of P-P value, I ^Ps(x) be the single big gun imaging of P-S value, all big guns records are offset and superpose, just obtain final elastic wave imaging results.

Images