CN104749623A - Seismic data imaging processing method - Google Patents

Abstract

The invention discloses a seismic data imaging processing method. The seismic data imaging processing method includes that selecting a reference frequency and a highest frequency according to the dominant frequency and frequency width of seismic data, and confirming the initial width of a vector Gaussian beam according to the longitudinal wave velocity field and transverse wave velocity field of an exploration area so as to give an initial ray parameter interval; directly emitting elastodynamic Gaussian beams at a focus and a receiving point based on earth surface inclination information, and performing kinematical ray tracing and dynamic ray tracing to calculate the travel time and amplitude information of the elastodynamic Gaussian beams; using cross-correlation imaging conditions to calculate corresponding single-channel PP wave and PS wave imaging values according to the travel time and amplitude information of the elastodynamic Gaussian beams and polarity information, circulating each channel to obtain single-shot PP wave and PS wave imaging values, and overlapping the PP wave and PS wave imaging values of all the shots to obtain a final PP wave and PS wave imaging section.

Description

A kind of imaging of seismic data disposal route

Technical field

The invention belongs to seism processing field, specifically a kind of imaging of seismic data disposal route.

Background technology

Along with deepening continuously of exploratory development, the center of gravity of oil-gas exploration is just turning to two complex regions with complicated landform and complicated underground geologic bodies, and geological data also develops into elastic wave data by P wave data gradually.But under the MODEL OVER COMPLEX TOPOGRAPHY of near surface elevation and speed acute variation, conventional static correction can have a strong impact on migration imaging quality to the distortion that wave field causes.Because multi-wave and multi-component data comprehensively and comparatively truly can reflect underground medium and formation lithology information, the method that the two complex condition of development directly carries out Elastic Wave Migration imaging becomes the focus in seismic data process field gradually.

For solving MODEL OVER COMPLEX TOPOGRAPHY to the impact of migration imaging, many scholars have carried out large quantifier elimination.The first kind is the complicated landform offset imaging method based on ACOUSTIC WAVE EQUATION.Wiggins (1984), based on the reference field continuation of geological data, propose the Kirchhoff offset method being applicable to MODEL OVER COMPLEX TOPOGRAPHY, and Alkhalifah and Bagaini carries out base-level correction by setting velocity equivalent; Gray (2005) then proposes a kind of local static correction Gaussian beam offset method being applicable to complicated earth surface.Berryhill (1979; 1984) method of wave equation datuming is then proposed.Subsequently, Yilmaz and Lucas (1986), Schneider and Pillip etc. (1995), Bevc (1997) and Yang etc. (1999) on this basis binding layer substitute thought and have further developed the method.Above method all " offsets " formation method indirectly based on base-level correction, although base-level correction can eliminate the distortion that relief surface causes lineups effectively, the horizontal change of near-surface velocity can have a strong impact on the precision of base-level correction.For the defect of base-level correction, Reshef (1991), Beasley and Lynn (1992), He Ying etc. (2002), Shragge (2005), Cheng Jiubing (2006), Ye Yueming etc. (2008) and Lv Bin etc. (2011) successively propose multiple " directly offseting " formation method based on wave equation and solve complicated earth surface problem.Xu Yi (2008) then achieves sound wave reverse-time migration in relief surface situation with Triangular meshes method.Liu Hongwei etc. (2011) adopt filter method to earth's surface free boundary, achieve the reverse-time migration of relief surface sound wave with GPU.But under complex near surface conditions, ray class methods have the advantages such as flexible, efficient concurrently.Gray etc. (1995) have developed the method for directly carrying out Kirchhoff skew at complicated earth surface, Jager etc. (2003) propose the guarantor's width Kirchhoff offset method under relief surface condition, and Dong Chunhui etc. (2009) and Liu Guofeng etc. (2010) successively propose flow process and the application of the direct time migration of relief surface based on Kirchhoff operator.Yue etc. (2010) propose a kind of method of directly carrying out protecting the skew of width Gaussian beam at relief surface.The Gaussian beam migration operator of the above routine all grows up according to local dip superposition theory, Huang Jianpings etc. (2014) have then developed a kind of pair of complex condition non-inclined superposition based on ACOUSTIC WAVE EQUATION accurately restraints offset imaging method, and obtains good imaging effect.Equations of The Second Kind is the complicated landform offset imaging method based on equations for elastic waves.In the process of traditional multi-component earthquake data process, often first by wave field separation, then separately imaging is carried out to the P ripple after separation and S ripple.But, in the process of wave field separation, often there is non-remaining type wave energy, the much noise in imaging results can be caused, have a strong impact on image quality.In recent years, Chinese scholars, based on flatly surface condition, proposes elastic wave reverse-time migration and elastic wave Kirchhoff skew.Yue Yubo (2011) then proposes elastic wave Gaussian beam offset method, and the method has taken into account imaging precision and counting yield.And be still in the starting stage based on the research of the Elastic Wave Migration formation method of MODEL OVER COMPLEX TOPOGRAPHY, Huang Jianpings etc. (2014), according to local dip superposition theory, propose the Gaussian beam offset method based on non-separation multi-component earthquake data under a kind of relief surface condition.But, owing to having horizontal range difference and difference of elevation between bundle center and acceptance point, the method needs to introduce phase correction factor with approximate sign Green function, when horizontal range difference is larger, this approximate amplitude error caused can not be ignored, especially when landform big rise and fall and near-surface velocity change violent time, above-mentioned approximate processing greatly can reduce the effect of migration imaging.

