CN103576192A - Method for determining seismic wavelets - Google Patents

Abstract

The invention relates to a method for determining seismic wavelets during petroleum exploration seismic data processing. A bound term for determining seismic wavelet energy to be minimum is additionally arranged in an original objective function, influences of noise in a seismic channel on seismic wavelet determination accuracy are considered, the best balance point of seismic channel residual error energy and seismic wavelet energy is searched to determine the best cut-off order of a reflection coefficient convolution matrix singular value, and the reflection coefficient convolution matrix singular value is cut off to weaken influences of the noise in the seismic channel to determine the seismic wavelets. The method has no influence on contributions of a large singular value on wavelet determination, and compared with a common method, the method has the advantages of being favorable for weakening the influences of the noise in the seismic channel on the seismic wavelet determination accuracy and favorable for improving the seismic wavelet determination accuracy.

Description

A kind of method of definite seismic wavelet

Technical field

The present invention relates to petroleum exploration and development technology, the method for a kind of definite seismic wavelet in specifically seismic data processing, explanation.

Background technology

The definite of seismic wavelet is a very crucial problem in seismic prospecting.In forward problem, need in conjunction with seismic wavelet, form forward simulation geological data by wave equation or convolution model.In inverting and deconvolution problem, also need to extract a seismic wavelet by seismic trace, different seismic wavelets obtains different inversion results.

The basis of seismic wavelet extraction is convolution model, and it comprises seismic wavelet, reflection coefficient, contains noisy seismic channel data.Methods of seismic wavelet extraction is to utilize well-log information to calculate reflection coefficient, then in conjunction with seismic trace, by convolution model, obtains seismic wavelet.At present conventional well shake associating wavelet determines that method is based on synthetic road and seismic trace residual energy (hereinafter to be referred as seismic trace residual energy) minimum principle:

$\underset{t=t_1}{\overset{t=t_2}{\mathrm{\Σ}}}{[x\left(t\right)-\underset{\mathrm{\τ}=-p}{\overset{\mathrm{\τ}=q}{\mathrm{\Σ}}}w\left(\mathrm{\τ}\right)r(t-\mathrm{\τ})]}^2\→\min ---\left(1\right)$

In formula (1), r (t) is the time domain reflection coefficient being calculated by well-log information, and window scope is at that time: t ₁≤ t≤t ₂, x (t) is seismic trace amplitude, window scope is at t at that time ₁+ p≤t≤t ₂-q, w (t) be will be definite seismic wavelet, its beginning and end time scope be-p is to q, formula (1) is converted into matrix operation form and can obtains following expression:

E＝||RW-X|| ₂→min (2)

By $\frac{\∂E}{\∂W}=0$ :

W＝R ^-1X （3）

Formula (3) is at present conventional well shake and combines the computing formula of determining seismic wavelet, and R is reflection coefficient convolution matrix.Setting seismic wavelet sampling point number is that n(n is odd number), reflection coefficient sampling point number is m, reflection coefficient convolution matrix R is (m-n+1) * n rank matrixes, specifically can be expressed as:

$R=\left|\begin{array}{ccccc}{r}_{n-1}& r_{n-2}& ...& r_1& r_0\\ r_n& r_{n-1}& ...& r_2& {\mathrm{<figure-callout\; id="1"\; label="r"\; filenames="HDA00001964806600021.png"\; state="\{\{state\}\}">r</figure-callout>}}_1\\ .& .& .& .& .\\ .& .& .& .& .\\ .& .& .& .& .\\ r_{m-1}& r_{m-2}& ...& r_{n-1}& r_{n-2}\\ r_m& r_{m-1}& ...& r_n& r_{n-1}\end{array}\right|---\left(4\right)$

R in formula (4) _ifor reflectance value corresponding to each reflection sampling point;

X is (m-n+1) rank column vector that seismic trace amplitude forms, and specifically can be expressed as:

$X=\left|\begin{array}{c}{x}_0\\ x_1\\ .\\ .\\ .\\ x_{m-n}\\ x_{m-n+1}\end{array}\right|---\left(5\right)$

X in formula (5) _iamplitude for each sampling point of seismic trace;

Matrix R is carried out to svd:

