CN103576192B - A kind of method determining seismic wavelet - Google Patents

Abstract

The present invention is the method determining seismic wavelet in petroleum exploration seism processing, the bound term determining that seismic wavelet energy is minimum is added in original object function, consider the impact on determining seismic wavelet precision of the noise in seismic channel, by finding optimal " equilibrium point " of seismic channel residual energy and seismic wavelet energy, that determines reflection coefficient convolution singular values of a matrix most preferably blocks exponent number, by reflection coefficient convolution singular values of a matrix being blocked the noise effect weakening seismic channel, determine seismic wavelet.The present invention does not interferes with the contribution that wavelet is determined by bigger singular value, is more beneficial for weakening noise in seismic channel wavelet determining, the impact of precision, beneficially seismic wavelet determine the raising of precision compared with common method.

Description

A kind of method determining seismic wavelet

Technical field

The present invention relates to a kind of method determining seismic wavelet in petroleum exploration and development technology, specifically seism processing, explanation.

Background technology

The determination of seismic wavelet is a problem the most crucial in seismic prospecting.In forward problem, need to combine seismic wavelet by wave equation or convolution model and form forward simulation geological data.In inverting and deconvolution problem, it is also desirable to extracting a seismic wavelet by seismic channel, different seismic wavelets obtains different inversion results.

The basis of seismic wavelet extraction is convolution model, it include seismic wavelet, reflection coefficient, containing noisy seismic channel data.Methods of seismic wavelet extraction is to utilize well-log information to calculate reflection coefficient, is obtained seismic wavelet then in conjunction with seismic channel by convolution model.The most conventional well shake associating wavelet determines that method is based on synthesis road and seismic channel residual energy (hereinafter referred to as seismic channel 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 that window scope is at that time: t by well-log information calculated time domain reflection coefficient₁≤t≤t₂, x (t) is seismic channel amplitude, and window scope is at t at that time₁+p≤t≤t₂-q, w (t) are seismic wavelet to be determined, its beginning and end time scope is-p to q, formula (1) is converted into matrix operations form and can get following expression:

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

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

W=R^-1X(3)

Formula (3) is the most conventional well shake and combines the computing formula determining seismic wavelet, and R is reflection coefficient convolution matrix.Enactive earthquake wavelet sampling point number be n(n be odd number), reflection coefficient sampling point number is m, and 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& r_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 the reflectance value that each reflection sampling point is corresponding；

X is (m-n+1) rank column vector that seismic channel amplitude is constituted, 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 channel；

Matrix R is carried out singular value decomposition:

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

∑=diag (δ in formula (6)₀,δ₁, δ_n-1) it is n × n diagonal matrix, δ_iFor the singular value of reflection coefficient convolution matrix R, and δ₀≥δ₁≥…≥δ_n-1With descending order；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.

Obtain seismic wavelet in conjunction with (6), (3) and determine formula:

$W_i=\underset{i=0}{\overset{n}{\mathrm{\Σ}}}\frac{u_i^T{x}_i}{{\mathrm{\δ}}_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；Remaining each parameter physical meaning is with aforementioned.

From formula (7) it can be seen that work as singular value δ_iTo become the least,To become very large, at this moment small in seismic channel noise, and the result determining seismic wavelet will be produced large effect, result can make the seismic wavelet determined and actual wavelet there is very large deviation.

Summary of the invention

A kind of method being to provide determination seismic wavelet being conducive to seismic wavelet precision to improve mesh of the present invention.

The present invention implements step:

1) log data and seismic channel data and priori wavelet are gathered；

Log data described in step 1) specifically includes that well-name, the geodetic coordinates of this well position, time domain reflection coefficient after composite traces accurate calibration.

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

Seismic channel data described in step 1), refers to the geological data in window during seismic wavelet target interval to be determined；The time crest of window top and bottom seismic channel to be deviateed and trough of seismic channel data, selects weak reflection position.

Priori wavelet described in step 1) is to be previously used for the known seismic wavelet of composite seismogram calibration or theoretical wavelet.

2) utilize step 1) priori seismic wavelet and time domain reflection coefficient to make composite traces, in the range of to timing window, time domain reflection coefficient is re-scaled；

Step 2) in re-scale be by synthesis road and each well lie data correlation analysis obtain reflection coefficient and each well lie best match position.

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

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

Seismic wavelet length described in step 3) determines according to actual seismic data, and shallow earthquake wavelet lengths is short, and deep seismic wavelet lengths is long；The needs of window scope when the time domain reflection coefficient that step 1) obtains can not meet step 3), reduce wavelet lengths, the minimum length of wavelet lengths, it is necessary to include the main lobe of wavelet and the integrity of secondary lobe.

