CN108363101A - A kind of inclined shaft crosshole seismic Gaussian beam pre-stack depth migration imaging method - Google Patents

Abstract

The invention discloses a kind of inclined shaft crosshole seismic Gaussian beam pre-stack depth migration imaging methods, include the following steps：(1) files such as well earthquake reflected wave script holder record are read in；(2) seismic wave field is resolved into a series of part plan portion, obtains corresponding data volume；(3) ray tracing is carried out from shot point, the time τ that arbitrary ray is propagated is calculated using eikonal equation；(4) central ray is obtained, Gaussian beam dynamics tracking is carried out, uses Runge Kutta equation solution kinetics ray-tracing equation；(5) seismic amplitude of grid node in beam coverage area is calculated；(6) kinematics and kinetics ray-tracing are carried out from geophone station, calculates and store the attribute information that every ray corresponds to grid node within the scope of beam；(7) shot point and the corresponding beam of geophone station are chosen to carrying out imaging calculating；(8) imaging results of cumulative all beams pair.The invention has the beneficial effects that：The imaging method has high efficiency and high-precision.

Description

A Prestack Depth Migration Imaging Method of Gaussian Beam Seismic Gaussian Beams in Deviated Wells

technical field

The invention relates to an imaging method, in particular to a Gaussian beam pre-stack depth migration imaging method for inter-well seismic in inclined wells, and belongs to the technical field of material or object detection.

Background technique

During cross-well seismic exploration, the seismic source and geophone are placed in the well, which can effectively avoid the absorption of high-frequency components of seismic signals by the low-velocity zone on the surface, and obtain seismic signals with a wide frequency band and extremely high resolution. Fine imaging of geological targets such as formations, structures, and reservoirs.

Conventional crosswell seismic reflection wave migration imaging methods are divided into wave equation type migration and ray type migration.

The wave equation migration method solves the numerical solution of the wave equation, using the recursive wave field for imaging, generally including: the common shot set migration algorithm based on the one-way wave equation and the offset domain migration algorithm and the inverse of the two-way wave equation time shift algorithm. These methods require precise imaging conditions and one-way wave equation migration cannot accurately image steep structures, and the calculation efficiency is low.

Migration algorithms based on geometric ray theory realize wavefield continuation imaging by calculating the amplitude and phase of seismic wavefields, such as the Kirchhoff migration method. However, the Kirchhoff migration method is based on zero-order ray theory, and cannot accurately image caustic regions, shadow regions, and multi-arrival regions. Therefore, this imaging method often cannot meet the needs of fine exploration in complex structural regions.

"The inverse Gaussian beam common-reflection-point stack imaging in crosswell seismic" was published at the 85th SEG Annual Meeting, which introduced the inverse Gaussian beam stack imaging scheme for crosswell seismic reflection waves in deviated wells. This scheme decomposes the energy inverse Gaussian at each detection point to the actual position of the reflection point, and performs homing imaging on the reflected wave data, which effectively improves the conventional VSP-CDP imaging scheme.

CN104391327A discloses "A Pre-stack Reverse Time Depth Migration Imaging Method for Offshore Deviated Well Interwell Seismic", which introduces a depth domain migration imaging scheme for interwell seismic data under the condition of inclined well. The scheme uses the TTI medium tomography inversion method to obtain the migration model, and divides it into grids, then uses the first-arrival ray tracing method to calculate the excitation time imaging conditions, and finally performs wave field time inversion on the interwell seismic reflection wave field The wave field values at different times are obtained by continuation, and imaging conditions are applied to each grid to obtain the cross-well seismic depth migration imaging section.

It can be seen that the existing imaging methods can solve the seismic reflection wavefield migration imaging between deviated wells to a certain extent, but the calculation accuracy of the ray-based imaging method is low, and the wave equation-based imaging method is complicated to implement and the calculation efficiency is relatively low. .

Contents of the invention

In order to solve the deficiencies of the prior art, the object of the present invention is to provide a high-efficiency and high-precision inter-well seismic Gaussian beam pre-stack depth migration imaging method.

A pre-stack depth migration imaging method for inter-well seismic Gaussian beams in deviated wells, characterized in that it comprises the following steps:

Step1: Read in the interwell seismic reflection wave field records, observation system files, migration velocity files and related parameter files, which include: the number of grid points involved in the calculation of the seismic work area, the grid spacing, the initial width of the Gaussian beam, and the reference Frequency, maximum frequency and number of seismic recording sampling points;

Step2: The wavefront of the Gaussian beam at the source position is characterized by plane waves, and the seismic wave field is decomposed into a series of local plane parts to obtain the corresponding data volume;