$R=\mathrm{U\&Sigma;}{V}^T=\underset{i=0}{\overset{n}{\mathrm{\&Sigma;}}}{u}_i{\mathrm{\&delta;}}_i{v}_i^T---\left(6\right)$

∑=diag (δ in formula (6) ₀, δ ₁, δ _n-1) be n * n diagonal matrix, δ _ifor the singular value of reflection coefficient convolution matrix R, and δ ₀>=δ ₁>=...>=δ _n-1to successively decrease, order is arranged; U is reflection coefficient convolution matrix left singular matrix; V ^tfor reflection coefficient convolution matrix right singular matrix transposition; u _icolumn vector for reflection coefficient convolution matrix left singular matrix;

column vector for reflection coefficient convolution matrix right singular matrix transposition.

In conjunction with (6), (3), obtain seismic wavelet and determine formula:

$W_i=\underset{i=0}{\overset{n}{\mathrm{\&Sigma;}}}\frac{u_i^T{x}_i}{{\mathrm{\&delta;}}_i}{v}_i---\left(7\right)$

In formula (7)

column vector transposition for reflection coefficient convolution matrix left singular matrix; v _icolumn vector for reflection coefficient convolution matrix right singular matrix; All the other each parameter physical meaning are with aforementioned.

From formula (7), can find out: as singular value δ _ito become very little,

to become very large, small noise in seismic trace at this moment, will be on determining that the result of seismic wavelet produces larger impact, and result can make definite seismic wavelet and actual wavelet have very large deviation.

Summary of the invention

Be to provide a kind of method that is conducive to definite seismic wavelet of seismic wavelet precision raising order of the present invention.

Specific implementation step of the present invention is:

1) gather log data and seismic channel data and priori wavelet;

Log data described in step 1) mainly comprises: the terrestrial coordinate of well-name, this well position, the time domain reflection coefficient after composite traces accurate calibration.

Seismic channel data described in step 1) is near the multichannel seismic data together or well point of earthquake-capturing.

Seismic channel data described in step 1), refers to the geological data in window in the time of will determining seismic wavelet target interval; The time window top and bottom of seismic channel data will depart from crest and the trough of seismic trace, selects weak reflection position.

Priori wavelet described in step 1) is previously for the synthesis of record, to demarcate known seismic wavelet or theoretical wavelet.

2) utilize step 1) priori seismic wavelet and time domain reflection coefficient to make composite traces, within the scope of timing window, time domain reflection coefficient is being demarcated again;

Step 2) demarcation again in is by synthetic road and each well lie data correlation analysis are obtained to reflection coefficient and each well lie best match position.

Step 2) the timing window scope of giving in refers to when given earthquake within the scope of window traveling time territory reflection coefficient within the scope of positive and negative 30ms.

3) window scope when window and reflection coefficient while determining seismic channel data according to given seismic wavelet length;

Seismic wavelet length described in step 3) is definite according to actual seismic data, and shallow earthquake wavelet length is short, and deep seismic wavelet length is long; When the time domain reflection coefficient obtaining when step 1) can not meet step 3), the needs of window scope, reduce wavelet length, and the minimum length of wavelet length must comprise the integrality of main lobe and the secondary lobe of wavelet.

Seismic wavelet length described in step 3) is between 60--120ms.

During seismic trace described in step 3), window is greater than 1～1.5 times of seismic wavelet length, and when during reflection coefficient, window is than seismic trace, window scope is respectively expanded seismic wavelet length up and down half.

4) build reflection coefficient convolution matrix;

5) build seismic trace column vector X;

6) according to following formula, build and determine seismic wavelet objective function;

$E=\{{\left|\right|R{W}^{\left(k\right)}-X\left|\right|}_2^2+{\left|\right|W^{\left(k\right)}\left|\right|}_2^2\}\&RightArrow;\min ---\left(9\right)$

In formula:

R is reflection coefficient convolution matrix;

W ^(k)expression is carried out k rank to reflection coefficient convolution matrix and is blocked seismic wavelet in situation;

X is the column vector of seismic trace;

K is the exponent number that blocks of reflection coefficient convolution singular values of a matrix;

7) according to following formula, calculate different rank k singular value and block the seismic wavelet W in situation ^(k):

