CN103984007A - Optimization design method for time delay parameters of directional seismic waves - Google Patents

Abstract

The invention relates to an optimization design method for time delay parameters of directional seismic waves. The optimization design method comprises the following steps: establishing an exploration area underground medium model, calculating through the model to obtain directional seismic waves in different directions, obtaining an optimal time delay parameter according to a principle of combination of the strongest energy of reflected wave with the equilibrium distribution of the energy in a receiving array. By adopting the optimization design method, the optimization selection of the time delay parameters of the directional seismic waves under the complex geological conditions is effectively achieved. Compared with the conventional directional seismic wave forming method, the optimization design method is effectively applicable to resource exploration under the complex geological conditions in strong noise environment, thereby well improving the quality of the seismic exploration data.

Description

Seismic wave delay parameter Optimization Design directionally

Technical field:

The present invention relates to the Optimization Design of optimum delay parameter under a kind of directionally seismic wave delay parameter Optimization Design, particularly complex dielectrics condition distributing based on balancing energy.

Background technology:

Directional beam formation technology comes from phased-array radar, by controlling the phase place of each excitation apparatus, ripple is strengthened in a direction, other direction weakens, form the directive wave beam of tool, target echo is strengthened, thereby reach high-quality and long-range detection object, this technology has been applied to field of seismic exploration.

CN1560651 discloses a kind of < < phased array seismic prospecting method > >, and it adopts a plurality of vibroseiss to arrange by proportional spacing conllinear mode.By controlling time delay or the phase differential of each vibroseis, make sweep signal that each vibroseis sends in underground a direction with superimposed, vibration signal is strengthened.By changing time delay and the phase differential of each focus, thereby realize the orientation of seismic beam, make directionally seismic exploration become possibility.Adopt phased array seismic prospecting method, for steep dip tiltedly the seismic prospecting of plastid a kind of effective ways are provided, realize directionally seismic exploration, can also improve in addition the seismic prospecting degree of depth and seismic prospecting resolution, solved the problem that in combination seism, focus number of units can not be too much.

CN101320095 discloses a kind of < < single-seismic source directional lighting seismic signal synthesizing method > >, on survey line, design more than 1 and even n the equidistant shot point of odd number, on this survey line, be also laid with more than 1 and even the track pitch geophone station such as N simultaneously, by vibroseis, obtain the single big gun geological data of shot point 1, mobile shot point, obtain the single big gun geological data of shot point 2, mobile shot point, obtain the single big gun geological data of shot point 3, analogize in this way and obtain and even single big gun geological data of n shot point, single big gun geological data after time delay is done to linear superposition, the synthetic directive directional lighting seismic signal of tool, the corresponding directed seismic wave field of this signal.Single-seismic source directional lighting seismic signal synthesizing method is compared with phased seismic signal, and the signal to noise ratio (S/N ratio) of signal improves, and has the focus of needs number few, and job costs are low, simple to operate, has overcome the problem of a plurality of vibration exciter inconsistencies.

The directional lighting method of seismic prospecting > > of CN101339252 < < based on single focus.By traditional geologic prospecting method, determine offset distance, track pitch, conllinear is arranged wave detector array, offset distance according to same group of interior each shot point is different, sequence number does not differ m the method that shot point offset distance is identical of n on the same group, gather single big gun geological data of the entire profile, to belonging to single big gun geological data of same group, sequentially do successively by number uniformly-spaced time delay, single big gun geological data after time delay is done to linear superposition, the synthetic directive directional lighting seismic signal of tool, the corresponding directed seismic wave field of this signal, utilize common seismic data processing method just can obtain the directional lighting seismic section based on single focus.

Directionally seismic wave formation technology, by Earthquake occurrence control beam direction, improves and receive target echo signal to noise ratio (S/N ratio) in data, and the control of beam direction is to realize by the control to delay time.Obtain optimum target echo signal to noise ratio (S/N ratio), need to determine optimum delay parameter.Descend simply under ambient condition, the method that delay parameter can provide by above-mentioned patent is chosen, and actual seismic exploration is faced, is often complex dielectrics.

Summary of the invention:

Object of the present invention is just for above-mentioned the deficiencies in the prior art, and a kind of source and receiving end beam-forming that is mainly used in field of seismic exploration is provided, for the directionally seismic wave delay parameter Optimization Design of complex dielectrics condition.

Main thought of the present invention is: by model, calculate, try to achieve the directed geological data of different directions, and principle that balancing energy distribution combine the strongest according to receiving array reflected waves energy, obtains optimum delay parameter.

