CN102590857A - True surface relief prestack depth domain two-way wave imaging method - Google Patents

Abstract

The invention discloses a true surface relief prestack depth domain two-way wave imaging method, which specifically comprises the following steps of: (1) acquiring the true surface elevation of a plane with seismic single-shot data by taking the form of a relief surface datum as a basis; (2) redefining a To value, and correcting a speed model; (3) constructing source wavelets, realizing the forward modeling of a source wave field by adopting an optimized difference method, and applying a true surface boundary condition; and (4) performing regular reverse time extrapolation on recorded data along a time direction, applying the relief surface boundary condition, and performing a seismic wave field correlation to acquire a depth domain seismic image. According to the method, the depth domain wave field extrapolation and imaging of the recorded data on a true surface are realized, so that the shortcomings of wave field distortion and untruthful structural form caused by a high surface elevation difference in the conventional technologies can be overcome, and an accurate seismic image and a true geological structural form can be obtained under a complex surface condition.

Description

True surface relief prestack depth territory round trip wave imaging method

Technical field

The present invention relates to the seismic exploration technique field, particularly a kind of true surface relief prestack depth territory round trip wave imaging method.

Background technology

The pre-stack seismic imaging is the seismic data processing technique that proposes the seventies in last century, comprises time domain imaging technique and Depth Domain imaging technique; In the Depth Domain formation method; The method that adopts at present comprises that these two class methods are difficult to handle the violent high steep dip structure imaging problem of seismic event lateral speed change based on the integral method of ray theory with based on the differential wave equation one way ripple solution of wave theory.

The development in recent years higher round trip wave imaging method of a kind of precision that gets up, this method adopt round trip ripple to find the solution the differential wave equation, have remedied the image defects that conventional seismic migration faced theoretically.But conventional prestack round trip ripple imaging technique is that wave field extrapolation is carried out on the basis with surface flatly, promptly supposes geophone station on identical plane, and the seismologic record of in fact gathering is often at irregular relief surface.In area, western mountain region, the face of land discrepancy in elevation adopts conventional reference field imaging technique to cause incorrect structural feature because of surface relief up to thousands of rice in this case.The problem of seismic imaging method aspect complicated earth surface, the main existence two in complex structure area: the one, surface relief is big, and shallow velocity structure is complicated; The 2nd, underground structure is complicated; This pair of labyrinth makes the seismic imaging method difficulty that accurately forms images in this type area very big.

The seismic data imaging of relief surface area is main adopts two kinds of schemes: 1. is first top layer wave field is proofreaied and correct, and migration imaging again, this is conventional reference field seismic imaging technology; 2. being directly to carry out the Depth Domain imaging from relief surface, is the relief surface imaging technique that will study.In fact, the seismic imaging technology of relief surface is applied in the past Kirchhoff integral method and one way wave method.The Kirchhoff integral method is mainly considered the kinematics character of seismic wave propagation, adapts to any observed pattern, can handle the surface conditions of fluctuating neatly, and the propagation of seismic wave time ratio of gained is [Wiggins (1984)] accurately.But caustic, the interference of integral method in can not the dealing with complicated wave field, tyrannical during to velocity variations imaging precision not high, and to complex geological tectonic the blind area can appear.One way ripple wave equation pre stack depth migration method is high to the complex structure imaging precision; Can handle many ripples to problem; For the inspection surface area that is applied to one way ripple wave equation migration to rise and fall; The thinking that Reshef can superpose based on wave field has proposed the thought of " progressively adding up " wave field extrapolation, can carry out the migration imaging of wave field extrapolation from any non-horizontal surface; The one way wave method is a kind of approximation method equally, and the imaging angle is limited in 90 degree.

To the problems referred to above, will face the two large problems that will solve: the one, complex area surface relief problem; The 2nd, complicated wave field angle of image degree restricted problem.

In sum, depth imaging technique is obtained under the relief surface, high steep dip (greater than 90 degree) is constructed accurately image and is become the technical barrier that current urgent need solves.

Summary of the invention