$W^{\left(k\right)}=\underset{i=0}{\overset{k}{\mathrm{\&Sigma;}}}\frac{u_i^T{x}_i}{{\mathrm{\&delta;}}_i}{v}_i---\left(10\right)$

In formula:

δ _isingular value for reflection coefficient convolution matrix R;

column vector transposition for reflection coefficient convolution matrix left singular matrix;

V _ifor reflection coefficient convolution matrix right singular matrix column vector;

X _ieach element for the column vector of seismic trace;

In step 7), the span of different rank k is 0 < k < n, and n is seismic wavelet sampling point number.

Seismic wavelet W in step 7) ^(k)formula is to obtain according to seismic wavelet objective function formula.

8) calculate energy and the seismic wavelet energy of seismic trace residual error, determine best k value:

Seismic trace residual energy $E_s^{\left(k\right)}={\left|\right|{\mathrm{RW}}^{\left(k\right)}-X\left|\right|}_2^2---\left(11\right),$

With seismic wavelet energy $E_w^{\left(k\right)}={\left|\right|W^{\left(k\right)}\left|\right|}_2^2;---\left(12\right)$

In formula:

for cutting the disconnected situation Xia Zhen of singular value road residual energy in k rank;

In formula:

for k rank, cut singular value and block situation seismic wavelet energy:

Seismic wavelet energy corresponding to different wavelets and seismic trace residual energy are done to X plot, choose seismic wavelet energy and seismic trace residual energy X plot flex point as best k value between the two;

9) utilize best k value and seismic wavelet W ^(k)formula determines that best seismic wavelet is definite;

10) to a plurality of seismic trace repeating steps 2)～9) obtain its corresponding best seismic wavelet, and then adopt multiple wavelet averaging method to determine a comprehensive seismic wavelet, be near seismic wavelet position, well point described in step 1).

The present invention and traditional well shake associating wavelet determine that method compares, and have increased the bound term of definite seismic wavelet energy minimum in original objective function, have formed new wavelet and have determined objective function.Compare with former objective function, considered that the noise in seismic trace, on determining the impact of seismic wavelet precision, is conducive to the raising that wavelet is determined precision.In specific implementation process, by X plot, obtain the best " equilibrium point " of seismic trace residual energy and seismic wavelet energy, determine the best of reflection coefficient convolution singular values of a matrix and block exponent number, by the method that reflection coefficient convolution singular values of a matrix is blocked, weaken the noise of seismic trace to determining the impact of seismic wavelet precision.

The present invention is by directly the less singular value of reflection coefficient convolution matrix is blocked to weaken seismic trace noise determines the impact of precision on seismic wavelet, do not have influence on the contribution definite to wavelet of larger singular value, add whitening factor method and compare with conventional, the present invention determines precision by more favourable raising wavelet.

Accompanying drawing explanation

Fig. 1 is through the calibrated time domain reflection coefficient of composite traces;

Fig. 2 is for determining the seismic channel data of seismic wavelet;

Fig. 3 is the seismic wavelet of previously demarcating for the synthesis of record;

Fig. 4 and Fig. 5 are previously for the synthesis of the spectral amplitude and the phase spectrum that record the seismic wavelet of demarcating;

Fig. 6 utilizes 105 road geological data to carry out wavelet while determining, the energy X plot (transverse axis represents seismic trace residual energy, and the longitudinal axis represents seismic wavelet energy) of synthetic road and well lie residual energy and seismic wavelet;

Fig. 7 is the seismic wavelet that utilizes the present invention definite;

Fig. 8 and Fig. 9 are respectively spectral amplitude and the phase spectrums of determined seismic wavelet;

Figure 10 finally determines that seismic wavelet makes synthetic road and seismic trace comparison diagram.

Specific embodiments

The present invention is on the minimum basis of original objective function seismic trace residual energy, has increased the bound term of a definite seismic wavelet energy minimum, has formed new seismic wavelet and has determined objective function.Compare with former objective function, considered that seismic trace noise determines the impact of precision on seismic wavelet, be conducive to the raising that seismic wavelet is determined precision.

Specific implementation process of the present invention is as follows:

1) gather log data and seismic channel data and priori wavelet;

Log data described in step 1) mainly comprises: the terrestrial coordinate of well-name, this well position, the time domain reflection coefficient after composite traces accurate calibration.