The present invention is achieved by the following technical solutions:

Seismic wave delay parameter Optimization Design directionally, comprises the following steps:

A, according to surveying district's geology, geophysical information, set up survey area underground medium model;

B, focal point and survey layout method of seismic prospecting design routinely, establishing focal point number is N, number order is respectively 1,2, ..., N, wherein N >=2, each focal point keeps equidistant d, wave detector road number is n, and track pitch is D, and source function is f (t), set up and solve equations for elastic waves, obtain N single shot record of a corresponding N focal point;

C, base area seismic wave kinematic principle, ask for the root-mean-square velocity v that covers medium on target reflection horizon;

D, establish seismic event main beam direction and normal angle is θ, from seismic wave directionally, form principle wherein delay time τ is the delay time parameter of directionally seismic wave formation method, then according to trigonometric function character-1≤sin θ≤1, has

$-\frac{d}{v}\&le;\mathrm{\&tau;}\&le;\frac{d}{v}$ (inequality 1)

It is leading that wherein delay parameter τ gets negative value interval scale signal,

Solve obtain the disaggregation t of inequality ₀≤ τ≤t ₁;

E, at [t ₀, t ₁] in, adopt linear interpolation uniformly-spaced to choose m delay time parameter, be designated as τ ₁, τ ₂..., τ _i..., τ _m, τ wherein ₁=t ₀, τ _m=t ₁, m>=20, for each delay parameter τ _i(i ∈ [1, m]), presses focal point numbering time delay 0,2 τ respectively by N single shot record _i, 3 τ _i..., (N-1) τ _i, then stack, obtains corresponding delay parameter τ _idirected seismologic record H _i;

F, with H _icentered by middle exploration targets reflected signal lineups, window while choosing isometric reflection wave, calculates directed seismologic record H _inei Ge road target echo ENERGY E _{i, 1}, E _{i, 2}..., E _i,j..., E _i,n, j ∈ [1, n] wherein, according to formula calculate receiving array Nei Ge road reflected energy density p _i,j;

G, basis calculate H _ithe object reflection wave gross energy h that interior receiving array receives _i; According to [(τ ₁, h ₁), (τ ₂, h ₂) ..., (τ _i, h _i) ..., (τ _m, h _m)] draw reflected energy with delay time parameter change curve, and therefrom find out front 5 maximum of points, establish delay parameter corresponding to these five maximum of points and be:

τ ₁',τ ₂',τ ₃',τ ₄',τ ₅'；

H, calculate corresponding τ ₁', τ ₂', τ ₃', τ ₄', τ ₅' the reflected energy density variance of 5 directed seismologic records, delay time corresponding to minimum variance is the directionally optimum delay time parameter of seismic wave.

Beneficial effect: through test, the directionally seismic wave delay parameter Optimization Design distributing based on balancing energy, effectively solved under complex geological condition, directionally the delay parameter of seismic wave is optimized On The Choice, than traditional directionally seismic wave formation method, effectively be applied to the resource exploration of complex geological condition under strong noise environment, improve better seismic exploration data quality.

Accompanying drawing explanation:

Fig. 1 dielectric model

Seismic beam after Fig. 2 delay parameter is optimized forms result (having removed direct wave) τ=1.875ms

The seismic beam that Fig. 3 does not carry out delay parameter optimization forms result (having removed direct wave) τ=3ms

Embodiment:

Below in conjunction with drawings and Examples, be described in further detail:

Somewhere, Liaoning inclination stratiform dielectric model of take is example, for first reflectance target layer, carried out the directionally seismic wave delay parameter optimal design based on balancing energy distribution, but delay parameter optimal design is not subject to the parameter limit providing in example.

The directionally seismic wave delay parameter optimal design distributing based on balancing energy, comprises the following steps:

A, according to surveying district's geology, geophysical information, set up survey area inclination stratiform dielectric model, as shown in Figure 1, the long 1220m of model area, dark 1220m, water intaking is flat is to the right the positive dirction of x axle, is z axle positive dirction straight down, the point of crossing of two directions is taken as true origin, sets up rectangular coordinate system;

B, focal point and survey layout method of seismic prospecting design routinely, focal point number is N=9, number order is respectively 1,2, ..., 9, focal point coordinate is followed successively by (192,20) (200,20), (208,20), (216,20), (224,20), (232,20) (240,20) (248,20) (252,20), focal point spacing d=8, wave detector is arranged to press track pitch D=4m conllinear within the scope of x positive dirction 240～1120m, and n=221 is counted in wave detector road, set up and solve and set up elastic wave one-order velocity-stress equation (a kind of form of equations for elastic waves)