The object of the invention is intended to overcome the wave field distortion that above-mentioned prior art causes more greatly because of the face of land discrepancy in elevation; The false defective of structural feature can be obtained the more high-precision seismic imaging and the true surface relief prestack depth territory round trip wave imaging method of structural feature more accurately and provide a kind of.

True surface relief prestack depth of the present invention territory round trip wave imaging method comprises the steps:

(1), the relief surface reference field is handled, obtain geological data CMP face height value, obtain shot point, geophone station true face of land height value and coordinate information, definition shot point locus; DATUM is a fixed datum, and CMP surface is CMP road collection and speed place face, and Turce Surface is a true face of land elevation face (like Fig. 1);

(2), the structure source wavelet adopts difference method realization source wavefield just to drill, and uses true face of land boundary condition;

$w\left(t\right)=\frac{1}{w_2-w_1}\left\{w_2\exp \right[-{\left(\mathrm{\π}{w}_{2}t\right)}^2]-w_1\exp [-{\left(\mathrm{\π}{w}_{1}t\right)}^2\left]\right\}---\left(1\right)$

Make up source wavelet according to earthquake recording frequency and formula (1); Utilize high-order limited difference method formula (2) just drilling source wavefield to the dominant record time; The big gun wave field just pushing away in the process whole writing time T be divided into plurality of sections; Every one step of continuation is used the relief surface free boundary condition, makes up filter function F _Il(x z) carries out filtering to the above wave field of relief surface; Store the two-layer wave field data of each section when initial simultaneously to the machine internal memory;

$P_{i,j,k}^{n+1}=2P_{i,j,k}^n-P_{i,j,k}^{n-1}+\frac{{\mathrm{<figure-callout\; id="2"\; label="c"\; filenames="HSA00000649366800032.png,HSA00000649366800041.png"\; state="\{\{state\}\}">c</figure-callout>}}^2\mathrm{\&Delta;}{\mathrm{<figure-callout\; id="2"\; label="t"\; filenames="HSA00000649366800032.png,HSA00000649366800041.png"\; state="\{\{state\}\}">t</figure-callout>}}^2}{\mathrm{\&Delta;}{x}^2}\underset{m=1}{\overset{M}{\mathrm{\&Sigma;}}}{a}_m(P_{i+m,j,k}^n+P_{i-m,j,k}^n-2P_{i,j,k}^n)$

$+\frac{{\mathrm{<figure-callout\; id="2"\; label="c"\; filenames="HSA00000649366800032.png,HSA00000649366800041.png"\; state="\{\{state\}\}">c</figure-callout>}}^2\mathrm{\&Delta;}{\mathrm{<figure-callout\; id="2"\; label="t"\; filenames="HSA00000649366800032.png,HSA00000649366800041.png"\; state="\{\{state\}\}">t</figure-callout>}}^2}{\mathrm{\&Delta;}{\mathrm{<figure-callout\; id="2"\; label="y"\; filenames="HSA00000649366800032.png,HSA00000649366800041.png"\; state="\{\{state\}\}">y</figure-callout>}}^2}\underset{m=1}{\overset{M}{\mathrm{\&Sigma;}}}{a}_m(P_{i,j+m,k}^n+P_{i,,j-m,k}^n-2P_{i,j,k}^n)$

$+\frac{{\mathrm{<figure-callout\; id="2"\; label="c"\; filenames="HSA00000649366800032.png,HSA00000649366800041.png"\; state="\{\{state\}\}">c</figure-callout>}}^2\mathrm{\&Delta;}{\mathrm{<figure-callout\; id="2"\; label="t"\; filenames="HSA00000649366800032.png,HSA00000649366800041.png"\; state="\{\{state\}\}">t</figure-callout>}}^2}{\mathrm{\&Delta;}{\mathrm{<figure-callout\; id="2"\; label="z"\; filenames="HSA00000649366800032.png,HSA00000649366800041.png"\; state="\{\{state\}\}">z</figure-callout>}}^2}\underset{m=1}{\overset{M}{\mathrm{\&Sigma;}}}{a}_m(P_{i,j,k+m}^n+P_{i,j,k-m}^n-2P_{i,j,k}^n)$