Seismic channel data described in step 1) is near the multichannel seismic data together or well point of earthquake-capturing.

Seismic channel data described in step 1), refers to the geological data in window in the time of will determining seismic wavelet target interval; The time window top and bottom of seismic channel data will depart from crest and the trough of seismic trace, selects weak reflection position.

Priori wavelet described in step 1) is previously for the synthesis of record, to demarcate known seismic wavelet or theoretical wavelet.

2) utilize step 1) priori seismic wavelet and time domain reflection coefficient to make composite traces, within the scope of timing window, time domain reflection coefficient is being demarcated again;

Step 2) demarcation again in is by synthetic road and each well lie data correlation analysis are obtained to reflection coefficient and each well lie best match position.

Step 2) the timing window scope of giving in refers to when given earthquake within the scope of window traveling time territory reflection coefficient within the scope of positive and negative 30ms.

3) window scope when window and reflection coefficient while determining seismic channel data according to given seismic wavelet length;

Seismic wavelet length described in step 3) is definite according to actual seismic data, and shallow earthquake wavelet length is short, and deep seismic wavelet length is long; When the time domain reflection coefficient obtaining when step 1) can not meet step 3), the needs of window scope, reduce wavelet length, and the minimum length of wavelet length must comprise the integrality of main lobe and the secondary lobe of wavelet.

Seismic wavelet length described in step 3) is between 60--120ms.

During seismic trace described in step 3), window is greater than 1～1.5 times of seismic wavelet length, and when during reflection coefficient, window is than seismic trace, window scope is respectively expanded seismic wavelet length up and down half.

4) build reflection coefficient convolution matrix;

The method that in step 4), reflection coefficient convolution matrix provides according to formula (4) builds, r in formula (4) _ifor reflectance value corresponding to each reflection sampling point, m, n are respectively the sampling point number of reflection coefficient and seismic wavelet;

5) build seismic trace column vector X;

The method that step 5) seismic trace column vector provides according to formula formula (5) builds, x in formula (5) _iamplitude for each sampling point of seismic trace; M, n(n are odd number) be respectively the sampling point number of reflection coefficient and seismic wavelet;

6) according to following formula, build and determine seismic wavelet objective function;

$E=\{{\left|\right|R{W}^{\left(k\right)}-X\left|\right|}_2^2+{\left|\right|W^{\left(k\right)}\left|\right|}_2^2\}\&RightArrow;\min ---\left(9\right)$

In formula:

R is reflection coefficient convolution matrix;

W ^(k)expression is carried out k rank to reflection coefficient convolution matrix and is blocked seismic wavelet in situation;

X is the column vector of seismic trace;

K is the exponent number that blocks of reflection coefficient convolution singular values of a matrix;

7) according to following formula, calculate different rank k singular value and block the seismic wavelet W in situation ^(k):

$W^{\left(k\right)}=\underset{i=0}{\overset{k}{\mathrm{\&Sigma;}}}\frac{u_i^T{x}_i}{{\mathrm{\&delta;}}_i}{v}_i---\left(10\right)$

In formula:

δ _isingular value for reflection coefficient convolution matrix R;

column vector transposition for reflection coefficient convolution matrix left singular matrix;

V _ifor reflection coefficient convolution matrix right singular matrix column vector;

X _ieach element for the column vector of seismic trace;

In step 7), the span of different rank k is 0 < k < n, and n is seismic wavelet sampling point number.

Seismic wavelet W in step 7) ^(k)formula is to obtain according to seismic wavelet objective function formula.

8) calculate energy and the seismic wavelet energy of seismic trace residual error, determine best k value:

Seismic trace residual energy $E_s^{\left(k\right)}={\left|\right|{\mathrm{RW}}^{\left(k\right)}-X\left|\right|}_2^2---\left(11\right),$

With seismic wavelet energy $E_w^{\left(k\right)}={\left|\right|W^{\left(k\right)}\left|\right|}_2^2;---\left(12\right)$

In formula:

for cutting the disconnected situation Xia Zhen of singular value road residual energy in k rank;

In formula:

for k rank, cut singular value and block situation seismic wavelet energy:

Seismic wavelet energy corresponding to different wavelets and seismic trace residual energy are done to X plot, choose seismic wavelet energy and seismic trace residual energy intersection flex point as best k value between the two;