$\left\{\begin{array}{c}\mathrm{\&rho;}\frac{{\&PartialD;v}_x}{\&PartialD;t}=\frac{{\&PartialD;\mathrm{\&tau;}}_{\mathrm{xx}}}{\&PartialD;x}+\frac{{\&PartialD;\mathrm{\&tau;}}_{\mathrm{xz}}}{\&PartialD;z}\\ \mathrm{\&rho;}\frac{\&PartialD;v_z}{\&PartialD;t}=\frac{{\&PartialD;\mathrm{\&tau;}}_{\mathrm{xz}}}{\&PartialD;x}+\frac{{\&PartialD;\mathrm{\&tau;}}_{\mathrm{zz}}}{\&PartialD;z}\\ \frac{{\&PartialD;\mathrm{\&tau;}}_{\mathrm{xx}}}{\&PartialD;t}=(\mathrm{\&lambda;}+2\mathrm{\&mu;})\frac{\&PartialD;v_x}{\&PartialD;x}+\mathrm{\&lambda;}\frac{\&PartialD;v_z}{\&PartialD;z}\\ \frac{{\&PartialD;\mathrm{\&tau;}}_{\mathrm{zz}}}{\&PartialD;t}=(\mathrm{\&lambda;}+2\mathrm{\&mu;})\frac{\&PartialD;v_z}{\&PartialD;z}+\mathrm{\&lambda;}\frac{{\&PartialD;v}_x}{\&PartialD;x}\\ \frac{{\&PartialD;\mathrm{\&tau;}}_{\mathrm{xz}}}{\&PartialD;t}=\mathrm{\&mu;}(\frac{\&PartialD;v_x}{\&PartialD;z}+\frac{{\&PartialD;v}_z}{\&PartialD;x})\end{array}\right.$

Wherein, v _x, v _zthe horizontal component and the vertical component that represent respectively speed, τ _xx, τ _zzrepresent respectively x, z direction normal stress, τ _xzfor shearing stress, λ, μ are Lame's constant, adopt staggered-mesh spatial domain quadravalence, time domain Using Second-Order Central Difference form,

$\left\{\begin{array}{c}{v}_{{\mathrm{xi},}^{\&prime;}{j}^{\&prime;}}^{l^{\&prime;}+1/2}=v_{{\mathrm{xi}}^{\&prime;},j^{\&prime;}}^{l^{\&prime;}-1/2}+\mathrm{\&rho;}(D_x{\mathrm{\&tau;}}_{\mathrm{xx}}+D_z{\mathrm{\&tau;}}_{\mathrm{xz}}){|}_{i^{\&prime;},j^{\&prime;}}^{l^{\&prime;}}\\ v_{{\mathrm{zi}}^{\&prime;}+1/2,j^{\&prime;}+1/2}^{l^{\&prime;}+1/2}=v_{{\mathrm{zi}}^{\&prime;}+1/2,j^{\&prime;}+1/2}^{l^{\&prime;}-1/2}+\mathrm{\&Delta;t\&rho;}(D_x{\mathrm{\&tau;}}_{\mathrm{xz}}+D_z{\mathrm{\&tau;}}_{\mathrm{zz}}){|}_{i^{\&prime;}+1/2,j^{\&prime;}+1/2}^{l^{\&prime;}}\\ {\mathrm{\&tau;}}_{{\mathrm{xxi}}^{\&prime;}+1/2,j^{\&prime;}}^{l^{\&prime;}+1}={\mathrm{\&tau;}}_{{\mathrm{xxi}}^{\&prime;}+1/2,j^{\&prime;}}^{l^{\&prime;}}+\mathrm{\&Delta;t}\left[(\mathrm{\&lambda;}+2\mathrm{\&mu;})\right]D_x{v}_x+\mathrm{\&lambda;}{D}_z{v}_z){|}_{i^{\&prime;}+1/2,j^{\&prime;}}^{l^{\&prime;}+1/2}\\ {\mathrm{\&tau;}}_{{\mathrm{zzi}}^{\&prime;}+1/2,j^{\&prime;}}^{l^{\&prime;}+1}={\mathrm{\&tau;}}_{{\mathrm{zzi}}^{\&prime;}+1/2,j^{\&prime;}}^{l^{\&prime;}}+\mathrm{\&Delta;t}\left[(\mathrm{\&lambda;}+2\mathrm{\&mu;})\right]D_z{v}_z+\mathrm{\&lambda;}{D}_x{v}_x){|}_{i^{\&prime;}+1/2,j^{\&prime;}}^{l^{\&prime;}+1/2}\\ {\mathrm{\&tau;}}_{{\mathrm{xzi}}^{\&prime;},j^{\&prime;}+1/2}^{l^{\&prime;}}={\mathrm{\&tau;}}_{{\mathrm{xzi}}^{\&prime;},j^{\&prime;}+1/2}^{l^{\&prime;}}+\mathrm{\&Delta;t}\left[\mathrm{\&mu;}(D_z{v}_x+D_x{v}_x)\right]{|}_{i^{\&prime;},j^{\&prime;}+1/2}^{l^{\&prime;}+1/2}\end{array}\right.$