Step3: Carry out ray tracing from the shot point along the downward direction of the excitation well to the upward direction of the excitation well, and calculate the propagation time of any ray (that is, the travel time of seismic waves) τ by using the equation of the equation:

Among them, v is the velocity value at the discrete point;

Step4: Gaussian beam dynamics tracking is performed after obtaining the center ray through kinematics tracking, and the Runge-Kutta equation is used to solve the dynamic ray tracing equation:

Among them, n is the vertical distance from the point on the adjacent ray to the central ray, and the initial values of p and q are respectively the component of the slowness vector in the direction perpendicular to the central ray and the distance from the adjacent ray to the central ray;

Step5: Calculate the seismic wave amplitude of the grid nodes within the coverage of the ray beam according to the central ray:

Among them, A ₀ is the amplitude of the shot point, q(R) is the distance between the adjacent ray at the receiver point and the central ray, N is the number of formations that the ray passes through when it reaches the receiver point, and R _i is the reflection coefficient of the i-th interface or transmission coefficient, α _i , β _i are the incident angle and transmission angle at the i-th interface, respectively, ρ _i (R), v _i (R) are the formation density and velocity before the ray passes through the i-th interface, respectively, Respectively, the formation density and velocity after the ray passes through the i-th interface;

Step6: Carry out kinematic and dynamic ray tracing from the receiver point along the downward direction of the receiving well to the upward direction of the receiving well, calculate and store the attribute information of the grid nodes within the range of each ray corresponding to the ray beam, the attribute information includes: seismic wave amplitude and travel time;

Step7: Select the ray beam pairs corresponding to the shot point and the receiver point for imaging calculation, and the calculation equation is as follows:

Among them, I _s (x) is the imaging value at point x, p _s is the slowness value of the ray from the shot point, p _bc is the slowness value of the ray from the receiver point, A is the weight function, D _s is the local plane wave decomposition As a result, L is the different windows divided in the single shot record, p' is the slowness parameter of the shot record for local oblique stacking, τ' is the travel time of the ray beam from the shot point to the imaging point and the arrival time of the ray beam from the receiver point to the imaging point The sum of travel times at the point;

Step8: Accumulate the imaging results of all ray beam pairs to obtain the final migration imaging result.

The benefits of the present invention are:

1. The cross-well seismic Gaussian ray beam propagation operator is used to solve the problem of low imaging efficiency of the wave equation method, and overcome the shortcomings of the ground seismic imaging method such as weak imaging of deep structures and insufficient resolution;

2. The use of dynamic ray tracing overcomes the difficulty that ray-like migration imaging methods cannot perform fine imaging in complex structural regions (such as caustic regions, singular regions, etc.);

3. The Gaussian beam has a certain effective width, and the superposition of the Gaussian beams that contribute to the calculation point can be selected to calculate the final wave field, which effectively improves the calculation efficiency and is an imaging method that takes into account both efficiency and precision.

Description of drawings

Fig. 1 is the flow chart of the inter-well seismic Gaussian beam pre-stack depth migration imaging method of deviated wells of the present invention;

Fig. 2(a) is a schematic diagram of the ray beam propagation range at the excitation point of cross-well seismic Gaussian beam migration;

Fig. 2(b) is a schematic diagram of the ray beam propagation range at the receiving point of cross-well seismic Gaussian beam migration;

Fig. 3 is a velocity model diagram;

Fig. 4 is the cross-well seismic Gaussian beam prestack depth migration imaging section;

Figure 5 is the surface seismic section through the well in the actual data imaging section of a certain block of CNOOC;

Figure 6 is the cross-well seismic Gaussian beam pre-stack depth migration imaging section of the actual data of a certain block of CNOOC.

Detailed ways

The present invention will be specifically introduced below in conjunction with the accompanying drawings and specific embodiments.

Referring to Fig. 1, the inter-well seismic Gaussian beam pre-stack depth migration imaging method of the present invention comprises the following steps:

Step1: read in the file

Read in the interwell seismic reflection wave field records, observation system files, migration velocity files and related parameter files, where the related parameter files include: the number of grid points involved in the calculation of the seismic work area, the grid spacing, the initial width of the Gaussian beam, the reference Frequency, maximum frequency and number of seismic record sampling points.

Step2: Decompose the seismic wave field

The wavefront of the Gaussian beam at the source position is characterized by plane waves, and the seismic wave field is decomposed into a series of local plane parts to obtain the corresponding data volume.

Step3: Carry out ray tracing from the shot point and calculate the seismic wave travel time

Ray tracing is carried out from the shot point along the downward direction of the excitation well to the upward direction of the excitation well, and the propagation time of any ray (that is, the travel time of seismic waves) τ is calculated by using the equation of the equation:

Among them, v is the velocity value at the discrete point.

