CN102313903B - Pre-stack time migration method in VTI medium based on wave equation extrapolation operator - Google Patents

Abstract

The invention discloses a pre-stack time migration method in a VTI medium based on wave equation extrapolation operator. The method comprises the following steps: (1) inputting CMP gather and time domain velocity field; (2) carrying out local slant stack to the CMP gather, mapping the CMP gather to a plane wave data body of intercept time-ray parameter field, sorting the plane wave data body according to the ray parameter so as to form a total ray parameter profile; (3) cycling the ray parameter; (4) determining whether a computer node is a master control node according to a process and entering into a master-slave parallel mode; (6) making all processes to be synchronous; (7) outputting the imaging data body corresponding to the current ray parameter by the master control node; (8) determining whether ray parameter cycle can be completed, and turning to a step 3 if the ray parameter cycle does not completed; (9) outputting a time offset profile and imaging gather. By using the method, the fidelity imaging gather can be provided. The method can provide better basic data for subsequent AVO/AVA analysis, which is beneficial to accurate oil and gas reservoir lithology prediction.

Description

Prestack time migration method in VTI medium based on wave equation extrapolation operator

Technical field

The present invention relates to a kind of prestack time migration technique, relate in particular to prestack time migration method in a kind of VTI medium.

Background technology

In the very not complicated region of tectonic structure, it is the gordian technique of oil and gas reservoir prediction that AVO (variation of amplitude the offset distance)/AVA (variation of amplitude incident angle) on the imaging road collection after pre-stack time migration analyzes.Prestack time migration technique is improving seismic imaging precision, for having brought into play vital role in the aspects such as the realization of prestack inversion.Current VTI (having the transverse anisotropy of vertical axis of symmetry) medium time migration is to use Kirchhoff integral method VTI pre-stack time migration algorithm substantially, that adopt is the Problem of Fourth-order Differential Equations (Alkhalifah with equivalent anisotropic Parametric Representation, T, 1997, Velocity analysis using nonhyperbolic moveout in transversely isotropic media, .Geophysics, 62,1839-1854), its in essence or based on layered medium hypothesis lower derive time apart from relation.Conventionally in the business software system adopting, be all to provide Offset imaging road collection with Kirchhoff integral method pre-stack time migration substantially, and carry out on this basis AVO analysis.In fact, Kirchhoff offset method use be root-mean-square velocity field, the ability of speed-adaptive horizontal change a little less than, when velocity field horizontal change, the amplitude weight factor that the Green function based in constant speed medium calculates is normally coarse.Therefore, imaging road collection amplitude infidelity, imaging road that Kirchhoff integral method provides are concentrated and are existed illusion to disturb, and such imaging road collection is unfavorable for that follow-up AVO analyzes and reservoir prediction.Attenuation by absorption and the anisotropic impact that imaging road collection and AVO are analyzed also will be eliminated in imaging process.Further, along with the continuous intensification of exploratory area, each oil field Songliao basin, more and more higher to the requirement of migration imaging precision, isotropy treatment technology has been difficult to meet the requirement of exploratory development.

Up to now, yet there are no the pre-stack time migration algorithm based on wave equation extrapolation operator.With Kirchhoff offset method Comparatively speaking, wave equation extrapolation operator is more natural to the processing of amplitude, and the adaptive faculty stronger to speed horizontal change.Therefore, be necessary to develop prestack time migration method in a set of VTI medium based on wave equation extrapolation operator.

Summary of the invention

For above-mentioned prior art, from adapting to slow lateral velocity variation, improve imaging resolution (attenuation by absorption compensation) and considering that the aspects such as anisotropic consider, and the invention provides prestack time migration method in a kind of VTI medium based on wave equation extrapolation operator.Be by offset distance plane wave wave equation prestack time migration method to VTI medium, thereby obtain the pre-stack time migration algorithm based on wave equation in VTI medium, the present invention has better guarantor's amplitude ability, the imaging road collection of fidelity can be provided, for follow-up AVO/AVA analyzes, provide better basic data, be conducive to more accurate oil and gas reservoir lithology prediction.The present invention is more accurately efficient wave equation prestack time migration method of one, can compensate the attenuation by absorption of the earth and process the imaging in anisotropic medium, good pre-stack time migration stacked section can be provided, and output ray parameter or angle domain imaging road collection are analyzed for follow-up residual velocity analysis and AVO.