Between a length of 60--120ms of seismic wavelet described in step 3).

During seismic channel described in step 3), window is greater than 1～1.5 times of seismic wavelet length, and during reflection coefficient, window respectively extends half seismic wavelet length than window scope during seismic channel.

4) reflection coefficient convolution matrix is built；

5) seismic channel column vector X is built；

6) seismic wavelet object function is determined according to below equation structure；

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

In formula:

R is reflection coefficient convolution matrix；

W^(k)Represent and reflection coefficient convolution matrix is carried out seismic wavelet in the case of k rank are blocked；

X is the column vector of seismic channel；

K be reflection coefficient convolution singular values of a matrix block exponent number；

7) the seismic wavelet W in the case of different rank k singular value is blocked is calculated according to below equation^(k):

$W^{\left(k\right)}=\underset{i=0}{\overset{k}{\mathrm{\Σ}}}\frac{u_i^T{x}_i}{{\mathrm{\δ}}_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 channel；

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 Equation.

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

Seismic channel 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:Cut in the case of singular value is broken for k rank and shake road residual energy；

In formula:Cut singular value for k rank and block situation seismic wavelet energy:

The seismic wavelet energy corresponding to different wavelets and seismic channel residual energy do cross plot, choose seismic wavelet energy and seismic channel residual energy cross plot flex point as optimal k value between the two；

9) optimal k value and seismic wavelet W are utilized^(k)Formula determines that optimal seismic wavelet determines；

10) multiple seismic channels are repeated step 2)～9) obtain its corresponding optimal seismic wavelet, use multiple wavelet averaging method to determine a comprehensive seismic wavelet the most again, for the seismic wavelet near position, well point described in step 1).

The present invention adds the bound term determining that seismic wavelet energy is minimum in original object function, defines new wavelet and determine object function compared with determining method with traditional well shake associating wavelet.Compared with former object function, it is contemplated that the impact on determining seismic wavelet precision of the noise in seismic channel, beneficially wavelet determine the raising of precision.Obtained optimal " equilibrium point " of seismic channel residual energy and seismic wavelet energy by cross plot during implementing, that determines reflection coefficient convolution singular values of a matrix most preferably blocks exponent number, weakens the impact on determining seismic wavelet precision of the noise of seismic channel by the method blocking reflection coefficient convolution singular values of a matrix.

The present invention is to weaken seismic channel noise seismic wavelet determines the impact of precision by directly blocking the reflection coefficient less singular value of convolution matrix, do not interfere with the contribution that wavelet is determined by bigger singular value, compared with commonly using and adding whitening factor method, the present invention will more favorably improve wavelet determine precision.

Accompanying drawing explanation

Fig. 1 time domain reflection coefficient after composite seismogram calibration；

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

Fig. 3 is the seismic wavelet being previously used for composite seismogram calibration；

Fig. 4 and Fig. 5 is amplitude spectrum and the phase spectrum of the seismic wavelet being previously used for composite seismogram calibration；

Fig. 6 is to utilize the 105th road geological data to carry out wavelet when determining, synthesis road and well lie residual energy and the energy cross plot (transverse axis represents seismic channel residual energy, and the longitudinal axis represents seismic wavelet energy) of seismic wavelet；

Fig. 7 is to utilize the seismic wavelet that present invention determine that；

Fig. 8 and Fig. 9 is amplitude spectrum and the phase spectrum of determined seismic wavelet respectively；

Figure 10 is finally to determine that seismic wavelet is made to synthesize road and seismic channel comparison diagram.

Specific embodiments

The present invention is on the basis of original object function seismic channel residual energy minimum, adds the bound term determining that seismic wavelet energy is minimum, defines new seismic wavelet and determine object function.Compared with former object function, it is contemplated that seismic channel noise determines that on seismic wavelet the impact of precision, beneficially seismic wavelet determine the raising of precision.

The present invention to implement process as follows:

1) log data and seismic channel data and priori wavelet are gathered；

Log data described in step 1) specifically includes that well-name, the geodetic coordinates of this well position, time domain reflection coefficient after composite traces accurate calibration.

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

Seismic channel data described in step 1), refers to the geological data in window during seismic wavelet target interval to be determined；The time crest of window top and bottom seismic channel to be deviateed and trough of seismic channel data, selects weak reflection position.

Priori wavelet described in step 1) is to be previously used for the known seismic wavelet of composite seismogram calibration or theoretical wavelet.

2) utilize step 1) priori seismic wavelet and time domain reflection coefficient to make composite traces, in the range of to timing window, time domain reflection coefficient is re-scaled；