$+{\mathrm{<figure-callout\; id="2"\; label="c"\; filenames="HSA00000649366800032.png,HSA00000649366800041.png"\; state="\{\{state\}\}">c</figure-callout>}}^2\mathrm{\&Delta;}{\mathrm{<figure-callout\; id="2"\; label="t"\; filenames="HSA00000649366800032.png,HSA00000649366800041.png"\; state="\{\{state\}\}">t</figure-callout>}}^2{\mathrm{\&delta;}}_{i,j,k}---\left(2\right)$

δ in the formula _{I, j}Be the focus item, f (n) is discrete source function, i _s, j _s, k _sBe respectively the corresponding x of focal point, y, the net point sequence number of z direction;

(3), be recorded as the input target, utilize finite difference method along the extrapolation of time orientation inverse time to the record wave field, formula (4) is used the relief surface free boundary condition simultaneously, formula (5) with the big gun collection; Obtain the inverse time extrapolation record wave field of relief surface;

$P_{i,j,k}^{n-1}=2P_{i,j,k}^n-P_{i,j,k}^{n+1}+\frac{{\mathrm{<figure-callout\; id="2"\; label="c"\; filenames="HSA00000649366800032.png,HSA00000649366800041.png"\; state="\{\{state\}\}">c</figure-callout>}}^2\mathrm{\&Delta;}{\mathrm{<figure-callout\; id="2"\; label="t"\; filenames="HSA00000649366800032.png,HSA00000649366800041.png"\; state="\{\{state\}\}">t</figure-callout>}}^2}{\mathrm{\&Delta;}{x}^2}\underset{m=1}{\overset{M}{\mathrm{\&Sigma;}}}{a}_m(P_{i+m,j,k}^n+P_{i-m,j,k}^n-2P_{i,j,k}^n)$

$+\frac{{\mathrm{<figure-callout\; id="2"\; label="c"\; filenames="HSA00000649366800032.png,HSA00000649366800041.png"\; state="\{\{state\}\}">c</figure-callout>}}^2\mathrm{\&Delta;}{\mathrm{<figure-callout\; id="2"\; label="t"\; filenames="HSA00000649366800032.png,HSA00000649366800041.png"\; state="\{\{state\}\}">t</figure-callout>}}^2}{\mathrm{\&Delta;}{\mathrm{<figure-callout\; id="2"\; label="y"\; filenames="HSA00000649366800032.png,HSA00000649366800041.png"\; state="\{\{state\}\}">y</figure-callout>}}^2}\underset{m=1}{\overset{M}{\mathrm{\&Sigma;}}}{a}_m(P_{i,j+m,k}^n+P_{i,,j-m,k}^n-2P_{i,j,k}^n)$

$+\frac{{\mathrm{<figure-callout\; id="2"\; label="c"\; filenames="HSA00000649366800032.png,HSA00000649366800041.png"\; state="\{\{state\}\}">c</figure-callout>}}^2\mathrm{\&Delta;}{\mathrm{<figure-callout\; id="2"\; label="t"\; filenames="HSA00000649366800032.png,HSA00000649366800041.png"\; state="\{\{state\}\}">t</figure-callout>}}^2}{\mathrm{\&Delta;}{\mathrm{<figure-callout\; id="2"\; label="z"\; filenames="HSA00000649366800032.png,HSA00000649366800041.png"\; state="\{\{state\}\}">z</figure-callout>}}^2}\underset{m=1}{\overset{M}{\mathrm{\&Sigma;}}}{a}_m(P_{i,j,k+m}^n+P_{i,j,k-m}^n-2P_{i,j,k}^n)$

$+{\mathrm{<figure-callout\; id="2"\; label="c"\; filenames="HSA00000649366800032.png,HSA00000649366800041.png"\; state="\{\{state\}\}">c</figure-callout>}}^2\mathrm{\&Delta;}{\mathrm{<figure-callout\; id="2"\; label="t"\; filenames="HSA00000649366800032.png,HSA00000649366800041.png"\; state="\{\{state\}\}">t</figure-callout>}}^2{\mathrm{\&delta;}}_{i,j,k}---\left(4\right)$

Relief surface free boundary condition: P _Ω(x, d)=0 (5)