In order to solve the problems of the technologies described above, the technical scheme that in a kind of VTI medium based on wave equation extrapolation operator of the present invention, prestack time migration method is achieved comprises the following steps:

Step 1: input CMP road collection and time domain velocity field;

Step 2: DuiCMP road collection carries out local dip stack, JiangCMP road collection is mapped to the plane wave data volume in intercept time-ray parameter territory; And according to ray parameter to the sorting of plane wave data volume, form cascode line parameter profile;

Step 3: ray parameter is circulated;

Step 4: process judges that whether place computer node is main controlled node, enters MS master-slave parallel schema thus; If that is: process is moved on main controlled node, carry out following (5-A) process, if process is not moved on main controlled node, but move on computing node, carry out following (5-B) process;

(5-A) first, read cascode line parameter profile corresponding to current ray parameter, then along time orientation, make one dimension fast fourier transform, the cascode line parameter plane wave datum body in generated frequency territory; Meanwhile, setting general assignment number is the minimum value among sum frequency number and computing node number; Distribution and the collection of parallel task in main controlled node control, and carry out following steps:

(5-A1) general assignment number is circulated;

(5-A2) the plane wave data volume of a frequency slice is issued to computing node, send the counter of task from increasing;

(5-A3) judge that whether task transmission finishes, and if not, turns to (5-A1); If so, turn to (5-A4);

(5-A4) sum frequency number is circulated;

(5-A5) receive the task identifier that computing node returns;

(5-A6) check and sent to the task number of computing node whether to be less than sum frequency number; If so, turn to (5-A7); If not, turn to (5-A8);

(5-A7) the plane wave data volume of current frequency slice is issued to computing node, send the counter of task from increasing, turn to (5-A4);

(5-A8) to all computing nodes, send task end identifier;

(5-B) first, judge whether this computing node number is greater than total frequency number; If so, turn to step 6; If not, turn to (5-B1);

(5-B1) process enters endless loop;

(5-B2) process is waited for and is received the data of sending from main controlled node;

(5-B3) process judges whether main controlled node sends end identifier; If so, turn to (5-B9); If not, turn to (5-B4);

(5-B4) along time orientation, carry out wave field extrapolation;

(5-B5) phase velocity of calculating offset distance plane wave propagation solves following implicit equation:

$f\left(\mathrm{\θ}\right)=2\frac{\sin \mathrm{\θ}}{v_p(\mathrm{\θ};v_{p0},\mathrm{\ϵ},\mathrm{\δ})}-p_h=0---\left(1\right)$

In formula (1), plane wave phase velocity v _p(θ; v _p0, ε, δ) can be expressed as:

$\frac{{V_{\mathrm{qP}}}^2\left(\mathrm{\θ}\right)}{{v_{p0}}^2}=\frac{1}{2}+\mathrm{\ϵ}{\sin }^2\mathrm{\θ}+\frac{1}{2}\sqrt{{(1+2\mathrm{\ϵ}{\sin }^2\mathrm{\θ})}^2-2(\mathrm{\ϵ}-\mathrm{\δ}){\sin }^{2}2\mathrm{\θ}}---\left(2\right)$

In formula (2), ε, δ are VTI medium anisotropy parameter, and θ is phase angle, v _p0for in VTI medium along axis of symmetry (Z axis) direction earthquake phase velocity of wave; By formula (1)～(2), can draw by ray parameter p _hthe phase velocity v of the offset distance plane wave propagation characterizing _p(p _h);

(5-B6) to the offset distance plane wave equation in VTI medium, adopt implicit expression finite difference scheme to carry out wave field extrapolation, realize offset distance plane wave pre-stack time migration; The offset equation adopting is:

$(\frac{4}{v_p^2\left(p_h\right)}\frac{{\∂}^2}{{\∂\mathrm{\τ}}^2}+\frac{{\∂}^2}{{\∂x}^2})\tilde{U}(x,\mathrm{\τ};\mathrm{\ω})+4\frac{{\mathrm{\ω}}^2}{v_p^2\left(p_h\right)}(1-\frac{v_p^2\left(p_h\right){p_h}^2}{4})\tilde{U}(x,\mathrm{\τ};\mathrm{\ω})=0---\left(3\right)$