Summary of the invention

The object of the invention is to the difficult problem for western part of China and southern mountain front exploratory area multiwave multicomponent earthquake data imaging processing, propose a kind of imaging of seismic data disposal route, at the elastic wave Gaussian beam Depth Domain imaging technique of two complex condition based on single track process, using non-separation multiwave multicomponent earthquake data and p-and s-wave velocity field as input, utilize earth's surface obliquity information, (vertical by the backward extension vector wave field characterizing decoupling zero at each acceptance point place outgoing elastokinetics Gaussian beam, shear wave field), and adopt cross-correlation image-forming condition to ask for imaging value.

The present invention includes following steps:

(1) reference frequency ω is chosen according to the dominant frequency of multiwave multicomponent earthquake data and frequency span _rand highest frequency ω _h;

(2) according to following criterion, calculating elastic dynamics Gaussian beam original width:

$w_0=\frac{2\mathrm{\π}{V}_{\mathrm{avg}}^v}{{\mathrm{\ω}}_r}$ (formula 1)

In formula: the geometrical mean of V-shaped ripple (compressional wave or shear wave) velocity field,

(3) V-shaped ripple initial ray parameter space is:

$\mathrm{\Δ}{p}^v=\frac{\mathrm{\π}}{4w_0\sqrt{{\mathrm{\ω}}_r{\mathrm{\ω}}_h}}$ (formula 2)

(4) based on earth's surface elevation, inclination angle and p-and s-wave velocity information, and in conjunction with the selected parameter also calculated of above-mentioned 1,2,3 steps, utilize single track processing mode, respectively at focus and acceptance point place direct outgoing elastokinetics Gaussian beam, and by solving kinematics ray tracing equation 3 and kinetics ray-tracing equation 4 and then calculating whilst on tour and the amplitude information of elastokinetics Gaussian beam, then the backward extension vector wave field of decoupling zero under complex near surface conditions is tried to achieve, as shown in Equation 5;

$\frac{{\mathrm{dx}}_i\left(s\right)}{\mathrm{d\τ}}=V_v^2\left(s\right)p_i\left(s\right)$

(formula 3)

$\frac{{\mathrm{dp}}_i\left(s\right)}{\mathrm{d\τ}}=-\frac{1}{V_v\left(s\right)}\frac{\∂V_v\left(s\right)}{\∂\left(x_i\right)}$

$\frac{\mathrm{dQ}}{\mathrm{ds}}=V_{v}P$

(formula 4)

$\frac{\mathrm{dP}}{\mathrm{ds}}=-V_v^{-2}\frac{{\∂}^2{V}_v}{\∂n^2}Q$

$\begin{array}{c}{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)]\end{array}$ (formula 5)

In formula, wherein, W ₁ ^v(x _r), W ₂ ^v(x _r) be the weights relevant with p-and s-wave velocity and polarization vector, for acceptance point x _rthe displacement at the x place that place causes, ρ (x _r) be x _rplace's Media density, S is the free earth's surface that rises and falls, u ₁(x _r, ω) and u ₂(x _r, ω) and be respectively horizontal and vertical component seismic record;

(5) based on the PP ripple of single track process and PS ripple imaging formula under the complex near surface conditions derived according to cross-correlation image-forming condition, the whilst on tour of elastokinetics Gaussian beam, amplitude and polarity information is utilized to calculate the imaging value of corresponding single track;