(4), be input with source wavefield and the record wave field of using the fluctuating boundary condition, use wave field dependent imaging condition, to obtain the imaging results in this moment.

Said is reference field with true ground surface, again definition of T o value and be specially according to true face of land elevation erection rate model: geological data is carried out conventional processing, Time Created territory, Depth Domain reference field aspect graph and define its two spatial relationship; Set up initial velocity model based on stationary plane, on the initial velocity model basis, initial velocity model is carried out iteration optimization, obtain velocity depth model accurately; Be the basis with the rough ground surface, definition rough ground surface in rate pattern, and between rough ground surface and reference field the filling replacement velocity, rate pattern is modified to the corresponding rate pattern of relief surface face.

The present invention also provides a kind of true surface relief prestack depth territory round trip wave imaging method, and treatment technology comprises following content:

Relief surface obtains, and is used for the reference field form of seismologic record is decomposed, and obtains the corresponding rough ground surface of seismologic record;

The rate pattern correction is used to set up the seismic velocity model, and definition time territory, Depth Domain reference field aspect graph time-space relationship, erection rate model obtain the corresponding rate pattern of relief surface face;

Source wavelet makes up with wave field and just drills, and is used to make up suitable source wavelet, on the source wavelet basis, carries out wave field and just drills, and obtains source wavefield;

Data ruleization is used for the seismologic record data processing, to obtain the basis of the regular geological data of said time domain as the extrapolation of record wave field inverse time;

The relief surface free boundary is handled, and is used to handle the relief surface boundary problem, makes up filter function F _Il(x z) carries out filtering to source wavefield and record wave field, and satisfying the relief surface boundary stress is zero boundary condition, to obtain the seismic wave field based on relief surface.

The present invention compared with prior art remarkable result is: the present invention carries out Depth Domain wave field extrapolation and imaging through record data on the true face of land; The wave field of having avoided causing more greatly because of the face of land discrepancy in elevation under the routine techniques distorts; The false defective of structural feature; Solve the deficiency of using round trip ripple imaging technique under the present surface relief condition, and can obtain accurate seismic imaging and more real Geological Structural Forms thereof under the complicated earth surface condition.

Description of drawings

Fig. 1 is the true surface relief prestack depth of a present invention territory round trip wave imaging method theory diagram

Fig. 2 is a relief surface reference field graph of a relation of the present invention

Fig. 3 a is a BLH zone time of the present invention territory floating datum aspect graph

Fig. 3 b is a BLH of the present invention area Depth Domain floating datum aspect graph

Fig. 4 is a BLH of the present invention area rate pattern diagrammatic cross-section

Fig. 5 a is that the present invention makes up source wavelet waveform synoptic diagram

Fig. 5 b is that the present invention makes up source wavelet wavelet spectrum synoptic diagram

Fig. 6 a is that theoretical model of the present invention is just being drilled the wave field synoptic diagram

Fig. 6 b is that theoretical model of the present invention is just being drilled single big gun synoptic diagram

Fig. 7 a is that the preceding synoptic diagram of making an uproar is pressed in relief surface round trip ripple imaging of the present invention

Fig. 7 b is that the back synoptic diagram of making an uproar is pressed in relief surface round trip ripple imaging of the present invention

Fig. 8 a is a relief surface one way ripple imaging synoptic diagram of the present invention

Fig. 8 b is a relief surface round trip ripple imaging synoptic diagram of the present invention

Fig. 9 a is a BLH of the present invention area one way ripple imaging results synoptic diagram

Fig. 9 b is a relief surface round trip ripple imaging synoptic diagram of the present invention

Embodiment

Below in conjunction with accompanying drawing and specific embodiment the present invention is further described

Embodiment 1

The seismologic record data that selection has complicated earth surface, complex structure and complicated wave field make an experiment; BLH area surface relief is bigger; Because the complicacy of the face of land and underground structure causes this area's seismic imaging problem serious; Particularly high steep dip can't form images greater than structure under stratum and the salt of 90 degree, receives the influence of surface relief and salt dome velocity sag, salt structure distortion down; Through repetition test repeatedly, its result adopts conventional stationary plane seismic imaging technology can't solve this area's accurately image problem.For this reason; We utilize a kind of true surface relief prestack depth of the present invention territory round trip wave imaging method; This method has been carried out true face of land seismic imaging and has been handled in the BLH basin, under relief surface, high steep dip is greater than having obtained very desirable imaging effect under the stratum of 90 degree and the salt.