In formula (3),

for frequency field plane wave data, ω is angular frequency, v _p(p _h) be the phase velocity of the plane wave propagation that calculates in (2);

Adopt zero-time image-forming condition, obtain cascode line parameter field imaging road collection;

$I\left({x,\mathrm{\τ};p}_h\right)=\underset{\mathrm{\ω}}{\∫}\tilde{U}(x,\mathrm{\τ};\mathrm{\ω})\·{\tilde{U}}^{*}(x,\mathrm{\τ};\mathrm{\ω})\mathrm{d\ω}---\left(4\right)$

(5-B7) judge whether wave field extrapolation finishes; If not, turn to (5-B4); If so, turn to (5-B8);

(5-B8) to main controlled node, send task end identifier, turn to (5-B1);

(5-B9) process on all computing nodes is to imaging data body corresponding to the current ray parameter of main controlled node stipulations, and stipulations functional symbol is summation;

Step 6: all Process Synchronizations;

Step 7: main controlled node is by imaging data body output corresponding current ray parameter;

Step 8: judge whether ray parameter circulation completes; If not, turn to step 3;

Step 9: output time migrated section and imaging road collection.

Compared with prior art, the invention has the beneficial effects as follows:

Prestack time method and technology in present business software system is all Kirchhoff integration method substantially.Mainly be to provide a kind of fast, adapt to the formation method of irregular system.It is unsuitable for meticulous fidelity imaging and follow-up AVO/AVA analyzes and reservoir prediction.The imaging road that it provides is concentrated and is existed illusion to disturb, and hi-fi of amplitude is poor, is unsuitable for follow-up AVO and analyzes.Kirchhoff integral method pre-stack time migration is difficult to combine with attenuation by absorption compensation.Some other wave equation migration method based on single square root or double square root equation all has huge calculated amount, can not be on a large scale for the production of reality, and this is also one of reason that wave equation pre-stack time migration can not be universal in actual applications.

The wave field extrapolation operator based on wave equation in VTI medium of the present invention in prestack time migration method can provide the imaging section of more protecting amplitude; At plane wave zone, realize skew and there is very high counting yield; Utilize the present invention in wave field extrapolation process, to utilize the present invention can adaptive Calculation Plane phase velocity of wave and without ray tracing; Can export ray parameter or angle domain imaging road collection and be used for doing AVO analysis or residual velocity analysis simultaneously.

To sum up, compared with prior art, in a kind of VTI medium of the present invention, the advantage having of prestack time migration method is:

(1) be the imaging of Offset plane wave data, can there is very high counting yield;

(2) be method of finite difference imaging, use be that time domain interval velocity is not root-mean-square velocity, can adapt to slow lateral velocity variation;

(3) be wave equation imaging, can obtain good hi-fi of amplitude imaging effect;

(4) can produce easily the angular-trace gather for AVO/AVA analysis and migration velocity analysis;

(5) can process easily attenuation by absorption compensation;

(6) can carry out the pre-stack time migration of VTI medium, real data migration imaging successful is better than similar Kirchhoff integral offset method.

Accompanying drawing explanation

Accompanying drawing is prestack time migration method process flow diagram in VTI medium of the present invention.

Embodiment

Below in conjunction with embodiment, the present invention is described in further detail.

As shown in drawings, prestack time migration method in a kind of VTI medium based on wave equation extrapolation operator of the present invention, comprises the following steps:

Step 1: input CMP road collection and time domain velocity field;

Step 2: DuiCMP road collection carries out local dip stack, JiangCMP road collection is mapped to the plane wave data volume in intercept time-ray parameter territory; And according to ray parameter to the sorting of plane wave data volume, form cascode line parameter profile;

Step 3: ray parameter is circulated;

Step 4: process judges that whether place computer node is main controlled node, enters MS master-slave parallel schema thus;

If that is: process is moved on main controlled node, carry out following (5-A) process, if process is not moved on main controlled node, but move on computing node, carry out following (5-B) process;

(5-A) first, read cascode line parameter profile corresponding to current ray parameter, then along time orientation, make one dimension fast fourier transform, the cascode line parameter plane wave datum body in generated frequency territory; Meanwhile, setting general assignment number is the minimum value among sum frequency number and computing node number; Distribution and the collection of parallel task in main controlled node control, and carry out following steps:

(5-A1) general assignment number is circulated;

(5-A2) the plane wave data volume of a frequency slice is issued to computing node, send the counter of task from increasing;

(5-A3) judge that whether task transmission finishes, and if not, turns to (5-A1); If so, turn to (5-A4);

(5-A4) sum frequency number is circulated;

(5-A5) receive the task identifier that computing node returns;

(5-A6) check and sent to the task number of computing node whether to be less than sum frequency number; If so, turn to (5-A7); If not, turn to (5-A8);

(5-A7) the plane wave data volume of current frequency slice is issued to computing node, send the counter of task

From increasing, turn to (5-A4);

(5-A8) to all computing nodes, send task end identifier;

(5-B) first, judge whether this computing node number is greater than total frequency number; If so, turn to step 6; If not, turn to (5-B1);

(5-B1) process enters endless loop;

(5-B2) process is waited for and is received the data of sending from main controlled node;

(5-B3) process judges whether main controlled node sends end identifier; If so, turn to (5-B9); If not, turn to (5-B4);

(5-B4) along time orientation, carry out wave field extrapolation;

(5-B5) phase velocity of calculating offset distance plane wave propagation solves following implicit equation:

$f\left(\mathrm{\θ}\right)=2\frac{\sin \mathrm{\θ}}{v_p(\mathrm{\θ};v_{p0},\mathrm{\ϵ},\mathrm{\δ})}-p_h=0---\left(1\right)$

In formula (1), plane wave phase velocity v _p(θ; v _p0, ε, δ) can be expressed as:

$\frac{{V_{\mathrm{qP}}}^2\left(\mathrm{\θ}\right)}{{v_{p0}}^2}=\frac{1}{2}+\mathrm{\ϵ}{\sin }^2\mathrm{\θ}+\frac{1}{2}\sqrt{{(1+2\mathrm{\ϵ}{\sin }^2\mathrm{\θ})}^2-2(\mathrm{\ϵ}-\mathrm{\δ}){\sin }^{2}2\mathrm{\θ}}---\left(2\right)$

In formula (2), ε, δ are VTI medium anisotropy parameter, and θ is phase angle, v _p0for in VTI medium along axis of symmetry (Z axis) direction earthquake phase velocity of wave; By formula (1)～(2), can draw by ray parameter p _hthe phase velocity v of the offset distance plane wave propagation characterizing _p(p _h);

(5-B6) to the offset distance plane wave equation in VTI medium, adopt implicit expression finite difference scheme to carry out wave field extrapolation, realize offset distance plane wave pre-stack time migration; The offset equation adopting is:

$(\frac{4}{v_p^2\left(p_h\right)}\frac{{\∂}^2}{{\∂\mathrm{\τ}}^2}+\frac{{\∂}^2}{{\∂x}^2})\tilde{U}(x,\mathrm{\τ};\mathrm{\ω})+4\frac{{\mathrm{\ω}}^2}{v_p^2\left(p_h\right)}(1-\frac{v_p^2\left(p_h\right){p_h}^2}{4})\tilde{U}(x,\mathrm{\τ};\mathrm{\ω})=0---\left(3\right)$

In formula (3),

for frequency field plane wave data, ω is angular frequency, v _p(p _h) be the phase velocity of the plane wave propagation that calculates in (2);

Adopt zero-time image-forming condition, obtain cascode line parameter field imaging road collection;

$I\left({x,\mathrm{\τ};p}_h\right)=\underset{\mathrm{\ω}}{\∫}\tilde{U}(x,\mathrm{\τ};\mathrm{\ω})\·{\tilde{U}}^{*}(x,\mathrm{\τ};\mathrm{\ω})\mathrm{d\ω}---\left(4\right)$

(5-B7) judge whether wave field extrapolation finishes; If not, turn to (5-B4); If so, turn to (5-B8);

(5-B8) to main controlled node, send task end identifier, turn to (5-B1);

(5-B9) process on all computing nodes is to imaging data body corresponding to the current ray parameter of main controlled node stipulations, and stipulations functional symbol is summation;

Step 6: all Process Synchronizations;

Step 7: main controlled node is by imaging data body output corresponding current ray parameter;