Wherein, i', j', l' is respectively x, z, t direction loop variable, D _x, D _zrepresent respectively x, the first order differential operator of z, Δ t represents time step 0.125ms, with variable, g represents speed v _x, v _zor stress tensor τ _xx, τ _zz, τ _xz,

$\left\{\begin{array}{c}{D}_{x}g(x,z)=\frac{1}{\mathrm{\&Delta;x}}\underset{\mathrm{ii}=1}{\overset{4}{\mathrm{\&Sigma;}}}{C}_{\mathrm{ii}}^4[g(x+{\mathrm{\&Delta;x}}_{2\mathrm{ii}-1},z)-g(x-{\mathrm{\&Delta;x}}_{2\mathrm{ii}-1},z)]+o\left({\mathrm{\&Delta;x}}^8\right)\\ D_{z}g(x,z)=\frac{1}{\mathrm{\&Delta;z}}=\frac{1}{\mathrm{\&Delta;z}}\underset{\mathrm{jj}=1}{\overset{4}{\mathrm{\&Sigma;}}}{C}_{\mathrm{jj}}^4[g(x,z+{\mathrm{\&Delta;z}}_{2\mathrm{jj}-1})-g(x,z-{\mathrm{\&Delta;z}}_{2\mathrm{jj}-1})]+o\left({\mathrm{\&Delta;z}}^8\right)\end{array}\right.,$

Stable condition is v _maxfor maximum velocity of longitudinal wave, boundary condition adopts PML perfect matching layer to absorb and processes, and matching layer equation is accordingly,

$\left\{\begin{array}{c}{v}_x=v_x^x+v_x^z\\ \frac{{\&PartialD;v}_x^x}{\&PartialD;t}+d\left(x\right)v_x^x=\frac{1}{\mathrm{\&rho;}}\frac{{\&PartialD;\mathrm{\&tau;}}_{\mathrm{xx}}}{\&PartialD;x}\\ \frac{{\&PartialD;v}_x^z}{\&PartialD;t}+d\left(z\right)v_x^z=\frac{1}{\mathrm{\&rho;}}\frac{{\&PartialD;\mathrm{\&tau;}}_{\mathrm{xz}}}{\&PartialD;z}\end{array}\right.$ $\left\{\begin{array}{c}{v}_z=v_z^x+v_z^z\\ \frac{{\&PartialD;v}_z^x}{\&PartialD;t}+d\left(x\right)v_z^x=\frac{1}{\mathrm{\&rho;}}\frac{{\&PartialD;\mathrm{\&tau;}}_{\mathrm{xz}}}{\&PartialD;x}\\ \frac{{\&PartialD;v}_z^z}{\&PartialD;t}+d\left(z\right)v_z^z=\frac{1}{\mathrm{\&rho;}}\frac{{\&PartialD;\mathrm{\&tau;}}_{\mathrm{zz}}}{\&PartialD;z}\end{array}\right.$ $\left\{\begin{array}{c}{\mathrm{\&tau;}}_{\mathrm{xx}}={\mathrm{\&tau;}}_{\mathrm{xx}}^x+{\mathrm{\&tau;}}_{\mathrm{xx}}^z\\ \frac{{\&PartialD;\mathrm{\&tau;}}_{\mathrm{xx}}^x}{\&PartialD;t}+d\left(x\right){\mathrm{\&tau;}}_{\mathrm{xx}}^x=(\mathrm{\&lambda;}+2\mathrm{\&mu;})\frac{{\&PartialD;v}_x}{\&PartialD;x}\\ \frac{{\&PartialD;\mathrm{\&tau;}}_{\mathrm{xx}}^z}{\&PartialD;t}+d\left(z\right){\mathrm{\&tau;}}_{\mathrm{xx}}^z=\mathrm{\&lambda;}\frac{{\&PartialD;v}_z}{\&PartialD;z}\end{array}\right.$

$\left\{\begin{array}{c}{\mathrm{\&tau;}}_{\mathrm{zz}}={\mathrm{\&tau;}}_{\mathrm{zz}}^x+{\mathrm{\&tau;}}_{\mathrm{zz}}^z\\ \frac{{\&PartialD;\mathrm{\&tau;}}_{\mathrm{zz}}^x}{\&PartialD;t}+d\left(x\right){\mathrm{\&tau;}}_{\mathrm{zz}}^x=\mathrm{\&lambda;}\frac{{\&PartialD;v}_x}{\&PartialD;x}\\ \frac{{\&PartialD;\mathrm{\&tau;}}_{\mathrm{zz}}^z}{\&PartialD;t}+d\left(z\right){\mathrm{\&tau;}}_{\mathrm{zz}}^z=(\mathrm{\&lambda;}+2\mathrm{\&mu;})\frac{{\&PartialD;v}_z}{\&PartialD;z}\end{array}\right.$ $\left\{\begin{array}{c}{\mathrm{\&tau;}}_{\mathrm{xz}}={\mathrm{\&tau;}}_{\mathrm{xz}}^x+{\mathrm{\&tau;}}_{\mathrm{xz}}^z\\ \frac{{\&PartialD;\mathrm{\&tau;}}_{\mathrm{xz}}^x}{\&PartialD;t}+d\left(x\right){\mathrm{\&tau;}}_{\mathrm{xz}}^x=\mathrm{\&mu;}\frac{{\&PartialD;v}_z}{\&PartialD;x}\\ \frac{{\&PartialD;\mathrm{\&tau;}}_{\mathrm{xz}}^z}{\&PartialD;t}+d\left(z\right){\mathrm{\&tau;}}_{\mathrm{xz}}^z=\mathrm{\&mu;}\frac{{\&PartialD;v}_x}{\&PartialD;z}\end{array}\right.$

Wherein, d (x), d (z) represent respectively the attenuation coefficient of x direction and z direction, adopt Collino attenuation model,

$\left\{\begin{array}{c}d\left(x\right)=\frac{{3V}_{\max }}{2\mathrm{\&delta;}}\log \left(\frac{1}{R}\right){\left(\frac{x}{\mathrm{\&delta;}}\right)}^2\\ d\left(z\right)=\frac{{3V}_{\max }}{2\mathrm{\&delta;}}\log \left(\frac{1}{R}\right){\left(\frac{z}{\mathrm{\&delta;}}\right)}^2\end{array}\right.$

δ=200m is matching layer thickness, and R is theoretical boundary reflection coefficient, gets 10 herein ^-6, source function f (t) adopts the autocorrelation wavelet of linear raising frequency Chirp signal, and the formula of Chirp signal is sweep signal initial frequency f ₁=10Hz, termination frequency f ₂=100Hz, sweep signal amplitude A=1, sweep signal duration T=2s, signal while taking from related function central point left and right 45ms in window is as source signal, get spatial sampling step-length 4m, calculate the seismic wave field that each focal point forms, and obtain N=9 single shot record of a corresponding N=9 focal point;

C, base area seismic wave kinematic principle, asking for the root-mean-square velocity of covering medium on target reflection horizon is v=2000m/s;

D, establish seismic event main beam direction and normal angle is θ, from seismic wave directionally, form principle wherein delay time τ is the delay time parameter of directionally seismic wave formation method, then according to trigonometric function character-1≤sin θ≤1, has

$-\frac{d}{v}\&le;\mathrm{\&tau;}\&le;\frac{d}{v}$ (inequality 1)

It is leading that wherein delay parameter τ gets negative value interval scale signal, solve inequality 1, obtain the disaggregation-4≤τ≤4ms of inequality, the model providing in the present embodiment, focus is in the left side of receiving array, belong to monolateral reception, by seismic reflection principle, known main beam direction necessarily meets 0≤θ≤90, delay parameter should be on the occasion of, therefore the scope that can further dwindle τ is 0 to 4ms, i.e. t ₀=0ms, t ₁=4ms;

E, at [t ₀, t ₁] in, adopt linear interpolation uniformly-spaced to choose m=33 delay time parameter, be designated as τ ₁, τ ₂..., τ _i..., τ ₃₃, τ wherein ₁=0ms, τ ₃₃=4ms, for each delay parameter τ _i(i ∈ [1,33]), presses focal point numbering time delay 0,2 τ respectively by N single shot record _i, 3 τ _i..., 32 τ _irear stack, obtains corresponding delay parameter τ _idirected seismologic record H _i;

F, with H _icentered by middle exploration targets reflected signal lineups, window while choosing the reflection wave of lineups center each 50ms of left and right, calculates directed seismologic record H _inei Ge road target echo energy:

E _{i, 1}, E _{i, 2}..., E _i,j..., E _{i, 221}, j ∈ [1,221] wherein,

According to formula calculate receiving array Nei Ge road reflected energy density p _i,j;

G, basis can calculate H _ithe object reflection wave gross energy h that interior receiving array receives _i; According to [(τ ₁, h ₁), (τ ₂, h ₂) ..., (τ _i, h _i) ..., (τ ₃₂, h ₃₂)] can draw reflected energy with delay time parameter change curve, and therefrom find out front 5 maximum of points, establishing delay parameter corresponding to these five maximum of points is τ ₁'=1.625ms, τ ₂'=1.750ms, τ ₃'=1.875ms, τ ₄'=2.000ms, τ ₅'=2.125ms;

H, calculate corresponding τ ₁'=1.625ms, τ ₂'=1.750ms, τ ₃'=1.875ms, τ ₄'=2.000ms, τ ₅the reflected energy density variance of 5 directed seismologic records of '=2.125ms is respectively 7.14 * 10 ^-4, 5.13 * 10 ^-4, 4.03 * 10 ^-4, 5.74 * 10 ^-4, 6.69 * 10 ^-4, minimum variance is 4.03 * 10 ^-4,, corresponding delay time 1.875ms, τ=1.875ms is the optimum delay time parameter of directionally seismic wave.

Claims (1)

1. a seismic wave delay parameter Optimization Design directionally, is characterized in that, comprises the following steps:

A, according to surveying district's geology, geophysical information, set up exploratory area underground medium model;

B, focal point and survey layout method of seismic prospecting design routinely, establishing focal point number is N, number order is respectively 1,2, ..., N, wherein N >=2, each focal point keeps equidistant d, wave detector road number is n, and track pitch is D, and source function is f (t), set up and solve equations for elastic waves, obtain N single shot record of a corresponding N focal point;

C, ask for the root-mean-square velocity v that covers medium on target reflection horizon;

D, establish seismic event main beam direction and normal angle is θ, wherein delay time τ is the delay time parameter of directionally seismic wave formation method, then according to trigonometric function character-1≤sin θ≤1, has

$-\frac{d}{v}\&le;\mathrm{\&tau;}\&le;\frac{d}{v}$ (inequality 1)

It is leading that wherein delay parameter τ gets negative value interval scale signal, solves inequality 1, obtains the disaggregation of inequality

t ₀≤τ≤t ₁；

E, at [t ₀, t ₁] in, adopt linear interpolation uniformly-spaced to choose m delay time parameter, be designated as τ ₁, τ ₂..., τ _i..., τ _m, τ wherein ₁=t ₀, τ _m=t ₁, m>=20, for each delay parameter τ _i(i ∈ [1, m]), presses focal point numbering time delay 0,2 τ respectively by N single shot record _i, 3 τ _i..., (N-1) τ _i, then stack, obtains corresponding delay parameter τ _idirected seismologic record H _i;

F, with H _icentered by middle exploration targets reflected signal lineups, window while choosing isometric reflection wave, calculates directed seismologic record H _inei Ge road target echo ENERGY E _{i, 1}, E _{i, 2}..., E _i,j..., E _i,n, j ∈ [1, n] wherein, according to formula calculate receiving array Nei Ge road reflected energy density p _i,j;

G, basis calculate the object reflection wave gross energy h that in Hi, receiving array receives _i; According to [(τ ₁, h ₁), (τ ₂, h ₂) ..., (τ _i, h _i) ..., (τ _m, h _m)] draw reflected energy with delay time parameter change curve, and therefrom find out front 5 maximum of points, establish delay parameter corresponding to these five maximum of points and be:

τ ₁',τ ₂',τ ₃',τ ₄',τ ₅'；

H, calculate corresponding τ ₁', τ ₂', τ ₃', τ ₄', τ ₅' the reflected energy density variance of 5 directed seismologic records, delay time corresponding to minimum variance is the directionally optimum delay time parameter of seismic wave.