True surface relief prestack depth territory round trip wave imaging method, concrete steps are following:

1), Fig. 1 provide for the embodiment of the invention 1 based on true face of land inverse time Depth Domain round trip wave imaging method process flow diagram; At first set up the floating surface form, confirm shot point position and imaging space, according to true face of land information erection rate model; Carry out that true face of land boundary condition source wavefield is just being drilled and the extrapolation of wave field inverse time, use the relevant output of wave field imaging results.

2), definition concerns between the rough ground surface, and is as shown in Figure 2; DATUM is a fixed datum, and CMP surface is CMP road collection and speed place face, and Turce Surface is a true face of land elevation face; Definition time territory, Depth Domain reference field aspect graph are shown in Fig. 3 a, Fig. 3 b; The relief surface reference field is handled, obtained geological data CMP face height value, obtain shot point, geophone station true face of land height value and coordinate information, definition shot point locus;

3), definition shot point position, according to true face of land information erection rate model, as shown in Figure 4;

4), make up source wavelet, carry out true face of land boundary condition source wavefield and just drill; On seismologic record, carry out spectrum analysis and confirm frequency range W ₁～W ₂, make up source wavelet according to earthquake recording frequency and formula (1) and carry out wave field and just drill, shown in Fig. 5 a, Fig. 5 b;

5), use the fluctuating boundary condition, wave field dependent imaging condition is obtained the seismic imaging result; Through from record data, obtaining true face of land elevation; According to relief surface definition shot point locus, and realize that true face of land free boundary source wavefield just drilling and noting down wave field along the extrapolation of time orientation inverse time, shown in Fig. 6 a, Fig. 6 b; Obtain Depth Domain seismic imaging data so that the earthquake field wave is relevant, shown in Fig. 7 a, Fig. 7 b, Fig. 8 a, Fig. 8 b;

6), the moving equation Depth Domain imaging technique of relief surface round trip wave-wave; Solve under the salt, the imaging problem of salt dome flank; Realized accurate playback to salt dome border and salt dome flank; Eliminated the salt dome velocity sag to the time domain structure distortion that underlying formation causes, the salt sub-surface can accurately be formed images in Depth Domain, as shown in Figure 9.This method solves the defective of using round trip ripple imaging technique under the present surface relief condition, to obtain more high-precision imaging and more real Geological Structural Forms.

The above embodiment is the part of instance of the present invention, and except that the foregoing description, the present invention can also have other embodiments; All employings are equal to the technical scheme of replacement or equivalent transformation form, all drop on the protection domain of requirement of the present invention.

Claims (2)

1. true surface relief prestack depth territory round trip wave imaging method, it is characterized in that: this method comprises the steps:

(1), the relief surface reference field is handled, obtain geological data CMP face height value, obtain shot point, geophone station true face of land height value and coordinate information, definition shot point locus; DATUM is a fixed datum, and CMP surface is CMP road collection and speed place face, and Turce Surface is a true face of land elevation face;

(2), the structure source wavelet adopts difference method realization source wavefield just to drill, and uses true face of land boundary condition; $w\left(t\right)=\frac{1}{w_2-w_1}\left\{w_2\mathrm{Exp}\right[-{\left(\mathrm{\&pi;}{w}_{2}t\right)}^2]-w_1\mathrm{Exp}[-{\left(\mathrm{\&pi;}{w}_{1}t\right)}^2\left]\right\}---\left(1\right)$

Make up source wavelet according to earthquake recording frequency and formula (1); Utilize high-order limited difference method formula (2) just drilling source wavefield to the dominant record time; The big gun wave field just pushing away in the process whole writing time T be divided into plurality of sections, every one step of continuation is used the relief surface free boundary condition.Make up filter function F _Il(x z) carries out filtering to the above wave field of relief surface; Store the two-layer wave field data of each section when initial simultaneously to the machine internal memory;