Step 8: judge whether ray parameter circulation completes; If not, turn to step 3;

Step 9: output time migrated section and imaging road collection.

Although in conjunction with figure, invention has been described above; but the present invention is not limited to above-mentioned embodiment; above-mentioned embodiment is only schematic; rather than restrictive; those of ordinary skill in the art is under enlightenment of the present invention; in the situation that not departing from aim of the present invention, can also make a lot of distortion, within these all belong to protection of the present invention.

Claims (1)

1. a prestack time migration method in the VTI medium based on wave equation extrapolation operator, is characterized in that: comprise the following steps:

Step 1: input CMP road collection and time domain velocity field;

Step 2: DuiCMP road collection carries out local dip stack, JiangCMP road collection is mapped to the plane wave data volume in intercept time-ray parameter territory; And according to ray parameter to the sorting of plane wave data volume, form cascode line parameter profile;

Step 3: ray parameter is circulated;

Step 4: process judges that whether place computer node is main controlled node, enters MS master-slave parallel schema thus;

If that is: process is moved on main controlled node, carry out following (5-A) process, if process is not moved on main controlled node, but move on computing node, carry out following (5-B) process;

(5-A) first, read cascode line parameter profile corresponding to current ray parameter, then along time orientation, make one dimension fast fourier transform, the cascode line parameter plane wave datum body in generated frequency territory; Meanwhile, setting general assignment number is the minimum value among sum frequency number and computing node number; Distribution and the collection of parallel task in main controlled node control, and carry out following steps:

(5-A1) general assignment number is circulated;

(5-A2) the plane wave data volume of a frequency slice is issued to computing node, send the counter of task from increasing;

(5-A3) judge that whether task transmission finishes, and if not, turns to (5-A1); If so, turn to (5-A4);

(5-A4) sum frequency number is circulated;

(5-A5) receive the task identifier that computing node returns;

(5-A6) check and sent to the task number of computing node whether to be less than sum frequency number; If so, turn to (5-A7); If not, turn to (5-A8);

(5-A7) the plane wave data volume of current frequency slice is issued to computing node, send the counter of task from increasing, turn to (5-A4);

(5-A8) to all computing nodes, send task end identifier;

(5-B) first, judge whether this computing node number is greater than total frequency number; If so, turn to step 6; If not, turn to (5-B1);

(5-B1) process enters endless loop;

(5-B2) process is waited for and is received the data of sending from main controlled node;

(5-B3) process judges whether main controlled node sends task end identifier; If so, turn to (5-B9); If not, turn to (5-B4);

(5-B4) along time orientation, carry out wave field extrapolation;

(5-B5) phase velocity of calculating offset distance plane wave propagation solves following implicit equation:

（1）

In formula (1), plane wave phase velocity

can be expressed as:

（2）

In formula (2),

,

for VTI medium anisotropy parameter,

for phase angle, for in VTI medium along axis of symmetry Z-direction earthquake phase velocity of wave; By formula (1) ~ (2), can draw by ray parameter

the phase velocity of the offset distance plane wave propagation characterizing

wherein,

;

(5-B6) to the offset distance plane wave equation in VTI medium, adopt implicit expression finite difference scheme to carry out wave field extrapolation, realize offset distance plane wave pre-stack time migration; The offset equation adopting is:

（3）

In formula (3),

for frequency field plane wave data,

for lateral coordinates, under two-dimensional case,

for big gun line direction,

for the vertical two way travel time of seismic event, represent pseudotime territory,

for angular frequency,

for the phase velocity of the plane wave propagation that calculates in (2);

Adopt zero-time image-forming condition, obtain cascode line parameter field imaging road collection;

（4）

(5-B7) judge whether wave field extrapolation finishes; If not, turn to (5-B4); If so, turn to (5-B8);

(5-B8) to main controlled node, send task end identifier, turn to (5-B1);

(5-B9) process on all computing nodes is to imaging data body corresponding to the current ray parameter of main controlled node stipulations, and stipulations functional symbol is summation;

Step 6: all Process Synchronizations;

Step 7: main controlled node is by imaging data body output corresponding current ray parameter;

Step 8: judge whether ray parameter circulation completes; If not, turn to step 3;

Step 9: output time migrated section and imaging road collection.

Images