Step 2) in re-scale be by synthesis road and each well lie data correlation analysis obtain reflection coefficient and each well lie best match position.

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

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

Seismic wavelet length described in step 3) determines according to actual seismic data, and shallow earthquake wavelet lengths is short, and deep seismic wavelet lengths is long；The needs of window scope when the time domain reflection coefficient that step 1) obtains can not meet step 3), reduce wavelet lengths, the minimum length of wavelet lengths, it is necessary to include the main lobe of wavelet and the integrity of secondary lobe.

Between a length of 60--120ms of seismic wavelet described in step 3).

During seismic channel described in step 3), window is greater than 1～1.5 times of seismic wavelet length, and during reflection coefficient, window respectively extends half seismic wavelet length than window scope during seismic channel.

4) reflection coefficient convolution matrix is built；

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

5) seismic channel column vector X is built；

The method that step 5) seismic channel column vector provides according to formula formula (5) builds, x in formula (5)_iAmplitude for each sampling point of seismic channel；M, n(n are odd number) it is respectively reflection coefficient and the sampling point number of seismic wavelet；

6) seismic wavelet object function is determined according to below equation structure；

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

In formula:

R is reflection coefficient convolution matrix；

W^(k)Represent and reflection coefficient convolution matrix is carried out seismic wavelet in the case of k rank are blocked；

X is the column vector of seismic channel；

K be reflection coefficient convolution singular values of a matrix block exponent number；

7) the seismic wavelet W in the case of different rank k singular value is blocked is calculated according to below equation^(k):

$W^{\left(k\right)}=\underset{i=0}{\overset{k}{\mathrm{\Σ}}}\frac{u_i^T{x}_i}{{\mathrm{\δ}}_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 channel；

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 Equation.

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

Seismic channel 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:Cut in the case of singular value is broken for k rank and shake road residual energy；

In formula:Cut singular value for k rank and block situation seismic wavelet energy:

The seismic wavelet energy corresponding to different wavelets and seismic channel residual energy do cross plot, choose seismic wavelet energy and seismic channel residual energy intersection flex point as optimal k value between the two；

9) optimal k value and seismic wavelet W are utilized^(k)Formula determines that optimal seismic wavelet determines；

10) multiple seismic channels are repeated step 2)～9) obtain its corresponding optimal seismic wavelet, use multiple wavelet averaging method to determine a comprehensive seismic wavelet the most again, for the seismic wavelet near position, well point described in step 1).

Experiment embodiment of the present invention such as figure.Time domain reflection coefficient after composite seismogram calibration is as it is shown in figure 1, this tests selected position, well point (the earthquake sampling interval is 4ms) near the 105th road as shown in Figure 2.Fig. 3 is to be previously used for the theoretical wavelet (crest frequency is the Ricker wavelet of 25hz) that synthetic seismogram carries out demarcating, Fig. 4 and Fig. 5 is amplitude spectrum and the phase spectrum of the seismic wavelet being previously used for composite seismogram calibration；

The given a length of 88ms determining wavelet, selecting the seismic channel determined for seismic wavelet is 104 as shown in Figure 4,105,106 3 road geological datas, the computation window scope of seismic profile is: 320ms-420ms, and time window length is 100ms, and time window length is slightly larger than given wavelet lengths；According to preceding method obtain for the time domain reflection coefficient that seismic wavelet determines time window scope be: 276ms-464ms, time window length is 188ms, time window length is slightly larger than 2 times of given wavelet lengths, and the wavelet that during seismic channel, window and reflection coefficient time window length meet the present invention and provide determines the requirement that imposes a condition.

Determining seismic wavelet with the 104th, 105,106 3 road geological datas respectively according to preceding method, three seismic wavelets using usual way to obtain three seismic channels the most again carry out " average " and obtain a comprehensive seismic wavelet.

The earthquake sampling interval is 4ms, a length of 88ms of seismic wavelet to be determined, seismic wavelet sampling point number is 23, i.e. reflection coefficient convolution singular values of a matrix maximum order is 23 (not blocking).Fig. 6 is to utilize the 105th road geological data to carry out wavelet when determining, reflection coefficient convolution singular values of a matrix carries out different rank when blocking, the energy of a series of seismic wavelets obtained and seismic channel residual energy cross plot As can be seen from Figure 6: along with K value increases (reduction of truncated singular value exponent number), determine seismic wavelet W^(k)EnergyBecome big, seismic channel residual energyDiminish；Along with K value reduces (increase of truncated singular value exponent number), determine seismic wavelet W^(k)EnergyDiminish, seismic channel residual energyBecome big.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 channel residual energy and seismic wavelet energy.For these data provided seismic channel residual energy and equilibrium point of seismic wavelet energy when truncated singular value exponent number is 18 (corner point) wavelet.