(6) the computing flow process of 4,5 steps is repeated, calculate single big gun imaging value, the imaging value of all big guns is carried out superposing and is obtained required imaging section the most at last, two complex condition based on the PP ripple of single track process and PS ripple list big gun migration imaging formula such as formula shown in 6a-6b:

$\begin{array}{c}{I}^{\mathrm{pp}}\left(x\right)=\∫U_2^{p^{*}}(x,x_s,\mathrm{\ω})u_2^p(x,x_r,\mathrm{\ω})\mathrm{d\ω}\\ =\frac{{\mathrm{\ω}}_r{w}_0^2}{16{\mathrm{\π}}^2}\∫\mathrm{d\ω}\frac{\mathrm{i\ω}}{V_p^2\left(x_s\right)V_p^2\left(x_r\right)}\sqrt{\frac{\mathrm{\ρ}\left(x_r\right)}{\mathrm{\ρ}\left(x_s\right)}}\\ \×\∫\∫\frac{{\mathrm{dp}}_1^p\left(x_s\right){\mathrm{dp}}_1^p\left(x_r\right)}{p_2^p\left(x_s\right)p_2^p\left(x_r\right)}{\hat{U}}_2^{p^{*}}(x,x_s,\mathrm{\ω}){\hat{U}}_2^{p^{*}}(x,x_r,\mathrm{\ω})\\ \×[W_1^p\left(x_r\right)u_1(x_r,\mathrm{\ω})+W_2^p\left(x_r\right)u_2(x_r,\mathrm{\ω})]\end{array}$ (formula 6a)

$\begin{array}{c}{I}^{\mathrm{ps}}\left(x\right)=\∫U_2^{p^{*}}(x,x_s,\mathrm{\ω})u_1^s(x,x_r,\mathrm{\ω})\mathrm{d\ω}\\ =\frac{{\mathrm{\ω}}_r{w}_0^2}{16{\mathrm{\π}}^2}\∫\mathrm{d\ω}\frac{\mathrm{i\ω}}{V_p^2\left(x_s\right)V_p^2\left(x_r\right)}\sqrt{\frac{\mathrm{\ρ}\left(x_r\right)}{\mathrm{\ρ}\left(x_s\right)}}\\ \×\∫\∫\frac{{\mathrm{dp}}_1^p\left(x_s\right){\mathrm{dp}}_1^s\left(x_r\right)}{p_2^p\left(x_s\right)p_2^s\left(x_r\right)}{\hat{U}}_2^{p^{*}}(x,x_s,\mathrm{\ω}){\hat{U}}_1^{s^{*}}(x,x_r,\mathrm{\ω})\\ \×[W_1^s\left(x_r\right)u_1(x_r,\mathrm{\ω})+W_2^s\left(x_r\right)u_2(x_r,\mathrm{\ω})]\mathrm{sgn}\left(\mathrm{\φ}\right)\end{array}$ (formula 6b)

In formula, for corresponding V-shaped ripple is at x _rthe initial slowness vector at place, φ is P ripple incident angle, and sgn (φ) is sign function;

Single track processing mode is utilized during the vector wave field of beneficial effect of the present invention backward extension decoupling zero under MODEL OVER COMPLEX TOPOGRAPHY, can effectively avoid the common elastic ripple Gaussian beam based on local dip superposition theory to offset amplitude error, Resolution Error and the noise problem brought, be a kind of high precision migration and imaging techniques based on non-separation multiwave multicomponent earthquake data.

Accompanying drawing explanation

Fig. 1 a and Fig. 1 b is respectively Zhongyuan Oil Field FAULT MODEL compressional wave and shear wave velocity field.This model has typical small hill cheuch complicated earth surface, and the medium and deep of model comprises well-developed fault tectonic, is the target area of oil-gas exploration.Model parameter is as follows: model meshes is 436 × 600, is respectively 4 meters and 10 meters in length and breadth to sampling interval;

Fig. 2 a and 2b is respectively corresponding single big gun z component and x component recording.Earthquake big gun data 157 big guns altogether, big gun is spaced apart 20 meters, and Mei Bao 121 road receives, per pass 4000 sampled points, and sampling rate is 0.5 millisecond, and road is spaced apart 10 meters;