9) utilize best k value and seismic wavelet W ^(k)formula determines that best seismic wavelet is definite;

10) to a plurality of seismic trace repeating steps 2)～9) obtain its corresponding best seismic wavelet, and then adopt multiple wavelet averaging method to determine a comprehensive seismic wavelet, be near seismic wavelet position, well point described in step 1).

Experiment embodiment of the present invention is as figure.Through composite traces calibrated time domain reflection coefficient as shown in Figure 1, this tests selected position, well point near 105 roads as shown in Figure 2 (earthquake sampling interval is 4ms).Fig. 3 is the theoretical wavelet (Ricker wavelet that crest frequency is 25hz) of previously demarcating for the synthesis of seismologic record, and Fig. 4 and Fig. 5 are previously for the synthesis of the spectral amplitude and the phase spectrum that record the seismic wavelet of demarcating;

The length of given definite wavelet is 88ms, selecting for the definite seismic trace of seismic wavelet is 104,105 ,106 tri-road geological datas as shown in Figure 4, the computation window scope of seismic section is: 320ms-420ms, and window length is 100ms, window length is slightly larger than given wavelet length; While obtaining for the definite time domain reflection coefficient of seismic wavelet according to preceding method, window scope is: 276ms-464ms, window length is 188ms, window length is slightly larger than 2 times of given wavelet length, and during seismic trace, window and reflection coefficient window length meet wavelet provided by the invention and determine the requirement that imposes a condition.

According to preceding method, with the 104th, 105 ,106 tri-road geological datas, determine seismic wavelet respectively, and then adopt three seismic wavelets that usual way obtains three seismic traces to carry out " comprehensively average " to obtain a comprehensive seismic wavelet.

Earthquake sampling interval is 4ms, and seismic wavelet length that be definite is 88ms, and seismic wavelet sampling point number is 23, and reflection coefficient convolution singular values of a matrix maximum order is 23 (not blocking).Fig. 6 utilizes 105 road geological data to carry out wavelet while determining, reflection coefficient convolution singular values of a matrix is carried out to different rank while blocking, the energy of a series of seismic wavelets that obtain and seismic trace residual energy X plot

as can be seen from Figure 6: along with K value increases (truncated singular value exponent number reduces), determine seismic wavelet W ^(k)energy

become large, seismic trace residual energy

diminish; Along with K value reduces (increase of truncated singular value exponent number), determine seismic wavelet W ^(k)energy diminish, seismic trace residual energy

become large.There is an obvious corner point in curve in the change procedure of K value, and this point is exactly the equilibrium point of seismic trace residual energy and seismic wavelet energy.The equilibrium point of the data that provide for this seismic trace residual energy and seismic wavelet energy when truncated singular value exponent number is 18 (corner point) wavelet.

Fig. 7 adopts wavelet provided by the invention to determine method, utilize time domain reflection coefficient as shown in Figure 1 and the 104th, 105,106 3 seismic wavelets that seismic trace is determined as shown in Figure 2, Fig. 8 and Fig. 9 are respectively spectral amplitude and the phase spectrums of determined seismic wavelet.From above three figure, can find out the seismic wavelet of determining by method provided by the invention, stable in time domain wave mode feature, there are a complete main lobe and two secondary lobes, its phase spectrum is vary stable in effectively frequent wide region.

Figure 10 is synthetic road and the seismic trace comparison diagram of making of the seismic wavelet shown in Fig. 7, and as can be seen from Figure 10 synthetic road and seismic trace have very strong correlativity, and related coefficient is 0.863.

Claims (12)

1. a method for definite seismic wavelet, feature is to realize by following steps:

1) gather log data and seismic channel data and priori seismic wavelet;

2) utilize the time domain reflection coefficient of priori seismic wavelet to make composite traces, time time domain reflection coefficient is demarcated again within the scope of window;

3) length of given seismic target earthquakes wavelet, window scope when window and reflection coefficient while determining seismic channel data;

4) build reflection coefficient convolution matrix;

5) build seismic trace column vector X;

6) according to following formula, determine seismic target earthquakes wavelet objective function;