Figure 2(a) shows the propagation range of the ray beam emitted by the shot point (excitation point) of the Gaussian beam prestack depth migration shot point (excitation point) at a certain angle in the inclined well interwell seismic Gaussian beam.

Step4: Solve the dynamic ray tracing equation

Gaussian beam dynamics tracing is performed after the central ray is obtained through kinematics tracing, and the Runge-Kutta equation is used to solve the dynamic ray tracing equation:

Among them, n is the vertical distance from the point on the adjacent ray to the central ray, and the initial values of p and q are respectively the component of the slowness vector in the direction perpendicular to the central ray and the distance from the adjacent ray to the central ray.

Step5: Obtain the amplitude information of the grid nodes within the coverage of the ray beam

Compute seismic wave amplitudes at grid nodes within the coverage of the ray beam from the central ray:

Among them, A ₀ is the amplitude of the shot point, q(R) is the distance between the adjacent ray at the receiver point and the central ray, N is the number of formations that the ray passes through when it reaches the receiver point, and R _i is the reflection coefficient of the i-th interface or transmission coefficient, α _i , β _i are the incident angle and transmission angle at the i-th interface, respectively, ρ _i (R), v _i (R) are the formation density and velocity before the ray passes through the i-th interface, respectively, (R), v _i (R) are the formation density and velocity after the ray passes through the i-th interface, respectively.

Step6: Perform ray tracing from the receiver point and calculate the attribute information of the grid node

Carry out kinematics and dynamics ray tracing from the receiver point along the downward direction of the receiving well to the upward direction of the receiving well, calculate and store the attribute information of the grid nodes within the range of each ray beam, the attribute information includes: seismic wave amplitude and travel time .

Figure 2(b) shows the propagation range of the ray beam emitted by the prestack depth migration receiver point (receiving point) of the cross-well seismic Gaussian beam in a deviated well at a certain angle.

Step7: Carry out imaging calculation

The ray beam pairs corresponding to the shot point and the receiver point are selected for cross-correlation imaging calculation, and the calculation equation is as follows:

Among them, I _s (x) is the imaging value at point x, p _s is the slowness value of the ray from the shot point, p _bc is the slowness value of the ray from the receiver point, A is the weight function, D _s is the local plane wave decomposition As a result, L is the different windows divided in the single shot record, p' is the slowness parameter of the shot record for local oblique stacking, τ' is the travel time of the ray beam from the shot point to the imaging point and the arrival time of the ray beam from the receiver point to the imaging point The sum of travel times at the point.

Step8: Accumulate imaging results

The imaging results of all ray beam pairs are accumulated to obtain the final migration imaging result.

In order to verify that the imaging method provided by the present invention has a good imaging effect, we conducted cross-well seismic data imaging tests in the forward modeling model and a block of CNOOC respectively.

1. Conduct cross-well seismic data imaging test on the forward modeling model

The velocity model (forward model) is shown in Figure 3.

It can be seen from Fig. 3 that there are reverse faults, pinch-out of sand bodies and formations with inclined occurrences between the two wells.

We adopt the method provided by the invention to carry out imaging, wherein:

(1) Parameters read in: crosswell seismic reflection wavefield record (reflectwavefield.cds), observation system file (geom.lge), smoothed velocity file (velosmooth.dat).

(2) Relevant parameter files set: grid point nx, grid spacing dx, reference frequency f, maximum frequency f _max , seismic record sampling point nt, and Gaussian beam initial width w involved in the calculation of the seismic work area.

(3) Ray tracing is carried out from the shot point along the downward direction of the excitation well to the upward direction of the excitation well, and the incident angle range is -20°～160°.

(4) Carry out kinematic and dynamic ray tracing from the receiver point along the downward direction of the receiving well to the upward direction of the receiving well, and the incident angle: 30°～165°.

The cross-well seismic Gaussian beam prestack depth migration imaging section of the deviated well is shown in Fig. 4.

It can be seen from FIG. 4 that the imaging method of the present invention corresponds very well to the model for the stratum imaging effect of the reverse fault and inclined occurrence shown in FIG. 3 .

It can be seen that the Gaussian beam prestack depth migration imaging has been successfully performed on a typical inter-well complex structure model by using the imaging method provided by the present invention, and the correctness, effectiveness and stability of the method have been verified.

2. Conducted cross-well seismic data imaging test in a certain block of CNOOC

Figure 5 shows the cross-section of the ground seismic exploration well in a certain block of CNOOC.

It can be seen from Fig. 5 that the target layer is between 0.4s and 0.6s, and micro-fault structures are developed in the target layer, and the resolution of microstructures in the target area by ground seismic imaging method is relatively low.

We adopt the method provided by the invention to carry out imaging, wherein:

(1) Files read in: reflected wavefield record (oriwavefield.cds), observation system file (geometry.lge), smoothed velocity file (velo_smooth.dat).

(2) Relevant parameter files set: grid point nx, grid spacing dx, reference frequency f, maximum frequency f _max , seismic record sampling point nt, and Gaussian beam initial width w involved in the calculation of the seismic work area.

(3) Ray tracing is carried out from the shot point along the downward direction of the excitation well to the upward direction of the excitation well, and the incident angle range is -5°～85°.

(4) Carry out kinematic and dynamic ray tracing from the receiver point along the downward direction of the receiving well to the upward direction of the receiving well, and the incident angle: -160°～20°.

The cross-well seismic Gaussian beam pre-stack depth migration imaging section of the deviated well interwell seismic data in this block of CNOOC is shown in Fig. 6.

It can be seen from Fig. 6 that the interlayer structural characteristics of the cross-well seismic Gaussian beam pre-stack depth migration imaging section of the deviated well are more clearly described, the stratum occurrence is basically consistent with the surface seismic imaging results, and the CDP interval in the surface seismic section is 12.5 m, and the CDP interval in the interwell seismic migration profile is 3.125m.

In addition, it can also be seen from Fig. 6 that the inter-well seismic Gaussian beam prestack depth migration imaging method for deviated wells provided by the present invention can finely describe the geological structure, and the resolution of the imaging result is higher.

It can be seen that the application of the imaging method provided by the present invention to Gaussian beam pre-stack depth migration imaging processing of the actual data of a CNOOC interwell seismic exploration block with a relatively complex structure has obtained better geological imaging effects.

It should be noted that the above embodiments do not limit the present invention in any form, and all technical solutions obtained by means of equivalent replacement or equivalent transformation fall within the protection scope of the present invention.

Claims (3)

1. a kind of inclined shaft crosshole seismic Gaussian beam pre-stack depth migration imaging method, which is characterized in that include the following steps：

Step1：Read in well earthquake reflected wave script holder record, observation system file, migration velocity file and relevant parameter file；

Step2：In hypocentral location, the wavefront features of Gaussian beam are plane wave, and seismic wave field is resolved into a series of part plan Portion obtains corresponding data volume；

Step3：Ray tracing in downward direction is carried out to excitation well upward direction from shot point along excitation well, utilizes eikonal equation meter Calculate the time τ that arbitrary ray is propagated：

Wherein, v is the velocity amplitude at discrete point；

Step4：It is tracked by kinematics and carries out Gaussian beam dynamics tracking after obtaining central ray, asked using Runge Kutta equation Solve kinematics ray tracing equation：

Wherein, n is vertical range of the point on adjacent ray to central ray, and the initial value of p, q are respectively perpendicular to central ray The component of slowness vector and adjacent ray leave the distance of central ray on direction；

Step5：The seismic amplitude of grid node in beam coverage area is calculated according to central ray：

Wherein, A₀For shot point amplitude, q (R) is the distance that adjacent ray leaves central ray at geophone station, and N is that ray reaches detection The ground number of plies passed through when point, R_iFor the reflectance factor or transmission coefficient at i-th of interface, α_i、β_iRespectively i-th interface enters Firing angle and angle of transmission, ρ_i(R)、v_i(R) it is respectively density of earth formations and speed of the ray before i-th of interface, Respectively density of earth formations and speed of the ray behind i-th of interface；

Step6：Ray tracing in downward direction is carried out to received well upward direction from geophone station along received well, calculates and stores every Ray corresponds to the attribute information of grid node within the scope of beam；

Step7：Shot point and the corresponding beam of geophone station are chosen to carrying out imaging calculating, accounting equation is as follows：

Wherein, I_s(x) be point x at picture value, p_sThe slowness value of ray, p are sent out for shot point_bcThe slow of ray is sent out for geophone station Angle value, A are weight function, D_sFor Local plane wave decomposition as a result, L is the different window divided in single shot record, p ' is that big gun record is used In the slowness parameter of local dip superposition, τ ' sends out for shot point and sends out ray with geophone station when walking at beam arrival imaging point Beam reach imaging point at the sum of when walking；

Step8：The imaging results of cumulative all beams pair, obtain final migration imaging result.

2. inclined shaft crosshole seismic Gaussian beam pre-stack depth migration imaging method according to claim 1, which is characterized in that In Step1, the relevant parameter file includes：Grid dimension, grid spacing, the Gaussian beam for participating in calculating earthquake work area are initially wide Degree, reference frequency, maximum frequency and earthquake record sampling number.

3. inclined shaft crosshole seismic Gaussian beam pre-stack depth migration imaging method according to claim 1, which is characterized in that In Step6, the attribute information includes：The amplitude of seismic wave and when walking.