Fig. 3 a and 3c is respectively PP ripple and the PS ripple imaging section of conventional method;

Fig. 3 b and 3d is respectively the PP ripple after application the present invention and PS ripple imaging section.

Embodiment

For making above and other object of the present invention, feature and advantage become apparent, cited below particularly go out preferred embodiment, and coordinate institute's accompanying drawings, be described in detail below:

(1) P-and S-wave velocity model (Fig. 1 a-1b) and corresponding unsegregated multi-wave and multi-component big gun record (Fig. 2 a-2b) conduct input is chosen.Rate pattern grid is 436 × 600, and transverse and longitudinal spacing distance is respectively 4 meters and 10 meters.Earthquake big gun data have 157 big guns, and big gun is spaced apart 20 meters, and Mei Bao 121 road receives, per pass 4000 sampled points, and sampling rate is 0.5 millisecond, and road is spaced apart 10 meters;

(2) choose suitable reference frequency 10Hz and highest frequency 50Hz according to the dominant frequency of earthquake big gun data and frequency span, and then calculate initial beamwidth and initial ray parameter space according to formula 1 and formula 2;

(3) according to the recording geometry determination focus described in the information such as earth's surface elevation, inclination angle and above-mentioned steps 1 and acceptance point in the position of complicated earth surface, then direct outgoing elastokinetics Gaussian beam calculate its whilst on tour and amplitude information;

(4) according to elastokinetics Gaussian beam whilst on tour, amplitude information and the polarity information thereof asked in step 3, application cross-correlation image-forming condition calculates corresponding single track PP ripple and PS ripple imaging value;

(5) repeat step 3 and 4, calculate single big gun PP ripple and PS ripple imaging value, the PP ripple of all big guns and PS ripple imaging value are carried out respectively superposing and are obtained required PP ripple and PS ripple imaging section Fig. 3 b and Fig. 3 d the most at last.

As can be seen from Fig. 3 b and 3d: no matter (1) is that the imaging of PP ripple or the imaging of PS ripple all clearly reflect the well-developed fault tectonic of relief surface and underground; (2) because conversion S wave velocity is less than p wave interval velocity, thus in PS ripple imaging results, lineups have the longitudinal frame higher than PP imaging results; (3) PS areas imaging is wider than PP areas imaging, and this is because conversion S wave reflection angle is less than P wave reflection angle, so wave detector can receive the conversion S ripple information from subsurface interface more broad range.(4) in model deep layer, PS ripple imaging energy Ratios PP ripple imaging energy is weak, causes this is because the conversion S wave energy produced under deep layer P ripple low-angle condition of incidence is more weak.

Comparison diagram 3a and 3b and Fig. 3 c and 3d respectively, can find: relative to conventional relief surface elastic wave Gaussian beam offset method, no matter be the imaging of PP ripple or the imaging of PS ripple, imaging effect of the present invention is all better, eliminate owing to restrainting the amplitude error existed between center and acceptance point caused by larger distance in conventional method, especially in the drawings the survey area of further object shown in dashed rectangle is shown, the imaging energy of context of methods is obviously better than conventional method, the more clearly well-developed fault tectonic of imaging, imaging results then embodies the superiority of the present invention compared to conventional method.

Claims (1)

1. the invention provides a kind of imaging of seismic data disposal route, it is characterized in that, comprise the following steps:

(1) reference frequency and highest frequency is chosen according to the dominant frequency of multiwave multicomponent earthquake data and frequency span;

(2) calculating elastic dynamics Gaussian beam original width;

(3) arranging initial ray parameter space is;

(4) according to earth's surface elevation, inclination angle and p-and s-wave velocity information, in conjunction with the selected parameter also calculated of above-mentioned 1,2,3 steps, utilize single track processing mode, respectively at focus and acceptance point place direct outgoing elastokinetics Gaussian beam, calculate whilst on tour and the amplitude information of elastokinetics Gaussian beam, then try to achieve the backward extension vector wave of decoupling zero under complex near surface conditions;

(5) based on the PP ripple of single track process and PS ripple imaging formula under the complex near surface conditions derived according to cross-correlation image-forming condition, the whilst on tour of elastokinetics Gaussian beam, amplitude and polarity information is utilized to calculate the imaging value of corresponding single track;

(6) repeat 4,5 steps, calculate single big gun imaging value, the imaging value of all big guns is carried out superposing and is obtained required imaging section the most at last.