Fig. 7 is that the wavelet using the present invention to provide determines method, utilizing time domain reflection coefficient as shown in Figure 1 and the 104th, 105,106 3 seismic wavelets that seismic channel is determined as shown in Figure 2, Fig. 8 and Fig. 9 is amplitude spectrum and the phase spectrum of determined seismic wavelet respectively.The seismic wavelet that the method provided by the present invention is determined is can be seen that from above three figures, stable at time domain Wave pattern, there are a complete main lobe and two secondary lobes, the change in the range of width the most again and again of its phase spectrum is stable.

Figure 10 is the synthesis road and seismic channel comparison diagram made of the seismic wavelet shown in Fig. 7, and as can be seen from Figure 10 synthesis road has the strongest dependency with seismic channel, and correlation coefficient is 0.863.

Claims (11)

1. the method determining seismic wavelet, feature is through the following steps that realize:

1) log data and seismic channel data and priori seismic wavelet are gathered；

2) the time domain reflection coefficient utilizing priori seismic wavelet makes composite traces, re-scales time domain reflection coefficient in the range of to timing window；

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

4) reflection coefficient convolution matrix is built；

5) seismic channel column vector X is built；

6) seismic target earthquakes wavelet object function is determined according to below equation；

$E=\left\{\right||{\mathrm{RW}}^{\left(k\right)}-X|{|}_2^2+\left|\right|W^{\left(k\right)}\left|{|}_2^2\right\}\→min---\left(9\right)$

In formula:

R is reflection coefficient convolution matrix；

W^(k)Represent and reflection coefficient convolution matrix is carried out seismic wavelet in the case of k rank are blocked；

X is the column vector of seismic channel；

K be reflection coefficient convolution singular values of a matrix block exponent number；

7) the seismic wavelet W in the case of different rank k singular value is blocked is calculated according to below equation^(k)；

$W^{\left(k\right)}=\underset{i=0}{\overset{k}{\Σ}}\frac{u_i^T{x}_i}{{\mathrm{\δ}}_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 channel；

8) energy and the seismic wavelet W of seismic channel residual error are calculated^(k)Energy, determines optimal k value；

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

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

In formula:Seismic channel residual energy in the case of blocking for k rank singular value；

In formula:Situation seismic wavelet energy is blocked for k rank singular value；

The seismic wavelet energy corresponding to different wavelets and seismic channel residual energy do cross plot, choose seismic wavelet energy and seismic channel residual energy intersection flex point as optimal k value between the two；

9) optimal k value and seismic wavelet W are utilized^(k)Formula determines optimal seismic wavelet；

10) multiple seismic channels are repeated step 2)～9) obtain its corresponding optimal seismic wavelet, use multiple wavelet averagely to obtain the seismic target earthquakes wavelet near position, well point the most again.

Method the most according to claim 1, feature is step 1) described in log data specifically include that well-name, the geodetic coordinates of this well position, time domain reflection coefficient after composite traces accurate calibration.

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

Method the most according to claim 1, feature is step 1) described in seismic channel data, refer to the seismic channel data in window during seismic wavelet target interval to be determined；The time crest of window top and bottom seismic channel to be deviateed and trough of described seismic channel data, selects weak reflection position.

Method the most according to claim 1, feature is step 1) described in priori seismic wavelet be to be previously used for the known seismic wavelet of composite seismogram calibration or theoretical wavelet.

Method the most according to claim 1, feature is step 2) in re-scale be by synthesis road and each well lie data correlation analysis are obtained reflection coefficient and each well lie best match position.

Method the most according to claim 1, feature is step 3) described in seismic target earthquakes wavelet lengths determine according to actual seismic data, shallow earthquake wavelet lengths is short, and deep seismic wavelet lengths is long；When time domain reflection coefficient can not meet step 3) reflection coefficient time window scope needs time, reduce wavelet lengths, the minimum length of wavelet lengths, it is necessary to include the main lobe of wavelet and the integrity of secondary lobe.

Method the most according to claim 1, feature is step 3) described in seismic target earthquakes wavelet lengths be between 60--120ms.

Method the most according to claim 1, feature is step 3) described in seismic channel data time window scope be greater than 1 times of seismic wavelet length, during reflection coefficient, window scope respectively extends half seismic wavelet length than window scope during seismic channel data.

Method the most according to claim 1, feature is step 7) in different k rank singular value to block span be 0 ＜ k ＜ n, n is seismic wavelet sampling point number.

11. methods according to claim 1, feature is step 7) in seismic wavelet W^(k)Energy theorem is to obtain according to seismic target earthquakes wavelet objective function Equation.