$E=\{{\left|\right|R{W}^{\left(k\right)}-X\left|\right|}_2^2+{\left|\right|W^{\left(k\right)}\left|\right|}_2^2\}\&RightArrow;\min ---\left(9\right)$

In formula:

R is reflection coefficient convolution matrix;

W ^(k)expression is carried out k rank to reflection coefficient convolution matrix and is blocked seismic wavelet in situation;

X is the column vector of seismic trace;

K is the exponent number that blocks of reflection coefficient convolution singular values of a matrix;

7) according to following formula, calculate different rank k singular value and block the seismic wavelet W in situation ^(k):

$W^{\left(k\right)}=\underset{i=0}{\overset{k}{\mathrm{\&Sigma;}}}\frac{u_i^T{x}_i}{{\mathrm{\&delta;}}_i}{v}_i---\left(10\right)$

In formula:

δ _isingular value for reflection coefficient convolution matrix R;

column vector transposition for reflection coefficient convolution matrix left singular matrix;

V _ifor reflection coefficient convolution matrix right singular matrix column vector;

X _ieach element for the column vector of seismic trace;

8) calculate energy and the seismic wavelet W of seismic trace residual error ^(k)energy, determine best k value:

Seismic trace residual energy $E_s^{\left(k\right)}={\left|\right|{\mathrm{RW}}^{\left(k\right)}-X\left|\right|}_2^2---\left(11\right),$

With seismic wavelet energy $E_w^{\left(k\right)}={\left|\right|W^{\left(k\right)}\left|\right|}_2^2;---\left(12\right)$

In formula:

for cutting the disconnected situation Xia Zhen of singular value road residual energy in k rank;

In formula:

for k rank, cut singular value and block situation seismic wavelet energy:

Seismic wavelet energy corresponding to different wavelets and seismic trace residual energy are done to X plot, choose seismic wavelet energy and seismic trace residual energy intersection flex point as best k value between the two;

9) utilize best k value and seismic wavelet W ^(k)formula is determined best seismic wavelet;

10) to a plurality of seismic trace repeating steps 2)～9) obtain its corresponding best seismic wavelet, and then adopt multiple wavelet on average to obtain near the seismic target earthquakes wavelet in position, well point.

2. according to the method for claim 1, feature is that the log data described in step 1) mainly comprises: the terrestrial coordinate of well-name, this well position, the time domain reflection coefficient after composite traces accurate calibration.

3. according to the method for claim 1, feature is the seismic channel data described in step 1), is near the multichannel seismic data together or well point of earthquake-capturing.

4. according to the method for claim 1, feature is the seismic channel data described in step 1), refers to the geological data in window in the time of will determining seismic wavelet target interval; The time window top and bottom of seismic channel data will depart from crest and the trough of seismic trace, selects weak reflection position.

5. according to the method for claim 1, feature is that the priori seismic wavelet described in step 1) is previously for the synthesis of record, to demarcate known seismic wavelet or theoretical wavelet.

6. according to the method for claim 1, feature is step 2) in demarcation be again by synthetic road and each well lie data correlation analysis are obtained to reflection coefficient and each well lie best match position.

7. according to the method for claim 1, feature is step 2) in time window scope refer to when given earthquake within the scope of window traveling time territory reflection coefficient within the scope of positive and negative 30ms.

8. according to the method for claim 1, feature is that the seismic target earthquakes wavelet length described in step 3) is determined according to actual seismic data, and shallow earthquake wavelet length is short, and deep seismic wavelet length is long; The needs of window scope when time domain reflection coefficient can not meet step 3), reduce wavelet length, and the minimum length of wavelet length must comprise the integrality of main lobe and the secondary lobe of wavelet.

9. according to the method for claim 1, feature is that the seismic target earthquakes wavelet length described in step 3) is between 60--120ms.

10. according to the method for claim 1, when feature is the seismic trace described in step 3), window is greater than 1～1.5 times of seismic wavelet length, and when during reflection coefficient, window is than seismic trace, window scope is respectively expanded seismic wavelet length up and down half.

11. according to the method for claim 1, feature be in step 7) different k to cut on rank singular value span be 0 < k < n, n is seismic wavelet sampling point number.

12. according to the method for claim 1, and feature is seismic wavelet W in step 7) ^(k)energy theorem is to obtain according to seismic target earthquakes wavelet objective function formula.

Images