$P_{i,j,k}^{n+1}=2P_{i,j,k}^n-P_{i,j,k}^{n-1}+\frac{c^2\mathrm{\&Delta;}{t}^2}{\mathrm{\&Delta;}{x}^2}\underset{m=1}{\overset{M}{\mathrm{\&Sigma;}}}{a}_m(P_{i+m,j,k}^n+P_{i-m,j,k}^n-2P_{i,j,k}^n)$

$+\frac{c^2\mathrm{\&Delta;}{t}^2}{\mathrm{\&Delta;}{y}^2}\underset{m=1}{\overset{M}{\mathrm{\&Sigma;}}}{a}_m(P_{i,j+m,k}^n+P_{i,,j-m,k}^n-2P_{i,j,k}^n)$

$+\frac{c^2\mathrm{\&Delta;}{t}^2}{\mathrm{\&Delta;}{y}^2}\underset{m=1}{\overset{M}{\mathrm{\&Sigma;}}}{a}_m(P_{i,j+m,k}^n+P_{i,,j-m,k}^n-2P_{i,j,k}^n)$

$+c^2\mathrm{\&Delta;}{t}^2{\mathrm{\&delta;}}_{i,j,k}---\left(2\right)$

δ in the formula _{I, j}Be the focus item, f (n) is discrete source function, i _s, j _s, k _sBe respectively the corresponding x of focal point, y, the net point sequence number of z direction;

(3), be recorded as the input target, utilize finite difference method along the extrapolation of time orientation inverse time to the record wave field, formula (4) is used the relief surface free boundary condition simultaneously, formula (5) with the big gun collection; Obtain the inverse time extrapolation record wave field of relief surface;

$P_{i,j,k}^{n-1}=2P_{i,j,k}^n-P_{i,j,k}^{n+1}+\frac{c^2\mathrm{\&Delta;}{t}^2}{\mathrm{\&Delta;}{x}^2}\underset{m=1}{\overset{M}{\mathrm{\&Sigma;}}}{a}_m(P_{i+m,j,k}^n+P_{i-m,j,k}^n-2P_{i,j,k}^n)$

$+\frac{c^2\mathrm{\&Delta;}{t}^2}{\mathrm{\&Delta;}{y}^2}\underset{m=1}{\overset{M}{\mathrm{\&Sigma;}}}{a}_m(P_{i,j+m,k}^n+P_{i,,j-m,k}^n-2P_{i,j,k}^n)$

$+\frac{c^2\mathrm{\&Delta;}{t}^2}{\mathrm{\&Delta;}{z}^2}\underset{m=1}{\overset{M}{\mathrm{\&Sigma;}}}{a}_m(P_{i,j,k+m}^n+P_{i,j,k-m}^n-2P_{i,j,k}^n)$

$+c^2\mathrm{\&Delta;}{t}^2{\mathrm{\&delta;}}_{i,j,k}---\left(4\right)$

Relief surface free boundary condition: P _Ω(x, d)=0 (5)

(4), be input with source wavefield and the record wave field of using the fluctuating boundary condition, use wave field dependent imaging condition, to obtain the imaging results in this moment.

2. true surface relief prestack depth according to claim 1 territory round trip wave imaging method is characterized in that, true surface relief prestack depth territory round trip wave imaging method, and treatment technology comprises:

A, relief surface obtain: the reference field form to seismologic record is decomposed, and obtains the corresponding rough ground surface of seismologic record;

B, rate pattern correction: set up the seismic velocity model, definition time territory, Depth Domain reference field aspect graph time-space relationship, erection rate model obtain the corresponding rate pattern of relief surface face;

C, source wavelet make up with wave field and just drill: make up suitable source wavelet, on the source wavelet basis, carry out wave field and just drill, obtain source wavefield;

D, data ruleization: the seismologic record data processing, to obtain the basis of the regular geological data of said time domain as the extrapolation of record wave field inverse time;

E, relief surface free boundary are handled: handle the relief surface boundary problem, make up filter function F _Il(x z) carries out filtering to source wavefield and record wave field, and satisfying the relief surface boundary stress is zero boundary condition, to obtain the seismic wave field based on relief surface.

Images